
close all
clear all

alpha=0.1;
min_npoints=15;


% A.-IMAGE PREPROCESSING

% OCULAR FUNDUS IMAGE
Ocular_Fundus_Image=imread('01_original.bmp');
Mascara=imread('01_Mask.gif');
% figure
% imshow(Ocular_Fundus_Image)
% Title('Ocular Fundus Image')


% GREEN CHANNEL
Green_Channel=double(Ocular_Fundus_Image(:,:,2))/255;
%Green_Channel=double(rgb2gray(Ocular_Fundus_Image))/255;
figure
imshow(Green_Channel)
Title('Green Channel')
[Filas,Columnas]=size(Green_Channel);
imwrite(Green_Channel,'Green_Channel.bmp')


% A.1.-BACKGROUND NORMALIZATION
Filtro_Media=fspecial('average',[31,31]);
Imagen_Filtrada=imfilter(Green_Channel,Filtro_Media,'conv');
% figure
% imshow(Imagen_Filtrada)
% Title('Imagen Filtrada')

Background_Normalization=Imagen_Filtrada-Green_Channel;
% Eliminamos el borde de la imagen a traves de la mascara
for i=1:Filas
    for j=1:Columnas
        if Mascara(i,j)==0
            Background_Normalization(i,j)=0;
        end
    end
end

figure
imshow(Background_Normalization)
Title('Background Normalization')
imwrite(Background_Normalization,'Background_Normalization.bmp')

    
% A.2.-THIN VESSEL ENHANCEMENT
LDF1=1/6*[-1 -1 -1; 2 2 2; -1 -1 -1];
LDF2=1/6*[2 -1 -1; -1 2 -1; -1 -1 2];
LDF3=1/6*[-1 2 -1; -1 2 -1; -1 2 -1];
LDF4=1/6*[-1 -1 2; -1 2 -1; 2 -1 -1];

Orientacion_0=imfilter(Background_Normalization,LDF1);
Orientacion_45=imfilter(Background_Normalization,LDF2);
Orientacion_90=imfilter(Background_Normalization,LDF3);
Orientacion_135=imfilter(Background_Normalization,LDF4);

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

Thin_Vessel_Enhancement=Background_Normalization+Realzado;

figure
imshow(Thin_Vessel_Enhancement)
Title('Thin Vessel Enhancement')            
imwrite(Thin_Vessel_Enhancement,'Thin_Vessel_Enhancement.bmp')


% B.-DETECTION OF CENTERLINE SEGMENT CANDIDATES

% B.1.-INITIAL SELECTION OF CANDIDATE POINTS

DoOG_Vertical=[-1 -2 0 2 1; -2 -4 0 4 2; -1 -2 0 2 1];
DoOG_Horizontal=[-1 -2 -1; -2 -4 -2; 0 0 0; 2 4 2; 1 2 1];
DoOG_45=[0 0 -1 0 0 0 0; 0 -2 0 -2 0 0 0; -1 0 -4 0 0 0 0; 0 -2 0 0 0 2 0; 0 0 0 0 4 0 1; 0 0 0 2 0 2 0; 0 0 0 0 1 0 0];
DoOG_135=[0 0 0 0 1 0 0; 0 0 0 2 0 2 0; 0 0 0 0 4 0 1; 0 -2 0 0 0 2 0; -1 0 -4 0 0 0 0; 0 -2 0 -2 0 0 0; 0 0 -1 0 0 0 0];

Vessel_Vertical=imfilter(Thin_Vessel_Enhancement,DoOG_Vertical);
Vessel_Horizontal=imfilter(Thin_Vessel_Enhancement,DoOG_Horizontal);
Vessel_45=imfilter(Thin_Vessel_Enhancement,DoOG_45);
Vessel_135=imfilter(Thin_Vessel_Enhancement,DoOG_135);

figure
imshow(Vessel_Vertical)
Title('Vessel Vertical')
imwrite(Vessel_Vertical,'Vessel_Vertical.bmp')

figure
imshow(Vessel_Horizontal)
Title('Vessel Horizontal')
imwrite(Vessel_Horizontal,'Vessel_Horizontal.bmp')

figure
imshow(Vessel_45)
Title('Vessel 45')
imwrite(Vessel_45,'Vessel_45.bmp')

figure
imshow(Vessel_135)
Title('Vessel 135')
imwrite(Vessel_135,'Vessel_135.bmp')

% Initial Set Candidate Points
ISCP_Vertical=zeros(Filas,Columnas);
ISCP_Horizontal=zeros(Filas,Columnas);
ISCP_45=zeros(Filas,Columnas);
ISCP_135=zeros(Filas,Columnas);

