function [PosWarp0, PosWarp1, Mu,Nu] = BeierNeelyaux(I0,I1,J0,J1,s); 
[N1,N2] = size(I0);
//positionnement des correspondances  
P0 = J0(1:2);  
Q0 = J0(3:4); 
V0 = Q0-P0; 
vec_t0 = V0/norm(V0,2); 
P1 = J1(1:2);
Q1 = J1(3:4);
V1 = Q1-P1;  
vec_t1 = V1/norm(V1,2);  
//calcul de t(s)
V = (1-s)*V0 + s*V1; 
NormeV = norm(V,2); 
vec_t = V/NormeV; 
//calcul des tableaux des composantes des xx 
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);
//calcul des tableaux des mu_s(x) et nu_s(x) 
Mu = (vec_t(1)*XX1 + vec_t(2)*XX2)/NormeV; 
Nu = -vec_t(2)*XX1 + vec_t(1)*XX2;
//Remplissage des 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); 
endfunction  
