Implement the 3D Structure Tensor in C/C++?

Question

I am trying to figure out the details on how to implement the 3D structure tensor in C/C++ in an easy but efficient way and need some advice!

For a discrete function $ I(x_i,y_j,z_k)$ the 3D structure tensor is given by: $$ S=\begin{pmatrix} W \ast I_x^2 & W \ast (I_xI_y) & W \ast (I_xI_z)\\ W \ast (I_xI_y) & W \ast I_y^2 & W \ast I_yI_z \\ W \ast (I_xI_z) & W \ast I_yI_z & W \ast I_z^2 \\ \end{pmatrix}$$ where W is a smoothing kernel, $\ \ast $ denotes convolution and subscript denotes partial derivative with respect to.

The calculation of the structure tensor have two main steps:

1. calculate the partial derivatives of the function in a way that is robust to noise over a window.
2. smooth products of the partial derivatives over another larger window.

I start by looking at step 2. I want to use a Gaussian as smoothing kernel. The normal distribution in 3 dimensions is separable: $$ g(x,y,z) = g(x)g(y)g(z) $$ where $$ g(x) = \frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{1}{2}(x/\sigma)^2) $$ etc.

The Fourier transform of the normal distribution in 3 dimensions is given by: $$ G(kx,ky,kz) = G(kx)G(ky)G(kz) $$ where $$ G(kx) = \frac{\sigma}{\sqrt{2\pi}}exp(-\frac{1}{2}(kx\sigma)^2) $$ etc.

How do I implement the Gaussian smoothing?

The simplest way to implement the Gaussian smoothing would be to loop over a 3D Gaussian kernel for each point in I.

Since the Gaussian is separable however it should be more efficent to perform a convolution with a 1D Gaussian in the x direction followed by the y direction and z direction.

Another even more efficient approach (?) would be to do the convolution in the Fourier space (where it becomes a multiplication): $$ g*a=\mathcal{F}^{-1}[GA] $$

Next I look at step 1.

The derivative of the Gaussian is given by: $$ G_x(x) = -x/\sigma^2G(x)G(y)G(z)$$ etc., where subscript denotes partial derivative with respect to.

The Gaussian derivative can be used to estimate the partial derivatives of I in a way that is robust to noise: $$ I_x=-x/\sigma^2G(x)*G(y)*G(z)*I $$

Here I am faced with the same decision: implement it without using separability, implement it using separability or implement it in Fourier space?

Any advice or comments?

Answer 1

There are many sides to this.

1. For typical kernel sizes, Fourier transform is not the most efficient method for convolution, if you are correctly utilizing the CPU for sliding window operations. That, of course, takes a good deal of software optimization. There are libraries to do that though (for example many of the deep learning libraries do 3D convolutions in CPU/GPU - I guess most of them do not use FFT). Given the specific symmetric nature of the Gaussian kernel, it is also very suitable for sliding operations. If the kernel size is very large, you could also try recursive convolutions.

2. The Fourier domain operations might be considered to be a better alternative if you have some luxury to pre-compute certain things. For example, if the image size is not changing and etc.

3. Probably the best answer would be to try and actually see it. I would prefer the separable convolutions. Vast amount of resources are also available for this.