% Buscamos los signos para los Vasos Verticales y sumamos mayor positivo y
% valor absoluto de mayor negativo
for i=3:(Filas-2)
    for j=3:(Columnas-2)
        if Vessel_Vertical(i,j)==0&Vessel_Vertical(i,j-1)>0&Vessel_Vertical(i,j+1)<0
            ISCP_Vertical(i,j)=Vessel_Vertical(i,j-1)+abs(Vessel_Vertical(i,j+1));
        elseif Vessel_Vertical(i,j-2)>0&Vessel_Vertical(i,j-1)>0&Vessel_Vertical(i,j)<0&Vessel_Vertical(i,j+1)<0
            m=1;
            X(m)=Vessel_Vertical(i,j-2);
            m=m+1;
            X(m)=Vessel_Vertical(i,j-1);
            m=m+1;
            X(m)=Vessel_Vertical(i,j);
            m=m+1;
            X(m)=Vessel_Vertical(i,j+1);
            Maximo=max(X);
            Minimo=min(X);
            ISCP_Vertical(i,j)=Maximo+abs(Minimo);
                    
        elseif Vessel_Vertical(i,j-2)>0&Vessel_Vertical(i,j-1)>0&Vessel_Vertical(i,j)<0
            if (Vessel_Vertical(i,j-2)+Vessel_Vertical(i,j-1)+Vessel_Vertical(i,j))>0
                m=1;
                X(m)=Vessel_Vertical(i,j-2);
                m=m+1;
                X(m)=Vessel_Vertical(i,j-1);
                m=m+1;
                X(m)=Vessel_Vertical(i,j);
                Maximo=max(X);
                Minimo=min(X);
                ISCP_Vertical(i,j)=Maximo+abs(Minimo);
                
            elseif Vessel_Vertical(i,j-1)>0&Vessel_Vertical(i,j)<0&Vessel_Vertical(i,j+1)<0
                if (Vessel_Vertical(i,j-1)+Vessel_Vertical(i,j)+Vessel_Vertical(i,j+1))<0
                    m=1;
                    X(m)=Vessel_Vertical(i,j-1);
                    m=m+1;
                    X(m)=Vessel_Vertical(i,j);
                    m=m+1;
                    X(m)=Vessel_Vertical(i,j+1);
                    Maximo=max(X);
                    Minimo=min(X);
                    ISCP_Vertical(i,j)=Maximo+abs(Minimo);
                    
                else
                    ISCP_Vertical(i,j)=0;
                end
            end
        end
    end
end

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

figure
imshow(ISCP_Vertical)
Title('ISCP Vertical')
imwrite(ISCP_Vertical,'ISCP_Vertical.bmp')

m=0;
for i=1:Filas
    for j=1:Columnas
        if ISCP_Vertical(i,j)~=0
            m=m+1;
            X_Vertical(m)=ISCP_Vertical(i,j);
        end
    end
end

[Histograma,N]=hist(ISCP_Vertical,60000);
Suma=sum(Histograma,2);
Suma(1)=0;
[M,I]=max(Suma);
Agregacion_V=N(I);
clear Histograma N Suma

Mean_V=mean(X_Vertical);
Std_V=std(X_Vertical);
T_Seed_V=Mean_V-alpha*Std_V;
Mascara_V=zeros(Filas,Columnas);
Sumareg_V=zeros(Filas,Columnas);


m=1;
for i=1:Filas
    for j=1:Columnas
        if ISCP_Vertical(i,j)>=T_Seed_V
            [Sumareg_V,Mascara_V,Longitud_V,Intensidad_Media_V,Coordenadas]=CrecimientoRegionVH(ISCP_Vertical,i,j,Mascara_V,Sumareg_V,Agregacion_V);
            if Longitud_V>=min_npoints 
                Candidates_V(m,1)=i;
                Candidates_V(m,2)=j;
                Candidates_V(m,3)=Longitud_V;
                p=1;
                for k=4:2:(2*Longitud_V+3)
                    Candidates_V(m,k)=Coordenadas(p,1);
                    Candidates_V(m,k+1)=Coordenadas(p,2);
                    p=p+1;
                end
                Intensidad_V(m)=Intensidad_Media_V;
                m=m+1;
            end
        end
    end
end


figure
imshow(Sumareg_V)
Title('Sumareg V')
imwrite(Sumareg_V,'Sumareg_V.bmp')

[Vessel_Centerline_Segments_V,Possible_Vessel_Centerline_Segments_V,Sub_Intensidad_V,Sub_Candidates_V]=VasosSemiDefinitivos(Intensidad_V,Candidates_V,Filas,Columnas,Sumareg_V);

figure
imshow(Vessel_Centerline_Segments_V)
Title('Vessel Centerline Segments_V')
imwrite(Vessel_Centerline_Segments_V,'Vessel_Centerline_Segments_V.bmp')

figure
imshow(Possible_Vessel_Centerline_Segments_V)
Title('Possible Vessel Centerline Segments V')
imwrite(Possible_Vessel_Centerline_Segments_V,'Possible_Vessel_Centerline_Segments_V.bmp')


% Buscamos los signos para los Vasos Horizontales y sumamos mayor positivo y
% valor absoluto de mayor negativo

for i=3:(Filas-2)
    for j=3:(Columnas-2)
        if Vessel_Horizontal(i,j)==0&Vessel_Horizontal(i-1,j)>0&Vessel_Horizontal(i+1,j)<0
            ISCP_Horizontal(i,j)=Vessel_Horizontal(i-1,j)+abs(Vessel_Horizontal(i+1,j));
        elseif Vessel_Horizontal(i-2,j)>0&Vessel_Horizontal(i-1,j)>0&Vessel_Horizontal(i,j)<0&Vessel_Horizontal(i+1,j)<0
            m=1;
            X(m)=Vessel_Horizontal(i-2,j);
            m=m+1;
            X(m)=Vessel_Horizontal(i-1,j);
            m=m+1;
            X(m)=Vessel_Horizontal(i,j);
            m=m+1;
            X(m)=Vessel_Horizontal(i+1,j);
            Maximo=max(max(X));
            Minimo=min(min(X));
            ISCP_Horizontal(i,j)=Maximo+abs(Minimo);
        elseif Vessel_Horizontal(i-2,j)>0&Vessel_Horizontal(i-1,j)>0&Vessel_Horizontal(i,j)<0
            if (Vessel_Horizontal(i-2,j)+Vessel_Horizontal(i-1,j)+Vessel_Horizontal(i,j))>0
                m=1;
                X(m)=Vessel_Horizontal(i-2,j);
                m=m+1;
                X(m)=Vessel_Horizontal(i-1,j);
                m=m+1;
                X(m)=Vessel_Horizontal(i,j);
                Maximo=max(max(X));
                Minimo=min(min(X));
                ISCP_Horizontal(i,j)=Maximo+abs(Minimo);
            elseif Vessel_Horizontal(i-1,j)>0&Vessel_Horizontal(i,j)<0&Vessel_Horizontal(i+1,j)<0
                if (Vessel_Horizontal(i-1,j)+Vessel_Horizontal(i,j)+Vessel_Horizontal(i+1,j))<0
                    m=1;
                    X(m)=Vessel_Horizontal(i-1,j);
                    m=m+1;
                    X(m)=Vessel_Horizontal(i,j);
                    m=m+1;
                    X(m)=Vessel_Horizontal(i+1,j);
                    Maximo=max(max(X));
                    Minimo=min(min(X));
                    ISCP_Horizontal(i,j)=Maximo+abs(Minimo);
                else
                    ISCP_Horizontal(i,j)=0;
                end
            end
        end
    end
