close all; clear all; clc;

%% Skript sloužící pro segmentaci obrazu na základě základních segmentačních metod

% Algoritmus slouží nejprve pro výpočet vzorové masky jako průměr 
% z mrtvých buněk v obrazu, následně je na celý snímek aplikována
% mikrostrukturní analýza a dále je snímek procházen a srovnáván se
% vzorovou maskou a na základě toho hledány a označeny mrtvé buňky.

%% První část algoritmu pro určení vzorové masky

% načtení jednoho z 25 snímků
obr = im2double(rgb2gray(imread('1.tif')));

% vykreslení originálního snímku
figure;
imshow(obr, []);
title('Originální snímek')

% zakomentovaná část kódu, která slouží pro vytvoření vzorové masky, na
% základě vytvoření hranic kolem mrtvé buňky

% mrtve = zeros(100,100,4);
% 
% for i = 1:2
%     kruh = roipoly(obr);
% 
%     hranice_kruh = bwmorph(kruh,'remove');
%     figure
%     imshow(hranice_kruh)
% 
%     SE = strel('disk', 3);
%     hranice = imdilate(hranice_kruh,SE); % hranice misto hranice_kruh
%     
%     figure
%     imshow(hranice)
%     
%     [vyrez,sour] = imcrop(hranice_kruh);
%     
%     mrtve(:,:,i) = imresize(vyrez,[100,100]);
% end

% načtení již předem vytvořené matice hranic mrtvých buněk
load('mrtve.mat');

% průěrování hranic mrtvých buněk
mrtve_bunky = (mrtve(:,:,1)+mrtve(:,:,2)+mrtve(:,:,3)+mrtve(:,:,4))/4;

% prahování
mrtve_bunky(mrtve_bunky >= 0.5) = 1;
mrtve_bunky(mrtve_bunky < 0.5) = 0;

% dilatace hranic pro jejich zvýraznění
SE2 = strel('disk', 3);
mrtve_bunky_dilate = imdilate(mrtve_bunky,SE2);

figure
imshow(mrtve_bunky_dilate)
title('Dilatovaná hranice průměru z mrtvých buněk')

% vyplnění kompletní vzorové masky
for i = 1:size(mrtve_bunky_dilate,1)
    radek = mrtve_bunky_dilate(i,:);
    pozice = find(radek==1);
    if isempty(pozice);
        continue;
    end
    for j = 1:length(pozice)-1
        if (pozice(j+1) - pozice(j)) > 1;
            radek(pozice(j):pozice(j+1)) = 1;
        end
    end
    mrtve_bunky_dilate(i,:) = radek;
end

% prahování
mrtve_bunky_dilate(mrtve_bunky_dilate >= 0.5) = 1;
mrtve_bunky_dilate(mrtve_bunky_dilate < 0.5) = 0;

% vykreslení vyplněné zprůměrované mrtvé buňky
figure
imshow(mrtve_bunky_dilate)
title('Vyplněná zprůměrovaná mrtvá buňka')

%% Mikrostrukturni analyza

% v této části algoritmu dochází k aplikování mikrostrukturní analýzy na
% celý původní snímek na základě konvoluce snímku se souborem 9 masek

tic
% načtení masek pro mikrostrukturní analýzu
load('Laws.mat'); 
[m n] = size(obr);

% konvoluce snímku s maskami
for i = 1:9
    konv(:,:,i) = abs(conv2(obr, imrotate(law(:,:,i),180),'same'));
end

% vytvoření parametrické mapy a třídění pomocí k-means
par_mapa = reshape(konv,m*n,9);
trideni = kmeans(abs(par_mapa),2);
zpet = reshape(trideni,m,n);

zpet = mat2gray(zpet);

% ošetření vytvoření negativu při mikrostrukturní analýze
if mode(mode(zpet)) == 1
    zpet = ones(size(zpet)) - zpet;
end

% vykreslení mikrotstrukturní analýzy
figure
imshow(zpet,[])
title('Mikrostrukturni analyza')

% morfologická operace eroze pro odstranění nežádoucích pixelů
SE3 = strel('disk', 1);
obr_erode = imerode(zpet,SE3);

% vykreslení erodovaného obrazu
figure
imshow(obr_erode,[])
title('Erodovaný obraz bez nežádoucích pixelů')

% vložení černého okraje (1px) pro následnou analýzu
obr_okraj = zeros(size(obr_erode,1)+2,size(obr_erode,2)+2);
obr_okraj(2:size(obr_erode,1)+1, 2:size(obr_erode,2)+1) = obr_erode;


%% Procházení obrazu oknem 100x100

% Poslední část algoritmu, kdy je obraz po mikrostrukturní analýze
% procházen nejprve oknem 100x100 px a následně oknem 70x70 px a
% analyzovány různé podmínky pro segmentaci mrtvých buněk na základě
% podobnosti se vzorovou maskou.

% suma pixelů vzorové masky a určení prahu
pixel = sum(sum(mrtve_bunky_dilate));
prah = pixel*0.2;

[a b] = size(mrtve_bunky_dilate);

% určení průměrné velikosti vzorové masky v ose X a v ose Y na základě
% pozice první a poslední jedničky
PrvniX = find(mrtve_bunky_dilate(a/2,:),1,'first');
PosledniX = find(mrtve_bunky_dilate(a/2,:),1,'last');
PrvniY = find(mrtve_bunky_dilate(:,b/2),1,'first');
PosledniY = find(mrtve_bunky_dilate(:,b/2),1,'last');

