Uma abordagem avançada para detectar limites usando as paredes dos navios como exemplo

Informação interessante

A figura abaixo mostra uma reconstrução tridimensional do coração, obtida como resultado do trabalho de um tomógrafo moderno:


Para escala, a espessura da lâmpada aórtica é de 3,2 cm, pense! No entanto, quando as pessoas têm problemas cardíacos devido a vasos sanguíneos, então, como regra, não estamos falando sobre esses grandes. A imagem mostra que o coração está cercado por vasos menores, e alguns deles se ramificam diretamente das grandes artérias. Estas são as chamadas artérias coronárias que alimentam o coração diretamente com sangue. Se ocorrer um estreitamento do lúmen (estenose), por exemplo, devido à formação de cálcio, o fluxo sanguíneo diminui. Quando a estenose é pronunciada, ocorre necrose tecidual, ou seja, um ataque cardíaco. A seguir, falarei sobre nossa abordagem para calcular os limites dos vasos sanguíneos, o que, como resultado, permite que você encontre automaticamente o estreitamento e faça uma estimativa.

Para entender o material, você precisa ter um entendimento superficial do volume, voxels e suas intensidades. Você pode descobrir lendo o começo deste artigo .

Para avaliar o estreitamento de uma embarcação, precisamos conhecer o lúmen da embarcação ou sua borda interna. Para fazer isso, pelo menos, detecte todo o cálcio. Também encontramos o limite externo, porque permite estimar a espessura da parede, o que também é útil. Para começar, vamos dar uma olhada no esquema completo de detecção de limites e, em seguida, analisaremos cada estágio em detalhes:

Traçando a linha central

O estágio mais difícil em termos de implementação (pelo menos levou mais tempo). O método baseia-se no uso de uma matriz Hessian (segmentação de navios em várias escalas usando aprimoramento da matriz Hessian). Mais detalhes no artigo já mencionado.

Fatiar

Temos apenas uma linha central e precisamos de intensidades espacialmente dependentes de voxels que possam ser convenientemente trabalhadas. Para obtê-los, uma “pilha” de fatias está indo. Para iniciantes, os pontos são definidos através de distâncias fixas na linha central. Então uma perpendicular é construída a partir de cada ponto . Depois é o segundo perpendicular . Onde - a direção da linha central em um ponto. As duas perpendiculares são normalizadas. Em todo ponto vetor de linha central formar um sistema de coordenadas 2D. Assim, as fatias são formadas:


A posição do voxel é definida como onde aqui São as coordenadas reais do voxel, k é o número da fatia. Fórmula inversa para coordenadas reais: . Ao mover para um novo sistema de coordenadas, o espaço formado pelas fatias é simplificado:


Do que precisamos!

Construção da borda externa do navio

Vamos dar uma olhada no diagrama:


Cortamos nossa pilha de fatias obtidas no estágio anterior em oito planos (semelhante ao corte de um bolo) e realizamos todos os cálculos no espaço dos planos:


Se você exibir os valores normalizados das intensidades dos voxels que atingem o plano de corte, obteremos a seguinte imagem:


Para detectar os limites do navio, a abordagem clássica (detecção de arestas por gradiente) é usada em conjunto com a busca do caminho. Esquema:


1. Aplique a suavização gaussiana com um valor pequeno para suprimir ruídos:
Para um ponto com coordenadas :
onde retorna o valor da intensidade em um ponto ; r normalmente leva o valor ( - arredondamento); - coeficiente de suavização.

2. Em cada ponto, encontramos o gradiente e seu valor, os cálculos são realizados com valores de intensidade suavizada:
,
onde - derivativos privados. Eles são encontrados pelo método das diferenças finitas:

,
onde - intensidade em um ponto depois de alisar.
Em seguida, precisamos da direção do gradiente ( É a normalização do vetor) e o valor do gradiente ( É o comprimento do vetor)

3. A direção do gradiente é traduzida em graus ou radianos:
(atan2 () é a função do arco tangente em C ++, que não deve ser confundida com atan ()), arredondamos o ângulo para que ele possa ter 4 valores em incrementos de 45 graus, ou seja, superior e inferior são considerados na mesma direção:

4. Execute a supressão de valores não máximos. Se o valor do gradiente pelo menos em um dos dois pontos vizinhos (de acordo com a direção do gradiente) for maior ou igual ao valor do gradiente no ponto atual, esse ponto não poderá pertencer ao limite:


5. Todos os voxels restantes são considerados limites. Com base no valor do gradiente, no limiar de cálcio (que não está disponível imediatamente) e na proximidade do centro "vertical", cada ponto recebe um determinado custo (quanto mais brilhante o voxel, maior a prioridade ao procurar um caminho):

Nesta forma, os limites da embarcação são quase exclusivamente definidos.

