How to "fill in" 3-dimensional Laplacian kernels

Question

The other day, I asked this question 3x3x3 Laplace Kernel?, regarding what the 3x3x3 kernel was for applying a Laplacian convolution.

On that page, it mentions the kernels were "deduced by using discrete differential quotients." Does anyone know how this is done? For example, what if I need a 5x5x5 kernel? How do I "solve" for it? I tried looking in Matlab but didn't find anything.

The reason I want to do this is because if I have a 100x100x10 pixel image I want to apply a Laplacian convolution to, I can convolve with my 3x3x3 kernel, OR I can apply a Fourier transform to both, and then use point-wise multiplication to get the same solution, then use the inverse Fourier to get my result back, which should match the result I got using the direct convolution method. As I mentioned, to do the Fourier approach, I need two images - my original, and the "Laplacian" one. I figured this Laplacian one needs to be generated in the same way as the 3x3x3 one my previous question identified, but needs to also be the same dimensions as my real image (100x100x10 in this example) in order for the Fourier transform approach to work. So my question above is - how do I 'fill in' the 3-dimensional matrix that is my discrete Laplacian operator?

Sorry if this question doesn't make sense. In a very dumbed down sense, I just want to know how to solve for the 3x3x3 kernel myself, so that I can solve for 'larger kernels' like a 100x100x10 one.

Answer 1

If you want the "5x5x5" kernel, then it is no longer the Laplacian in the usual sense.

Let me quickly describe where that kernel you saw on the Wikipedia page comes from. The Laplacian, as a differential operator, is $\sum_i (\partial_i)^2$. Now if we discretize the space into a grid, we can approximate the partial derivative in the $x$ direction by the finite difference of the function at a point $f(x_0,y_0,z_0)$ and $f(x_0 + 1, y_0,z_0)$. Here we use the usual formula $$ \frac{\partial}{\partial x}f(x_0,y_0,z_0) = \lim_{h\to 0} \frac{f(x_0 + h, y_0, z_0) - f(x_0,y_0,z_0)}{h} $$ And since that the space has discrete steps, the smallest non-vanishing interval you can take is of size step 1. So morally speaking $f(x_0 + 1 , y_0,z_0) - f(x_0,y_0,z_0)$ is the derivative in the $+x$ direction of the function $f$ defined on the grid.

The second derivative can be expressed as $$ \frac{\partial^2}{\partial x^2} f(x_0,y_0,z_0) = \lim_{h\to 0} \frac{\partial_xf(x_0 + h, y_0, z_0) - \partial_xf(x_0,y_0,z_0)}{h}$$ $$ = \lim_{h\to 0} \frac{ f(x_0 + h, y_0, z_0) - f(x_0,y_0,z_0) - f(x_0,y_0,z_0) + f(x_0 - h, y_0, z_0)}{h^2} $$ So the discrete version is simply the difference between the discrete derivative in the $+x$ direction and that in the $-x$ direction. (Notice that since we discretize, the forward and backward derivatives are not necessary the same.) In other words, the discrete second derivative is $$ f(x_0 + 1, y_0, z_0) + f(x_0-1,y_0,z_0) - 2 f(x_0,y_0,z_0) $$

Now you just need to add it over the three spatial directions.

So naturally the discrete Laplace kernel has support with in a cube that is 3 on a side. So what you want to do, I think, is probably just trivially extend the filter: make the filter 100 on a side and put the 3x3x3 kernel smack in the center (or as close to it as you can), and pad the rest with zero. Then you can take the Fourier transform and multiply against the Fourier transform of your original image. Note that I suspect whatever software solution you have will take the Fourier transform by assuming the data is periodic, so if you do this the edges will wrap around and bleed into the opposite sides a little. So the two results probably will not match on a two pixels wide band around the boundary. But they should match in the interior.

Answer 2

You could just use the 3x3x3 kernel from last time and ignore the extra data, that would be fine, but very wasteful. The payoff for using a big kernel should be a high degree of approximation. You may observe that a degree $n$ scheme (error of order $n+1$), if one exists, will give an exact answer on all polynomials of degree $\leq n$. This gives us a linear system to solve. Proceeding naively, we see that the number of polynomials of degree $\leq i$ in $3$ dimensions is given by $3 + i \choose i$. In particular, lets try for an order 7 scheme, in which case there are ${3+7 \choose 7} = 120$ numbers to match exactly and five free parameters. Assuming that there is no degeneracy, there ought to be a five-dimensional space of choices, now just pick one. For a thorough discussion, I recommend Strikwerda's book on finite difference schemes (http://www.amazon.com/Difference-Schemes-Partial-Differential-Equations/dp/0898715679), where you can find a lot of heuristics for generating these kernels as well as some rigorous analysis.

Answer 3

Not certain if this will help, but Mathematica has a LaplacianFilter[] function with adjustable dimensions in its ImageProcessing package:

LaplacianFilter[n1,n2,n3,...] // MatrixForm

There is a description here. The command above returns the $n_1 \times n_2 \times n_3 \times \cdots$ filter matrix. You might have to attend to which version of Mathematica you use, as the latest seems to have optimized the filter to directly apply to an image, rather than return the matrix explicitly.

Matlab has a similar function, described here:

fspecial('laplacian')

Not sure you can adjust the dimensions for fspecial() as in Mathematica.