end

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

figure
imshow(ISCP_Horizontal)
Title('ISCP Horizontal')
imwrite(ISCP_Horizontal,'ISCP_Horizontal.bmp')

m=0;
for i=1:Filas
    for j=1:Columnas
        if ISCP_Horizontal(i,j)~=0
            m=m+1;
            X_Horizontal(m)=ISCP_Horizontal(i,j);
        end
    end
end

[Histograma,N]=hist(ISCP_Horizontal,60000);
Suma=sum(Histograma,2);
Suma(1)=0;
[M,I]=max(Suma);
Agregacion_H=N(I);
clear Histograma N Suma

Mean_H=mean(X_Horizontal);
Std_H=std(X_Horizontal);
T_Seed_H=Mean_H-alpha*Std_H;
Mascara_H=zeros(Filas,Columnas);
Sumareg_H=zeros(Filas,Columnas);
T_Seed_H=T_Seed_H;

m=1;
for i=1:Filas
    for j=1:Columnas
        if ISCP_Horizontal(i,j)>=T_Seed_H
            [Sumareg_H,Mascara_H,Longitud_H,Intensidad_Media_H,Coordenadas]=CrecimientoRegionVH(ISCP_Horizontal,i,j,Mascara_H,Sumareg_H,Agregacion_H);
            if Longitud_H>=min_npoints 
                Candidates_H(m,1)=i;
                Candidates_H(m,2)=j;
                Candidates_H(m,3)=Longitud_H;
                p=1;
                for k=4:2:(2*Longitud_H+3)
                    Candidates_H(m,k)=Coordenadas(p,1);
                    Candidates_H(m,k+1)=Coordenadas(p,2);
                    p=p+1;
                end
                Intensidad_H(m)=Intensidad_Media_H;
                m=m+1;
            end
            
        end
    end
end


figure
imshow(Sumareg_H)
Title('Sumareg H')
imwrite(Sumareg_H,'Sumareg_H.bmp')

[Vessel_Centerline_Segments_H,Possible_Vessel_Centerline_Segments_H,Sub_Intensidad_H,Sub_Candidates_H]=VasosSemiDefinitivos(Intensidad_H,Candidates_H,Filas,Columnas,Sumareg_H);

figure
imshow(Vessel_Centerline_Segments_H)
Title('Vessel Centerline Segments H')
imwrite(Vessel_Centerline_Segments_H,'Vessel_Centerline_Segments_H.bmp')

figure
imshow(Possible_Vessel_Centerline_Segments_H)
Title('Possible Vessel Centerline Segments H')
imwrite(Possible_Vessel_Centerline_Segments_H,'Possible_Vessel_Centerline_Segments_H.bmp')

       
% Buscamos los signos para los Vasos de 45º y sumamos mayor positivo y
% valor absoluto de mayor negativo

for i=3:(Filas-2)
    for j=3:(Columnas-2)
        if Vessel_45(i,j)==0&Vessel_45(i-1,j-1)>0&Vessel_45(i+1,j+1)<0
            ISCP_45(i,j)=Vessel_45(i-1,j-1)+abs(Vessel_45(i+1,j+1));
        elseif Vessel_45(i-2,j-2)>0&Vessel_45(i-1,j-1)>0&Vessel_45(i,j)<0&Vessel_45(i+1,j+1)<0
            m=1;
            X(m)=Vessel_45(i-2,j-2);
            m=m+1;
            X(m)=Vessel_45(i-1,j-1);
            m=m+1;
            X(m)=Vessel_45(i,j);
            m=m+1;
            X(m)=Vessel_45(i+1,j+1);
            Maximo=max(max(X));
            Minimo=min(min(X));
            ISCP_45(i,j)=Maximo+abs(Minimo);
        elseif Vessel_45(i-2,j-2)>0&Vessel_45(i-1,j-1)>0&Vessel_45(i,j)<0
            if (Vessel_45(i-2,j-2)+Vessel_45(i-1,j-1)+Vessel_45(i,j))>0
                m=1;
                X(m)=Vessel_45(i-2,j-2);
                m=m+1;
                X(m)=Vessel_45(i-1,j-1);
                m=m+1;
                X(m)=Vessel_45(i,j);
                Maximo=max(max(X));
                Minimo=min(min(X));
                ISCP_45(i,j)=Maximo+abs(Minimo);
            elseif Vessel_45(i-1,j-1)>0&Vessel_45(i,j)<0&Vessel_45(i+1,j+1)<0
                if (Vessel_45(i-1,j-1)+Vessel_45(i,j)+Vessel_45(i+1,j+1))<0
                    m=1;
                    X(m)=Vessel_45(i-1,j-1);
                    m=m+1;
                    X(m)=Vessel_45(i,j);
                    m=m+1;
                    X(m)=Vessel_45(i+1,j+1);
                    Maximo=max(max(X));
                    Minimo=min(min(X));
                    ISCP_45(i,j)=Maximo+abs(Minimo);
                else
                    ISCP_45(i,j)=0;
                end
            end
        end
    end
