Eigenvalues and Eigenvectors of a 3D Image Laplacian

Question

I need to compute the eigenvalues and eigenvectors of a 3D image Laplacian. I'm trying to evaluate the heat kernel on the 3D uniform grid (the uniform structure generated by the voxelized image) at different time values, to implement a Volumetric Heat Kernel Signature (please see the "Numerical computation" section). The domain I'm working on is not rectangular, so I have 1s on some grid points and 0s on other points, producing a binary 3D shape. I've been reading about this for 2D: in 2D, the laplacian of every grid point is calculated and the eigenvalues and eigenvectors of the resulting matrix are then calculated. How can I do this for 3D?All the information and examples I have read are for 2D images. The specific questions are:

1. How can I obtain a matrix from the 3D image I obtain after applying it a (3x3x3) discrete Laplace operator?
2. If there is no way to obtain a matrix, then I should compute the eigenvalues and eigenvectors of a 3D structured grid. What is the simplest way to do this? Could you recommend me any soft for doing this?

I would really appreciate any hint.

Thanks in advance.

Answer 1

This answer comes a little late, but I think that it's necessary to clear up some of the confusion about what the eigenstructure of a Laplacian is and how it is calculated.

First of all, it's important to stress that this is not about properties of the local kernel used for calculating discrete derivatives. Instead you have to understand the Laplacian as linear operator on a vector space whose elements are in your case the volumetric data sets.

That means applying the Laplace operator is a linear map that maps one vector (data set) to another vector (data set). In this context, the eigenvectors of the Laplacian are again vectors (data sets) in that very same vector space. So the answer to your question would be a set of volumetric data sets, in fact as many as the dimension of your vector space (i.e the number of independent voxels).

Let's consider a very simple example. Take a 1-dimensional image, i.e. a single row of pixels, and let's also only use very few pixels, namely 4. Then the two center pixels have to direct neighbours, while the first and last pixel only have one neighbour each.

$$1 - 2 - 3 - 4$$

With this pixel geometry, we can give the discrete Laplace operator's result for the two center pixels 2 and 3 as linear functions of the pixels' values:

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

The other two pixels 1 and 4 don't have enough direct neighbours to calculate the discrete second derivative. We could fix that by assuming that pixels 1 and 4 are direct neighbours, closing the topology to a circle and imposing what is called a circular boundary condition. Or we just take the second derivatives at the boundaries to vanish. Let's do both, but start with the cyclic boundary condition. So we have:

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

This map is linear and we can write it as a matrix equation, mapping the column vector $p := (p[1],p[2],p[3],p[4])$ to $l = (L[1],L[2],L[3],L[4])$ by multiplication with the matrix $M$.

$$ L = M \cdot p $$

We call this matrix the discrete representation of the Laplace operator, and for our case it is $$M= \left( \begin{array}{cccc} -2 & 1 & 0 & 1 \\ 1 & -2 & 1 & 0 \\ 0 & 1 & -2 & 1 \\ 1 & 0 & 1 & -2 \\ \end{array} \right) $$

The eigenvectors of this matrix are $$ v_1 = (-1,1,-1,1) $$ $$ v_2 = (0,-1,0,1) $$ $$ v_3 = (-1,0,1,0) $$ $$ v_4 = (1,1,1,1) $$ with the associated eigenvalues $$ \lambda_1 = 4 , \lambda_1 = 2 , \lambda_1 = 2 , \lambda_1 = 0 $$

You may recognise these vectors as the basis vector of the discrete Fourier transform on this vector, and the eigenvalues as their discrete frequencies. This is true in general, and in fact the decomposition of a vector (or more generally, a function) into the eigenspectrum of a Laplace operator generalises the idea of the Fourier transform.

Now let's investigate what happens if we use the alternative boundary conditions where $l[1] = 0$ and $l[4]=0$. The Matrix $M$ is then

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

Cearly, this matrix has a different set of eigenvectors and eigenvalues. They're not as intuitive and insteresting, so I won't list them explicitly. It's worth noting however that we now get the eigenvalue $0$ twice, that means the eigensubspace where the laplacian vanishes is two dimensional.

