structure_tensor

% determine structure tensor from image patch
 
 im    = window
 tau   = differentiation scale
 sigma = integration scale

 T     = structure tensor 
       = \overline(g'_x g_x)
       = G(\sigma) * (G(\tau)*\nabla_g'  G(\tau)*\nabla_g )

 T     = structure_tensor(im,tau,sigma)

 Observe: 
 The structure tensor is not the sum ...
                         but the average of the squared gradient

 Wolfgang Förstner 4/2018
 wfoerstn@uni-bonn.de

0001 %% determine structure tensor from image patch
0002 %
0003 % im    = window
0004 % tau   = differentiation scale
0005 % sigma = integration scale
0006 %
0007 % T     = structure tensor
0008 %       = \overline(g'_x g_x)
0009 %       = G(\sigma) * (G(\tau)*\nabla_g'  G(\tau)*\nabla_g )
0010 %
0011 % T     = structure_tensor(im,tau,sigma)
0012 %
0013 % Observe:
0014 % The structure tensor is not the sum ...
0015 %                         but the average of the squared gradient
0016 %
0017 % Wolfgang Förstner 4/2018
0018 % wfoerstn@uni-bonn.de
0019 
0020 function T  = structure_tensor(im,tau,sigma)
0021 
0022 % size of image, number of channels
0023 [N,M,K] = size(im);
0024 
0025 % initiate T
0026 T = zeros(2);
0027 
0028 % sum up over all channels
0029 for k = 1:K
0030     % gradient at scale tau
0031     im_x = gaussFFT(double(im(:,:,k))/256,tau,'Gx'); 
0032     im_y = gaussFFT(double(im(:,:,k))/256,tau,'Gy'); 
0033     % squared gradient
0034     imxx = im_x.*im_x;
0035     imyy = im_y.*im_y;
0036     imxy = im_x.*im_y;
0037     % mean squared gradient
0038     txx = gaussFFT(imxx,sigma,'G');
0039     tyy = gaussFFT(imyy,sigma,'G');
0040     txy = gaussFFT(imxy,sigma,'G');
0041     % sum up channels
0042     T = T+[txx((N+1)/2,(M+1)/2) , txy((N+1)/2,(M+1)/2); ...
0043            txy((N+1)/2,(M+1)/2) , tyy((N+1)/2,(M+1)/2)];
0044 end