function [Iinit0,Iinit1,Warp] = MultiBeierNeely(I0,I1,J0,J1,a,b,p,s); 
[m,n] = size(I0); 
//encadrement prealable des tableaux 
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)];
//initialisation des tableaux Alpha (poids) et Warp(tableaux)
[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
//calcul des tableaux  Warpk, Muk,Nuk 
    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);  
    [Warpk0,Warpk1,Warpk,Muk,Nuk] = BeierNeely(I0,I1,J0(k,:),J1(k,:),s); 
//calcul du tableau des distance des pixels x a P_ks
    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); 
//calcul du tableau des distances des pixels x a Q_ks    
    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); 
//calcul du tableau des distances au segment [P_ks,Q_ks]    
    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 du tableau des poids de Warpk   
    alpha = (Lk^p./(a+Dk)).^b; 
//incrementation de Alpha et de Warp      
    Alpha = Alpha +  alpha;   
    WARP = WARP + alpha.*Warpk;  
end 
    Warp = WARP./Alpha ;     
endfunction  
