close all
clear all

% Cargo la Imagen Inicial
Imagen_Inicial=imread('01_original.bmp'); 
figure(1)
imshow(Imagen_Inicial)

% Obtenemos la componente verde
Verde=double(Imagen_Inicial(:,:,2))/255;
figure(2)
imshow(Verde)

%Obtenemos el filtro de media aritmetica grande y se lo aplicamos a la
%componente verde
Filtro_Media=fspecial('average',[60,60]);
Imagen_Filtrada=imfilter(Verde,Filtro_Media);
figure(3)
imshow(Imagen_Filtrada)

Imagen_Normalizada=Imagen_Filtrada-Verde;
figure(4)
imshow(Imagen_Normalizada)

[Filas,Columnas]=size(Imagen_Normalizada);

%Eliminamos el borde de la imagen a traves de la mascara
Mascara=imread('01_Mask.gif');

for i=1:Filas
    for j=1:Columnas
        if Mascara(i,j)==0
            Imagen_Normalizada(i,j)=0;
        end
    end
end

figure(5)
imshow(Imagen_Normalizada)
            
%Thin Vessel enhancement
H1=1/6*[-1 -1 -1; 2 2 2; -1 -1 -1];
H2=1/6*[2 -1 -1; -1 2 -1; -1 -1 2];
H3=1/6*[-1 2 -1; -1 2 -1; -1 2 -1];
H4=1/6*[-1 -1 2; -1 2 -1; 2 -1 -1];

Orientacion_0=imfilter(Imagen_Normalizada,H1);
Orientacion_45=imfilter(Imagen_Normalizada,H2);
Orientacion_90=imfilter(Imagen_Normalizada,H3);
Orientacion_135=imfilter(Imagen_Normalizada,H4);

for i=1:Filas
    for j=1:Columnas
        m=1;
        X(m)=Orientacion_0(i,j);
        m=m+1;
        X(m)=Orientacion_45(i,j);
        m=m+1;
        X(m)=Orientacion_90(i,j);
        m=m+1;
        X(m)=Orientacion_135(i,j);
        Realzado(i,j)=max(X);
    end
end

Imagen_Realzada=Imagen_Normalizada+Realzado;

So8=strel('disk',8);
Sc=strel('disk',1);
TopHat8=Imagen_Realzada-min(imopen(imclose(Imagen_Realzada,Sc),So8),Imagen_Realzada);

figure(6)
imshow(Imagen_Realzada)
        
% Obtenemos gradientes en direccion x e y 
[GradienteX,GradienteY]=gradient(TopHat8);

% Obtenemos la magnitud del gradiente 
Mag_Gradiente=sqrt(GradienteX.^2+GradienteY.^2);

% Obtenemos la orientacion del gradiente forzando que atan(0/0)=0 y asi no
% de error
for f=1:Filas
    for c=1:Columnas
        if GradienteX(f,c)==0 && GradienteY(f,c)==0
            Orient_Gradiente(f,c)=0;
        else
            Orient_Gradiente(f,c)=atan(GradienteY(f,c)/GradienteX(f,c));
        end
    end
end

% N*N tamaño de la vecindad alrededor de (i,j)
N=31;

% Obtenemos el valor M que se utilizara posteriormente para hacer los
% sumatorios y delimitarlos
M=(N-1)/2;

% Pasamos a Calcular el angulo de anisotropia, utilizando variables
% desconocidas hasta ahora que se explicaran a medida que vayan apareciendo
Angulo_Anisotropia=[]; 

% Analizamos primero las M primeras filas
Contador_F=0;
for i=1:M
    Contador_C=0;
    for j=1:M
        resultado_num=0;
        resultado_den=0;
        for m=(M-Contador_F):(N-1)   
            for n=(M-Contador_C):(N-1)   
                resultado_provisional_num=(Mag_Gradiente(i-M+m,j-M+n)^2)*(sin(2*Orient_Gradiente(i-M+m,j-M+n)));
                resultado_provisional_den=(Mag_Gradiente(i-M+m,j-M+n)^2)*(cos(2*Orient_Gradiente(i-M+m,j-M+n)));
                
                resultado_num=resultado_num+resultado_provisional_num;
                resultado_den=resultado_den+resultado_provisional_den;
            end
        end
        if resultado_den==0 && resultado_num==0
            Angulo_Anisotropia(i,j)=0;
        else
            Angulo_Anisotropia(i,j)=0.5*atan(resultado_num/resultado_den);
        end
        Contador_C=Contador_C+1;
    end
    Contador_F=Contador_F+1;
