%   KOMPLEXNI    [obr_vysl,obr_vysl_g,souradnice3] = komplexni(obr_barevny)
%   Funkce Komplexni provede globalni a lokalni dekteci krizeni cev
%   (bifurkací) na sítnici.
%
%   obr_vysl - vysledný obraz s vyznačenými detekovanými bifurkacemi
%   obr_vysl_g - obraz s vyznačenými detekovanými bifurkacemi(g-složka)
%   souradnice3 - matice se souřadnicemi nalezených bifurkací
%   obr_barevny - vstupni obraz sítnice (z fundus kamery nebo SLO)

function [obr_vysl,obr_vysl_g,souradnice3] = komplexni(obr_barevny)

%% upraveni vstupniho obrazu, aby program pracoval i se SLO obrazy
if numel(size(obr_barevny)) == 2
    disp('SLO')
    pom = obr_barevny;
    obr_barevny = [];
    obr_barevny(:,:,1) = pom;
    obr_barevny(:,:,2) = pom;
    obr_barevny(:,:,3) = pom;
    clear pom
end
%% zmena rozliseni obrazu pokud je vetsi nez 1500 x 1000 px
vel = size(obr_barevny);
if vel(1) > 1000 && vel(2) > 1500
    obr_barevny = imresize(obr_barevny,[1000 1500]);
    disp('Resize')
    clear vel
end
%%  lokální a globální detekce bifurkací

[obr2 kolo] = normal2(obr_barevny); % normalizace obrazu pro potlaceni nestejnorodého osvětlení
obr_vysl = obr2;
vel = size(obr2);   % zjištění rozllišení obrazu

pocet_x= 17; %pocet segmentu v rovine x
pocet_y = 13; %pocet segmentu v rovine y

% rozmer_x = floor(vel(1) / pocet_x);     %vypocet poctu pixelu v kazdem segmentu (pro zobrazení skeletonu)
% rozmer_y = floor(vel(2) / pocet_y);     %vypocet poctu pixelu v kazdem segmentu (pro zobrazení skeletonu)

sourad_global = globalni2(obr2);     %globalni detekce
disp('Konec globalni detekce');
[sourad_lokal1 d]= lokalni3(obr2,pocet_x,pocet_y);  % lokální detekce
disp('Konec lokalni detekce c.1');
pocet_x= 14; %pocet segmentu v rovine x
pocet_y = 10; %pocet segmentu v rovine y
sourad_lokal2 = lokalni3(obr2,pocet_x,pocet_y); % lokální detekce
disp('Konec lokalni detekce c.2');
%----------------------------------------------------------------------
%% sloceni souradnic z globalniho a lokalniho vyhledavani
souradnice = zeros(2,length(sourad_lokal1)+length(sourad_lokal2)+length(sourad_global));
souradnice(:,1:length(sourad_lokal1))=sourad_lokal1;
souradnice(:,length(sourad_lokal1)+1:end-length(sourad_global))=sourad_lokal2;
souradnice(:,end-length(sourad_global)+1:end)=sourad_global;

%----------------------------------------------------------------------
%% smazani bifurkaci mimo uzitecnou oblast
% figure(30)
% imshow(kolo)
kolo = bwmorph(kolo,'erode',20);
kolo(1:20,1:end) = 0;
kolo(end-20:end,1:end) = 0;
% figure(31)
% imshow(kolo)
for i = 1:vel(1)
    for j = 1:vel(2)
        if kolo(i,j) == 0
            souradnice(:,find (souradnice(1,:) == i & souradnice(2,:) == j))=[];
        end
    end
end

%----------------------------------------------------------------------
%% krizky
souradnice(:,find (souradnice(1,:)==0))=[]; % smazaní nulových souřadnic
% obr_kriz = obr2;
% for i = 1:length(souradnice)
%     obr_kriz(souradnice(1,i)-3:souradnice(1,i)+3,souradnice(2,i)) = 1;
%     obr_kriz(souradnice(1,i),souradnice(2,i)-3:souradnice(2,i)+3 ) = 1;
% end

% figure(3)
% imshow(obr_kriz,[])
% title('Vsechny nalezene bifurkace')
%----------------------------------------------------------------------

