% Tato funkce slouzi k slicovani dvou snimku.
%
% Vstupni promenne:
%   snimek1 - referencni, snimek2 - registrovany, hann=1 - eliminace
%   okrajoveho efektu Hannovym oknem, pocet iteraci - cislo
%
% Potrebne m-fily:
%   rozliseni.m, iluminace.m, transform.m

function [snimek1, snimek2] = fazova_korelace (snimek1, snimek2, hann, pocet_iteraci)

%% ÚPRAVA NA STEJNOU VELIKOST SNIMKU
[snimek1,snimek2]=rozliseni(snimek1,snimek2);
snimek_puvodni=snimek2;

%% INICIALIZACE PARAMETRU A GENEROVANI VEKTORU

[a,b,e]=size(snimek1);
y_stred=floor(a/2+1);       %pro prevod do log-polarnich souradnic   
x_stred=floor(b/2+1);       %pro prevod do log-polarnich souradnic
f=size(snimek2,3);          %pro vybrani vhodne slozky
matice=[1 0 0; 0 1 0; 0 0 1];   %inicializace transformacni matice
krok=0.25;                  %krok pro vektor theta
zaklad_logaritmu=1.01;      %zaklad logaritmu pro prevod do log. souradnic
theta=0:krok:180;           %generovani vektoru theta

%generování vektoru rho
i=1;                                    
pomocny=1;                              
while (pomocny<=(floor(min(a, b)/2))-2) 
    log_rho(i)=pomocny;                 
    pomocny=zaklad_logaritmu^i;         
    i=i+1;
end

%vahovani polarniho spektra
for i=1:length(log_rho)
    vaha(i)=i/length(log_rho);
end

%okenkova funkce pro amplitudova spektra v log-polarnich souradnicich
for i=1:length(theta)
    for j=1:length(log_rho)
        H1(i,j)=0.5*(1-cos(2*pi()*i/(a-1)))*0.5*(1-cos(2*pi()*j/(b-1)));
    end
end


%% -----------------------------ITERACE------------------------------------

for iterace=1:pocet_iteraci
    
% vybrani vhodne slozky
if e==3
    f1=(snimek1(:,:,2));
    f1_okoli=f1;
else
    f1=(snimek1(:,:,1));
    f1_okoli=f1;
end

if f==3
    f2=(snimek2(:,:,2));
    f2_okoli=f2;
else
    f2=(snimek2(:,:,1));
    f2_okoli=f2;
end

%uprava nerovnomerne iluminace
[f1, maska1]=iluminace(f1);
[f2, maska2]=iluminace(f2);

f1=double(f1);
f2=double(f2);

%zvyrazneni hran
H2=[-1 -1 -1 ;-1 9 -1 ;-1 -1 -1];
f1=conv2(f1,H2,'same');
f2=conv2(f2,H2,'same');

%eliminace vlivu masky
maska1=bwmorph(maska1,'erode',10);
maska2=bwmorph(maska2,'erode',10);
f1(maska1==0)=mean(mean(f1(maska1~=0)));
stredni=mean(mean(f2(maska2~=0)));
f2(maska2==0)=stredni;

%eliminace okrajoveho efektu
okraj=ones(a,b);
okraj=bwmorph(okraj,'erode',10);
if std(double(f1_okoli(okraj==0)))>1 || hann==1
    f1=okno(f1);
end

if std(double(f2_okoli(okraj==0)))>1 || hann==1
    f2=okno(f2);
end