end

Contador_F=0;
for i=1:M
    for j=(M+1):(Columnas-M)
        resultado_num=0;
        resultado_den=0;
        for m=(M-Contador_F):(N-1)   
            for n=0:(N-1)   
                resultado_provisional_num=(Mag_Gradiente(i-M+m,j-M+n)^2)*(sin(2*Orient_Gradiente(i-M+m,j-M+n)));
                resultado_provisional_den=(Mag_Gradiente(i-M+m,j-M+n)^2)*(cos(2*Orient_Gradiente(i-M+m,j-M+n)));
                
                resultado_num=resultado_num+resultado_provisional_num;
                resultado_den=resultado_den+resultado_provisional_den;
            end
        end
        if resultado_den==0 && resultado_num==0
            Angulo_Anisotropia(i,j)=0;
        else
            Angulo_Anisotropia(i,j)=0.5*atan(resultado_num/resultado_den);
        end
    end
    Contador_F=Contador_F+1;
end

Contador_F=0;
for i=1:M
    Contador_C=0;
    for j=(Columnas-M+1):Columnas
        resultado_num=0;
        resultado_den=0;
        for m=(M-Contador_F):(N-1)   
            for n=0:(N-2-Contador_C)   
                resultado_provisional_num=(Mag_Gradiente(i-M+m,j-M+n)^2)*(sin(2*Orient_Gradiente(i-M+m,j-M+n)));
                resultado_provisional_den=(Mag_Gradiente(i-M+m,j-M+n)^2)*(cos(2*Orient_Gradiente(i-M+m,j-M+n)));
                
                resultado_num=resultado_num+resultado_provisional_num;
                resultado_den=resultado_den+resultado_provisional_den;
            end
        end
        if resultado_den==0 && resultado_num==0
            Angulo_Anisotropia(i,j)=0;
        else
            Angulo_Anisotropia(i,j)=0.5*atan(resultado_num/resultado_den);
        end
        Contador_C=Contador_C+1;
    end
    Contador_F=Contador_F+1;
end

%Analizamos desde la fila M+1 hasta la Fila-M
for i=(M+1):(Filas-M)
    Contador_C=0;
    for j=1:M
        resultado_num=0;
        resultado_den=0;
        for m=0:(N-1)   % Recorremos todas las filas
            for n=(M-Contador_C):(N-1)   % Recorremos todas las columnas dentro de cada fila
                resultado_provisional_num=(Mag_Gradiente(i-M+m,j-M+n)^2)*(sin(2*Orient_Gradiente(i-M+m,j-M+n)));
                resultado_provisional_den=(Mag_Gradiente(i-M+m,j-M+n)^2)*(cos(2*Orient_Gradiente(i-M+m,j-M+n)));
                
                resultado_num=resultado_num+resultado_provisional_num;
                resultado_den=resultado_den+resultado_provisional_den;
            end
        end
        if resultado_den==0 && resultado_num==0
            Angulo_Anisotropia(i,j)=0;
        else
            Angulo_Anisotropia(i,j)=0.5*atan(resultado_num/resultado_den);
        end
        Contador_C=Contador_C+1;
    end
    
    for j=(M+1):(Columnas-M)
        
        resultado_num=0;
        resultado_den=0;
        
        for m=0:(N-1)   % Recorremos todas las filas
            for n=0:(N-1)   % Recorremos todas las columnas dentro de cada fila
                resultado_provisional_num=(Mag_Gradiente(i-M+m,j-M+n)^2)*(sin(2*Orient_Gradiente(i-M+m,j-M+n)));
                resultado_provisional_den=(Mag_Gradiente(i-M+m,j-M+n)^2)*(cos(2*Orient_Gradiente(i-M+m,j-M+n)));
                
                resultado_num=resultado_num+resultado_provisional_num;
                resultado_den=resultado_den+resultado_provisional_den;
            end
        end
       
        if resultado_den==0 && resultado_num==0
            Angulo_Anisotropia(i,j)=0;
        else
            Angulo_Anisotropia(i,j)=0.5*atan(resultado_num/resultado_den);
        end
    end
    
    Contador_C=0;
     
    for j=(Columnas-M+1):Columnas
        resultado_num=0;
        resultado_den=0;
        for m=0:(N-1)   
            for n=0:(N-2-Contador_C)   
                resultado_provisional_num=(Mag_Gradiente(i-M+m,j-M+n)^2)*(sin(2*Orient_Gradiente(i-M+m,j-M+n)));
                resultado_provisional_den=(Mag_Gradiente(i-M+m,j-M+n)^2)*(cos(2*Orient_Gradiente(i-M+m,j-M+n)));
                
                resultado_num=resultado_num+resultado_provisional_num;
                resultado_den=resultado_den+resultado_provisional_den;
            end
        end
        if resultado_den==0 && resultado_num==0
            Angulo_Anisotropia(i,j)=0;
        else
            Angulo_Anisotropia(i,j)=0.5*atan(resultado_num/resultado_den);
        end
        Contador_C=Contador_C+1;
    end
