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

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?

• It's not clear what your question is. – Jason R Aug 8 '12 at 21:46
• Hi Jason R! I edited my question. Hopefully it is clearer now. – Andy Aug 9 '12 at 9:25
• I think it depends on size of your kernel. For smaller kernels, doing three convolutions is more efficient. For large kernels, FT approach is more efficient. If you don't need accuracy, only smooth the data somehow, it is also possible to perform two or three passes of box filter (also separable) with sliding sum. Such filter approximates Gaussian and have very low complexity. – Libor Aug 9 '12 at 14:11