%% rotace + meritko
    
    F1=fft2(f1);                                                                       
    F1_shift=fftshift(abs(fft2(f1)));               
    F2_shift=fftshift(abs(fft2(f2)));               
    
    %prevod amplitudoveho spektra do polarnich souradnic
    F1_pol=zeros(length(theta),length(log_rho));    %inicializace matice pro amplitudove spektrum v log-pol. s.
    F2_pol=zeros(length(theta),length(log_rho));    %inicializace matice pro amplitudove spektrum v log-pol. s.
    for l=1:length(log_rho)
        for m=1:length(theta)
            x=log_rho(l)*cos(deg2rad(theta(m)));    %vypocet souradnice x v originalu
            y=log_rho(l)*sin(deg2rad(theta(m)));    %vypocet souradnice y v originalu
            x_dolu=floor(x);                        
            y_dolu=floor(y);                        
            x_nahoru=x_dolu+1;                      
            y_nahoru=y_dolu+1;
            F1_pol(m,l)=(F1_shift(y_stred+x_dolu,x_stred+y_dolu)*(x_nahoru-x)*(y_nahoru-y)+F1_shift(y_stred+x_nahoru,x_stred+y_dolu)*(x-x_dolu)*(y_nahoru-y)+F1_shift(y_stred+x_dolu,x_stred+y_nahoru)*(x_nahoru-x)*(y-y_dolu)+F1_shift(y_stred+x_nahoru,x_stred+y_nahoru)*(x-x_dolu)*(y-y_dolu))/((x_nahoru-x_dolu)*(y_nahoru-y_dolu))*vaha(l);    % pouziti bilinearni transformace pro vypocet hodnoty v log-polarnich souradnicich
            F2_pol(m,l)=(F2_shift(y_stred+x_dolu,x_stred+y_dolu)*(x_nahoru-x)*(y_nahoru-y)+F2_shift(y_stred+x_nahoru,x_stred+y_dolu)*(x-x_dolu)*(y_nahoru-y)+F2_shift(y_stred+x_dolu,x_stred+y_nahoru)*(x_nahoru-x)*(y-y_dolu)+F2_shift(y_stred+x_nahoru,x_stred+y_nahoru)*(x-x_dolu)*(y-y_dolu))/((x_nahoru-x_dolu)*(y_nahoru-y_dolu))*vaha(l);    % pouziti bilinearni transformace pro vypocet hodnoty v log-polarnich souradnicich            
        end
    end
        
    spektrum_F1_pol=fft2(F1_pol.*H1);       %spektrum vahovaneho amplitudoveho spektra v polarnich souradnicich
    spektrum_F2_pol=fft2(F2_pol.*H1);       %spektrum vahovaneho amplitudoveho spektra v polarnich souradnicich

    %fazova korelace
    H3=[1 1 1; 1 2 1; 1 1 1]/10;
    vypocet=conv2(ifft2(spektrum_F1_pol.*conj(spektrum_F2_pol)./abs(spektrum_F1_pol.*spektrum_F2_pol)),H3,'same');
    [rotace,log_meritko]=find(vypocet==max(max(vypocet)));  %zjisteni souradnic maxima

    %zjisteni spravneho meritka
    
    meritko1=zaklad_logaritmu^(length(log_rho)-(log_meritko-1));
    meritko2=zaklad_logaritmu^(-(log_meritko-1));
    rotace1=(rotace-1)*krok;
    rotace2=rotace1+180;
    
    if log_meritko==1 && (rotace1==0 || rotace1==360 || rotace1 == 180)
        
    vypocet=conv2(ifft2(F1.*conj(fft2(f2))./abs(F1.*fft2(f2))),H3,'same');    
    
    else
    matice_pomocna=zeros(3,3,4);
    f2_pomocny=zeros(a,b,4);
    vypocet_pomocny=zeros(a,b,4);
    maximum=zeros(1,4);
    
    matice_pomocna(:,:,1)=[1/meritko1 0 0; 0 1/meritko1 0; 0 0 1]*[1 0 b/2; 0 1 a/2; 0 0 1]*[cos(pi*rotace1/180) -sin(pi*rotace1/180) 0; sin(pi*rotace1/180) cos(pi*rotace1/180) 0; 0 0 1]*[1 0 -b/2; 0 1 -a/2; 0 0 1];
    matice_pomocna(:,:,2)=[1/meritko1 0 0; 0 1/meritko1 0; 0 0 1]*[1 0 b/2; 0 1 a/2; 0 0 1]*[cos(pi*rotace2/180) -sin(pi*rotace2/180) 0; sin(pi*rotace2/180) cos(pi*rotace2/180) 0; 0 0 1]*[1 0 -b/2; 0 1 -a/2; 0 0 1];
    matice_pomocna(:,:,3)=[1/meritko2 0 0; 0 1/meritko2 0; 0 0 1]*[1 0 b/2; 0 1 a/2; 0 0 1]*[cos(pi*rotace1/180) -sin(pi*rotace1/180) 0; sin(pi*rotace1/180) cos(pi*rotace1/180) 0; 0 0 1]*[1 0 -b/2; 0 1 -a/2; 0 0 1];
    matice_pomocna(:,:,4)=[1/meritko2 0 0; 0 1/meritko2 0; 0 0 1]*[1 0 b/2; 0 1 a/2; 0 0 1]*[cos(pi*rotace2/180) -sin(pi*rotace2/180) 0; sin(pi*rotace2/180) cos(pi*rotace2/180) 0; 0 0 1]*[1 0 -b/2; 0 1 -a/2; 0 0 1];
    
    f2_pomocny(:,:,1)=transform(f2, matice_pomocna(:,:,1), 'linear', stredni);
    f2_pomocny(:,:,2)=transform(f2, matice_pomocna(:,:,2), 'linear', stredni);
    f2_pomocny(:,:,3)=transform(f2, matice_pomocna(:,:,3), 'linear', stredni);
    f2_pomocny(:,:,4)=transform(f2, matice_pomocna(:,:,4), 'linear', stredni);
    
    vypocet_pomocny(:,:,1)=conv2(ifft2(F1.*conj(fft2(f2_pomocny(:,:,1)))./abs(F1.*fft2(f2_pomocny(:,:,1)))),H3,'same');
    vypocet_pomocny(:,:,2)=conv2(ifft2(F1.*conj(fft2(f2_pomocny(:,:,2)))./abs(F1.*fft2(f2_pomocny(:,:,2)))),H3,'same');
    vypocet_pomocny(:,:,3)=conv2(ifft2(F1.*conj(fft2(f2_pomocny(:,:,3)))./abs(F1.*fft2(f2_pomocny(:,:,3)))),H3,'same');
    vypocet_pomocny(:,:,4)=conv2(ifft2(F1.*conj(fft2(f2_pomocny(:,:,4)))./abs(F1.*fft2(f2_pomocny(:,:,4)))),H3,'same');
    
    maximum(1)=abs(max(max(vypocet_pomocny(:,:,1))));
    maximum(2)=abs(max(max(vypocet_pomocny(:,:,2))));
    maximum(3)=abs(max(max(vypocet_pomocny(:,:,3))));
    maximum(4)=abs(max(max(vypocet_pomocny(:,:,4))));
    
    matice=matice*matice_pomocna(:,:,find(maximum==max(maximum),1));
    f2=f2_pomocny(:,:,find(maximum==max(maximum),1));
    vypocet=vypocet_pomocny(:,:,find(maximum==max(maximum),1));
    
    end

%% posun
  
    [m,n]=find(vypocet==max(max(vypocet)));     %zjisteni souradnice maxima        
    if (m<a/2) 
        dy=-m+1;
    else
        dy=a-m;
    end
    
    if (n<b/2)
        dx=-n+1;
    else
        dx=b-n;
    end
    
    matice1=[1 0 dx; 0 1 dy; 0 0 1];
    matice=matice*matice1;        
    snimek2=transform(snimek_puvodni, matice, 'linear', 0);
end

snimek1=uint8(snimek1);
snimek2=uint8(snimek2);

end