end

% Analizamos las ultimas M Filas

Contador_F=0;
for i=(Filas-M+1):Filas
    Contador_C=0;
    for j=1:M
        resultado_num=0;
        resultado_den=0;
        for m=0:(N-2-Contador_F)   % Recorremos todas las filas
            for n=(M-Contador_C):(N-1)   % Recorremos todas las columnas dentro de cada fila
                resultado_provisional_num=(Mag_Gradiente(i-M+m,j-M+n)^2)*(sin(2*Orient_Gradiente(i-M+m,j-M+n)));
                resultado_provisional_den=(Mag_Gradiente(i-M+m,j-M+n)^2)*(cos(2*Orient_Gradiente(i-M+m,j-M+n)));
                
                resultado_num=resultado_num+resultado_provisional_num;
                resultado_den=resultado_den+resultado_provisional_den;
            end
        end
        if resultado_den==0 && resultado_num==0
            Angulo_Anisotropia(i,j)=0;
        else
            Angulo_Anisotropia(i,j)=0.5*atan(resultado_num/resultado_den);
        end
        Contador_C=Contador_C+1;
    end
    Contador_F=Contador_F+1;
end

Contador_F=0;
for i=(Filas-M+1):Filas
    for j=(M+1):(Columnas-M)
        resultado_num=0;
        resultado_den=0;
        for m=0:(N-2-Contador_F)   
            for n=0:(N-1)   
                resultado_provisional_num=(Mag_Gradiente(i-M+m,j-M+n)^2)*(sin(2*Orient_Gradiente(i-M+m,j-M+n)));
                resultado_provisional_den=(Mag_Gradiente(i-M+m,j-M+n)^2)*(cos(2*Orient_Gradiente(i-M+m,j-M+n)));
                
                resultado_num=resultado_num+resultado_provisional_num;
                resultado_den=resultado_den+resultado_provisional_den;
            end
        end
        if resultado_den==0 && resultado_num==0
            Angulo_Anisotropia(i,j)=0;
        else
            Angulo_Anisotropia(i,j)=0.5*atan(resultado_num/resultado_den);
        end
    end
    Contador_F=Contador_F+1;
end

Contador_F=0;
for i=(Filas-M+1):Filas
    Contador_C=0;
    for j=(Columnas-M+1):Columnas
        resultado_num=0;
        resultado_den=0;
        for m=0:(N-2-Contador_F)   
            for n=0:(N-2-Contador_C)   
                resultado_provisional_num=(Mag_Gradiente(i-M+m,j-M+n)^2)*(sin(2*Orient_Gradiente(i-M+m,j-M+n)));
                resultado_provisional_den=(Mag_Gradiente(i-M+m,j-M+n)^2)*(cos(2*Orient_Gradiente(i-M+m,j-M+n)));
                
                resultado_num=resultado_num+resultado_provisional_num;
                resultado_den=resultado_den+resultado_provisional_den;
            end
        end
        if resultado_den==0 && resultado_num==0
            Angulo_Anisotropia(i,j)=0;
        else
            Angulo_Anisotropia(i,j)=0.5*atan(resultado_num/resultado_den);
        end
        Contador_C=Contador_C+1;
    end
    Contador_F=Contador_F+1;
end

figure(7)
imshow(Angulo_Anisotropia)
Angulo_Orientacion=Angulo_Anisotropia+(pi/2); % Deben estar todos comprendidos entre 0 y pi
Angulo_Anisotropia=abs(Angulo_Anisotropia);


figure(8)
imshow(Angulo_Anisotropia)
figure(9)
imshow(Angulo_Orientacion)

% Calculamos la coherencia

Coherencia=[];

