% Tento program je součástí Semestrální práce: 
% Metody pro zpracování termovizních snímků s detekcí stanovené oblasti obličeje          
% Autor: Dana Kolářová
% ID: 124347
% Škola: VUT Brno
% Fakulta: FEKT                                                                    
% Obor: Biomedicínská technika a bioinformatika
%************************************************
% Aplikace je spustitelná v prostředi Matlab.
% Použitá verze : Matlab 7.9.0 (R2009b)
% tento skript slouží k segmentaci obličeje pomocí k-means shlukování.
% Program se volá jako:
% [X,Y,line,polovina1,polovina2,img_newsize,hranice]=cluster_analysis(he,Newpath1)
%
% Vstup:
% he - segmentovaný snímek
% Newpath1 - cesta ke snímku
%
% Výstup:
% [X,Y] - souradnice segmentované oblasti
% line - souřadnice obličejové osy
% polovina1 - souradnice pp rozdelene osou
% polovina2 - souradnice lp obliceje
% img_newsize - segmentovaný snímek v menším rozlišení
% hranice - souradnice kontury segmentované oblasti

function [X,Y,line,polovina1,polovina2,img_newsize,hranice]=cluster_analysis(he,Newpath1)
 c=[];
 he=imread([Newpath1,he]); % nastení snímku
img_grey=im2double(rgb2gray(he)); % převod na šedotonový obraz 
entropy=entropyfilt(img_grey);  % výpočet lokální entropie
cform = makecform('srgb2lab'); 
lab_he = applycform(he,cform); % převod do barevného prostoru L*a*b*
ab = double(lab_he(:,:,2:3)); % převod na foramát double 
[nrows,ncols] = size(ab(:,:,1));  % výpečet rozměru obrázku
ab = reshape(ab,nrows*ncols,2); % přeskupení matice zpět do formátu mxn
en= reshape(entropy,nrows*ncols,1); % přidání informace o entropii k jednotlivým pixelům
ab(:,3)=en; % umocnění entropie na třetí, ab hodnoty měly stejnou váhu jako umístění barev v barevných osách
A=zeros(7,3,3); % nastavení matice pro nastavení počátečních centroidu shluku

A2=[191	147	5;169	72	3;173	175	5;148	91	3;149	192	4;163	99	4;128	205	4]; % nastavení souřadnic centroidů -  na základě předchozího tetsovaní
A3=[191	147	5;169	72	3;173	165	5;148	91	3;149	172	4;153	108	4;128	205	4];
A4=[192.26	132.17	4.77;179.02	67.90	3.02;171.04	195.74	4.87;161.03	75.41	3.39;118.57	188.17	3.76;128.80	127.31	3.72;130.80	209.45	4.47];
A(:,:,1)=A4;
A(:,:,2)=A2;
A(:,:,3)=A3;

nColors = 7; % definice počtu shluků
% % nastavení shlukování, euklidovská vzdálenost, počet shluků,opakování shlukování 4 čtyřikrát abychom se vyhnuli lokálním minimům
[cluster_idx cluster_center] = kmeans(ab,nColors,'distance','sqEuclidean', ...
    'emptyaction','singleton','Replicates',3,'Start',A);

pixel_labels = reshape(cluster_idx,nrows,ncols); % přeskupení matice vystupu shlukovaní zpět do velikosti obrázku mxn
segmented_images = cell(1,3); % uložení pixelů z jednotlivých shluků do oddělaených obrázků
rgb_label = repmat(pixel_labels,[1 1 3]); 
%%
for k = 1:nColors   % vykreslení
    color = he;
    color(rgb_label ~= k) = 0;
    segmented_images{k} = color;
end

%% nalezení shluků představujících obličej
number=[];
for ind=1:nColors 
    x=find((pixel_labels==ind)); % najdi počet pixelů v jednotlivých shlucích
    num=length(x);
    if num>ncols*nrows/nColors % pokud je více než 10000 ( hodnota stanovena empiricky) uloží se číslo shluku
        number(2,end+1)=num;
        number(1,end)=ind;
    end      
