image processing - How do I implement the Structure tensor (in Matlab)?

I am trying to figure out the details on how to implement the structure tensor in Matlab and need some advice!

For an image $\ I(x,y) $ the structure tensor S is given by:

$$ S=\begin{pmatrix} W \ast I_x^2 & W \ast (I_xI_y)\\ W \ast (I_xI_y) & W \ast I_y^2 \end{pmatrix}$$ where W is a smoothing kernel, $\ \ast $ denotes convolution and subscript denotes partial derivative with respect to.

I want to use a Gaussian as smoothing kernel. This is separable: $$ W = G_{\sigma}*(G_{\sigma})^T$$ The constant $\ \sigma $ controls the outer/integration scale. E.g. if I truncate it at $\ +/- 2\sigma$ and use $\ \sigma = 1$ I get:

$$ G_{\sigma} = [0.0103,0.2076,0.5642,0.2076,0.0103] $$

However the input to the Gaussian might have a mean different from 0.

How do I compensate for this?

I was thinking substracting the mean first.

I could use central differences to estimate the partial derivatives: $$ I_x \approx [-0.5 \hspace{2 pt} 0 \hspace{2 pt} 0.5]*I $$ However since the image is noisy it is probably better to use a more robust differentian operator. The Sobel operator is such an operator. However according to wikipedia it is not perfectly rotational symmetric. The Scharr operators are supposed to be perfectly rotational symmetric and robust approximations of the partial derivatives:

$$ I_x \approx \frac{1}{32}\begin{pmatrix} -3 & 0 & 3\\ -10 & 0 & 10\\ -3 & 0 & 3 \end{pmatrix} * I$$

$$ I_y \approx \frac{1}{32}\begin{pmatrix} -3 & -10 & -3\\ 0 & 0 & 0\\ 3 & 10 & 3 \end{pmatrix} * I$$

The Scharr operators are separable: $$ \begin{pmatrix} -3 & 0 & 3\\ -10 & 0 & 10\\ -3 & 0 & 3 \end{pmatrix} = [-1,0,1]*[3,10,3]^T $$

$$ \begin{pmatrix} -3 & -10 & -3\\ 0 & 0 & 0\\ 3 & 10 & 3 \end{pmatrix} = [3,10,3]*[-1,0,1]^T$$

The elements of the structure tensor now becomes series of consecutive horizontal and vertical convolutions with the original image.

EDIT: Another option (thanks to nikie!) is the Gaussian derivative. Introducing $$ GD_{\sigma}=\frac{-x}{\sigma}G{\sigma}$$

Is this then correct? $$ Ix = G{\sigma}^T*(GD_{\sigma}*I) $$ $$ Iy = G{\sigma}*(GD_{\sigma}^T*I) $$

How do I Implement this easily in Matlab?

How about in C++?

Any comments regarding choice of smoothing kernel and differentiation operator?

EDIT 2: My Matlab code so far:

% I: image
% si: inner (differetiation) scale,
% so: outer (integration) scale

function [Sxx, Sxy, Syy] = structureTensor(I,si,so)
[m n] = size(I);

Sxx = NaN(m,n);
Sxy = NaN(m,n);
Syy = NaN(m,n);

% Robust differentiation by convolution with derivative of Gaussian:
x  = -2*si:2*si;
g  = exp(-0.5*(x/si).^2);

g  = g/sum(g);
gd = -x.*g/si; % is this normalized?

Ix = conv2( g',conv2(gd,I) );
Iy = conv2( g,conv2(gd',I) );

Ixx = Ix.^2;
Ixy = Ix.*Iy;
Iyy = Iy.^2;



% Smoothing:
x  = -2*so:2*so;
g  = exp(-0.5*(x/so).^2);
Sxx = conv2( g',conv2(g,Ixx) ); 
Sxy = conv2( g',conv2(g,Ixy) );
Syy = conv2( g',conv2(g,Iyy) );

How do I normalize the Gaussian derivative?

Answer

However the input to the Gaussian might have a mean different from 0.

Unless your image is completely constant, $I_x^2$ and $I_y^2$ are guaranteed to have a mean > 0 (assuming the input pixels are all reals, of course).

How do I compensate for this?

Why would you want to do that?

I was thinking substracting the mean first.

Subtracting the mean before applying a derivative operator like Sobel or Scharr? That shouldn't make any difference.

The Scharr operators are supposed to be perfectly rotational symmetric

Perfect rotational symmetry? Is that even possible with a finite kernel? I think the Scharr operator is just optimal (in a least-squares sense) for that filter size. With larger filter sizes, you get better rotational symmetry.

How do I Implement this easily in Matlab?

You can use imfilter to apply the Sobel or Scharr filters, then .* to multiply/square the results pixelwise, then imfilter again to apply the gaussian smoothing kernel.

How about in C++?

Any decent image processing library should contain functions for linear filters and pixel-wise multiplication. e.g. for the IPP, the functions would be called ippiFilterRow/Column..., ippiMul... and ippiSqr...

Any comments regarding choice of smoothing kernel and differentiation operator?

I think the gaussian smoothing kernel is pretty standard. Sobel and Scharr are both common choices for the differentiation operator. Another good choice would be a gaussian derivative filter. Which one's best depends on your application: how noisy are your images? How sharp are the contrasts you're looking for? (If there's lots of noise, and the contrasts are sharp, a derivative filter with a larger aperture like the gaussian derivative might be better). How important is rotational symmetry?

How do I normalize the Gaussian derivative?

The gaussian is: $\frac{1}{2 \pi \sigma ^2}e^{-\frac{x^2+y^2}{2 \sigma ^2}}$

Deriving by x yields:

$\partial _x\left(\frac{1}{2 \pi \sigma ^2}e^{-\frac{x^2+y^2}{2 \sigma ^2}}\right) = -\frac{e^{-\frac{x^2+y^2}{2 \sigma ^2}} x}{2 \pi \sigma ^4}$

So I think you should divide by gd = -x.*g/(si^2); because sum(g)=2*pi*(si^2)