
function Icont = contornos(I,umbral)
% CONTORNOS Detecta contornos en imag. de niv. de gris empleando gradiente   
% Obtencion de los contornos usando un filtro de Sobel                 %
% Se buscan los máximos en la dirección del gradiente                  %                                

%______________________________________________________________________%
%                                                                      %
% NOMBRE DEL FICHERO:   Icont = contornos (Ient,umbral)                %
%                                                                      %
% RECIBE: Imagen de niveles de gris y un umbral                        %
% DEVUELVE:  La misma imagen con los pixels de contorno en blanco      %
% ALGORITMO EMPLEADO: Tecnica basada en el gradiente                   %
% VARIABLES:   Ient = imagen introducida                               %
%              Icont = imagen de contornos resultante.Binaria          %
%                      Fondo negro y contornos en blanco               %
%              umbral entre 0 y 1. Cuanto más alto menos contornos     %
%              Ix = Imagen de la derivada respecto a x                 %
%              Iy = Imagen de la derivada respecto a y                 %
%              M = Imagen Modulo del gradiente                         %
%              A = Imagen Argumento del gradiente                      %
% LIBRO: VISION ARTIFICIAL INDUSTRIAL. UNIVERSIDAD DE VALLADOLID 2012  %   
% TEMA 7 LA IMAGEN DIGITAL pág.142 Algoritmo 7.1                       %
% AUTORES: Eusebio de la Fuente y Félix Miguel Trespaderne             %
%______________________________________________________________________%
%  Ejemplo de ejecucion:                                               %
% >>I=imread('cameraman.tif');                                         %
% >> figure,imshow(I)                                                  %
% >> Icont = contornos(I,0.3);
% >> figure,imshow(Icont); 
%______________________________________________________________________%
% Con funciones de matlab
% >>I=imread('cameraman.tif');
% >>C=edge(I,'sobel');
% >>imshow(C);
%______________________________________________________________________%


[filas,columnas] = size(I);
I = double(I);
M = zeros(filas,columnas);          %Imagen Modulo del gradiente
A = zeros(filas,columnas);          %Imagen Argumento del gradiente
Icont = zeros(filas,columnas);      %Imagen con el contorno resultante

%% Filtrado Sobel para obtener Ix e Iy
Sx=[-1, 0, 1 ;            %otra forma: S=fspecial('sobel'); Sx=S'; Sy=S;
    -2, 0, 2;
    -1, 0, 1];
Sy=	[-1,-2,-1;  
      0, 0, 0;  
      1, 2, 1];
Ix=imfilter(I,Sx);
Iy=imfilter(I,Sy);

%% Modulo y argumento grad 
for i=2:filas-1
    for j=2:columnas-1
        M(i,j) = abs(Ix(i,j))+abs(Iy(i,j));
        if Ix(i,j)~=0        %si grad no es vertical
            A(i,j) = atan(Iy(i,j)/Ix(i,j));
        else
            A(i,j)=-1.5708;  % 90º, gradiente vertical 
        end
    end
end

M=mat2gray(M);        %reescala el modulo entre 0 y 1

%% Eliminacion puntos con gradiente bajo
for i=2:filas-1
    for j=2:columnas-1
        if M(i,j)<umbral
            M(i,j)=0.0;
        end
    end
end


%% Cuantifica direc. gradiente a cuatro direcciones
for i=2:filas-1
    for j=2:columnas-1
        if  M(i,j)==0 
            A(i,j)=0;  %si el gradiente es cero lo marcamos a 0 como fondo
        elseif A(i,j)>=-1.5708 && A(i,j)<-1.1781 	%entre -90º y -67,5º
            A(i,j)=1;
            
        elseif A(i,j)>=-1.1781 && A(i,j)<-0.3927 %entre -67,5º y -22,5º
            A(i,j)=2;
            
        elseif A(i,j)>=-0.3927 && A(i,j)<0.3927 	%entre -22,5º y 22,5º
            A(i,j)=3;
            
        elseif A(i,j)>=0.3927 && A(i,j)<1.1781	%entre 22,5º y 67,5º
            A(i,j)=4;
            
        elseif A(i,j)>=1.1781 && A(i,j)<1.5708	%entre 67,5º y 90º
            A(i,j)=1;
        end
    end
end


%% Maximos en direc. gradiente
for i=2:filas-1
    for j=2:columnas-1
        switch A(i,j)
            case 0
                %Es el fondo. No comparamos nada
            case 1
                if M(i-1,j)<=M(i,j) && M(i,j)>M(i+1,j)
                    Icont(i,j)=1;
                end
            case 2
                if M(i+1,j-1)<=M(i,j) && M(i,j)>M(i-1,j+1)
                    Icont(i,j)=1;
                end
            case 3
                if M(i,j-1)<=M(i,j) && M(i,j)>M(i,j+1)
                    Icont(i,j)=1;
                end
            case 4
                if M(i-1,j-1)<=M(i,j) && M(i,j)>M(i+1,j+1)
                    Icont(i,j)=1;
                end
        end
    end
end