%% slouceni obrazu po skeletizaci
% obr_skel = double(zeros(vel(1),vel(2)));
% for i = 1:pocet_x
%     for j = 1:pocet_y  
%         obr_skel(rozmer_x*(i-1)+1:i*rozmer_x,rozmer_y*(j-1)+1:j*rozmer_y) = d(:,:,i,j) ;
%     end
% end
%----------------------------------------------------------------------

%zobrazeni obrazku po skeletizaci a nalezene bifurkace
%----------------------------------------------------------------------
%% znovuprojiti nalezenych bifurkaci:

sirka = 20;        %urcuje velikost okoli nalezene bifurkace
souradnice2 = souradnice;
clear sourad
sourad = zeros(2,300,length(souradnice2));

for i = 1:length(souradnice2)
    if souradnice2(1,i)<sirka %pokud je x-ova souradnice mensi nez zvolena sirka okoli bifurkace
        obr_p = obr2(1:souradnice2(1,i)+sirka,souradnice2(2,i)-sirka:souradnice2(2,i)+sirka);   %vyber okoli bifurkace
        stred = [souradnice2(1,i) sirka];   %pozice puvodni bifurkace
    elseif souradnice2(1,i)+sirka > 1000    %pokud by okoli bifurkace zasahovalo mimo puvodni obraz
        obr_p = obr2(souradnice2(1,i)-sirka:1000,souradnice2(2,i)-sirka:souradnice2(2,i)+sirka);    %vyber okoli bifurkace
        stred = [souradnice2(1,i) sirka];   %pozice puvodni bifurkace
    else        
        obr_p = obr2(souradnice2(1,i)-sirka:souradnice2(1,i)+sirka,souradnice2(2,i)-sirka:souradnice2(2,i)+sirka);  %vyber okoli bifurkace
        stred = [sirka sirka];   %pozice puvodni bifurkace     
    end
    
    obr_p=obr_p*1/max(max(obr_p));  
    his = imhist(obr_p);
    h = find(his>3);
    if (h(1)/256)>0.8   %pokud vsechny intenzity lezi nad 0.8 -> smaze bifurkaci
        souradnice2(1,i) = 1;
        souradnice2(2,i) = 1;
    else    
%         figure(7)
%         subplot(2,3,1)
%         imshow(obr_p)
%         subplot(2,3,2)
%         imhist(obr_p)   
%         hold on
        
%         lvl2 = graythresh(obr_p); 
%         lvl = prah_kittler2(obr_p);
        [lvl T] = prah_entropy_kapur(obr_p);    % zjištění prahu Kapůrovou metodou
        
%         plot([lvl lvl],[0 1200],'Color','red')
%         plot([lvl2 lvl2],[0 1200],'Color','green')
%         hold off
        
%         obr_p2 = im2bw(obr_p,lvl2);
        obr_p = im2bw(obr_p,lvl);   % prahování
     
%         subplot(2,3,3)
%         imshow(obr_p)
%         title('Kittler')
%         subplot(2,3,4)
%         imshow(obr_p2)
%         title('graytresh')
        
%         obr_p = obr_p + obr_p2;
%         obr_p(obr_p == 2) = 1;
        
%         subplot(2,3,5)
%         imshow(obr_p2)
%         title('Slouceny')
        
        if sum(sum(obr_p))>0.7*sum(h)   % pokud je poč. 1 v prahovanem obraze vetší než 70 % všech px
             souradnice2(1,i) = 1;
            souradnice2(2,i) = 1;
            continue
        end
        pocet2 = 6;     %počáteční nastavení, aby se spustil cyklus while
        dil = 0;
        watchd = 0; 
        while pocet2 > 5    %opakovani, dokud pocet nalezenych bifurkaci presahuje 5
            watchd = watchd + 1;
            if watchd > 20  % brání zacyklení programu
                obr_p = zeros(size(obr_p));
            else
                
                obr_p = bwmorph(obr_p,'bridge');  % homomorfická operace propojování
                obr_p = bwmorph(obr_p,'clean'); % homomorfická operace čištění
                obr_p = bwmorph(obr_p,'thin',30);   % homomorfická operace  ztenčování
                obr_p = bwmorph(obr_p,'dilate',dil);    % homomorfická operace dilatace
                obr_p = bwmorph(obr_p,'thin',10);   % homomorfická operace ztenčování
                obr_p = bwmorph(obr_p,'clean'); % čištění
                dil = dil + 1;  % zvyšení počtu dilatací
            end
        