end

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

figure
imshow(ISCP_45)
Title('ISCP 45')
imwrite(ISCP_45,'ISCP_45.bmp')

m=0;
for i=1:Filas
    for j=1:Columnas
        if ISCP_45(i,j)~=0
            m=m+1;
            X_45(m)=ISCP_45(i,j);
        end
    end
end

[Histograma,N]=hist(ISCP_45,60000);
Suma=sum(Histograma,2);
Suma(1)=0;
[M,I]=max(Suma);
Agregacion_45=N(I);
clear Histograma N Suma

Mean_45=mean(X_45);
Std_45=std(X_45);
T_Seed_45=Mean_45-alpha*Std_45;
Mascara_45=zeros(Filas,Columnas);
Sumareg_45=zeros(Filas,Columnas);

m=1;
for i=1:Filas
    for j=1:Columnas
        if ISCP_45(i,j)>=T_Seed_45
            [Sumareg_45,Mascara_45,Longitud_45,Intensidad_Media_45,Coordenadas]=CrecimientoRegion(ISCP_45,i,j,Mascara_45,Sumareg_45,Agregacion_45);
            if Longitud_45>=min_npoints 
                Candidates_45(m,1)=i;
                Candidates_45(m,2)=j;
                Candidates_45(m,3)=Longitud_45;
                p=1;
                for k=4:2:(2*Longitud_45+3)
                    Candidates_45(m,k)=Coordenadas(p,1);
                    Candidates_45(m,k+1)=Coordenadas(p,2);
                    p=p+1;
                end
                Intensidad_45(m)=Intensidad_Media_45;
                m=m+1;
            end
            
        end
    end
end


figure
imshow(Sumareg_45)
Title('Sumareg 45')
imwrite(Sumareg_45,'Sumareg_45.bmp')

[Vessel_Centerline_Segments_45,Possible_Vessel_Centerline_Segments_45,Sub_Intensidad_45,Sub_Candidates_45]=VasosSemiDefinitivos(Intensidad_45,Candidates_45,Filas,Columnas,Sumareg_45);

figure
imshow(Vessel_Centerline_Segments_45)
Title('Vessel Centerline Segments 45')
imwrite(Vessel_Centerline_Segments_45,'Vessel_Centerline_Segments_45.bmp')

figure
imshow(Possible_Vessel_Centerline_Segments_45)
Title('Possible Vessel Centerline Segments 45')
imwrite(Possible_Vessel_Centerline_Segments_45,'Possible_Vessel_Centerline_Segments_45.bmp')


% Buscamos los signos para los Vasos de 45º y sumamos mayor positivo y
% valor absoluto de mayor negativo

for i=3:(Filas-2)
    for j=3:(Columnas-2)
        if Vessel_135(i,j)==0&Vessel_135(i+1,j-1)>0&Vessel_135(i-1,j+1)<0
            ISCP_135(i,j)=Vessel_135(i+1,j-1)+abs(Vessel_135(i-1,j+1));
        elseif Vessel_135(i+2,j-2)>0&Vessel_135(i+1,j-1)>0&Vessel_135(i,j)<0&Vessel_135(i-1,j+1)<0
            m=1;
            X(m)=Vessel_135(i+2,j-2);
            m=m+1;
            X(m)=Vessel_135(i+1,j-1);
            m=m+1;
            X(m)=Vessel_135(i,j);
            m=m+1;
            X(m)=Vessel_135(i-1,j+1);
            Maximo=max(max(X));
            Minimo=min(min(X));
            ISCP_135(i,j)=Maximo+abs(Minimo);
        elseif Vessel_135(i+2,j-2)>0&Vessel_135(i+1,j-1)>0&Vessel_135(i,j)<0
            if (Vessel_135(i+2,j-2)+Vessel_135(i+1,j-1)+Vessel_135(i,j))<0
                m=1;
                X(m)=Vessel_135(i+2,j-2);
                m=m+1;
                X(m)=Vessel_135(i+1,j-1);
                m=m+1;
                X(m)=Vessel_135(i,j);
                Maximo=max(max(X));
                Minimo=min(min(X));
                ISCP_135(i,j)=Maximo+abs(Minimo);
            elseif Vessel_135(i+1,j-1)>0&Vessel_135(i,j)<0&Vessel_135(i-1,j+1)<0
                if (Vessel_135(i+1,j-1)+Vessel_135(i,j)+Vessel_135(i-1,j+1))<0
                    m=1;
                    X(m)=Vessel_135(i+1,j-1);
                    m=m+1;
                    X(m)=Vessel_135(i,j);
                    m=m+1;
                    X(m)=Vessel_135(i-1,j+1);
                    Maximo=max(max(X));
                    Minimo=min(min(X));
                    ISCP_135(i,j)=Maximo+abs(Minimo);
                else
                    ISCP_135(i,j)=0;
                end
            end
        end
    end
end

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

figure
imshow(ISCP_135)
Title('ISCP 135')
imwrite(ISCP_135,'ISCP_135.bmp')

m=0;
for i=1:Filas
    for j=1:Columnas
        if ISCP_135(i,j)~=0
            m=m+1;
            X_135(m)=ISCP_135(i,j);
        end
    end