end
dist=[[1:nColors]' , sqrt((nrows/2-cluster_center(:,1)).^2+(ncols/2-cluster_center(:,2)).^2)]; % výpočet vzdálenosti centroidů od středu obrázku
dist = sortrows(dist,2); % seřazení těchto vzdáleností
to=[dist(1:3,1),cluster_center(dist(1:3,1))]; % nalezení 3 shluků s nejmenší vzdáleností od středu
if ismember(to(1,1),number)==1 % pokud shluk s nejmenší vzdáleností má více než 10 000 pixelů je označen jako obličejová část
top=to(1,1);
elseif ismember(to(2,1),number)==1
    top=to(2,1); % pokud ne, je jako nejlepší označet clust s druhou nejmenší vzdáleností ( jelikož v prvním mohou byýt pouze oční koutky)
else c=1; 
    top=to(1,1);
end

dist2=[[1:nColors]' , sqrt((cluster_center(top,1)-cluster_center(:,1)).^2+(cluster_center(top,2)-cluster_center(:,2)).^2)]; % k nejlepšímu shluku se hledá ten, který má centroid nejblíže
cluster=find(dist2(:,2)<26); % rozdíl musí být menší než 26( opět stanoveno empiricky)
if isempty(c)==0 
    cluster=to(:,1);
end
C=cluster_center;
C(:,end+1)=1:length(C(:,1));
C(cluster,2)=0;
C = sortrows(C,2);
if C(end,end)==top
Cen=find(C(:,2)<cluster_center(top,2));
contoure=segmented_images{C(Cen(end),end)};
C=C(Cen(end),end);
else C=C(end,end);
end

contoure=segmented_images{C};

%% vykresleni a spojení obličejových shluků

 x=[]; 
 y=[];
for l=1:length(cluster)
 [x1,y1]=find((segmented_images{cluster(l)}~=0)); % postupne nalezení souřadnici ve shlucích reprezentující obličej
 x=[x;x1]; 
 y=[y;y1];
end

img=uint8(zeros(nrows,ncols,3)); % Vytvoření proměnné pro vytvoření výslednho obrazu
 for p=1:length(x)
     img(x(p),y(p))=he(x(p),y(p)); % vykreslení pixelů
 end
img_segmented=img;
%% oddělání tech kousíčku nespojitých s hlavní obličejovou částí
img(end+1,:,:)=0;
help_i=ceil(im2double(rgb2gray(img))); % převod na binární obraz
L=bwlabel(help_i,4); % označení bílých segmentů
STATS = regionprops(L,'Area'); % výpočet plochy segmentů
[~,m]=max([STATS(:,1).Area]); % nalezení největší ploch
[x,y]=find(L~=m); % nalezení indexů malých ploch

for p=1:length(x)  % začernění nalezených pixelů
     img(x(p),y(p),:)=0;
     help_i(x(p),y(p))=0;
end
img(end,:,:)=[];

 %% horni slupka + kontura toho dosavadniho
 i=1;
 help_i=ceil(im2double(rgb2gray(img))); % převedení na binární obraz
while 1<2                               % nekonečný cyklus
 x1=nrows+1-i;                          % postupuje se od spodu...
y1=min([find(help_i(x1,:)~=0)]);        % pokud jsou nenulové pixely v prohledávané řadě najde se minimum 
y2=max([find(help_i(x1,:)~=0)]);        % .. a maximum
if isempty(y1)==0                       % pokud byly body nalezeny
    break                               % cyklus se zastaví
end
i=i+1;
end
B = bwtraceboundary(help_i,[x1,y1],'N',8); % nalezení kontury
%%

 [x,y,obrazek,koo]=contoure_track(contoure,B,[x1,y1,y2],help_i); %kontura sloučených shluku-B, contoure =svrchni shluk

for p=1:length(x)
     img(x(p),y(p))=he(x(p),y(p)); % vykreslení pixelů
end
  hranice=bwtraceboundary(ceil(im2double(rgb2gray(img))),[x1,y1],'N',8);
%% oddělání tech kousíčku
img(end+1,:,:)=0;
help_i=ceil(im2double(rgb2gray(img)));      % převod na binární obraz
L=bwlabel(help_i,4);                        % označení bílých segmentů
STATS = regionprops(L,'Area');              % výpočet plochy segmentů
[~,m]=max([STATS(:,1).Area]);               % nalezení největší ploch
[x,y]=find(L~=m);                           % nalezení indexů malých ploch
for p=1:length(x)                           % začernění nalezených pixelů
     img(x(p),y(p),:)=0;
     help_i(x(p),y(p))=0;
end
x=find(help_i(nrows,:)==1);
if ~isempty(x);
img(end,:,:)=[];
img(end,x(1):x(end),:)=he(end,x(1):x(end),:);
help_i=ceil(im2double(rgb2gray(img))); % binární obraz
end

%% vyplnění vnitřních oblastí
help_i2=1-help_i; % převrácený binární obraz
L=bwlabel(help_i2,4); % nalezení bilých ploch
STATS = regionprops(L,'Area','PixelList','Centroid');  % vypočet velikosti ploch
[~,m]=max([STATS(:,1).Area]); % nalezení maxima
[x,y]=find(L~=m); % největší je pozadí, hledá se vše kromě největší plochy
for p=1:length(x)
     img(x(p),y(p),:)=he(x(p),y(p),:); % vyplnění chybějícíhc částí
end
[img_final,line]=koutky(img,STATS,img_segmented); % nalezení očí nosu úst osy hlavy

%% vykleslení
polovina1=img_final;  % nalezení 1. poloviny podle osy
for i=1:nrows
 polovina1(i,line(i,2)+1:end,:)=0;  % nalezení 1. poloviny obličeje podle osy
end
polovina2=img_final;    % nalezení 2. poloviny obličeje podle osy
for i=1:nrows
 polovina2(i,1:line(i,2),:)=0;  % nalezení 2. poloviny podle osy
end
h=waitbar(1/1);         % poslední waitbar
 set(findobj(h,'type','patch'), ...  % změna barvy waitbaru
            'edgecolor','b','facecolor','b')
 img_newsize=imresize(ceil(im2double(rgb2gray(img_final))),[150,200]); % změna velikosti na velikost metadat       % změna velikosti, velikost orig. snímků z termokamery
 polovina1_newsize=imresize(ceil(im2double(rgb2gray(polovina1))),[150,200]);  % změna velikosti, velikost orig. snímků z termokamery
 polovina2_newsize=imresize(ceil(im2double(rgb2gray(polovina2))),[150,200]); % změna velikosti, velikost orig. snímků z termokamery
% figure(3),imshow(img_newsize)
[X,Y]=find(polovina1_newsize);         % nalezení pixelů představujících segmentovanou polovinu obličeje
polovina1=[X,Y];                       % uložení pixelů do matice
[X,Y]=find(polovina2_newsize);         % nalezení pixelů představujících segmentovanou polovinu obličeje
polovina2=[X,Y];                       % uložení pixelů do matice
% hodnoty=1-(img_newsize);
% img_newsize(end-5:end,:)=0;
% figure(6);imshow(img_newsize)
[X,Y]=find(img_newsize);              % nalezení pixelů představující segmentovanou část
close(h)                              % zavření waitbaru
