%% image read and prep
image = imread(""); %Path ke snimku
s = size(image);
%vytvoření pracovních matic a vektorů
imin = zeros(s(1),s(2)); %image do double array pro zpracování
imout = imin;
imageout = image; %výstupní image
bot = zeros(1,s(2)); %vector pro ukládání horní hrany
top = zeros(1,s(2)); %vector pro ukládání dolní hrany

%load vstupní data do matic
for i=1:s(1)
    for j=1:s(2)
        imin(i,j)=image(i,j);
        imout(i,j)=imin(i,j);
        imageout(i,j)=255;
    end
end

% Processing
limdark = 0; %limit černá
limbright = 255; %limit - bílá
tolerance = 1.35; %Nastavení tolerance doporučeno 1,1 - 1,5
flip = false; %Přechod od detekce dolní hrany k horní
row = 1; %row - řádek
col = 1; %column - sloupec

for j=1:s(2) %nový sloupec
    col = j;
    row = s(1);
    if(imin(row,col) ~= 0)
        limdark = imin(row,col)*tolerance;
    else
        limdark = 20; %pokud je na první pozici černá.
    end

    
    for i=1:s(1) %nový řádek
        
        if (i > s(1)-1) %kontrola dojetí sloupce bez detekce
            flip = false;
            top(col)=0;
            bot(col)=0;
            break
        end
        
        if ((imin(row,col) < limdark) && (flip == false))
            %imout(row,col) = 0;

        elseif (flip == false)
            flip = true;
            bot(col) = s(1)-row;
            limbright = imin(row-1,col)/tolerance;
            
        end

        if((imin(row,col) < limbright) && (flip == true))
            flip = false;
            top(col)=s(1)-row;
            %imout(row,col) = 0;
            break
        
        elseif (flip==true)
            %imout(row+1,col) = 255;
        end


        row = row-1;
    end
end

figure;
tiledlayout(4,1);
nexttile;
plot(bot);
hold on;
plot(top);

% %odstranění podezřelých hodnot
avbot = floor(mean(bot));
for i=1:s(2)
    if (top(i) < avbot) %pokud je snímek nakřivo, je lepší tuto if() zakomentovat
        top(i) = 0;
    end
    if (bot(i) <= 2)
        bot(i) = avbot; %interpolate bot
    end
end


%lineární aproximace a doplnění pro top 
for i=1:s(2)
    forward = 0;
    forward_index = 0;
    back = 0;
    back_index = 0;
    if ((top(i) == 0) && (i ~= 1)) %doplnění

        for j=1:i-1
            if (i+j) > s(2) %Ochrana proti překročení velk. vektoru
                top(i) = ceil(mean(top));
                break
            end
            if ((top(i-j) ~= 0) && (back_index == 0))
                back = top(i-j);
                back_index = i-j;
            end
            if ((top(i+j) ~= 0) && (forward_index == 0))
                forward = top(i+j);
                forward_index = i+j;
            end
            if ((forward ~= 0) && (back ~= 0))

               top(i) = ceil(back+(i-back_index)*((forward-back)/(forward_index-back_index)));
               if (isnan(top(i)))
                   top(i)=ceil(top(i));
               end

                break
            end
        end
    elseif (i ~= 1)
        %NOP
    else 
        top(i) = ceil(mean(top));

    end
end

%Filtrování - vyhlazení 
tophamp = hampel(top,5);
bothamp = hampel(bot,5);

%výsledky do výstupní image
for i=1:s(1)
    for j=1:s(2)
        a=floor(s(1)-tophamp(j)+1);
        b=floor(s(1)-bothamp(j)+1);
        imout(a,j) = 255;
        imout(b,j) = 0;
        imageout(i,j)=imout(i,j);
    end
end

thickness = mean(top) - mean(bot);
%thickness = mean(top(1:end)) - mean(bot(1:end));%Možné vyřazení špatných výsledků

%Vykreslení a výpis

nexttile;
plot(bothamp);
hold on;
plot(top,"--r")
plot(tophamp);
nexttile;
imshow(image);
nexttile;
imshow(imageout);

scale = 95.23; %měřítko nm/pixel.
%47,39 pro 2048 snímky//95,23 pro 1024 snímky
disp("Tloušťka (nm):");
disp(thickness*scale); 