end

[Histograma,N]=hist(ISCP_135,60000);
Suma=sum(Histograma,2);
Suma(1)=0;
[M,I]=max(Suma);
Agregacion_135=N(I);
clear Histograma N Suma

Mean_135=mean(X_135);
Std_135=std(X_135);
T_Seed_135=Mean_135-alpha*Std_135;
Mascara_135=zeros(Filas,Columnas);
Sumareg_135=zeros(Filas,Columnas);

m=1;
for i=1:Filas
    for j=1:Columnas
        if ISCP_135(i,j)>=T_Seed_135
            [Sumareg_135,Mascara_135,Longitud_135,Intensidad_Media_135,Coordenadas]=CrecimientoRegion(ISCP_135,i,j,Mascara_135,Sumareg_135,Agregacion_135);
            if Longitud_135>=min_npoints 
                Candidates_135(m,1)=i;
                Candidates_135(m,2)=j;
                Candidates_135(m,3)=Longitud_135;
                p=1;
                for k=4:2:(2*Longitud_135+3)
                    Candidates_135(m,k)=Coordenadas(p,1);
                    Candidates_135(m,k+1)=Coordenadas(p,2);
                    p=p+1;
                end
                Intensidad_135(m)=Intensidad_Media_135;
                m=m+1;
            end
            
        end
    end
end

figure
imshow(Sumareg_135)
Title('Sumareg 135')
imwrite(Sumareg_135,'Sumareg_135.bmp')

[Vessel_Centerline_Segments_135,Possible_Vessel_Centerline_Segments_135,Sub_Intensidad_135,Sub_Candidates_135]=VasosSemiDefinitivos(Intensidad_135,Candidates_135,Filas,Columnas,Sumareg_135);

figure
imshow(Vessel_Centerline_Segments_135)
Title('Vessel Centerline Segments 135')
imwrite(Vessel_Centerline_Segments_135,'Vessel_Centerline_Segments_135.bmp')

figure
imshow(Possible_Vessel_Centerline_Segments_135)
Title('Possible Vessel Centerline Segments 135')
imwrite(Possible_Vessel_Centerline_Segments_135,'Possible_Vessel_Centerline_Segments_135.bmp')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Media=mean(Sub_Candidates_V,1);
Longitud_Media_V=Media(3);
Intensidad_Media_V=mean(Sub_Intensidad_V)-std(Sub_Intensidad_V);

Media=mean(Sub_Candidates_H,1);
Longitud_Media_H=Media(3);
Intensidad_Media_H=mean(Sub_Intensidad_H)-std(Sub_Intensidad_H);
        
Media=mean(Sub_Candidates_45,1);
Longitud_Media_45=Media(3)/sqrt(2);
Intensidad_Media_45=mean(Sub_Intensidad_45)-std(Sub_Intensidad_45);

Media=mean(Sub_Candidates_135,1);
Longitud_Media_135=Media(3)/sqrt(2);
Intensidad_Media_135=mean(Sub_Intensidad_135)-std(Sub_Intensidad_135);

Longitud_Media_G=(Longitud_Media_V+Longitud_Media_H+Longitud_Media_45+Longitud_Media_135)/4;
Intensidad_Media_G=AnalisisImagenGlobal(Candidates_V,Candidates_H,Candidates_45,Candidates_135,Intensidad_V,Intensidad_H,Intensidad_45,Intensidad_135);



% BUSCAMOS LOS VALORES DE REFERENCIA DE INTENSIDAD
if Intensidad_Media_V>=Intensidad_Media_G
    Intensidad_Ref_V=Intensidad_Media_V;
else
    Intensidad_Ref_V=Intensidad_Media_G;
end
if Intensidad_Media_H>=Intensidad_Media_G
    Intensidad_Ref_H=Intensidad_Media_H;
else
    Intensidad_Ref_H=Intensidad_Media_G;
end
if Intensidad_Media_45>=Intensidad_Media_G
    Intensidad_Ref_45=Intensidad_Media_45;
else
    Intensidad_Ref_45=Intensidad_Media_G;
end
if Intensidad_Media_135>=Intensidad_Media_G
    Intensidad_Ref_135=Intensidad_Media_135;
else
    Intensidad_Ref_135=Intensidad_Media_G;
end

%BUSCAMOS LOS VALORES DE REFERENCIA DE LONGITUD
if Longitud_Media_V>=Longitud_Media_G
    Longitud_Ref_V=Longitud_Media_V;
else
    Longitud_Ref_V=Longitud_Media_G;
end
if Longitud_Media_H>=Longitud_Media_G
    Longitud_Ref_H=Longitud_Media_H;
else
    Longitud_Ref_H=Longitud_Media_G;
end
if Longitud_Media_45>=Longitud_Media_G
    Longitud_Ref_45=Longitud_Media_45;
else
    Longitud_Ref_45=Longitud_Media_G;
end
if Longitud_Media_135>=Longitud_Media_G
    Longitud_Ref_135=Longitud_Media_135;
else
    Longitud_Ref_135=Longitud_Media_G;
end

Valor_Ref_V=Intensidad_Ref_V*sqrt(Longitud_Ref_V);
Valor_Ref_H=Intensidad_Ref_H*sqrt(Longitud_Ref_H);
Valor_Ref_45=Intensidad_Ref_45*sqrt(Longitud_Ref_45);
Valor_Ref_135=Intensidad_Ref_135*sqrt(Longitud_Ref_135);


%NOS QUEDAMOS CON LOS VASOS DEFINITIVOS

