%
% This source code is freely distributed from the "colormonogenic" website:
% http://xlim-sic.labo.univ-poitiers.fr/projets/colormonogenic/
% published in 2018,
% which presents the main research results by 
% Raphaël Soulard & Philippe Carré,
% from the XLIM Laboratory (UMR CNRS 7252),
% University of Poitiers, France.
%
% Author: R. Soulard.
%

% Extraction of structure tensor features
% from the marginal gradient of a matrix

% (No dependencies)


function [N,lambdaP,coher,orient]=tensor(Ix,Iy,sigm)
% -------------------------------------------------------------- %
% Ix,Iy: marginal gradient
% sigm: width of gaussian filtering, if 0 : no filtering
% N2: squared norm of the gradient
% lambdaP: higher eigen value of tensor
% coher: coherency in [0;1]
% orient: orientation of eigenvector for higher eigenvalue
% -------------------------------------------------------------- %
% tensor components ( T=[a,b;b,c] ):
a = sum(Ix.^2  , 3);
b = sum(Ix.*Iy , 3);
c = sum(Iy.^2  , 3);

% Smoothing of tensor components:
if sigm>0 
  % Create space domain Gaussian mask:
  [M,N]=size(a);
  xx = [0:ceil(M/2)-1 , -floor(M/2):-1]' * ones(1,N);
  yy = ones(M,1) * [0:ceil(N/2)-1 , -floor(N/2):-1];
  N2 = xx.^2+yy.^2;
  h = (1/(2*pi*sigm*sigm)) * exp( -N2/(2*sigm*sigm) );
  % Filtering of a, b and c:
  a = real(ifft2(fft2( a ).*fft2(h))); 
  b = real(ifft2(fft2( b ).*fft2(h)));
  c = real(ifft2(fft2( c ).*fft2(h))); 
end

N2 = a+c;

idx = N2<10^(-16); % get singular cases where Gradient is null
a(idx) = 0; % there, set the gradient to actual zero to avoid unstabilities
b(idx) = 0; % there, set the gradient to actual zero to avoid unstabilities
c(idx) = 0; % there, set the gradient to actual zero to avoid unstabilities
N2(idx) = 0;% there, set the gradient to actual zero to avoid unstabilities

N = sqrt(N2);
tmp =  (a-c) + 1i*(2*b) ;
r        = abs(   tmp );
twotheta = angle( tmp );
orient = twotheta/2; % in [-pi/2;pi/2]
lambdaP = 0.5*(N2+r);
%lambdaM = 0.5*(N2-r);
coher = r; % init
coher(~idx) = r(~idx)./N2(~idx); % coherency = r/N2 = (lambdaP-lambdaM)/(lambdaP+lambdaM)