Contador_F=0;
for i=1:M
    Contador_C=0;
    for j=1:M
        resultado_num=0;
        resultado_den=0;
        for m=(M-Contador_F):(N-1)   % Recorremos todas las filas
            for n=(M-Contador_C):(N-1)   % Recorremos todas las columnas dentro de cada fila
                resultado_provisional_num=Mag_Gradiente(i-M+m,j-M+n)*cos(Orient_Gradiente(i-M+m,j-M+n)-Angulo_Anisotropia(i,j));
                resultado_provisional_den=Mag_Gradiente(i-M+m,j-M+n);
                
                resultado_num=resultado_num+resultado_provisional_num;
                resultado_den=resultado_den+resultado_provisional_den;
            end
        end
        if resultado_num==0&resultado_den==0
            Coherencia(i,j)=0;
        else        
            Coherencia(i,j)=Mag_Gradiente(i,j)*(resultado_num/resultado_den);
        end
        Contador_C=Contador_C+1;
    end
    Contador_F=Contador_F+1;
end

Contador_F=0;
for i=1:M
    for j=(M+1):(Columnas-M)
        resultado_num=0;
        resultado_den=0;
        for m=(M-Contador_F):(N-1)   
            for n=0:(N-1)   
                resultado_provisional_num=Mag_Gradiente(i-M+m,j-M+n)*cos(Orient_Gradiente(i-M+m,j-M+n)-Angulo_Anisotropia(i,j));
                resultado_provisional_den=Mag_Gradiente(i-M+m,j-M+n);
                
                resultado_num=resultado_num+resultado_provisional_num;
                resultado_den=resultado_den+resultado_provisional_den;
            end
        end
        if resultado_num==0&resultado_den==0
            Coherencia(i,j)=0;
        else        
            Coherencia(i,j)=Mag_Gradiente(i,j)*(resultado_num/resultado_den);
        end
        
    end
    Contador_F=Contador_F+1;
end

Contador_F=0;
for i=1:M
    Contador_C=0;
    for j=(Columnas-M+1):Columnas
        resultado_num=0;
        resultado_den=0;
        for m=(M-Contador_F):(N-1)   
            for n=0:(N-2-Contador_C)   
                resultado_provisional_num=Mag_Gradiente(i-M+m,j-M+n)*cos(Orient_Gradiente(i-M+m,j-M+n)-Angulo_Anisotropia(i,j));
                resultado_provisional_den=Mag_Gradiente(i-M+m,j-M+n);
                
                resultado_num=resultado_num+resultado_provisional_num;
                resultado_den=resultado_den+resultado_provisional_den;
            end
        end
        if resultado_num==0&resultado_den==0
            Coherencia(i,j)=0;
        else        
            Coherencia(i,j)=Mag_Gradiente(i,j)*(resultado_num/resultado_den);
        end
        
        Contador_C=Contador_C+1;
    end
    Contador_F=Contador_F+1;
end

%Analizamos desde la fila M+1 hasta la Fila-M
for i=(M+1):(Filas-M)
    Contador_C=0;
    for j=1:M
        resultado_num=0;
        resultado_den=0;
        for m=0:(N-1)   % Recorremos todas las filas
            for n=(M-Contador_C):(N-1)   % Recorremos todas las columnas dentro de cada fila
                resultado_provisional_num=Mag_Gradiente(i-M+m,j-M+n)*cos(Orient_Gradiente(i-M+m,j-M+n)-Angulo_Anisotropia(i,j));
                resultado_provisional_den=Mag_Gradiente(i-M+m,j-M+n);
                
                resultado_num=resultado_num+resultado_provisional_num;
                resultado_den=resultado_den+resultado_provisional_den;
            end
        end
        if resultado_num==0&resultado_den==0
            Coherencia(i,j)=0;
        else        
            Coherencia(i,j)=Mag_Gradiente(i,j)*(resultado_num/resultado_den);
        end
        
        Contador_C=Contador_C+1;
    end
    
    for j=(M+1):(Columnas-M)
        
        resultado_num=0;
        resultado_den=0;
        
        for m=0:(N-1)   % Recorremos todas las filas
            for n=0:(N-1)   % Recorremos todas las columnas dentro de cada fila
                resultado_provisional_num=Mag_Gradiente(i-M+m,j-M+n)*cos(Orient_Gradiente(i-M+m,j-M+n)-Angulo_Anisotropia(i,j));
                resultado_provisional_den=Mag_Gradiente(i-M+m,j-M+n);
                
                resultado_num=resultado_num+resultado_provisional_num;
                resultado_den=resultado_den+resultado_provisional_den;
            end
        end
        if resultado_num==0&resultado_den==0
            Coherencia(i,j)=0;
        else        
            Coherencia(i,j)=Mag_Gradiente(i,j)*(resultado_num/resultado_den);
        end      
        
    end

    Contador_C=0; 
    
    for j=(Columnas-M+1):Columnas
        resultado_num=0;
        resultado_den=0;
        for m=0:(N-1)   
            for n=0:(N-2-Contador_C)   
                resultado_provisional_num=Mag_Gradiente(i-M+m,j-M+n)*cos(Orient_Gradiente(i-M+m,j-M+n)-Angulo_Anisotropia(i,j));
                resultado_provisional_den=Mag_Gradiente(i-M+m,j-M+n);
                
                resultado_num=resultado_num+resultado_provisional_num;
                resultado_den=resultado_den+resultado_provisional_den;
            end
        end
        if resultado_num==0&resultado_den==0
            Coherencia(i,j)=0;
        else        
            Coherencia(i,j)=Mag_Gradiente(i,j)*(resultado_num/resultado_den);
        end
        Contador_C=Contador_C+1;
    end