Prumer_templatu_X = PosledniX - PrvniX;
Prumer_templatu_Y = PosledniY - PrvniY;

% definování výstupu
vystup = zeros(size(obr_okraj,1),size(obr_okraj,2),3);
vystup(:,:,1) = obr_okraj;
vystup(:,:,2) = obr_okraj;
vystup(:,:,3) = obr_okraj;

for i = 1:size(obr_okraj,1)-99
    for j = 1:size(obr_okraj,2)-99     
        okno = obr_okraj(i:(i+a-1),j:(j+b-1));  % vytvoření okna 100x100 px
        
        % ošetření podmínky jestli je suma pixelů okna větší než prah
        % pokud ano pokračuje se dál v cyklu, pokud ne přejde se na další
        % krok (j+1)
        pom = prah(prah<sum(sum(okno)));
        if isempty(pom)
            continue;
        end
        
        % ošetření toho, že v okně se nachází celá samostatná buňka, musí
        % být suma prvního a posledního řádku a prvního a posledního
        % sloupce v okně rovna 0, pokud není přejde se na další krok cyklu
        % (j+1)
        prvni_radek = sum(okno(1,:));
        posledni_radek = sum(okno(end,:));
        prvni_sloupec = sum(okno(:,1));
        posledni_sloupec = sum(okno(:,end));
        
        if prvni_radek > 0 || posledni_radek > 0 || prvni_sloupec > 0 || posledni_sloupec > 0
            continue;
        else
            pixel_okno = sum(sum(okno));
            
            
           % srovnání velikosti buňky v okně s průměrnou velikostí vzorové
           % masky, pokud nejsou splněny podmínky pokračuje se v dalším
           % kroku cyklu (j+1)
            FirstX = find(okno(a/2,:),1,'first');
            LastX = find(okno(a/2,:),1,'last');
            FirstY = find(okno(:,b/2),1,'first');
            LastY = find(okno(:,b/2),1,'last');
            
            Prumer_okna_X = LastX - FirstX;
            Prumer_okna_Y = LastY - FirstY;
            
            if (Prumer_okna_X > Prumer_templatu_X) | (Prumer_okna_Y > Prumer_templatu_Y)
                continue;
            else
                
                % srovnání sumy pixelů v okně s prahem
                [X,Y] = find(okno==1);
                if pixel_okno >= prah
                    for n = 1:length(X)
                        vystup(i+X(n)-1,j+Y(n)-1,1) = 1;  % obarvení mrtvé buňky na červeno
                        vystup(i+X(n)-1,j+Y(n)-1,2) = 0;  
                        vystup(i+X(n)-1,j+Y(n)-1,3) = 0;
                    end
                end
            end
        end
    end
end

figure
imshow(vystup)
title('Výstup po průchodu oknem 100x100 px')

%%  Procházení obrazu oknem 70x70 px

% Postup je stejný jako u procházení obrazu oknem 100x100 px akorát je
% vyřezané okno, ve kterém probíhá analýza, zmenšeno na velikost 70x70 px
% stejně jako vzorová maska. 

% změna velikosti vzrové masky
mrtve_bunky_dilate = imresize(mrtve_bunky_dilate,[70,70]);

pixel = sum(sum(mrtve_bunky_dilate));
prah = (pixel*0.15:pixel*0.2);
[a b] = size(mrtve_bunky_dilate);

PrvniX = find(mrtve_bunky_dilate(a/2,:),1,'first');
PosledniX = find(mrtve_bunky_dilate(a/2,:),1,'last');
PrvniY = find(mrtve_bunky_dilate(:,b/2),1,'first');
PosledniY = find(mrtve_bunky_dilate(:,b/2),1,'last');

Prumer_templatu_X = PosledniX - PrvniX;
Prumer_templatu_Y = PosledniY - PrvniY;

for i = 1:size(obr_okraj,1)-69
    for j = 1:size(obr_okraj,2)-69     
        okno = obr_okraj(i:(i+a-1),j:(j+b-1));
        
        pom = prah(prah<sum(sum(okno)));
        if isempty(pom)
            continue;
        end
        
        prvni_radek = sum(okno(1,:));
        posledni_radek = sum(okno(end,:));
        prvni_sloupec = sum(okno(:,1));
        posledni_sloupec = sum(okno(:,end));

        if prvni_radek > 0 || posledni_radek > 0 || prvni_sloupec > 0 || posledni_sloupec > 0
            continue;
        else
            pixel_okno_pom = prah(prah>=sum(sum(okno)));
            
            FirstX = find(okno(a/2,:),1,'first');
            LastX = find(okno(a/2,:),1,'last');
            FirstY = find(okno(:,b/2),1,'first');
            LastY = find(okno(:,b/2),1,'last');
            
            Prumer_okna_X = LastX - FirstX;
            Prumer_okna_Y = LastY - FirstY;
            
            if (Prumer_okna_X > Prumer_templatu_X) | (Prumer_okna_Y > Prumer_templatu_Y)
                continue;
            else
                [X,Y] = find(okno==1);
                if isempty(pixel_okno_pom)
                    for n = 1:length(X)
                        vystup(i+X(n)-1,j+Y(n)-1,1) = 1;  
                        vystup(i+X(n)-1,j+Y(n)-1,2) = 0;  
                        vystup(i+X(n)-1,j+Y(n)-1,3) = 0;
                    end
                end
            end
        end
    end
end

figure
imshow(vystup)
title('Výstup po průchodu oknem 100x100 a oknem 70x70')
toc