function [Warp0,Warp1,Warp,Mu,Nu] = BeierNeely(I0,I1,J0,J1,s); 

% 
% function [Warp0,Warp1,Warp,Mu,Nu] = BeierNeely(I0,I1,J0,J1,s) ; 
% 
% Calcule, etant donnees deux images I0 et I1 de meme taille, les 
% images transformees par morphing Warp0 et Warp1 dans l'algorithme de Beier-Neely 
% et leur interpolation Warp. Le morphing est ici effectue a partir de la 
% correspondance d'une paire de points (P0,Q0) dans I0 et de la paire
% correspondante (P1,Q1) dans I1. La routine calcule aussi les nouvelles 
% coordonnees Mu et Nu permettant d'exprimer la position du pixel reaffecte. 
% Les matrices J0 et J1 sont des matrices a 1 ligne et 4 colonnes figurant les pixels 
% des points P0, Q0 (pour J0) et ceux des points (P1,Q1) pour J1. Le
% parametre s est a prendre entre 0 et 1.  


% On met l'image dans un cadre de zeros. 
[m,n] = size(I0); 
D=sqrt(m^2+n^2)/2; 
MM = ceil(D-m/2); 
NN = ceil(D-n/2); 
Iinit0= [zeros(MM,n);I0;zeros(MM,n)]; 
Iinit0 = [zeros(m+2*MM,NN) Iinit0 zeros(m+2*MM,NN)];
Iinit1= [zeros(MM,n);I1;zeros(MM,n)]; 
Iinit1 = [zeros(m+2*MM,NN) Iinit1 zeros(m+2*MM,NN)];
[N1,N2] = size(Iinit0); 

% Il faut donc modifier la position des extremites des segments pour les
% considerer cette fois dans Iinit. 
J0init = J0 + [MM,NN,MM,NN]; 
J1init = J1 + [MM,NN,MM,NN]; 

%On fait le calcul des positions des images I0 et I1 reassignees
%(interpolation bilineaire) 
[WarpPos0,WarpPos1,Mu,Nu] = BeierNeelyaux(Iinit0,Iinit1,J0init,J1init,s);
Warp0=zeros(N1,N2);
Warp1=zeros(N1,N2); 
for j=1: N1
    for k=1:N2  
        j0reass = WarpPos0(j,k,1); 
        k0reass = WarpPos0(j,k,2); 
        j1reass = WarpPos1(j,k,1); 
        k1reass = WarpPos1(j,k,2); 
   
        % On remplit Warp0 
        if j0reass*(j0reass-(N1-1))>0 
             Warp0(j,k)=0; 
        elseif k0reass*(k0reass-(N2-1))>0 
                Warp0(j,k)=0; 
        else
            tt=j0reass-floor(j0reass); 
            ss=k0reass-floor(k0reass); 
            indexXinf = floor(j0reass)+1; 
            indexXsup = ceil(j0reass)+1; 
            indexYinf = floor(k0reass)+1; 
            indexYsup = ceil(k0reass) +1; 
            Iaux1 = (1-tt)*Iinit0(indexXinf,indexYinf) + tt*Iinit0(indexXinf,indexYsup); 
            Iaux2 = (1-tt)*Iinit0(indexXsup,indexYinf) + tt*Iinit0(indexXsup,indexYsup); 
            Warp0(j,k) = (1-ss)*Iaux1 + ss*Iaux2;
        end 
        % On remplit Warp1 
        if j1reass*(j1reass-(N1-1))>0 
             Warp1(j,k)=0; 
        elseif k1reass*(k1reass-(N2-1))>0 
                Warp1(j,k)=0; 
        else
            tt=j1reass-floor(j1reass); 
            ss=k1reass-floor(k1reass); 
            indexXinf = floor(j1reass)+1; 
            indexXsup = ceil(j1reass)+1; 
            indexYinf = floor(k1reass)+1; 
            indexYsup = ceil(k1reass) +1; 
            Iaux1 = (1-tt)*Iinit1(indexXinf,indexYinf) + tt*Iinit1(indexXsup,indexYinf); 
            Iaux2 = (1-tt)*Iinit1(indexXsup,indexYinf) + tt*Iinit1(indexXsup,indexYsup); 
            Warp1(j,k) = (1-ss)*Iaux1 + ss*Iaux2;
        end           
    end
end 
Warp = (1-s)*Warp0 + s*Warp1; 
end 