So how do things change if we have a proper image instead of a single row of pixels? Not much. We only need to write down the Laplacian for each single pixel while considering the direct neighbour relationship, or topology, of the image. To make things a little trickier, let's go with an irregularly shaped two dimensional image. $$ \begin{array}{ccccccccc} && 1 &-& 2 &-& 3 && \\ && | & & | & & | \\ 4 &-& 5 &-& 6 &-& 7 &-& 8 \\ | & & | & & | & & | & & | \\ 9 &-& 10 &-& 11 &-& 12 &-& 13 \\ && | & & | && | & & \\ && 14 &-& 15 &-& 16 && \\ \end{array} $$

Obviously now we have to take a 2-dimensional Laplacian by summing the second partial derivatives in horizontal and vertical direction. For that we require the two direct neighbours in each direction. The inner points $5,6,7,10,11,12$ therefore have a full 2-d laplace expansion. For point $5$ it looks like this for example:

$$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]$$

For the corner points $1,3,4,8,9,13,14,16$ we cannot construct a discrete second derivative so we use a boundary condition, e.g. $l[1]=0$

There are two points left, $2$ and $15$. Both have two direct neighbours in horizontal direction, but not in vertical direction. We can therefore apply a boundary condition that only affects the vertical direction, by setting the vertical second derivative to zero, while we evaluate the discrete second derivative horizontally and get $l[2]=p[1]-2p[2]+p[3]$ and likewise for $l[15]$.

Following this construction we get a linear equation for every point that relates it to the pixel values. Again, we write it as a matrix equation $L=Mp$, where the matrix in this case has $16\times 16$ entries. To be specific, it's $$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)$$

And again we can solve for the eigensystem of this matrix. A nice physical interpretation of the images that you get as eigenvectors is that they represent the vibrational modes of a membrane shaped like your image, with frequencies given by the eigenvalues.

You can easily step this game up to any number of dimensions, as long as you know the neighbourhood relations between your voxels. Simply formulate the individual linear equation like above, construct a matrix, find the eigensystem.

With the information gained from the eigenvectors of the Laplacian solving difference equations involving the discrete Laplacian can be greatly simplified. Once the eigenstructure is found, depending only on the geometry of the region, all data sets can be easily decomposed into the eigenbasis and the difference equations become trivial.

Answer 2

Maybe this helps (?).

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

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

where ¤ is the Kronecker product (it's not the proper symbol, didn't know how to get it here). I think the eigenvalues and eigenvectors are computed from the resulting matrix per voxel.

The tensor-concept as described in the following paper Section 4 seems related, though it's not exactly the same. It is, however, closely related to the Hessian matrix as mentioned @nikie's answer.

Carsten Steger, Subpixel-Precise Extraction of Lines and Edges, International Archives of Photogrammetry and Remote Sensing (2000).

The paper can be found here: http://ias.in.tum.de/people/steger/publications

Answer 3

It's not clear what you're trying to do: The Laplacian is a scalar value, not a matrix. So it has no eigenvalues. Maybe you mean the eigenvalues/-vectors of the Hessian?

To calculate the eigenvalues and eigenvector of the Hessian, you would first calculate the Hessian (a symmetric 3x3 matrix, containing the second derivatives in each of the 3 directions) for each pixel. A simple way to do this is to apply three gradient filters (in x,y,z direction) to your 3d image. Now you have 3 images. To each of those images, you again apply three gradient filters, giving 3x3 results images, containing the 3x3 entries of the Hessian matrix for each pixel. (Since the gradient filters commute, I * gx * gy = I * gy * gx, you really only have to calculate 6 of them.)

There are two ways to get the eigenvalues and -vectors of this 3x3 matrix: Either find the roots of the characteristic polynomial or use an iterative method. The characteristic polynomial is cubic, so I'm guessing iterative methods would be faster in practice (because cube roots take far longer to calculate the simple addition/multiplications).