end

% Analizamos las ultimas M Filas

Contador_F=0;
for i=(Filas-M+1):Filas
    Contador_C=0;
    for j=1:M
        resultado_num=0;
        resultado_den=0;
        for m=0:(N-2-Contador_F)   % Recorremos todas las filas
            for n=(M-Contador_C):(N-1)   % Recorremos todas las columnas dentro de cada fila
                resultado_provisional_num=Mag_Gradiente(i-M+m,j-M+n)*cos(Orient_Gradiente(i-M+m,j-M+n)-Angulo_Anisotropia(i,j));
                resultado_provisional_den=Mag_Gradiente(i-M+m,j-M+n);
                
                resultado_num=resultado_num+resultado_provisional_num;
                resultado_den=resultado_den+resultado_provisional_den;
            end
        end
        if resultado_num==0&resultado_den==0
            Coherencia(i,j)=0;
        else        
            Coherencia(i,j)=Mag_Gradiente(i,j)*(resultado_num/resultado_den);
        end
        Contador_C=Contador_C+1;
    end
    Contador_F=Contador_F+1;
end

Contador_F=0;
for i=(Filas-M+1):Filas
    for j=(M+1):(Columnas-M)
        resultado_num=0;
        resultado_den=0;
        for m=0:(N-2-Contador_F)   
            for n=0:(N-1)   
                resultado_provisional_num=Mag_Gradiente(i-M+m,j-M+n)*cos(Orient_Gradiente(i-M+m,j-M+n)-Angulo_Anisotropia(i,j));
                resultado_provisional_den=Mag_Gradiente(i-M+m,j-M+n);
                
                resultado_num=resultado_num+resultado_provisional_num;
                resultado_den=resultado_den+resultado_provisional_den;
            end
        end
        if resultado_num==0&resultado_den==0
            Coherencia(i,j)=0;
        else        
            Coherencia(i,j)=Mag_Gradiente(i,j)*(resultado_num/resultado_den);
        end
    end
    Contador_F=Contador_F+1;
end

Contador_F=0;

for i=(Filas-M+1):Filas
    Contador_C=0;
    for j=(Columnas-M+1):Columnas
        resultado_num=0;
        resultado_den=0;
        for m=0:(N-2-Contador_F)   
            for n=0:(N-2-Contador_C)   
                resultado_provisional_num=Mag_Gradiente(i-M+m,j-M+n)*cos(Orient_Gradiente(i-M+m,j-M+n)-Angulo_Anisotropia(i,j));
                resultado_provisional_den=Mag_Gradiente(i-M+m,j-M+n);
                
                resultado_num=resultado_num+resultado_provisional_num;
                resultado_den=resultado_den+resultado_provisional_den;
            end
        end
        if resultado_num==0&resultado_den==0
            Coherencia(i,j)=0;
        else        
            Coherencia(i,j)=Mag_Gradiente(i,j)*(resultado_num/resultado_den);
        end
        Contador_C=Contador_C+1;
    end
    Contador_F=Contador_F+1;
end

figure(10)
imshow(Coherencia)

Maximo_Valor=max(max(Coherencia));

Coherencia_Maxima=Coherencia/Maximo_Valor;
figure(11)
imshow(Coherencia_Maxima)