Vessel_Centerline_Segments_V=VasosDefinitivos(Sub_Intensidad_V,Sub_Candidates_V,Valor_Ref_V,Vessel_Centerline_Segments_V);
Vessel_Centerline_Segments_H=VasosDefinitivos(Sub_Intensidad_H,Sub_Candidates_H,Valor_Ref_H,Vessel_Centerline_Segments_H);
Vessel_Centerline_Segments_45=VasosDefinitivos(Sub_Intensidad_45,Sub_Candidates_45,Valor_Ref_45,Vessel_Centerline_Segments_45);
Vessel_Centerline_Segments_135=VasosDefinitivos(Sub_Intensidad_135,Sub_Candidates_135,Valor_Ref_135,Vessel_Centerline_Segments_135);

figure
imshow(Vessel_Centerline_Segments_V)
Title('Vessel Centerline Segments V Definitivo')
imwrite(Vessel_Centerline_Segments_V,'Vessel_Centerline_Segments_V_Definitivo.bmp')
figure
imshow(Vessel_Centerline_Segments_H)
Title('Vessel Centerline Segments H Definitivo')
imwrite(Vessel_Centerline_Segments_H,'Vessel_Centerline_Segments_H_Definitivo.bmp')
figure
imshow(Vessel_Centerline_Segments_45)
Title('Vessel Centerline Segments 45 Definitivo')
imwrite(Vessel_Centerline_Segments_45,'Vessel_Centerline_Segments_45_Definitivo.bmp')
figure
imshow(Vessel_Centerline_Segments_135)
Title('Vessel Centerline Segments V Definitivo')
imwrite(Vessel_Centerline_Segments_135,'Vessel_Centerline_Segments_135_Definitivo.bmp')


Resultado_Provisional=Vessel_Centerline_Segments_V+Vessel_Centerline_Segments_H+Vessel_Centerline_Segments_45+Vessel_Centerline_Segments_135;

figure
imshow(Resultado_Provisional)
Title('Resultado Provisional')
imwrite(Resultado_Provisional,'Resultado_Provisional.bmp')


% VESSEL ENHANCEMENT

clear Orientacion_0 Orientacion_45 Orientacion_90 Orientacion_135 ISCP_Vertical ISCP_Horizontal ISCP_45 ISCP_135
clear Mascara_V Mascara_H Mascara_45 Mascara_135 Sumareg_H Sumareg_V Sumareg_45 Sumareg_135
clear Possible_Vessel_Centerline_Segments_V Possible_Vessel_Centerline_Segments_H Possible_Vessel_Centerline_Segments_45 Possible_Vessel_Centerline_Segments_135
clear Vessel_Centerline_Segments_V Vessel_Centerline_Segments_H Vessel_Centerline_Segments_45 Vessel_Centerline_Segments_135
clear Vessel_Vertical Vessel_Horizontal Vessel_45 Vessel_135 

So1=strel('disk',1);
So2=strel('disk',2);
So3=strel('disk',3);
So4=strel('disk',4);
So5=strel('disk',5);
So6=strel('disk',6);
So7=strel('disk',7);
So8=strel('disk',8);

Sc=strel('disk',1);

TopHat1=Background_Normalization-min(imopen(imclose(Background_Normalization,Sc),So1),Background_Normalization);
TopHat2=Background_Normalization-min(imopen(imclose(Background_Normalization,Sc),So2),Background_Normalization);
TopHat3=Background_Normalization-min(imopen(imclose(Background_Normalization,Sc),So3),Background_Normalization);
TopHat4=Background_Normalization-min(imopen(imclose(Background_Normalization,Sc),So4),Background_Normalization);
TopHat5=Background_Normalization-min(imopen(imclose(Background_Normalization,Sc),So5),Background_Normalization);
TopHat6=Background_Normalization-min(imopen(imclose(Background_Normalization,Sc),So6),Background_Normalization);
TopHat7=Background_Normalization-min(imopen(imclose(Background_Normalization,Sc),So7),Background_Normalization);
TopHat8=Background_Normalization-min(imopen(imclose(Background_Normalization,Sc),So8),Background_Normalization);

TopHat_1=(TopHat1+TopHat2)/2;
TopHat_2=(TopHat3+TopHat4)/2;
TopHat_3=(TopHat5+TopHat6)/2;
TopHat_4=(TopHat7+TopHat8)/2;

clear TopHat1 TopHat2 TopHat3 TopHat4 TopHat5 TopHat6 TopHat7 TopHat8

figure
imshow(TopHat_1)
Title('TopHat 1')
imwrite(TopHat_1,'TopHat_1.bmp')

figure
imshow(TopHat_2)
Title('TopHat 2')  
imwrite(TopHat_2,'TopHat_2.bmp')

figure
imshow(TopHat_3)
Title('TopHat 3')  
imwrite(TopHat_3,'TopHat_3.bmp')

figure
imshow(TopHat_4)
Title('TopHat 4')            
imwrite(TopHat_4,'TopHat_4.bmp')

% VESSEL SEGMENT RECONSTRUCTION

Tamano_1=0;
Tamano_2=0;
Tamano_3=0;
Tamano_4=0;
for i=1:Filas
    for j=1:Columnas
        if TopHat_1(i,j)~=0
            Tamano_1=Tamano_1+1;
        end
        if TopHat_2(i,j)~=0
            Tamano_2=Tamano_2+1;
        end
        if TopHat_3(i,j)~=0
            Tamano_3=Tamano_3+1;
        end
        if TopHat_4(i,j)~=0
            Tamano_4=Tamano_4+1;
        end
    end
end


