boundtriang.m

数学优化工具箱

function [hh, lowerbound, objval] = boundtriang(nbrimages, numvar, inrect, xm, camlist)

b=sparse(zeros(numvar,1));
At_l=sparse(zeros(0,numvar));
c_l=sparse(zeros(0,1));
At=sparse(zeros(0,numvar));
A_f = sparse(zeros(0,numvar));
c_f=sparse(zeros(0,1));
c=sparse(zeros(0,1));

clear K;
K.l = 0;
K.q = [];
K.f = 0;

% Point parametrization : [x1, x2, x3, 1]
% Variable order : x1, x2, x3, {t, s, sp, r, rp}_repeat

lower_vars = 7:5:7+5*(nbrimages-1);
numerator_vars = 4:5:4+5*(nbrimages-1);
denominator_vars = 5:5:5+5*(nbrimages-1);

for i = 1:nbrimages
    
    %% Rectangle bounds
    tL = inrect.tllist(i);  tU = inrect.tulist(i);
    sL = inrect.sllist(i);  sU = inrect.sulist(i);
    
    shiftpos = 3+(i-1)*5;
    cshift = numvar*(i-1);
    
    b(shiftpos+4) = -1;
    
    p = camlist{i};

    u = xm(1,i);
    v = xm(2,i);


    %% SOCP constraints

    %% Eqn 15

    ctmp = sparse(3,1);  Atmp = sparse(3,numvar);

    ctmp(2,1) = u*p(3,4)-p(1,4);
    ctmp(3,1) = v*p(3,4)-p(2,4);

    Atmp(1,3+(i-1)*5+1) = -1;
    Atmp(2,1) = -u*p(3,1)+p(1,1);
    Atmp(2,2) = -u*p(3,2)+p(1,2);
    Atmp(2,3) = -u*p(3,3)+p(1,3);
    Atmp(3,1) = -v*p(3,1)+p(2,1);
    Atmp(3,2) = -v*p(3,2)+p(2,2);
    Atmp(3,3) = -v*p(3,3)+p(2,3);
    
    c = [c; ctmp];
    At = [At; Atmp];
    K.q = [K.q, 3];


    %% Eqn 1

    ctmp = []; Atmp = [];
    ctmp = sparse(3,1);  Atmp = sparse(3,numvar);

    ctmp(2,1) = 2*sqrt(tL)*tU/(tU-tL);
    
    Atmp(1,shiftpos+3) = -1;
    Atmp(1,shiftpos+5) = -1;
    Atmp(2,shiftpos+1) = 2*sqrt(tL)/(tU-tL);
    Atmp(3, shiftpos+3) = 1;
    Atmp(3,shiftpos+5) = -1;
    
    c = [c; ctmp];
    At = [At; Atmp];
    K.q = [K.q, 3];


    %% Eqn 2

    ctmp = []; Atmp = [];
    ctmp = sparse(3,1);  Atmp = sparse(3,numvar);

    ctmp(2,1) = 2*sqrt(tU)*tL/(tU-tL);
    
    Atmp(1,shiftpos+2) = -1;
    Atmp(1,shiftpos+3) = 1;
    Atmp(1,shiftpos+4) = -1;
    Atmp(1,shiftpos+5) = 1;
    Atmp(2,shiftpos+1) = 2*sqrt(tU)/(tU-tL);
    Atmp(3,shiftpos+2) = 1;
    Atmp(3,shiftpos+3) = -1;
    Atmp(3,shiftpos+4) = -1;
    Atmp(3,shiftpos+5) = 1;
    
    c = [c; ctmp];
    At = [At; Atmp];
    K.q = [K.q, 3];


    %% Linear constraints
    
    ctmp = []; Atmp = [];
    ctmp = sparse(12,1);  Atmp = sparse(12,numvar);
    
    ctmp(1,1) = -sL*tU;
    ctmp(2,1) = -sU*tL;
    ctmp(3,1) = sU*tU;
    ctmp(4,1) = sL*tL;
    ctmp(9,1) = -tL;
    ctmp(10,1) = tU;
    ctmp(11,1) = -sL;
    ctmp(12,1) = sU;
%    ctmp(13,1) = p(3,4);
    c_f(i,1) = p(3,4);

    Atmp(1,shiftpos+1) = -sL;
    Atmp(1,shiftpos+3) = -(tU-tL);
    Atmp(2,shiftpos+1) = -sU;
    Atmp(2,shiftpos+2) = tU-tL;
    Atmp(2,shiftpos+3) = -(tU-tL);
    Atmp(3,shiftpos+1) = sU;
    Atmp(3,shiftpos+3) = tU-tL;
    Atmp(4,shiftpos+1) = sL;
    Atmp(4,shiftpos+2) = -(tU-tL);
    Atmp(4,shiftpos+3) = tU-tL;
    Atmp(5,shiftpos+5) = -1;
    Atmp(6,shiftpos+4) = -1;
    Atmp(6,shiftpos+5) = 1;
    Atmp(7,shiftpos+2) = -1;
    Atmp(7,shiftpos+3) = 1;
    Atmp(8,shiftpos+3) = -1;
    Atmp(9,shiftpos+1) = -1;
    Atmp(10,shiftpos+1) = 1;
    Atmp(11,shiftpos+2) = -1;
    Atmp(12,shiftpos+2) = 1;    
    %Atmp(13,1) = -p(3,1);  Atmp(13,2) = -p(3,2);
    %Atmp(13,3) = -p(3,3);  Atmp(13,shiftpos+2) = 1;
    
    A_f(i,1) = -p(3,1);  A_f(i,2) = -p(3,2);
    A_f(i,3) = -p(3,3);  A_f(i,shiftpos+2) = 1;
    K.f = K.f+1;
    
    c_l = [c_l; ctmp];
    At_l = [At_l; Atmp];
    K.l = K.l+12;

end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Call to SeDuMi
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pars=[];
pars.fid=0;
pars.eps=0;
[x,y,info]=sedumi([A_f;At_l;At],b,[c_f(:);c_l;c],K,pars);

if info.dinf==1, % If dual infeasible
    disp('Infeasible. Setting lower bound to infinity.');
    res=Inf*ones(1,nbrimages);
    lowerbound=Inf;
    objval = Inf;
    hh.y = zeros(numvar,1);
else
    lowerbound = sum(y(lower_vars));
    objval = sum(y(numerator_vars) ./ y(denominator_vars));	
    hh.y=y;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%