lin_sym.m: Computes the linear symmetry, written by Stéphane Michel

function [orientation]=lin_sym(pict,sigma_d,sigma_i,filter_size,aff);
%function [orientation]=lin_sym(pict,sigma_d,sigma_i,filter_size,aff)
% OUTPUT: Complex image containing x et y component of linear symmetry
%         if aff=1 display different images else no display else put 0
% Summary of processing:
%	
%-subsampling with a factor 2
%-computation of a gaussian filter and its derivative of size filter_size
%-real and imaginary part of linear symmetry computation
%

%
% ADDITION OF BORDERS
%
[m n]=size(pict);
border=pict(2,2)*ones(m+40, n+40);
border(21:20+m, 21:20+n) = pict;
picture=border;
clear border

if aff==1
 gcf=figure;
 set(gcf,'Name','Image initiale')
 set(gcf,'Position',[485 850 taille(1,2) taille(1,1)]);
 image(pict)
 colormap(gray(256))
 truesize;
 axis off 
end

[m n] = size( picture );

%
% COMPUTATION OF GAUSSIAN AND DERIVATIVE OF GAUSSIAN
% COEFFICIENTS
%
demi_size=(filter_size-1)/2;

x=-demi_size:1:demi_size;

% les constantes 1/3 et 2 permettent de normaliser
% la reponse en frequences des filtres

filt_gauss=(exp((-x.^2)/(2*sigma_d.^2))/3);
filt_derivgauss=(x.*filt_gauss)*2;

if aff==1
 gcf=figure;
 set(gcf,'Name','Filtres')
 subplot(2,1,1), plot(filt_gauss);
 subplot(2,1,2), plot(filt_derivgauss);
end

dx=conv2(picture,filt_gauss','same');
dx=conv2(dx,filt_derivgauss,'same');

dy=conv2(picture,filt_derivgauss','same');
dy=conv2(dy,filt_gauss,'same');

if aff==1
 gcf= figure;
 axes('Position',[0 0 1 1]);
 set(gcf,'Name','Derivee selon x');
 set(gcf,'Position',[ 20 485 n m]);
 imagesc(dx)
 colormap(gray(64))
 axis off;
 gcf=figure;
 axes('Position',[0 0 1 1]);
 set(gcf,'Name','Derivee selon y');
 set(gcf,'Position',[20 80 n m]);

 imagesc(dy);
 colormap(gray(64))
 axis off;
end

img=2*(dx.*dy);
rl=dy.^2 - dx.^2;

filtrage=1.333*fspecial('gaussian',10,sigma_i);
rl=filter2(filtrage,rl);
img=filter2(filtrage,img);

if aff==1
 gcf=figure;
 axes('Position',[0 0 1 1]);
 set(gcf,'Name','Composante X');
 set(gcf,'Position',[450 485 n m]);
 imagesc(rl);
 colormap(gray(64))
 axis off;
 gcf=figure;
 axes('Position',[0 0 1 1]);
 set(gcf,'Name','Composante Y');
 set(gcf,'Position',[450 80 n m]);
 imagesc(img)
 colormap(gray(64))
 axis off;
end

orientation=rl+i.*img;
maxi=abs(orientation);

maxi=max(maxi(:));
orientation=orientation/(maxi);
orientation=orientation( 21:m-20, 21:n-20);

if aff==1
 gcf=figure;
 axes('Position',[0 0 1 1]);
 set(gcf,'Name','Norme');
 set(gcf,'Position',[850 485 n m]);
 temp=abs(orientation);
 imagesc(temp);
 colormap(gray(64))
 axis off;
 gcf=figure;
 axes('Position',[0 0 1 1]);
 set(gcf,'Name','Phase');
 set(gcf,'Position',[850 80 n m]);
 temp=angle(orientation); 
 imagesc(temp)
 axis off;
end
clear rl img dx dy