function imcanny(Img,sigma,th,tl);
% SINTASSI: imcanny(I,sigma,th,tl)--> Estrae i bordi di un'immagine in ingresso "I" 
% mediante la regola di Canny, il parametro "sigma" fornisce la deviazione standard del 
% filtro gaussiano, "th" e "tl" sono i parametri per effettuare la sogliatura con isteresi

if th<tl 
    error('Th deve essere maggiore di Tl');
end    	

I = double(Img);               	% trasformo l'immagine in double
immostra(I,gray(256));          	% Immagine 1: Immagine originale
title('Immagine originale');

x	= [-(5*sigma):(5*sigma)];    	% Il vettore contiene circa il 97% dell'area sottesa dalla Gaussiana
G	= normpdf(x,0,sigma);        	% Filtro Gaussiano monodimensionale
dG	= diff(G);                   		% Derivata del Gaussiano monodimensionale

figure;
plot(G);				% Immagine 2: Filtro Gaussiano
title('Filtro Gaussiano');
figure;
plot(dG);				% Immagine 3: Filtro derivata del Gaussiano
title('Derivata del Gaussiano');

H = conv(G,dG);
Ix = conv2(I,H,'same');           	% Convoluzione dell'immagine con H lungo le righe
Iy = conv2(I',H,'same')';   	% Convoluzione dell'immagine con H lungo le colonne

immostra(Ix+128,gray(256));         % Immagine 4: Gradiente lungo X
title('Gradiente lungo X');
immostra(Iy+128,gray(256));         % Immagine 5: Gradiente lungo Y
title('Gradiente lungo Y');

Imod =sqrt(Ix.^2 + Iy.^2);      
Imod=(Imod./max(max(Imod)))*255;	% Normalizzazione, stretching  del modulo del gradiente
immostra(Imod,gray(256));      	% Immagine 6: Modulo del gradiente
title('Modulo del Gradiente');

% Ricerca dei massimi locali
%
%       1         2
%   X--------X--------X
%   |                 |
% 4 |    -Iy ^        | 3
%   |        |        |
%   |        |   Ix   |
%   X        P---->   X  
%   |                 |
% 3 |                 | 4
%   |                 |
%   |                 |
%   X--------X--------X
%       2        1  




Imod2=Imod>=tl; 
[r,c]=size(Imod);
for i=2:r-1
    for j=2:c-1
        if Imod2(i,j)
           switch (direzione(Ix(i,j),Iy(i,j)))
            case 1
                d=abs(Ix(i,j)/Iy(i,j));
                Imod2(i,j)=Imod(i,j)>(1-d)*Imod(i+1,j)+d*Imod(i+1,j+1) & Imod(i,j)>(1-d)*Imod(i-1,j)+d*Imod(i-1,j-1);
            case 2
                d=abs(Ix(i,j)/Iy(i,j));
                Imod2(i,j)=Imod(i,j)>(1-d)*Imod(i-1,j)+d*Imod(i-1,j+1) & Imod(i,j)>(1-d)*Imod(i+1,j)+d*Imod(i+1,j-1);
            case 3
                d=abs(Iy(i,j)/Ix(i,j));
                Imod2(i,j)=Imod(i,j)>(1-d)*Imod(i,j+1)+d*Imod(i-1,j+1) & Imod(i,j)>(1-d)*Imod(i,j-1)+d*Imod(i+1,j-1);
            case 4
                d=abs(Iy(i,j)/Ix(i,j));
                Imod2(i,j)=Imod(i,j)>(1-d)*Imod(i,j+1)+d*Imod(i+1,j+1) & Imod(i,j)>(1-d)*Imod(i,j-1)+d*Imod(i-1,j-1);
            end
        end
   end 
end

Imod=Imod.*Imod2;
immostra(Imod,gray(256));   % Immagine 7: Modulo del gradiente dopo nonmaximum suppression
title('Modulo del Gradiente dopo N.M.S.');

%---      Sogliatura con isteresi      ---%
Imodsog=Imod<th;
for i=2:r-1
    for j=2:c-1
        if not(Imodsog(i,j))
            for h=i-1:3
                for k=j-1:3
    Imodsog(h,k)=Imod(h,k)<tl;
                end
            end                    
        end        
    end
end
immostra(Imodsog*255,gray(256));    % Immagine 8: Bordi dopo la sogliatura (binaria)
title('Bordo');

%---     Sottofunzione usata per calcolare la direzione del gradiente     ---%
function dir = direzione(gx,gy);
if abs(gy)>abs(gx)
    if(gx*gy)>0 
        dir=1;
    else
        dir=2;          % Contiene caso gx=0
    end
else
    if(gx*gy)>0
        dir=4;          % Contiene caso gx=gy >0
    else
        dir=3;          % Contiene caso gx=0 e gx=gy<0
    end
end
