% % % % % elastix_raw_data
% Tento program je součástí diplomové práce 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é inženýrství a bioinformatika
%************************************************
% Aplikace je spustitelná v prostředi Matlab.
% Použitá verze : Matlab 7.9.0 (R2009b)
% Tato funkce slouží k  identifikování a vybrání oblastí očí nosu úst a
% dalších oblastí (př. vousatá brada).
%
% Tato funkce se volá jako:
% [img_final,line]=koutky(img,STATS,img_segmented)
%
% Vstupy:
% img - segmentovaný obrázek celkový
% img_segmented - snímek s chybějícími oblastmi,
% STATS -informace o techto chybějících oblastech
%
% Výstupy: 
% img_final  - segmentovaný snímek z vybranými částmi
% line - osa obličeje

function [img_final,line]=koutky(img,STATS,img_segmented)

[nrows,ncols,c]=size(img); % zjištění velikosti snímku
[eyebrow_position,radius]=HT(img); % nelezení očí
eyebrow_position=[sortrows([eyebrow_position]',2)]'; % seřadit tváře pravá a levá
points1=sort([eyebrow_position(1,1),eyebrow_position(1,2)]);
point2=points1(1)-radius(max(eyebrow_position(3,:))); % nalzenení pozic pro výřez oblastí
img0=im2double(rgb2gray(img)); % převod na šedotónový snímek
[xx,yy]=find(img0==0); % nalezení černých pixelů
for i=1:length(xx) 
    img0(xx(i),yy(i),:)=0.8; % převod na šedé aby nezasahovali do dalších detekcí
end

img0([1:point2+5-subplus((max(eyebrow_position(3,:))-10)),points1(2)+1:end],:)=0; % výřez oblasti na základě detekovaného obočí
prah=[0.98,0.98]; % nastavení prahů
middle_point=mean([eyebrow_position(2,1),eyebrow_position(2,2)]); % nalezení středu oblasti

%% levá strana
ll=length(eyebrow_position(2,1)-radius(eyebrow_position(3,1)):middle_point); % delka
img1b=img0(:,eyebrow_position(2,1)-radius(eyebrow_position(3,1))+round(ll)/2:eyebrow_position(2,1)-radius(eyebrow_position(3,1))+round(ll));% výpočet oblasti s vnějším koutkem
img1a=img0(:,eyebrow_position(2,1)-(radius(eyebrow_position(3,1))):eyebrow_position(2,1)-radius(eyebrow_position(3,1))+round(ll/2)); % výpočet oblasti s vnitřním koutkem

img1a(:,1)=0;
img1b(:,end)=0;

while 1<2
    [x11,y12]=find(img1b>=prah(1)*max(max(img1b))); % nalezení pixelů nad prahem v 1. oblasti
    [x21,y22]=find(img1a>=prah(2)*max(max(img1a))); % nalezení pixelů nad prahem v 2. oblasti
    %% vykreslení
    for i=1:length(x11)
        img1b(x11(i),y12(i),:)=0;
    end
    for i=1:length(x21)
        img1a(x21(i),y22(i),:)=0;
    end

    help_i=1-ceil(img1b); % hledání černých oblastí v 1 vyřezu
    help_i2=1-ceil(img1a); % hledání černých oblastí v 2 vyřezu
    STATS = regionprops(logical(help_i),'Area','PixelList','Centroid'); % výpočet plochy segmentů
    STATS2 = regionprops(logical(help_i2),'Area','PixelList','Centroid'); % výpočet plochy segmentů
    
    if (length(STATS)==1||length(STATS2)==1) && img0(1,1)==0 % pokud se nenašli 2 oblasti s nejvyšší teplotou
        [~,a]=min([length(STATS),length(STATS2)]);
        prah(a)=prah(a)-0.0005; % snížení prahu
        
        if prah(1)<=0.95 || prah(2)<=0.95 % pokud už jepráh moc malý,a přeto se nenašla vhodná oblast, znamená to že sledovaná oblast je moc velká a splynula s okolím, nutno ohraničit šedým polema  hledat znovu
            img0([1:point2+8,points1(2)+1:end],:)=0.8; % šedý rámec
            prah=[0.98,0.98]; % nastavení původního prahu
            middle_point=mean([eyebrow_position(2,1),eyebrow_position(2,2)]); % vypočtení středového bodu
            
            %% levá strana
            % znovu provedení výřezu
            ll=length(eyebrow_position(2,1)-radius(eyebrow_position(3,1)):middle_point);
            img1b=img0(:,eyebrow_position(2,1)-radius(eyebrow_position(3,1))+round(ll)/2:eyebrow_position(2,1)-radius(eyebrow_position(3,1))+round(ll));
            img1a=img0(:,eyebrow_position(2,1)-radius(eyebrow_position(3,1)):eyebrow_position(2,1)-radius(eyebrow_position(3,1))+round(ll/2));
            
            img1a(1,:)=0;
            img1b(:,end)=0;
            
        end
        
    else break
    end
end
img0(1,1)=0;
prah(:)=0.98;
[~,m]=max([STATS(:,1).Area]); % výpočet největší oblast (okolí)
STATS(m,1).Area=0;  % vynulováni
[~,m]=max([STATS(:,1).Area]); % skutečně nejětší oblast

[~,m1]=max([STATS2(:,1).Area]); % výpočet největší oblast (okolí)
STATS2(m1,1).Area=0; % vynulováni
[~,m1]=max([STATS2(:,1).Area]); % skutečně nejětší oblast

points=[STATS(m,1).Centroid;STATS2(m1,1).Centroid]; % nalezení centroidů
points(1,1)=points(1,1)+eyebrow_position(2,1)-radius(eyebrow_position(3,1))+round(ll)/2; % vypočtení skutečných souřadnic v původním obraze
points(2,1)=points(2,1)+eyebrow_position(2,1)-radius(eyebrow_position(3,1))-1; % vypočtení skutečných souřadnic v původním obraze

bod1a=(round(max(STATS(m,1).PixelList(:,1)))+eyebrow_position(2,1)-radius(eyebrow_position(3,1))+round(ll)/2+1); % nalezení souřadnic středu oka
bod1=[bod1a,floor(mean([min(STATS2(m1,1).PixelList(:,2)),max(STATS(m,1).PixelList(:,2))]))]; % nalezení souřadnic středu oka
% img(round(mean(points(:,2))),round(mean(points(:,1))),:)=[0,255,0]; % vykreslení

%% pravé oko -  analogicky
ll=length(middle_point:eyebrow_position(2,2)+2*radius(eyebrow_position(3,2))/3);
img2a=img0(:,middle_point:middle_point+round(ll/2));
img2b=img0(:,middle_point+round(ll/2)+1:eyebrow_position(2,2)+radius(eyebrow_position(3,2))/2);

img2b(:,1)=0;
img2a(:,1)=0;

while 1<2
    [x31,y32]=find(img2b>=prah(1)*max(max(img2b)));
    [x41,y42]=find(img2a>=prah(2)*max(max(img2a)));
    for i=1:length(x31)
        img2b(x31(i),y32(i),:)=0;
    end
    for i=1:length(x41)
        img2a(x41(i),y42(i),:)=0;
    end

    help_i=1-ceil(img2b);
    help_i2=1-ceil(img2a);
    STATS = regionprops(logical(help_i),'Area','Centroid','PixelList'); % výpočet plochy segmentů
    STATS2 = regionprops(logical(help_i2),'Area','Centroid','PixelList'); % výpočet plochy segmentů
    
    if (length(STATS)==1||length(STATS2)==1) && img0(1,1)==0 % pokud se nenašli 2 oblasti s nejvyšší teplotou
        [~,a]=min([length(STATS),length(STATS2)]);
        prah(a)=prah(a)-0.001;
        
        if prah(1)<=0.97 || prah(2)<=0.97
            img0([1:point2+8,points1(2)+1:end],:)=0.8;
            prah(:)=0.98;
            ll=length(eyebrow_position(2,1)-radius(eyebrow_position(3,1)):middle_point);
            img2a=img0(:,middle_point:middle_point+round(ll/2));
            img2b=img0(:,middle_point+round(ll/2)+1:eyebrow_position(2,2)+radius(eyebrow_position(3,2))/2);
            
            img2a(:,1)=0;
            img2b(:,end)=0;
            
        end
        
    else break
    end
end

[~,m]=max([STATS(:,1).Area]); % nalezení největší ploch
STATS(m,1).Area=0;
[~,m]=max([STATS(:,1).Area]);
[~,m1]=max([STATS2(:,1).Area]); % nalezení největší ploch
STATS2(m1,1).Area=0;
[~,m1]=max([STATS2(:,1).Area]);

bod2a=(floor(min(STATS2(m1,1).PixelList(:,1)))+middle_point);
bod2=[bod2a,floor(mean([min(STATS2(m1,1).PixelList(:,2)),max(STATS(m,1).PixelList(:,2))]))];
points=[points;STATS(m,1).Centroid;STATS2(m1,1).Centroid];
points(3,1)=points(3,1)+middle_point+round(ll/2);
points(4,1)=points(4,1)+middle_point;
point=round([mean(points(1:2,:));mean(points(3:4,:))]);
img(round(mean(points(3:4,2))),round(mean(points(3:4,1))),:)=[0,1,0];

%% vypocet primky aprox oci
p1=polyfit(point(:,1),point(:,2),1); % proložení bodů představujících oči
y=polyval(p1,1:ncols); % výpočet y-ových hodnot
line1=y; % osa appox oči horizontálně

%% vypocet kolmice

x=[bod1(1),bod2(1)]; % xove souřadnice
xm0 = mean(x); % nalezení bdu na kořenu nosu
shift=round(abs(x(2)-x(1))/2*0.8); % rozšíření oblasti ve které se budou hledat potenciálne chladnější oblasti
y=round(polyval(p1,[bod1(1),bod2(1)])); % výpočet y-obých hodnot
if diff(y)==0 % pokud neni rozdíl mezi y-vými souřadnciemi jsou konstanty nulové
    PLm=0;
    PLc=round(xm0);
else
    m = ((diff(y))/diff(x)); % jiank se vypočíta sklon
    ym0 = mean(point(:,2)); % a průměrná hodnota
    PLm = tan( atan(m)+pi/2); % otočení
    PLc = ym0 - PLm*xm0; % nová rovnice přímky
end

p1(2)=PLc; 
p1(1)=PLm;
y=polyval(p1,1:nrows); % výpošet nových souřadnic přímky
[~,x]=find(y>0);
[~,x1]=find(y<nrows); % nepočítá se pro souřadnice mimo rozlišení snímku, je nutné ouřadnice intepolovat do celých čísel - jednotlivých pixelů
x=[x(1):x1(end)];
if p1(1)~=0
    p1=polyfit([y(x(1)),y(x(2))],[x(1),x(2)],1); % výběr dvou bodu 
    y=polyval(p1,[1:nrows]); % výpočet již celočíselných souřadnic
end

for i=1:length(y)
    img(i,round(y(i)),:)=0; % zobrazení osy
end
line=[[1:i]',round(y)'];

%% hledání oblastí nosu a úst
save=[];
y_v=unique(round(y)); % hledají se neopakující se souřadnice
help_i2=1-ceil(im2double(rgb2gray(img_segmented))); % převrácený binární obraz
L=bwlabel(help_i2,4); % nalezení bilých ploch
STATS2 = regionprops(L,'Area','PixelList','Centroid');  % vypočet velikosti ploch
[~,m]=max([STATS2(:,1).Area]); % nalezení maxima
STATS2(m,1).Area=0;
STATS2(m,1).PixelList=[0,0];
x=[];y=[];
%% celá oblast
for ii=1:length(STATS2) % postupné procházení oblasti podél nalezené osy  v rozšířeé oblasti +-shift
    j=1;
    while j<STATS2(ii,1).Area
        for i=1:length(line(:,1))
            if (line(i,2)-shift)<=STATS2(ii,1).PixelList(j,1) && (line(i,2)+shift)>=STATS2(ii,1).PixelList(j,1)  && line(i,1)==STATS2(ii,1).PixelList(j,2) % pokud osa protíná potenciální oblast s jinou teplotou
                save=[save,ii]; % uložení
                x=[x;STATS2(ii,1).PixelList(:,1)]; % zařazení souřadnic
                y=[y;STATS2(ii,1).PixelList(:,2)];
                j=STATS2(ii,1).Area;
                break
            end
        end
        j=j+1;
    end
end

%%
for i=1:length(x) % odstranění chladných oblastí
    img(y(i),x(i),:)=0;
end
x=[];y=[];
%% elipsa levé oko
R1=1.3*(abs(point(1,1)-bod1(1))); % výpočet hlavní poloosy jednoho oka
R2=(point(1,2)-(eyebrow_position(1,1)-radius(eyebrow_position(3,1)))); % výpočet vedlejší poloosy jednoho oka
R3=1.3*(abs(point(2,1)-bod2(1))); % výpočet hlavní poloosy druhého oka
R4=(point(2,2)-(eyebrow_position(1,2)-radius(eyebrow_position(3,2)))); % výpočet vedlejší poloosy druhého oka
if abs(R2-R4)>=3 % oči by mely být symetické takže pokud se liší...
    [~,m]=max([R2,R4]); % veme se vyšší vlikost vedlejíš poloosy -  odchylka zpusobená HT obočí -  obočí je více asymetrické
    if m==1 % upravení menší poloosy
        R4=R2;
    else R2=R4;
    end
end
if abs(R1-R3)>=3 % oči by mely být symetické takže pokud se liší...
    [~,m]=max([R1,R3]); % veme se vyšší vlikost vedlejíš poloosy -  odchylka zpusobená HT obočí -  obočí je více asymetrické
    if m==1 % upraví se menší poolosa
        R3=R1;
    else R1=R3;
    end
end
if R1/R2>2 || R3/R4>2 % pokud je moc velký pomer mezi velikosti jednotlivých poolos
    R2=R2*1.5; % tak se upraví osy a zvýší se
    R4=R2;
end
if mean([R1,R3])<mean([R2,R4]) % pokud je větší vedlejší poloosa tak se musí zvětšit
    R1=R1+abs(mean([R1,R3])-mean([R2,R4]));
    R3=R3+abs(mean([R1,R3])-mean([R2,R4]));
end

cc=1:ncols;
rr=[1:nrows]';
cy1=round(bod1(1)-4*R1/5); % y-ova souradnice
cx1=line1(cy1)-R2/3; % x ove souřadnice

f=@(xx,yy) ((xx-cx1).^2)./(R2^2)+((yy-cy1).^2)./(R1^2)<1; % vytvoření elipsy
Mask = bsxfun(f,rr,cc);

[e1,e2]=find(Mask); % nalezení souřadnic elipsy
for i=1:length(e1) % vykreslení elipsy
    img(e1(i),e2(i),:)=0;
end

%% elipsa pravé oko
cy2=round(bod2(1)+4*R3/5); % vypočet y-ove osy
cx2=line1(cy2)-R4/3;  % vypčet x-ove osy
f=@(xx,yy) ((xx-cx2).^2)./(R4^2)+((yy-cy2).^2)./(R3^2)<1; % vytvoření elips o daných rozměrech
Mask = bsxfun(f,rr,cc);

[e1,e2]=find(Mask); % nalezení souřadnic elipsy
for i=1:length(e1)
    img(e1(i),e2(i),:)=0; % zamaskování pozic ve snímku
end

for i=1:length(x)
    img(y(i),x(i),:)=0;
end
img_final=img;
img_final(end-10:end,:)=0; % zkrácení snímku o 10 px ze zdola, kvuli deform. poli lícování
