%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Appariement de deux images - G. Besnard
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
close all
nbpoints=15; % nombre de points qu'on veut cliquer
choix=0; %%% 0 : cliquer les points, 1 : faire seulement la transformation sur un set de points déjà cliqués
dir='/Users/florentmathieu/Travail/Jonjonc/Essai_1/'; % Endroit où sont les images
nomg='refg';
nomd='refd';
correc=2;%1 rigide, 2 affine et 3 parabolique

%%%%%%%%%%%%%%%%%%%%%%%%%%%


img = double(imread(sprintf('%s%s.tif',dir,nomg))); % Image de référence
imd = double(imread(sprintf('%s%s.tif',dir,nomd))); % Image à corriger
imdd=imd;
imcor=sprintf('%scor',nomd); % Nom de l'image corrigée pour sauvegarde
D=100;
dec=10;
%%%%%%%%%%%%%%%%%%%%%%%%
if choix==0
    cx21g=zeros(nbpoints,1);
    cy21g=zeros(nbpoints,1);
    cx21g(1,1)=NaN;
    cy21g(1,1)=NaN;
    [cx21g,cy21g]=clique_points(nbpoints,img,cx21g,cy21g,D);
    close all
    affichage_image(img,cx21g,cy21g,dec,nbpoints);
    %%%%droite%%%%
    cx21d=zeros(nbpoints,1);
    cy21d=zeros(nbpoints,1);
    cx21d(1,1)=NaN;
    cy21d(1,1)=NaN;
    [cx21d,cy21d]=clique_points(nbpoints,imd,cx21d,cy21d,D);
    %%sauvegarde%%
    fid = fopen('points_transformation_gd.txt','wt');
    for ni=1:nbpoints
        fprintf(fid,'%8.3f %8.3f %8.3f %8.3f\n',cx21g(ni),cy21g(ni),cx21d(ni),cy21d(ni));
    end
    fclose(fid);
    X=cx21g;
    Y=cy21g;
    xx=cx21d;
    yy=cy21d;
    [haut long] = size(img);
    [m n] = meshgrid(1:long,1:haut);
    %%%estimation de la transformation%%%
    [M2,N2]=transformation_image(X,Y,xx,yy,img,imd,nbpoints,m,n,haut,long,correc);
    imdtrans=interp2(m,n,imdd,M2,N2);
    dep_long=M2-m;
    dep_haut=N2-n;
else
    load points_transformation_gd.txt;
    cx21g=points_transformation_gd(:,1);
    cy21g=points_transformation_gd(:,2);
    cx21d=points_transformation_gd(:,3);
    cy21d=points_transformation_gd(:,4);
    X=cx21g;
    Y=cy21g;
    xx=cx21d;
    yy=cy21d;
    [haut long] = size(img);
    [m n] = meshgrid(1:long,1:haut);
    %%%estimation de la transformation%%%
    [M2,N2]=transformation_image(X,Y,xx,yy,img,imd,nbpoints,m,n,haut,long,correc);
    imdtrans=interp2(m,n,imdd,M2,N2);
    dep_long=M2-m;
    dep_haut=N2-n;
    %%affichage final%%
end
figure
subplot(222)
imagesc(imd)
hold on
plot(cx21d,cy21d,'+r','markersize',10,'linewidth',2)
colormap(gray)
subplot(224)
imagesc(imdtrans)
hold on
plot(cx21g,cy21g,'+r','markersize',10,'linewidth',2)
colormap(gray)
subplot(221)
imagesc(img)
hold on
plot(cx21g,cy21g,'+r','markersize',10,'linewidth',2)
subplot(223)
imagesc(img-imdtrans)
imwrite(uint16(imdtrans),sprintf('%s%s.tif',dir,imcor),'tif');
save(sprintf('%sdeplongi.mat',dir),'dep_long')
save(sprintf('%sdeptrans.mat',dir),'dep_haut')