T_Marker_1=0.1*Tamano_1;
T_Mask_1=0.25*Tamano_1;
T_Marker_2=0.1*Tamano_2;
T_Mask_2=0.2*Tamano_2;
T_Marker_3=0.05*Tamano_3;
T_Mask_3=0.15*Tamano_3;
T_Marker_4=0.01*Tamano_4;
T_Mask_4=0.1*Tamano_4;

%CALCULO DE MARKER IMAGE Y MASK IMAGE DE TOPHAT_1

Bins=60000;
[Hist,N]=hist(TopHat_1,Bins);
Suma=sum(Hist,2);
clear Hist

Contador=0;
C_Marker_1=Bins;

while Contador<=T_Marker_1,
    C_Marker_1=C_Marker_1-1;
    Contador=Contador+Suma(C_Marker_1);
end

Umbral_Marker_1=N(C_Marker_1)    
        
Marker_Image_1=TopHat_1;

for i=1:Filas
    for j=1:Columnas
        if Marker_Image_1(i,j)<Umbral_Marker_1
            Marker_Image_1(i,j)=0;
        else
            Marker_Image_1(i,j)=1;
        end
    end
end

figure
imshow(Marker_Image_1)
Title('Marker Image 1')
imwrite(Marker_Image_1,'Marker_Image_1.bmp')


Contador=0;
C_Mask_1=Bins;

while Contador<=T_Mask_1,
    C_Mask_1=C_Mask_1-1;
    Contador=Contador+Suma(C_Mask_1);
end

Umbral_Mask_1=N(C_Mask_1);       
        
Mask_Image_1=TopHat_1;

for i=1:Filas
    for j=1:Columnas
        if Mask_Image_1(i,j)<Umbral_Mask_1
            Mask_Image_1(i,j)=0;
        else
            Mask_Image_1(i,j)=1;
        end
    end
end

figure
imshow(Mask_Image_1)
Title('Mask Image 1') 
imwrite(Mask_Image_1,'Mask_Image_1.bmp')

%%CALCULO DE MARKER IMAGE Y MASK IMAGE DE TOPHAT_2

[Hist,N]=hist(TopHat_2,Bins);
Suma=sum(Hist,2);
clear Hist

Contador=0;
C_Marker_2=Bins;

while Contador<=T_Marker_2,
    C_Marker_2=C_Marker_2-1;
    Contador=Contador+Suma(C_Marker_2);
end

Umbral_Marker_2=N(C_Marker_2);       
        
Marker_Image_2=TopHat_2;

for i=1:Filas
    for j=1:Columnas
        if Marker_Image_2(i,j)<Umbral_Marker_2
            Marker_Image_2(i,j)=0;
        else
            Marker_Image_2(i,j)=1;
        end
    end
end

figure
imshow(Marker_Image_2)
Title('Marker Image 2')  
imwrite(Marker_Image_2,'Marker_Image_2.bmp')

Contador=0;
C_Mask_2=Bins;

while Contador<=T_Mask_2,
    C_Mask_2=C_Mask_2-1;
    Contador=Contador+Suma(C_Mask_2);
end

Umbral_Mask_2=N(C_Mask_2);       
        
Mask_Image_2=TopHat_2;

for i=1:Filas
    for j=1:Columnas
        if Mask_Image_2(i,j)<Umbral_Mask_2
            Mask_Image_2(i,j)=0;
        else
            Mask_Image_2(i,j)=1;
        end
    end
end

figure
imshow(Mask_Image_2)
Title('Mask Image 2') 
imwrite(Mask_Image_2,'Mask_Image_2.bmp')

%%CALCULO DE MARKER IMAGE Y MASK IMAGE DE TOPHAT_3

[Hist,N]=hist(TopHat_3,Bins);
Suma=sum(Hist,2);
clear Hist

Contador=0;
C_Marker_3=Bins;

while Contador<=T_Marker_3,
    C_Marker_3=C_Marker_3-1;
    Contador=Contador+Suma(C_Marker_3);
end

Umbral_Marker_3=N(C_Marker_3);       
        
Marker_Image_3=TopHat_3;

for i=1:Filas
    for j=1:Columnas
        if Marker_Image_3(i,j)<Umbral_Marker_3
            Marker_Image_3(i,j)=0;
        else
            Marker_Image_3(i,j)=1;
        end
    end
end

figure
imshow(Marker_Image_3)
Title('Marker Image 3')  
imwrite(Marker_Image_3,'Marker_Image_3.bmp')

Contador=0;
C_Mask_3=Bins;

while Contador<=T_Mask_3,
    C_Mask_3=C_Mask_3-1;
    Contador=Contador+Suma(C_Mask_3);
end

Umbral_Mask_3=N(C_Mask_3);       
        
Mask_Image_3=TopHat_3;

for i=1:Filas
    for j=1:Columnas
        if Mask_Image_3(i,j)<Umbral_Mask_3
            Mask_Image_3(i,j)=0;
        else
            Mask_Image_3(i,j)=1;
        end
    end
end

figure
imshow(Mask_Image_3)
Title('Mask Image 3') 
imwrite(Mask_Image_3,'Mask_Image_3.bmp')

%%CALCULO DE MARKER IMAGE Y MASK IMAGE DE TOPHAT_4

[Hist,N]=hist(TopHat_4,Bins);
Suma=sum(Hist,2);
clear Hist

Contador=0;
C_Marker_4=Bins;

while Contador<=T_Marker_4,
    C_Marker_4=C_Marker_4-1;
    Contador=Contador+Suma(C_Marker_4);
end

Umbral_Marker_4=N(C_Marker_4);       
        
Marker_Image_4=TopHat_4;

