How do I implement the Structure tensor (in Matlab)?
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)