function [PDelta1,PDelta2]=RPCBlockAdjust(X3D, r1_a,c1_a,r2_a,c2_a, scale_offsets_1,scale_offsets_2,p1_1, p2_1, p3_1, p4_1, p1_2, p2_2, p3_2, p4_2)


    %LevMar BA refinement
 
npts=length(r1_a);
x1(1,:)=c1_a;
x1(2,:)=r1_a;
x2(1,:)=c2_a;
x2(2,:)=r2_a;

matches_2D_Inliers = [x1(2,:)', x1(1,:)', x2(2,:)', x2(1,:)'];
param_val = X3D;
X3Dinit=reshape(param_val,[3*size(X3D,2) 1]);

initial_guess=[0;0;0;0;X3Dinit];
reference=matches_2D_Inliers;
f = @(param)metricdist(param, reference, X3D,...
    scale_offsets_1,scale_offsets_2,p1_1,p2_1,p3_1,p4_1,p1_2,p2_2,p3_2,p4_2);
options = optimoptions(@lsqnonlin,'Display','iter', 'Algorithm', ...
    'levenberg-marquardt', 'TolX', 1e-13, 'TolFun', 1e-13);
options.MaxIter = 10;
[solution,~,~,~] = lsqnonlin(f,initial_guess, [], [], options);
    
PDelta1(1)=solution(1);PDelta1(2)=solution(2);
PDelta2(1)=solution(3);PDelta2(2)=solution(4);
end

function F = metricdist(param, matches_2D_Inliers,X3D,scale_offsets_1,scale_offsets_2, ...
    p1_1,p2_1,p3_1,p4_1,p1_2,p2_2,p3_2,p4_2)
%        normalizeSpaceCoordinates(rpc, lon,lat,alt)
    Xs_1=scale_offsets_1(3);Xo_1=scale_offsets_1(4);
    Ys_1=scale_offsets_1(1);Yo_1=scale_offsets_1(2);
    Zs_1=scale_offsets_1(5);Zo_1=scale_offsets_1(6);
    Xs_2=scale_offsets_2(3);Xo_2=scale_offsets_2(4);
    Ys_2=scale_offsets_2(1);Yo_2=scale_offsets_2(2);
    Zs_2=scale_offsets_2(5);Zo_2=scale_offsets_2(6);
  
    for ij=1:size(matches_2D_Inliers,1)
        X1 = RPCnormalization(param(4+(ij-1)*3+1),Xo_1,Xs_1);
        Y1 = RPCnormalization(param(4+(ij-1)*3+2),Yo_1,Ys_1);
        Z1 = RPCnormalization(param(4+(ij-1)*3+3),Zo_1,Zs_1);
        X2 = RPCnormalization(param(4+(ij-1)*3+1),Xo_2,Xs_2);
        Y2 = RPCnormalization(param(4+(ij-1)*3+2),Yo_2,Ys_2);
        Z2 = RPCnormalization(param(4+(ij-1)*3+3),Zo_2,Zs_2);

        numr1(ij,1)=RPCforwardform(p1_1,p2_1,X1,Y1,Z1)+param(1);
        numc1(ij,1)=RPCforwardform(p3_1,p4_1,X1,Y1,Z1)+param(2);
        numr2(ij,1)=RPCforwardform(p1_2,p2_2,X2,Y2,Z2)+param(3);
        numc2(ij,1)=RPCforwardform(p3_2,p4_2,X2,Y2,Z2)+param(4); 
         
        numr1(ij,1) = RPCunnormalization(numr1(ij,1),scale_offsets_1(7),scale_offsets_1(9));
        numc1(ij,1) = RPCunnormalization(numc1(ij,1),scale_offsets_1(8),scale_offsets_1(10));
        numr2(ij,1) = RPCunnormalization(numr2(ij,1),scale_offsets_2(7),scale_offsets_2(9));
        numc2(ij,1) = RPCunnormalization(numc2(ij,1),scale_offsets_2(8),scale_offsets_2(10));
        
    end
    
    F = abs(matches_2D_Inliers - [numr1, numc1, numr2, numc2] ) ;
end