function [Iinit0,Iinit1,Warp] = MultiBeierNeely(I0,I1,J0,J1,a,b,p,s); 
% 
% function [Iinit0,Iinit1,Warp] = MultiBeierNeely(I0,I1,J0,J1,a,b,p,s);
%
% Calcule l'image Warp obtenue par morphing entre deux images I0 et I1 a partir 
% des correspondances entre plusieurs segments, respectivement entre les
% points [(j_k0,k_k0) (jj_k0,kk_k0)] et les points [(j_k1,k_k1) (jj_k1,kk_k1)]; 
% les coordonnees de ces points en correspondance sont figurees comme les
% lignes des tableaux J0 et J1. Les parametres a,b,p sont a prendre
% respectivement dans ]0,1], [1/2,2], [0,1]. Les deux images I0 et I1 ont 
% ete au prealable transformees en les images encadrees Iinit0 et Iinit1.

[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); 
X1 = [1:N1]'*ones(1,N2); 
X2 = ones(N1,1)*[1:N2]; 

K = size(J0,1); 
Alpha = zeros(N1,N2); 
WARP  = zeros(N1,N2);  
for k=1:K
    % Declaration des segments de la k-eme correspondance 
    Pk0 = J0(k,1:2)+[MM,NN]; 
    Qk0 = J0(k,3:4)+[MM,NN]; 
    Pk1 = J1(k,1:2)+[MM,NN];  
    Qk1 = J1(k,3:4)+[MM,NN];  
    Vk  = (1-s)*(Qk0-Pk0) + s*(Qk1-Pk1); 
    Lk = norm(Vk,2); 
    % Calcul du Warpk et de Muk et Nuk pour la k-eme correspondance 
    [Warpk0,Warpk1,Warpk,Muk,Nuk] = BeierNeely(I0,I1,J0(k,:),J1(k,:),s); 
    % Calcul de la matrice des distances entre les pixels x et le point Pks
    % (Pks :=(1-s)*Pk0 + s*Pk1) 
    XX1 = X1 - (1-s)*Pk0(1) - s*Pk1(1); 
    XX2 = X2 - (1-s)*Pk0(2) - s*Pk1(2);
    distancesPks =(XX1.^2 + XX2.^2).^(1/2); 
    
    %Calcule de la matrice des distances entre les pixels x et le point Qks 
    % (Qks :=(1-s)*Qk0 + s*Qk1) 
    XXX1 = X1 - (1-s)*Qk0(1) - s*Qk1(1); 
    XXX2 = X2 - (1-s)*Qk0(2) - s*Qk1(2);
    distancesQks =(XXX1.^2 + XXX2.^2).^(1/2); 
    
    % Construction de la matrices des distances des pixels au segment 
    % [Pks,Qks]
    Dk=zeros(N1,N2); 
    for jj=1:N1; 
        for kk=1:N2; 
            if Muk(jj,kk)>1 
                Dk(jj,kk) = distancesQks(jj,kk); 
            elseif Muk(jj,kk)<0 
                Dk(jj,kk) = distancesPks(jj,kk); 
            else 
                Dk(jj,kk)= abs(Nuk(jj,kk)); 
            end 
        end 
    end 
    
    % Calcul de la matrices des coefficients ponderateurs correspondant a
    % la k-eme paire de segments en correspondance
    alpha = (Lk^p./(a+Dk)).^b; 
    
    % Incrementation de la somme des coefficients ponderateurs 
    Alpha = Alpha +  alpha;  
    
    % Incrementation de la matrice WARP 
    WARP = WARP + alpha.*Warpk;  
end 

% Calcul du Warp obtenu comme moyenne ponderee des Warpk avec les
% coefficients Alphak.  
    Warp = WARP./Alpha ;     
end 