%                   DIPLOMOVÁ PRÁCE
%                                   ADAM BŘÍNEK

% pro použití tohoto programu stačí zadat název programu v místě kde se
% soubor otevírá, poté se projde soubory vybrané ve for cyklu. Pro ukazku
% je vybrán 1 řez. Vykresléné obrazky ukazují hranici obrázku, hranici
% kontaktu a původní obrázek.
%%
clear all
close all

tic                                       % Měření doby výpočtu

zub1=[];                                    % Proměná do které se bude ukládat obrázek
M = 1;

prah2 = 0.6;                               % Ruční volba prahu určující výplňový materiál
prah3 = 0.37;                               % Ruční volba prahu určující vzduchové bubliny

Maska1 = [-1 -1 -1;-1 9 -1; -1 -1 -1]./9;   % Ostřící maska
row=755;  col=915;                          % Velikost obrazu

celypovrchA = [];
castecnypovrchA = [];
celypovrchB = [];
castecnypovrchB = [];
poz1s = [(col*1/3), (row/2)];
poz2s = [(col*2/3), (row/2)];
%%
% Načítání dat řez po řezu
for k = 1000:1:1000
    if k >0 && k < 10
        l = ['000' num2str(k)];
    elseif k > 9 && k < 100
        l = ['00' num2str(k)];
    elseif k > 99 && k < 1000
        l = ['0' num2str(k)];
    else
        l = num2str(k);
    end
    