%             subplot(2,3,6)
%             imshow(obr_p)
%             pause()
            
            [sourad(:,:,i) pocet2] = krizeni3(obr_p); %detekce krizeni
    
        end
%-----------------------------------------------------------------
%% oprava souradnic
        sourad2 = sourad(:,:,i);
        sourad2(:,find (sourad2(1,:)==0))=[]; %smazani pozic 0,0
        if pocet2 == 0 %pokud je pocet nalezenych birkaci=0 -> smaze puvodni bifurkaci
            souradnice2(1,i) = 1;
            souradnice2(2,i) = 1;  
            
        else
            e = zeros(size(sourad2(1,:)));
            for j = 1:length(sourad2(1,:))  %vypocet vzdalenosti nalezenych bifurkaci k puvodni
                e(j) = sqrt((sourad2(1,j)-stred(1))^2 + (sourad2(2,j)-stred(2))^2);
            end
            e2 =sort(e);
            s = find(e==e2(1));       %nalezeni souradnice s nejmensi vzdalenosti od puvodni
            s = s(1);
            if e2(1) > 11   %pokud nejkratsi vzdalenost presahuje 11, nejspise se jedna o chybne nalezenou bifurkaci
                souradnice2(1,i) = 1;
                souradnice2(2,i) = 1;  
            else
                souradnice2(1,i) = souradnice2(1,i) + sourad2(1,s)-stred(1);    %oprava souradnic bifurkace
                souradnice2(2,i) = souradnice2(2,i) + sourad2(2,s)-stred(2);
            end
        end
        %------------------------------------------------------------------
        
    end
    
end
souradnice2(:,find (souradnice2(1,:)==1))=[]; %smazani pozic 1,1
%-------------------------------------------------------------------------
%% odstraneni blizkych pozic
vzdalenost = 12;     %mezni vzdalenost
l=0;
for i = 1:length(souradnice2(1,:))
    stred = souradnice2(:,i);
    e = zeros(size(souradnice2(1,:)));
    for j = i+1:length(souradnice2(1,:))
        if i ~= j
            e(j) = sqrt((souradnice2(1,j)-stred(1))^2 + (souradnice2(2,j)-stred(2))^2); 
        end
    end
    e2 =sort(e);
    for k=1:length(e2)
            if (e2(k)-vzdalenost)<= 0
                    l=l+1;
                    s = find(e==e2(k));           
                    s2(l) = s(1);   %souradnice, ktere jsou blizko k jine
            end
        
    end
end

souradnice3 = souradnice2;
souradnice3(:,s2)=[];

%% odstraneni duplicity
t1 = souradnice3(1,:)*10000 + souradnice3(2,:);
[s i] = unique(t1,'first');
t1 = t1(sort(i));

souradnice3(:,find (souradnice3(1,:)==0))=[];

souradnice3 = zeros(2,length(t1));
souradnice3(1,:) = floor(t1./10000); 
souradnice3(2,:) = t1 - souradnice3(1,:)*10000;
%-------------------------------------------------------------------------
%% krizky
% obr_kriz3 = obr_vysl;
% obr_kriz2 = obr_vysl;
% obr_kriz0 = obr_vysl;
% 
% for i = 1:length(souradnice3)
%     obr_kriz3(souradnice3(1,i)-3:souradnice3(1,i)+3,souradnice3(2,i)) = 255;
%      obr_kriz3(souradnice3(1,i),souradnice3(2,i)-3:souradnice3(2,i)+3 ) = 255;
% end
% for i = 1:length(souradnice3)
%     obr_kriz2(souradnice2(1,i)-3:souradnice2(1,i)+3,souradnice2(2,i)) = 255;
%      obr_kriz2(souradnice2(1,i),souradnice2(2,i)-3:souradnice2(2,i)+3 ) = 255;
% end
% 
% 
% sourad_global(:,find (sourad_global(1,:)==0))=[];
% sourad_global(:,find (sourad_global(2,:)==0))=[];
% for i = 1:length(sourad_global)
%     obr_kriz0(sourad_global(1,i)-3:sourad_global(1,i)+3,sourad_global(2,i)) = 255;
%      obr_kriz0(sourad_global(1,i),sourad_global(2,i)-3:sourad_global(2,i)+3 ) = 255;
% end