6, 7. Para construir limites, usamos a pesquisa do caminho. Os pontos extremos mais próximos e com o menor custo são tomados como inicial e final. Para procurar um caminho, use uma pesquisa simples da primeira largura, que seleciona os pontos de limite de menor custo. Saltos também estão disponíveis, mas seu preço é alto. Os limites superior e inferior do vaso são pesquisados ​​separadamente e a suavização é aplicada a eles:


Este procedimento é realizado para cada plano, o que resulta em anéis de dezesseis segmentos para cada fatia na pilha. Esses anéis formam as bordas externas do navio.

Como você pode ver na imagem, há áreas nas quais as bordas são detectadas incorretamente. Isso ocorre devido à presença de cálcio, que resulta na detecção de limites de cálcio em vez de limites de vasos. Para impedir que isso aconteça após a primeira detecção de limites, é necessário determinar o limiar de cálcio (mais sobre isso posteriormente) e, em seguida, executar a segunda detecção de limites, ignorando os voxels relacionados ao cálcio. Então temos:

Detecção de limiar

Agora considere a maximização da expectativa do algoritmo de cluster (a seguir denominada EM). Vamos começar com a função de distribuição normal: é caracterizada por uma expectativa matemática e desvio padrão . Veja como eles afetam o tipo de distribuição:



Suponha que tenhamos dados (pontos) provenientes de uma fonte "amarela" e de uma fonte "azul":


Em seguida, usando fórmulas padrão, encontramos facilmente a média e o desvio padrão para cada fonte. Fórmulas para a fonte "a":




Mas e se soubermos o número de fontes, mas não soubermos quais pontos pertencem a qual fonte? E se tivermos uma foto dessas?


Se alguém viesse nos contar os parâmetros de distribuição, poderíamos calcular facilmente a probabilidade de cada ponto pertencente a cada uma das fontes. A probabilidade de o ponto pertencer à fonte "a":

onde




E se precisarmos calcular os parâmetros de origem, conhecendo as probabilidades:



Assim, é obtido um círculo vicioso: se conhecêssemos os parâmetros das fontes, calcularíamos qual ponto pertence a qual fonte, mas não conhecemos os parâmetros. E se soubéssemos qual ponto pertence a qual fonte, calcularíamos seus parâmetros, mas não sabemos qual ponto pertence a qual fonte. Equilibrar esses fatos é exatamente o que o algoritmo EM faz.

Na inicialização, o EM recebe alguns parâmetros predefinidos para fontes, que podem ser simplesmente selecionados aleatoriamente. Obviamente, se os parâmetros forem conhecidos, podemos calcular a probabilidade de cada ponto pertencente a cada uma das fontes. Agora que as probabilidades são conhecidas, podemos calcular novos parâmetros mais precisos. Então, tudo começa novamente, mas com novos parâmetros. Após cada ciclo, os parâmetros das fontes estão se tornando mais precisos.

Como podemos usar esse conhecimento em relação aos navios? Vamos dar uma olhada na estrutura de um deles:

Nos diagramas, esse ponto geralmente é omitido, mas formações de gordura e cálcio também podem estar presentes no vaso. Assim, cada voxel pertencerá a um dos tecidos. Após os experimentos, verificou-se que basta fazer as seguintes divisões:

- gordo
- parede # 1
- parede # 2
- meio de contraste
- cálcio

A distribuição das intensidades do voxel em cada caso é normal. I.e. temos tudo o que precisamos para usar o EM para encontrar os parâmetros de cada fonte.


Agora que conhecemos os parâmetros de cada fonte, podemos calcular os limites - os valores das intensidades, na interseção das quais, a associação ao voxel muda de uma fonte para outra. Estamos interessados ​​em:

1. O limiar do limite externo do navio. Se a intensidade do voxel estiver abaixo desse valor, considera-se que ele não pertence ao vaso;

2. O limiar do limite interno do vaso. Se a intensidade do voxel for maior que esse valor, então
refere-se ao lúmen do vaso, isto é, a uma mistura de sangue e meio de contraste;

3. O limiar de cálcio. Se o valor da intensidade do voxel for maior que esse valor, ele se refere ao cálcio.

Construção da borda interna do navio

Como sempre, vamos começar com o diagrama; os cálculos são realizados para cada fatia.



Se você exibir visualmente os dados de acordo com os limites obtidos na etapa anterior, obteremos a seguinte imagem:



A cor vermelha é a parede do vaso. Cor verde - folga. Branco é cálcio.

A primeira coisa que chama sua atenção é o cálcio "pendurado", que não fica ao lado de nenhuma das paredes. Isso é considerado normal e surge devido à suavização aplicada pelo próprio tomógrafo.

Primeiro, você precisa obter os limites de acordo com os limites e, para isso, o algoritmo de quadrados de marcha é usado. Pode-se dizer, passa em duas etapas. Primeiro, a área é dividida por uma grade discreta e os quadrados em que os valores de intensidade são maiores ou iguais ao limiar são considerados "positivos", o restante é considerado "negativo".

