function [Hl,Hr,F,rn,rns,Kl,Kr] = viprectification_homographies(Pleft,Pright,...
    w,h,n,lambda)
%Compute rectifying homographies.
%  [HL,HR,F,RESNORM,RESNORMS,KL,KR] = VIPRECTIFICATION_HOMOGRAPHIES(...
%  Pleft,Pright,W,H,N,LAMBDA) accepts corresponding [x;y] points Pleft and
%  Pright, and image width and height, W and H, and it calls the
%  VIPRECTIFICATION_LSQNL function. It does N optimizations with randomized
%  starts and takes the best result. It returns left and right homographies
%  HL and HR, the fundamental matrix F, minimization outputs RESNORM and
%  RESNORMS, and derived left and right intrinsic matrices KL and KR.
%  LAMBDA is the normalization factor to feed to VIPRECTIFICATION_LSQNL.
%  The default N is 1; the default LAMBDA uses VIPRECTIFICATION_LSQNL's
%  default value.

    if nargin < 5 || isempty(n)
        n = 1;
    end
    if nargin < 6 || isempty(lambda)
        lambda = [];
    end

    % Run optimization.
    options = optimset('Jacobian','off', 'Display','off');
    options.MaxFunEvals = 10000;
    options.MaxIter = 1000;
    bnds = [1; pi/2*ones(5,1)];
    best_m = zeros(6,1);
    best_rn = inf;
    rns = zeros(1,n);
    for i=1:length(rns)
        start = [2*rand-1.0; 0.2*rand(5,1)-0.1];
        [m,rn] = lsqnonlin(@(x)viprectification_lsqnl(x,...
            double(Pleft),double(Pright),w,h,lambda),start,-bnds,bnds,options);
        rns(i) = rn;
        if rn < best_rn
            best_rn = rn;
            best_m = m;
        end
    end
    m = best_m;
    rn = best_rn;
    
    a = (3^m(1))*(w+h)';
    K = [a 0 w/2; 0 a h/2; 0 0 1];
    Rr = viprectification_rotmat(m(2:4), 'YZX');
    Rl = viprectification_rotmat(m(5:6), 'YZ');
    F = inv(K)'*(Rr'*[0 0  0
                      0 0 -1
                      0 1  0]*Rl)/K;
    Hr = (K*Rr)/K;
    Hl = (K*Rl)/K;
    
    % Recenter projections
    wh2 = [w/2; h/2];
    l_cent = Hl*[wh2;1]; % Reprojected left center.
    del_l = wh2 - l_cent(1:2)/l_cent(3);
    r_cent = Hr*[wh2;1]; % Reprojected right center.
    del_r = wh2 - r_cent(1:2)/r_cent(3);
    % Set vertical displacements same, favoring the left.
    del_r(2) = del_l(2);
    % Modify the old intrinsics to construct different new ones.
    Kr = K;
    Kr(1:2,3) = Kr(1:2,3) + del_r;
    Kl = K;
    Kl(1:2,3) = Kl(1:2,3) + del_l;
    % Recompute the homographies.
    Hr = (Kr*Rr)/K;
    Hl = (Kl*Rl)/K;
end


