Eigenvalues and eigenvectors of a 3D image Laplacian
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 &amp; 1 &amp; 0 &amp; 1 \\
 1 &amp; -2 &amp; 1 &amp; 0 \\
 0 &amp; 1 &amp; -2 &amp; 1 \\
 1 &amp; 0 &amp; 1 &amp; -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 &amp; 0 &amp; 0 &amp; 0 \\
 1 &amp; -2 &amp; 1 &amp; 0 \\
 0 &amp; 1 &amp; -2 &amp; 1 \\
 0 &amp; 0 &amp; 0 &amp; 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}
    &amp;&amp; 1 &amp;-&amp; 2 &amp;-&amp; 3 &amp;&amp;  \\
    &amp;&amp; | &amp; &amp; | &amp; &amp; | \\
 4 &amp;-&amp; 5 &amp;-&amp; 6 &amp;-&amp; 7 &amp;-&amp; 8 \\
 | &amp; &amp; | &amp; &amp; | &amp; &amp; | &amp; &amp; | \\
 9 &amp;-&amp; 10 &amp;-&amp; 11 &amp;-&amp; 12 &amp;-&amp; 13 \\
    &amp;&amp; |  &amp; &amp; |  &amp;&amp;  | &amp; &amp;  \\
    &amp;&amp; 14 &amp;-&amp; 15 &amp;-&amp; 16 &amp;&amp; \\
\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 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 \\
 1 &amp; -2 &amp; 1 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 \\
 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 \\
 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 \\
 1 &amp; 0 &amp; 0 &amp; 1 &amp; -4 &amp; 1 &amp; 0 &amp; 0 &amp; 0 &amp; 1 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 \\
 0 &amp; 1 &amp; 0 &amp; 0 &amp; 1 &amp; -4 &amp; 1 &amp; 0 &amp; 0 &amp; 0 &amp; 1 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 \\
 0 &amp; 0 &amp; 1 &amp; 0 &amp; 0 &amp; 1 &amp; -4 &amp; 1 &amp; 0 &amp; 0 &amp; 0 &amp; 1 &amp; 0 &amp; 0 &amp; 0 &amp; 0 \\
 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 \\
 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 \\
 0 &amp; 0 &amp; 0 &amp; 0 &amp; 1 &amp; 0 &amp; 0 &amp; 0 &amp; 1 &amp; -4 &amp; 1 &amp; 0 &amp; 0 &amp; 1 &amp; 0 &amp; 0 \\
 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 1 &amp; 0 &amp; 0 &amp; 0 &amp; 1 &amp; -4 &amp; 1 &amp; 0 &amp; 0 &amp; 1 &amp; 0 \\
 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 1 &amp; 0 &amp; 0 &amp; 0 &amp; 1 &amp; -4 &amp; 1 &amp; 0 &amp; 0 &amp; 1 \\
 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 \\
 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 \\
 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 1 &amp; -2 &amp; 1 \\
 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 0 &amp; 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.