Cada vez que estamos em algum tipo de nó, precisamos desenhar um contorno em torno dos quadrados "positivos". Para tomar uma decisão, consideraremos os sinais de quatro quadrados vizinhos: superior esquerdo, inferior esquerdo, superior direito, inferior direito. Excluindo a simetria, estamos interessados ​​em três casos.

1. Três quadrados de um sinal e um oposto, o movimento do contorno ocorre na diagonal:


2. Dois quadrados do mesmo sinal e dois opostos, e os quadrados com o mesmo sinal estão em um lado, o movimento do contorno é vertical ou horizontal:


3. Dois quadrados do mesmo sinal e dois opostos, quadrados com os mesmos sinais são colocados em lados opostos:


Este é um caso excepcional, para tomar uma decisão, o valor médio da intensidade em todos os quatro quadrados é obtido e, se for maior ou igual ao limiar, o centro é positivo; em outros casos, é negativo. Também é importante qual nó está atual no momento:


O algoritmo dos quadrados de marchas constrói um contorno de maneira precisa e inequívoca. No exemplo abaixo, eu deliberadamente mudei a linha do centro do lado para que cada passo fosse claramente visível.


Para cada seção do navio, encontramos dois contornos principais - este é o contorno da borda externa e o contorno da borda interna. Imediatamente "cortamos" o contorno externo com o outro contorno externo, encontrado no início do artigo, procurando caminhos. Isso é feito para ignorar os galhos da embarcação. Ignoramos os contornos do cálcio que estão muito distantes da parede interna, como se eles não existissem, encontramos o resto e os usamos no futuro. Se o centro do vaso estiver dentro do cálcio, então o deslocamos na direção mais próxima do circuito do cálcio até que o centro esteja no lúmen (na região verde). Um centro tão atualizado, vou chamar a posição inicial.

De acordo com o esquema, pode haver dois casos: simples e complexo.


Se a posição inicial estiver dentro de um loop fechado de cálcio (por exemplo, se um stent estiver instalado), imediatamente equipararemos a borda interna a esse loop. As coisas ficam mais complicadas quando o centro está fora do cálcio. Nesse caso, precisamos expandir o circuito de partida para que ele flua suavemente ao redor do cálcio e da borda interna:


Para alcançar o resultado desejado, um algoritmo especial foi desenvolvido com base nas idéias usadas nos motores 2d físicos, em particular na resolução de colisões de polígonos e no teorema do eixo separado.

Dois conceitos básicos que não podem ser dispensados: para a interseção de polígonos convexos, o vetor mtv (vetor mínimo de conversão) é o menor deslocamento de um dos polígonos, após o qual a interseção para.

Também precisamos das normais poligonais - em 2D elas são perpendiculares às faces e indicam "fora":

Para não prolongar o artigo, omitirei o restante dos detalhes sobre os motores físicos. A única coisa que observo é que, durante cada iteração, cada ponto do contorno acumula a influência de forças sobre si mesmo na forma de um vetor e, no final de cada iteração, ele muda pelo comprimento desse vetor em sua direção. Considere as forças:

1. Duas forças agem em cada vértice do contorno inicial na direção dos vértices vizinhos, e as forças são diretamente proporcionais à distância desses vértices. Isso faz o contorno encolher e se esforça para manter uma forma arredondada;

2. Se o vértice do contorno lateral cair dentro do contorno inicial, o deslocamento proporcional ao vetor vértice mtv do vértice será aplicado à face mais próxima do contorno inicial;

3. Se a parte superior do contorno inicial estiver dentro do contorno lateral, será aplicado um deslocamento proporcional ao vetor de vértice mtv do vértice na parte superior do contorno inicial. Isso, juntamente com o parágrafo anterior, não permite que o circuito ultrapasse os limites de outros circuitos;

4. Se os casos 2 e 3 não funcionaram para o vértice, então uma força é aplicada na direção média das duas normais das faces adjacentes. Isso garante o crescimento do contorno "em largura".

Importante: você não pode mudar completamente um vértice ou uma face pelo comprimento do vetor mtv - ele é multiplicado por um certo coeficiente no intervalo de 0,2 a 0,8. Os coeficientes para cada força nos demais casos são selecionados experimentalmente.

Graças a essa abordagem, encontramos o lúmen do vaso, levando em consideração o fato de que o cálcio não se aproxima das paredes. Agora basta combinar os resultados de todas as seções e obter as bordas interna e externa da embarcação:


A aparente imprecisão da borda após o stent é causada por bifurcação anômala:


A área da borda interna no corte caracterizará a mesma folga em que estamos interessados. Uso adicional desses dados considero trivial e não requer consideração. Por fim, deixarei a imagem da borda interna exportada em 3D: