Autovalores e autovetores de uma imagem 3D Laplaciana

Preciso calcular os autovalores e autovetores de uma imagem em 3D da Lapônia. Estou tentando avaliar o núcleo de calor na grade uniforme 3D (a estrutura uniforme gerada pela imagem voxelizada) em diferentes valores de tempo, para implementar uma assinatura volumétrica de núcleo de calor (consulte a seção "Computação numérica"). O domínio em que estou trabalhando não é retangular, por isso tenho 1s em alguns pontos da grade e 0s em outros pontos, produzindo uma forma 3D binária. Eu tenho lido sobre isso no 2D: no 2D, o lapso de cada ponto da grade é calculado e os autovalores e autovetores da matriz resultante são calculados. Como posso fazer isso em 3D? Todas as informações e exemplos que li são para imagens em 2D. As perguntas específicas são:

1. Como posso obter uma matriz a partir da imagem 3D obtida após a aplicação de um operador Laplace discreto (3x3x3)?
2. Se não há como obter uma matriz, devo calcular os autovalores e autovetores de uma grade estruturada em 3D. Qual é a maneira mais simples de fazer isso? Você poderia me recomendar algum software para fazer isso?

Eu realmente aprecio qualquer dica.

Desde já, obrigado.

Essa resposta chega um pouco tarde, mas acho que é necessário esclarecer um pouco da confusão sobre o que é a estrutura própria de um Laplaciano e como é calculada.

Antes de tudo, é importante enfatizar que não se trata de propriedades do kernel local usado para calcular derivadas discretas. Em vez disso, você precisa entender o Laplaciano como operador linear em um espaço vetorial cujos elementos são, no seu caso, os conjuntos de dados volumétricos.

Isso significa que a aplicação do operador Laplace é um mapa linear que mapeia um vetor (conjunto de dados) para outro vetor (conjunto de dados). Nesse contexto, os vetores próprios do Laplaciano são novamente vetores (conjuntos de dados) no mesmo espaço vetorial. Portanto, a resposta para sua pergunta seria um conjunto de conjuntos de dados volumétricos, na verdade até a dimensão do seu espaço vetorial (ou seja, o número de voxels independentes ).

Vamos considerar um exemplo muito simples. Tire uma imagem unidimensional, ou seja, uma única linha de pixels, e vamos usar apenas muito poucos pixels, ou seja, 4. Então, os dois pixels centrais precisam direcionar os vizinhos, enquanto o primeiro e o último pixel têm apenas um vizinho cada.

$$1 - 2 - 3 - 4$$

Com essa geometria de pixels, podemos fornecer o resultado discreto do operador Laplace para os dois pixels centrais 2 e 3 como funções lineares dos valores dos pixels:

$$l[2] = p[1] - 2*p[2] + p[3]$$ $$l[3] = p[2] - 2*p[3] + p[4]$$

Os outros dois pixels 1 e 4 não têm vizinhos diretos suficientes para calcular a segunda derivada discreta. Poderíamos consertar isso assumindo que os pixels 1 e 4 são vizinhos diretos, fechando a topologia em um círculo e impondo o que é chamado de condição de contorno circular. Ou simplesmente pegamos os segundos derivativos nos limites para desaparecer. Vamos fazer as duas coisas, mas comece com a condição de contorno cíclico. Então nós temos:

$$l[1] = p[4] - 2*p[1] + p[2]$$ $$l[4] = p[3] - 2*p[4] + p[1]$$

Este mapa é linear e podemos escrevê-lo como uma equação matricial, mapeando o vetor da coluna $p := (p[1],p[2],p[3],p[4])$ para $l = (L[1],L[2],L[3],L[4])$ por multiplicação com a matriz $M$.

$$L = M \cdot p$$

Chamamos essa matriz de representação discreta do operador Laplace e, para o nosso caso, é

$$M= \left( \begin{array}{cccc} -2 & 1 & 0 & 1 \\ 1 & -2 & 1 & 0 \\ 0 & 1 & -2 & 1 \\ 1 & 0 & 1 & -2 \\ \end{array} \right)$$

Os autovetores dessa matriz são

$$v_1 = (-1,1,-1,1)$$ $$v_2 = (0,-1,0,1)$$ $$v_3 = (-1,0,1,0)$$ $$v_4 = (1,1,1,1)$$ com os autovalores associados $$\lambda_1 = 4 , \lambda_1 = 2 , \lambda_1 = 2 , \lambda_1 = 0$$

Você pode reconhecer esses vetores como o vetor base da transformada de Fourier discreta nesse vetor e os valores próprios como suas frequências discretas. Isso é verdade em geral e, de fato, a decomposição de um vetor (ou mais geralmente, uma função) no espectro eigens de um operador de Laplace generaliza a idéia da transformação de Fourier.

Agora vamos investigar o que acontece se usarmos as condições de contorno alternativas em que $l[1] = 0$ e $l[4]=0$. O Matrix$M$ é então

$$M=\left( \begin{array}{cccc} 0 & 0 & 0 & 0 \\ 1 & -2 & 1 & 0 \\ 0 & 1 & -2 & 1 \\ 0 & 0 & 0 & 0 \\ \end{array} \right)$$