for i=1:Filas
    for j=1:Columnas
        if Marker_Image_4(i,j)<Umbral_Marker_4
            Marker_Image_4(i,j)=0;
        else
            Marker_Image_4(i,j)=1;
        end
    end
end

figure
imshow(Marker_Image_4)
Title('Marker Image 4')  
imwrite(Marker_Image_4,'Marker_Image_4.bmp')

Contador=0;
C_Mask_4=Bins;

while Contador<=T_Mask_4,
    C_Mask_4=C_Mask_4-1;
    Contador=Contador+Suma(C_Mask_4);
end

Umbral_Mask_4=N(C_Mask_4);       
        
Mask_Image_4=TopHat_4;

for i=1:Filas
    for j=1:Columnas
        if Mask_Image_4(i,j)<Umbral_Mask_4
            Mask_Image_4(i,j)=0;
        else
            Mask_Image_4(i,j)=1;
        end
    end
end

figure
imshow(Mask_Image_4)
Title('Mask Image 4') 
imwrite(Mask_Image_4,'Mask_Image_4.bmp')

%%MORFOLOGICAL RECONSTRUCTION OPERATOR

Mascara_Resultado_1=zeros(Filas,Columnas);

for i=1:Filas
    for j=1:Columnas
        if (Marker_Image_1(i,j)==1)&(Mascara_Resultado_1(i,j)==0)
            Mascara_Resultado_1=ObtencionOperadorMorfologico(i,j,Marker_Image_1,Mask_Image_1,Mascara_Resultado_1);
        end
    end
end

figure
imshow(Mascara_Resultado_1)
Title('Mascara Resultado 1') 
imwrite(Mascara_Resultado_1,'Mascara_Resultado_1.bmp')

Mascara_Resultado_2=zeros(Filas,Columnas);

for i=1:Filas
    for j=1:Columnas
        if (Marker_Image_2(i,j)==1)&(Mascara_Resultado_2(i,j)==0)
            Mascara_Resultado_2=ObtencionOperadorMorfologico(i,j,Marker_Image_2,Mask_Image_2,Mascara_Resultado_2);
        end
    end
end

figure
imshow(Mascara_Resultado_2)
Title('Mascara Resultado 2') 
imwrite(Mascara_Resultado_2,'Mascara_Resultado_2.bmp')

Mascara_Resultado_3=zeros(Filas,Columnas);

for i=1:Filas
    for j=1:Columnas
        if (Marker_Image_3(i,j)==1)&(Mascara_Resultado_3(i,j)==0)
            Mascara_Resultado_3=ObtencionOperadorMorfologico(i,j,Marker_Image_3,Mask_Image_3,Mascara_Resultado_3);
        end
    end
end

figure
imshow(Mascara_Resultado_3)
Title('Mascara Resultado 3') 
imwrite(Mascara_Resultado_3,'Mascara_Resultado_3.bmp')

Mascara_Resultado_4=zeros(Filas,Columnas);

for i=1:Filas
    for j=1:Columnas
        if (Marker_Image_4(i,j)==1)&(Mascara_Resultado_4(i,j)==0)
            Mascara_Resultado_4=ObtencionOperadorMorfologico(i,j,Marker_Image_4,Mask_Image_4,Mascara_Resultado_4);
        end
    end
end

figure
imshow(Mascara_Resultado_4)
Title('Mascara Resultado 4') 
imwrite(Mascara_Resultado_4,'Mascara_Resultado_4.bmp')

%CALCULO DE RESULTADO DEFINITIVO

Semilla_1=zeros(Filas,Columnas);
Semilla_2=zeros(Filas,Columnas);
Semilla_3=zeros(Filas,Columnas);
Resultado_Final=zeros(Filas,Columnas);

for i=1:Filas
    for j=1:Columnas
        if (Resultado_Provisional(i,j)==1)&(Semilla_1(i,j)==0)
            Semilla_1=ReconstruccionMorfologica(i,j,Resultado_Provisional,Mascara_Resultado_1,Semilla_1);
        end
    end
end

figure
imshow(Semilla_1)
Title('Semilla 1') 
imwrite(Semilla_1,'Semilla_1.bmp')


for i=1:Filas
    for j=1:Columnas
        if (Semilla_1(i,j)==1)&(Semilla_2(i,j)==0)
            Semilla_2=ReconstruccionMorfologica(i,j,Semilla_1,Mascara_Resultado_2,Semilla_2);
        end
    end
end

figure
imshow(Semilla_2)
Title('Semilla 2')
imwrite(Semilla_2,'Semilla_2.bmp')

for i=1:Filas
    for j=1:Columnas
        if (Semilla_2(i,j)==1)&(Semilla_3(i,j)==0)
            Semilla_3=ReconstruccionMorfologica(i,j,Semilla_2,Mascara_Resultado_3,Semilla_3);
        end
    end
end

figure
imshow(Semilla_3)
Title('Semilla 3')
imwrite(Semilla_3,'Semilla_3.bmp')

for i=1:Filas
    for j=1:Columnas
        if (Semilla_3(i,j)==1)&(Resultado_Final(i,j)==0)
            Resultado_Final=ReconstruccionMorfologica(i,j,Semilla_3,Mascara_Resultado_4,Resultado_Final);
        end
    end
end

figure
imshow(Resultado_Final)
Title('Resultado Final')
imwrite(Resultado_Final,'Resultado_Final_Inicial.bmp')

Resultado_Final=LimpiezaPixeles(Resultado_Final);

figure
imshow(Resultado_Final)
Title('Resultado Final')
imwrite(Resultado_Final,'Resultado_Final.bmp')

