%Programa: Phase Only Correlation
%Descripción: Calcula el desplazamiento óptimo entre un par de imágenes 
%usando la transformada de Fourier 
%Autor: FJHL
%Fecha Inicial: 31/Agosto/2015
%Fecha Actualización: 02/Septiembre/2015

clear all;
close all;

%leer imagen
I_k0=imread('../images/Triangulo_k0.bmp');
I_k1=imread('../images/Triangulo_k1.bmp');

[nFilas,nColumnas,channels]=size(I_k0);
if(channels>1)
    I_k0=double(rgb2gray(I_k0));
    I_k1=double(rgb2gray(I_k1));
else
    I_k0=double(I_k0);
    I_k1=double(I_k1);
end

%Transformar la I_k0 para obtener una I_k1
%Crear la matriz de trasformación
%Affine=maketform('affine',[1 0 0;0 1 0; 15 5 1]);
%Aplicar la matriz de transformación a la imagen I_k0
%I_k1=imtransform(I_k0,Affine,'XData',[1 size(I_k0,2)],'YData',[1 size(I_k0,1)]);
%con rotación
theta=5;
Affine=maketform('affine',[cosd(theta) -sind(theta) 0;sind(theta) cosd(theta) 0; 15 5 1]);
I_k1=imtransform(I_k0,Affine,'XData',[1 size(I_k0,2)],'YData',[1 size(I_k0,1)]);


r=POC(I_k0,I_k1);
valor_maximo=max(r(:));
[peak_row,peak_col]=find(r==valor_maximo);

%Considerando que los desplazamientos calculados no pueden ser mayores que 
%la mitad del tamaño de la imagen, entonces
if(peak_row>nFilas/2),
    dx2_min=nFilas-peak_row +1;
else
    dx2_min=-peak_row+1;
end
if(peak_col>nColumnas/2),
    dx1_min=nColumnas-peak_col +1;
else
    dx1_min=-peak_col+1;
end
%óptimo via POC
dx1_min
dx2_min

residual_d0=(I_k0-I_k1).^2;%Con Desplazamiento d=(0,0)

%Calcular el error al aplicar d_min
%Transformar la imagen I_k0 con la d_min
%Crear la matriz de trasformación
Affine=maketform('affine',[1 0 0;0 1 0; dx1_min dx2_min 1]);
%Aplicar la matriz de transformación a la imagen I_k0
I_new=imtransform(I_k0,Affine,'XData',[1 size(I_k0,2)],'YData',[1 size(I_k0,1)]);

residual_doptimo=(I_k1-I_new).^2;

figure;
imagesc(I_k0);
colormap gray;
figure;
imagesc(I_k1);
colormap gray;
figure;
mesh(r);
figure;
imagesc(I_new);
colormap gray;

figure;
imagesc(residual_d0);
colormap gray;
figure;
imagesc(residual_doptimo);
colormap gray;