Certamente, essa matriz tem um conjunto diferente de vetores próprios e valores próprios. Eles não são tão intuitivos e interessantes, então não os listarei explicitamente. No entanto, vale a pena notar que agora obtemos o valor próprio$0$ duas vezes, isso significa que o espaço próprio onde o laplaciano desaparece é bidimensional.

Então, como as coisas mudam se tivermos uma imagem adequada em vez de uma única linha de pixels? Não muito. Precisamos apenas escrever o Laplaciano para cada pixel único, considerando a relação direta com o vizinho, ou topologia, da imagem. Para tornar as coisas um pouco mais complicadas, vamos com uma imagem bidimensional de formato irregular.

$$\begin{array}{ccccccccc} && 1 &-& 2 &-& 3 && \\ && | & & | & & | \\ 4 &-& 5 &-& 6 &-& 7 &-& 8 \\ | & & | & & | & & | & & | \\ 9 &-& 10 &-& 11 &-& 12 &-& 13 \\ && | & & | && | & & \\ && 14 &-& 15 &-& 16 && \\ \end{array}$$

Obviamente agora temos que pegar um laplaciano bidimensional somando as segundas derivadas parciais nas direções horizontal e vertical. Para isso, exigimos os dois vizinhos diretos em cada direção. Os pontos internos$5,6,7,10,11,12$portanto, possui uma expansão laplace 2-d completa. Para ponto$5$ fica assim por exemplo:

$$l[5] = p[4] - 2 p[5] + p[6] + p[1] - 2 p[5] + p[10] \\= p[4] + p[6] + p[1] + p[10] - 4 p[5]$$

Para os pontos de canto $1,3,4,8,9,13,14,16$ como não podemos construir uma segunda derivada discreta, usamos uma condição de contorno, por exemplo $l[1]=0$

Restam dois pontos, $2$ e $15$. Ambos têm dois vizinhos diretos na direção horizontal, mas não na direção vertical. Podemos, portanto, aplicar uma condição de contorno que afeta apenas a direção vertical, definindo a segunda derivada vertical como zero, enquanto avaliamos a segunda derivada discreta horizontalmente e obtemos$l[2]=p[1]-2p[2]+p[3]$ e da mesma forma para $l[15]$.

Após essa construção, obtemos uma equação linear para cada ponto que a relaciona aos valores de pixel. Novamente, escrevemos como uma equação matricial$L=Mp$, onde a matriz neste caso tem $16\times 16$entradas. Para ser específico, é

$$M=\left( \begin{array}{cccccccccccccccc} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & -2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \end{array} \right)$$

E, novamente, podemos resolver o sistema próprio dessa matriz. Uma boa interpretação física das imagens que você obtém como autovetores é que elas representam os modos vibracionais de uma membrana em forma de imagem, com frequências dadas pelos valores próprios.

Você pode facilmente escalar este jogo para qualquer número de dimensões, desde que conheça as relações de vizinhança entre seus voxels. Simplesmente formule a equação linear individual como acima, construa uma matriz, encontre o sistema próprio.

Com as informações obtidas dos vetores próprios do Laplaciano, as equações das diferenças envolvendo o Laplaciano discreto podem ser bastante simplificadas. Uma vez que a estrutura própria é encontrada, dependendo apenas da geometria da região, todos os conjuntos de dados podem ser facilmente decompostos na base própria e as equações das diferenças se tornam triviais.

Não está claro o que você está tentando fazer: o Laplaciano é um valor escalar, não uma matriz. Portanto, não possui valores próprios. Talvez você queira dizer os autovalores / vetores do Hessiano?

Para calcular os valores próprios e o vetor próprio do Hessian, você deve primeiro calcular o Hessian (uma matriz 3x3 simétrica, contendo as segundas derivadas em cada uma das três direções) para cada pixel. Uma maneira simples de fazer isso é aplicar três filtros de gradiente (na direção x, y, z) à sua imagem 3d. Agora você tem 3 imagens. Para cada uma dessas imagens, você aplica novamente três filtros de gradiente, fornecendo imagens de resultados 3x3, contendo as entradas 3x3 da matriz Hessian para cada pixel. (Como os filtros de gradiente são comutados, I * gx * gy = I * gy * gx, você realmente só precisa calcular 6 deles.)

Há duas maneiras de obter os autovalores e os vetores dessa matriz 3x3: Encontre as raízes do polinômio característico ou use um método iterativo. O polinômio característico é cúbico, portanto, acho que os métodos iterativos seriam mais rápidos na prática (porque as raízes do cubo demoram muito mais para calcular a adição / multiplicação simples).

Talvez isso ajude (?).

http://en.m.wikipedia.org/wiki/Kronecker_sum_of_discrete_Laplacians

L = Dxx ¤ I ¤ I + I ¤ Dyy ¤ I + I ¤ I ¤ Dzz

onde ¤ é o produto Kronecker (não é o símbolo adequado, não sabia como obtê-lo aqui). Penso que os autovalores e autovetores são calculados a partir da matriz resultante por voxel.

O conceito de tensor, conforme descrito na seção 4 do artigo seguinte, parece relacionado, embora não seja exatamente o mesmo. Está, no entanto, intimamente relacionado à matriz de Hessian, conforme mencionado na resposta de @ nikie.

Carsten Steger, extração sub-precisa de linhas e arestas, arquivos internacionais de fotogrametria e sensoriamento remoto (2000).

O artigo pode ser encontrado aqui: http://ias.in.tum.de/people/steger/publications