% % % % %HT%
% Tento program je součástí Diplomové 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)
% tento skript slouží k hledání obočí v termosnímcích.
% program se spouští jako: [eyebrow_position,radius]=HT(img)
%
% Vstup:
% img - segmentovaný snímek
%
% Výstup:
% eyebrow_position - informace o nalezených kružnicích approximujícíc obočí
% radius - polomer nastavených polomerů pro hledání

 function [eyebrow_position,radius]=HT(img)

img_grey = im2double(rgb2gray(img)); % change to type double and grey scale
[M,N] = size(img_grey); % size of image
 wide=0; % prednastavení širky
 help_i=ceil(img_grey);
  [s2,s1]=find(help_i==1); % souřadnice obličeje
  souradnice=[s1,s2];
  
for h=0.4*M:0.6*M % částečné omezení prohledávané oblasti

    row=find(h==souradnice(:,2)); % nalezení aktivních pixelů v prohledávaném řádku 
    if length(row)>=2 % 
        if souradnice(row(end))-souradnice(row(1))>wide % pokud je vetši než dosavadni šířka uloží se souřadnice
            wide=souradnice(row(end))-souradnice(row(1));
            xy=[souradnice(row(1)),souradnice(row(end))]; % jen x,  nepotřebuju
        end
    row=[];
    end
end

%% 
help_i=ceil(im2double(rgb2gray(img))); % prevedeni na biární obraz
i=1;

% hledání výchozích bodu pro nalezení kontury 
while 1<2 
 x1=M+1-i;
y1=min(find(help_i(x1,:)~=0)); % nalezení minima v řádku
y2=max(find(help_i(x1,:)~=0)); % nalezení maxima v řádku
i=1+i;
if isempty(y1)==0
    break
end
end
B = bwtraceboundary(help_i,[x1,y1],'N',8); % nalezení kontury
top=min(B(:,1));
%%
threshold=1;
eyebrow_position=[]; % prednastaveni matice pro ukladaní detekovaných pozc
while size(eyebrow_position,2)<2

% img_edge = mat2gray(edge(img_grey,'sobel',0.025)); %0,15 2 pro canny
 img_edge = mat2gray(edge(img_grey,'canny',[0.08,0.1]*threshold,2)); %0,15 2 pro canny
% img_edge = mat2gray(edge(img_grey,'prewitt',0.025)); %0,15 2 pro canny
% img_edge = mat2gray(edge(img_grey,'log',0.002)); %0,15 2 pro canny
% img_edge = mat2gray(edge(img_grey,'roberts',0.02)); %0,15 2 pro canny
% img_edge = mat2gray(edge(img_grey,'zerocross')); %0,15 2 pro canny

%%
radius = (N/8-15:N/8+10); % velikost polomeru hledanych kruznic
akumulator = zeros(M,N,length(radius)); % prednastaveni prostoru
hough_slice = zeros(M,N); % nastaveni Houghovy roviny pro vypocet

[points_x points_y] = find(img_edge); % hledani pozic hran
circle_x=cos(-pi/2:pi/360:pi/2); % souradnice x jednotkove kruznice pro pozdejsi vypocet
circle_y=sin(-pi/2:pi/360:pi/2); % souradnice y jednotkove kruznice pro pozdejsi vypocet
for r=radius
    for j=1:length(points_x)
        a = round(points_x(j) + r*circle_x); %vypocet kruznic
        b = round(points_y(j) + r*circle_y);
        for i=1:length(circle_x)
            if a(i)>0 && a(i)<M && b(i)>0 && b(i)<N % kotrola zda souradnice bodu na kruznici lezi uvnitr obrazu
                hough_slice(a(i),b(i))=hough_slice(a(i),b(i))+1; % inkrmentujeme hodnotu v akumulatoru
            end
        end
    end
    akumulator(:,:,r)=hough_slice(:,:); % do matice akumulator je ulozen Houghuv prostor pro dany polomer
    hough_slice(:,:)=zeros(M,N); % matici pro konkretni polomer vynulujeme pro dalsi vypocet
end
akumulator = akumulator(:,:,radius); %uprava akumulatoru
hough_space = mat2gray(akumulator); %normalizace Houghova prostoru
% 
% figure(3)
% imshow(hough_space(:,:,1))
%% Vyhodnoceni nalzených maxim
[R C V] = ind2sub(size(hough_space), find(hough_space>0.6)); % nalezeni pozic kde je hodnota px vetsi nez zvoleny prah
max_position = [R C V]'; % ulozeni do matice
for i=1:length(C)
    max_position(4,i) = hough_space(R(i),C(i),V(i));
end