% otevirani souboru
    fin=fopen(['final', l, '.raw'],'r');
    I=fread(fin, [col row],'uint16=>uint16'); 
    zub=im2double(I');
    I = [];
    fclose(fin);

    % v matlabu 2016 se může segmentovat pomocí prahy = multithresh(zub,4);
    
% histogram našeho jednoho obrázku
% figure
[x,h] = hist(zub(:),512);
% bar(h(2:end),x(2:end))

    % probehne pouze v pripade, ze je v kanalku vyplnujici material
    if sum(x(1,200:end)) > 100
        zub = zub.^0.8;
        Vysledek1 = conv2(zub,Maska1,'same');
        zub = zub + Vysledek1;

% figure
% [x,h] = hist(zub(:),512);
% bar(h(2:end),x(2:end))
        
% prahovani obrazu
        zub_bin = zub>prah2;
        
%         figure
%         imshow(zub_bin)

% Uzavření binárního obrazu      
        kol = imclose(zub_bin,strel('disk',30,4));
%         kol = imerode(kol,strel('disk',2,0));

        kol2 = kol.*zub;
        
% zjistovani okoli vyplnoveho materialu
        
        I2 = imdilate(kol,strel('disk',8,0));
        I2 = ~I2.*zub;      
        konecna = zub - I2 - kol2;
        
% figure
% [x,h] = hist(konecna(:),512);
% bar(h(2:end),x(2:end))

% prahování oblasti vzduchovách bublin
        konecna(konecna>prah3)=0;
        konecna(konecna>0)=1;
        konecna = imclose(konecna,strel('disk',8,0));
        konecna = medfilt2(konecna,[2 2]);

        konecna(konecna>0)=1;               % konecna je oblast kde je vedle okoli vyplnoveho materialu bublina
        
%figure
%imshow(konecna)

        konecna_bin = konecna + kol;
        pom1 = imclose(konecna_bin,strel('disk',30,4));
       
        pom1 = pom1.*zub;
        pom1 = zub - pom1;
        pom1 = zub - pom1;
%         pom1 = medfilt2(pom1,[5 5]);
        pom1(pom1>0)=1;                     % Pom1 je obraz segmentovaného zubního kanálku       

% urceni centroidů pro rozdělení dat na dvě části, podle zubních kanálků
        Ilabel = bwlabel(pom1);
        s  = regionprops(Ilabel, 'centroid');
        centroids = cat(1, s.Centroid);

% figure
% imshow(pom1,[])
% hold on
% plot(centroids(:,1), centroids(:,2), 'r*')
% hold off

% určení kolik míst se nachází v binárním obraze zubního kanálku.
        [pom2, ~] = size(s);
        if pom2 == 1
            [poz11] = s.Centroid;
        elseif pom2 == 2
            [poz1, poz2] = s.Centroid;
        else 
            pom2 = 2;
            poz1 = poz1s;
            poz2 = poz2s;
        end
        
        %% pokud jsou dva kanalky, tak se udelaji dva vzorky
        if (pom2 == 2) && (abs(poz1(1,2)-poz2(1,2)) < 75)
            
            if poz1(1,1) > poz2(1,1)            
                poz = poz1;
                poz1 = poz2;
                poz2 = poz;
            end
            
% rozdeleni obrazku pro dva vzorky ze vzdalenosti obou kanalku
            pom3 = (round(abs((poz1(1,1) - poz2(1,1))/2) + poz1(1,1)));
           
            pom1A = pom1;
            pom1A(:,pom3:end) = 0;
            
% potvrzení, že se v pom1A nachází zubní kanálek
            pom1A = jednotaS(pom1A,poz1,row,col);   
            
% vytvoření masky pro aplikaci do VG studia 
%             o = ['maskaA', num2str(k),'.tif'];
%             imwrite(pom1A,o); 
            
            konecnaA = konecna;
            konecnaA(:,pom3:end) = 0;
            
            pom1B = pom1;
            pom1B(:,1:pom3) = 0;
            
% potvrzení, že se v pom1B nachází zubní kanálek
            pom1B = jednotaS(pom1B,poz2,row,col);
            
% vytvoření masky pro aplikaci do VG studia             
%             o = ['maskaB', num2str(k),'.tif'];
%             imwrite(pom1B,o); 
            
            konecnaB = konecna;
            konecnaB(:,1:pom3) = 0;
            
% výpočet hranice kanálku A v proměné výsledekA a hranice kontaktu v porměné obA
            [obA, vysledekA] = hranice(pom1A, konecnaA, row, col);

% zapamatování oblasti kanálku A pro spodní oblast zubu
            pom1As = pom1A;
            
% výpočet hranice kanálku B v proměné výsledekB a hranice kontaktu v porměné obB
            [obB, vysledekB] = hranice(pom1B, konecnaB, row, col);
            
% zapamatování oblasti kanálku B pro spodní oblast zubu
            pom1Bs = pom1B;            
            
% výpočet délky hranice kanálku A a 
            [uA,~] = find(vysledekA);
            [wA,~] = find(obA);
            celypovrchA = [celypovrchA, length(uA)];
            castecnypovrchA = [castecnypovrchA, length(wA)];
            
            [uB,~] = find(vysledekB);
            [wB,~] = find(obB);
            celypovrchB = [celypovrchB, length(uB)];
            castecnypovrchB = [castecnypovrchB, length(wB)];
            poz1s = poz1;
            poz2s = poz2;
        elseif pom2 == 1 
            
        % rozdeleni podle vzorku A nebo B
            N11 = abs(poz11(1,1) - poz1s(1,1));
            N12 = abs(poz11(1,2) - poz1s(1,2));
            N21 = abs(poz11(1,1) - poz2s(1,1));
            N22 = abs(poz11(1,2) - poz2s(1,2));
            %% v oblasti kde se v jednom z kanálků nenachází výplňový materiál
            if (N11 < 75) && (N12 < 75)

% V kanálku A se nachází výplňový materiál, provede se výpočet jako v předchozí části 
                pom1A = pom1;
                pom1A(:,pom3:end) = 0;
                pom1A = jednotaS(pom1A,poz11,row,col);
                
%                 o = ['maskaA', num2str(k),'.tif'];
%                 imwrite(pom1A,o); 
                
                konecnaA = konecna;
                konecnaA(:,pom3:end) = 0;
                       
                [obA, vysledekA] = hranice(pom1A, konecnaA, row, col);
                [uA,~] = find(vysledekA);
                [wA,~] = find(obA);
                
                celypovrchA = [celypovrchA, length(uA)];
                castecnypovrchA = [castecnypovrchA, length(wA)];

% V kanálku B se nenachází výplňový materiál
% zvolíme oblast z předešlého řezu, kde se v kanálku nacházel výplňový
% materiál
                pom2B = imdilate(pom1Bs,strel('disk',16,4));
                
                pom2B = pom2B.*zub;

% provedeme prahování oblasti
                pom2B(pom2B > prah2) = 1;
                pom2B(pom2B < prah3 )= 1;
                pom2B(pom2B ~= 1) = 0;
                pom2B = imclose(pom2B,strel('disk',16,4));

%                 
                pom2B = jednotaS(pom2B,poz2,row,col);
                
%                 o = ['maskaB', num2str(k),'.tif'];
%                 imwrite(pom1Bs,o); 
                
                pom1Bs = pom2B;
                
                [vysledekB] = hranice1(pom2B, row, col);
                obB = vysledekB;
                [uB,vB] = find(vysledekB);
                celypovrchB = [celypovrchB, length(uB)];
                castecnypovrchB = [castecnypovrchB, 0];
                
                poz1s = poz11;
            elseif (N21 < 75) && (N22 < 75)
                %%  předešlý výpočet jako v předchozí části, ale pro oblast kanálku A           
                pom1B = pom1;
                pom1B(:,1:pom3) = 0;
                pom1B = jednotaS(pom1B,poz11,row,col);
                
                o = ['maskaB', num2str(k),'.tif'];
%                 imwrite(pom1B,o); 
                 
                konecnaB = konecna;
                konecnaB(:,1:end) = 0;
                       
                [obBA, vysledekB] = hranice(pom1B, konecnaB, row, col);
                [uB,vB] = find(vysledekB);
                [wB,zB] = find(obB);
                
                celypovrchB = [celypovrchB, length(uB)];
                castecnypovrchB = [castecnypovrchB, length(wB)];
                
                pom2A = imdilate(pom1As,strel('disk',16,4));
                pom2A = pom2A.*zub;
                pom2A(pom2A > prah2) = 1;
                pom2A(pom2A < prah3 )= 1;
                pom2A(pom2A ~= 1) = 0;
                pom2A = imclose(pom2A,strel('disk',16,4));
                
                pom2A = jednotaS(pom2A,poz1,row,col);
                
                o = ['maskaA', num2str(k),'.tif'];
%                 imwrite(pom2A,o); 
                
                pom1As = pom2A;
                
                [vysledekA] = hranice1(pom2A, row, col);
                obA = vysledekA;
                [uA,vA] = find(vysledekA);
                celypovrchA = [celypovrchA, length(uA)];
                castecnypovrchA = [castecnypovrchA, 0];
                
                poz2s(1,1) = poz3(1,1);
            
            else %% v oblasti se nenašli centroidy na pozici zubních kanálků
                        % využili jsme během tvoření algoritmu jako
                        % kontrolu jeho práce
                         
                A = 'Chyba rozdeleno do vzorku';
                disp(A)
                
%                 pom1 = zeros(row,col);
   
%                 o = ['maskaA', num2str(k),'.tif'];
%                 imwrite(pom1,o);
%         
%                  o = ['maskaB', num2str(k),'.tif'];
%                  imwrite(pom1,o);
            end
        
        end


    else %% pro případ, že se v obraze nenachází zubní kanálek
        
        A = ['neni kanal' ,num2str(M)];
        disp(A)
        
%         
        celypovrchA = [celypovrchA, 0];
        castecnypovrchA = [castecnypovrchA, 0];
        
        celypovrchB = [celypovrchB, 0];
        castecnypovrchB = [castecnypovrchB, 0];

% tvorba masek pro VG studio 

        pom1 = zeros(row,col);
       
%         o = ['maskaA', num2str(k),'.tif'];
%         imwrite(pom1,o);
        
%         o = ['maskaB', num2str(k),'.tif'];
%         imwrite(pom1,o);
        
    end

% moje(:,:,M) = zub + ob(:,:,M);
M = M+1;

% nulování hodnot
Vysledek2=[]; zub_bin=[]; kol=[]; poh1=[]; poh2=[]; poh3=[]; 
poh4=[]; poh5=[]; poh6=[]; poh7=[]; poh8=[];
u=[];v=[];w=[];z=[];
poh11=[]; poh12=[]; poh13=[]; poh14=[]; poh15=[]; poh16=[]; 
poh17=[]; poh18=[]; vysl=[];kol2=[];I2=[];moje=[];ob=[];vysledek=[];
se=[];konecna=[];konecna_bin=[];prom=[];pom=[];Xx=[];Yy=[];zub1=[];pom1=[];
pom1A=[];pom1B=[];konecnaA=[];konecnaB=[];
% obA=[];obB=[];
pom1A=[];pom1B=[];
% vysledekA=[];vysledekB=[];
% o = ['hranice', num2str(k),'.tif'];
% imwrite((zub + vysledek),o);    
end

cas = toc;                                  % konec časové analýzy algoritmu
% Původní obrázek
figure
imshow(zub)

% obazek 
figure
imshow(zub + vysledekA + vysledekB,[])

figure
imshow(zub + obA + obB,[])

% výpočet celkového procentuální kontaktu pro jednotlivé vzorky 
povrchA = sum(castecnypovrchA)/sum(celypovrchA);
A = ['Procentuální kotakt vzorek A' ,num2str(povrchA)];
disp(A)

povrchB = sum(castecnypovrchB)/sum(celypovrchB);
B = ['Procentuální kontakt vzorek B' ,num2str(povrchB)];
disp(B)





