function [PosWarp0, PosWarp1, Mu,Nu] = BeierNeelyaux(I0,I1,J0,J1,s) ; 
% 
% function [PosWarp0,PosWarp1,Mu,Nu] = BeierNeelyaux(I0,I1,J0,J1,s) ; 
% 
% Calcule, etant donnee deux images I0 et I1 de meme taille la position des pixels reassignes dans 
% l'algorithme pour les images intermediaires Warp0 et Warp1 dans l'algorithme de Beier-Neely de morphing  
% a  partir d'une paire de points (P0,Q0) dans une image I0 et de la paire
% correspondante (P1,Q1) dans une image 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. 

[N1,N2] = size(I0);
% coordonnees de l'origine P0 et de l'extremite Q0 du segment de la
% premiere image 
P0 = J0(1:2);  
Q0 = J0(3:4); 
V0 = Q0-P0; 
vec_t0 = V0/norm(V0,2); 
% coordonnees de l'origne P1 et de l'extremite Q1 du segment de la seconde
% image 
P1 = J1(1:2);
Q1 = J1(3:4);
V1 = Q1-P1;  
vec_t1 = V1/norm(V1,2);  
% coordonnees du vecteur Q_s-P_s (barycentre de V0 et V1)
V = (1-s)*V0 + s*V1; 
NormeV = norm(V,2); 
vec_t = V/NormeV; 

% On calcule le tableau dont les entrees sont les premieres composantes des XX et celle 
% dont les entrees sont les secondes composantes des XX = X- P_s.
X1 = [1:N1]'*ones(1,N2); 
X2 = ones(N1,1)*[1:N2]; 
XX1 = X1 - (1-s)*P0(1) - s*P1(1); 
XX2 = X2 - (1-s)*P0(2) - s*P1(2);
% On calcule des coordonnes mu_s et nu_s des pixels de I
Mu = (vec_t(1)*XX1 + vec_t(2)*XX2)/NormeV; 
Nu = -vec_t(2)*XX1 + vec_t(1)*XX2;
% On declare les tableaux PosWarp0 et PosWarp1 
PosWarp0 = zeros(N1,N2,2); 
PosWarp1 = zeros(N1,N2,2); 
PosWarp0(:,:,1) = P0(1) + Mu*(V0(1)) - Nu*vec_t0(2); 
PosWarp0(:,:,2) = P0(2) + Mu*(V0(2)) + Nu*vec_t0(1); 
PosWarp1(:,:,1) = P1(1) + Mu*(V1(1)) - Nu*vec_t1(2); 
PosWarp1(:,:,2) = P1(2) + Mu*(V1(2)) + Nu*vec_t1(1); 
end 