sort_position=(sortrows(max_position',4))'; % serazeni max pozic podle velikosti

s2=size(eyebrow_position,2); % prednastaveni znalosti velikosti matice eyebrow_position
c=0; % konstanta pro regulaci cyklu
i=length(sort_position(1,:)); % pocet opakování cyklu
while i>0
%     if length(sort_position(1~=0)
    if sort_position(1,i)>1.3*(top+radius(sort_position(3,i))) && sort_position(2,i)<0.8*N && sort_position(1,i)-radius(sort_position(3,i))<M/2 && sort_position(1,i)-radius(sort_position(3,i))>0.2*M && s2<2 && sort_position(2,i)>0.3*N% ...
        % sort_position(2,i)>0.2*N
        % podninka zda nalezeny bod lezi v oblasti obrazu, která by mohla
        % obsahovat oboci
        eyebrow_position(:,end+1)=sort_position(:,i); % ulozeni pozice daného bodu
    end
    if s2==2 && (sqrt((eyebrow_position(1,2)-eyebrow_position(1,1)).^2+ (eyebrow_position(2,2)-eyebrow_position(2,1)).^2)<N/8 || abs(eyebrow_position(1,1)-eyebrow_position(1,2))>(radius(eyebrow_position(3,1))/8)||abs(eyebrow_position(3,1)-eyebrow_position(3,2))>length(radius)/3 || abs(eyebrow_position(1,1)-eyebrow_position(1,2))>15 || sort_position(2,i)<0.3*N)% ||  abs(eyebrow_position(2,1)-eyebrow_position(2,2))>60)
        % testování zda nalezené pozice nelezi moc blízko u sebe
        eyebrow_position=eyebrow_position(:,end-1); % v případě že ano danou pozici smaze

    elseif s2==2 && sqrt((eyebrow_position(1,2)-eyebrow_position(1,1)).^2+ (eyebrow_position(2,2)-eyebrow_position(2,1)).^2)>N/8 && abs(eyebrow_position(3,1)-eyebrow_position(3,2))<length(radius)/3 && abs(eyebrow_position(1,1)-eyebrow_position(1,2))<15% && abs(eyebrow_position(2,1)-eyebrow_position(2,2))<60
        % nalezene pozice nelezi moc blizko u sebe, ale zaroven se nachazí
        % proporcionálně na spravnych pozicich
     i=1;
        break % zastaveni
    end
    i=i-1;
    s2=size(eyebrow_position,2); % prezkoumani poctu nalezenych bodu
    if i<=1 && s2<2
        % v připade, že k danemu maximu nebylo nalezeno žádne jiné do páru
        % se nalezený bod smaže a vyřadí se z dalšího prohledávání pomocí
        % navýšení konstanty
        c=c+1;
        eyebrow_position=[];
        i=length(sort_position(1,:))-c;
    elseif c==length(sort_position(1,:)) % v případě že nebyly nalezeny žádné pozice, které odpovídají nastaveným podmínkám
        disp('neuspesna detekce')
        break
    end
end
  threshold=threshold-threshold/10; % snížení prahu
end

%% vykreslení
% for g=1:length(eyebrow_position(1,:))
%     img(eyebrow_position(1,g),eyebrow_position(2,g),:)=0;
%     pom1=eyebrow_position(1,g);
%     pom2=eyebrow_position(2,g);
%     if pom1<M-2 && pom1>2
%         if pom2<N-2 && pom2>2
%             img(eyebrow_position(1,g)-1,eyebrow_position(2,g),:)=0;
%             img(eyebrow_position(1,g)-2,eyebrow_position(2,g),:)=0;
%             img(eyebrow_position(1,g)+1,eyebrow_position(2,g),:)=0;
%             img(eyebrow_position(1,g)+2,eyebrow_position(2,g),:)=0;
%             img(eyebrow_position(1,g),eyebrow_position(2,g)+1,:)=0;
%             img(eyebrow_position(1,g),eyebrow_position(2,g)+2,:)=0;
%             img(eyebrow_position(1,g),eyebrow_position(2,g)-1,:)=0;
%             img(eyebrow_position(1,g),eyebrow_position(2,g)-2,:)=0;
%         end
%     end
%     
%     a = round(eyebrow_position(1,g) - radius(eyebrow_position(3,g))*circle_x); % vypocet kruznic
%     b = round(eyebrow_position(2,g) - radius(eyebrow_position(3,g))*circle_y);
%     for i=1:length(circle_x)
%         if a(i)>0 && a(i)<M && b(i)>0 && b(i)<N % pokud souradnice bodu na kruznici lezi uvnitr obrazu
%             img(a(i),b(i),:)=0; % inkrmentujeme hodnotu v akumulatoru
%         end
%     end
% end

% figure(4)
% imshow(img)