% obr_kriz2_r = obr_barevny(:,:,1);
% obr_kriz2_g = obr_barevny(:,:,2);
% obr_kriz2_b = obr_barevny(:,:,3);
% 
% obr_kriz3_r = obr_barevny(:,:,1);
% obr_kriz3_g = obr_barevny(:,:,2);
% obr_kriz3_b = obr_barevny(:,:,3);

% for i = 1:length(souradnice3)
%     obr_kriz3_r(souradnice3(1,i)-3:souradnice3(1,i)+3,souradnice3(2,i)) = 255;
%     obr_kriz3_g(souradnice3(1,i)-3:souradnice3(1,i)+3,souradnice3(2,i)) = 255;
%     obr_kriz3_b(souradnice3(1,i)-3:souradnice3(1,i)+3,souradnice3(2,i)) = 255;
%     
%     obr_kriz3_r(souradnice3(1,i),souradnice3(2,i)-3:souradnice3(2,i)+3 ) = 255;
%     obr_kriz3_g(souradnice3(1,i),souradnice3(2,i)-3:souradnice3(2,i)+3 ) = 255;
%     obr_kriz3_b(souradnice3(1,i),souradnice3(2,i)-3:souradnice3(2,i)+3 ) = 255;
% end

for i = 1:length(souradnice3)
    obr_vysl(souradnice3(1,i)-3:souradnice3(1,i)+3,souradnice3(2,i)) = 255;
    
    obr_vysl(souradnice3(1,i),souradnice3(2,i)-3:souradnice3(2,i)+3 ) = 255;

end

obr_vysl_g = obr_barevny(:,:,2);

for i = 1:length(souradnice3)
    obr_vysl_g(souradnice3(1,i)-3:souradnice3(1,i)+3,souradnice3(2,i)) = 255;
    
    obr_vysl_g(souradnice3(1,i),souradnice3(2,i)-3:souradnice3(2,i)+3 ) = 255;

end

% for i = 1:length(souradnice2)
%     obr_kriz2_r(souradnice2(1,i)-3:souradnice2(1,i)+3,souradnice2(2,i)) = 255;
%     obr_kriz2_g(souradnice2(1,i)-3:souradnice2(1,i)+3,souradnice2(2,i)) = 255;
%     obr_kriz2_b(souradnice2(1,i)-3:souradnice2(1,i)+3,souradnice2(2,i)) = 255;
%     
%     obr_kriz2_r(souradnice2(1,i),souradnice2(2,i)-3:souradnice2(2,i)+3 ) = 255;
%     obr_kriz2_g(souradnice2(1,i),souradnice2(2,i)-3:souradnice2(2,i)+3 ) = 255;
%     obr_kriz2_b(souradnice2(1,i),souradnice2(2,i)-3:souradnice2(2,i)+3 ) = 255;
% end

% obr_kriz3(:,:,1) = obr_kriz3_r;
% obr_kriz3(:,:,2) = obr_kriz3_g;
% obr_kriz3(:,:,3) = obr_kriz3_b;
% 
% obr_kriz2(:,:,1) = obr_kriz2_r;
% obr_kriz2(:,:,2) = obr_kriz2_g;
% obr_kriz2(:,:,3) = obr_kriz2_b;





% figure(9)
% imshow(obr_kriz2)
% title('Pred odstranenim blizkych bifurkaci')



% figure(10)
% imshow(obr_kriz3)
% title('Po odstraneni blizkych bifurkaci')
% 
% figure(11)
% imshow(obr_kriz0)
% title('Bifurkace po globalni detekci')