# Eigenvalues and eigenvectors of a 3D image Laplacian

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.

• As pointed out by @nikie, please, clarify you're question. Do you want to compute 3D derivatives for the heat-distribution , i.e. generate a 3D vector-field over the grid ? – Dr.D. Jun 23 '13 at 10:49
• Thank you Dr.D. I've added more information to the question and a reference to the paper I'm trying to implement. – Federico Jun 23 '13 at 15:13
• I have no idea what a volumetric hear signature is, so I may be wrong, but I think the article you linked to already is about 3D shapes, isn't it? The way I understand it, the Laplacian can be seen as a matrix multiplication (since it's linear), and you're supposed to find the eigenvalues of that matrix (i.e. a matrix with one row for every voxel and one column for every voxel). – Niki Estner Jun 23 '13 at 19:04
• I think he is trying to find the eigenvalues of the laplacian in 3D. (but I may be wrong) – Beni Bogosel Dec 2 '13 at 12:47

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.

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).

• I'll state my problem hoping that it clarifies the question. I'm trying to evaluate the heat kernel on a 3D uniform grid at different time values. The volume 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 in 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? – Federico Jun 20 '13 at 21:56
• I see. Please update your question with this information, so I can delete my answer (it obviously doesn't answer the question). When you update your question, please note that a "2D matrix" usually means a matrix with 2 rows and 2 columns; a 3d matrix is a matrix with 3 rows and 3 columns. I think what you mean is not a 3d matrix, but a 3rd order tensor, which is not a matrix at all. – Niki Estner Jun 21 '13 at 11:09
• Actually, I think, you're answer may well be relevant, it's the same notion that I pointed to with that publication. It's just not quite clear, how it is applicable to the OPs 3D grid related problem. – Dr.D. Jun 23 '13 at 10:43

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