function [I1,I2]=sepRefLOQO(I,xf,yf,xb,yb)
%Computes a 2 layers additive decomposition of an input image, 
%given a sparse set of points whose gradients are labeled as belonging 
%to one of the layers.
%
%[I1,I2]=sepRefLOQO(I,xf,yf,xb,yb);
%    I - input superposition image.
%    xf,yf- coordinates of labeled foreground gradients.
%    xb,yb- coordinates of labeled background gradients.
%    I1,I2- output layers.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




itrN=3;


sigv=[0.0100    0.0520];
piw=[0.4000    0.6000];

w1=2;
w2=1;


[h,w]=size(I);
imgSize=h*w;


f_inds=imIndexToVect(yf,xf,h);
b_inds=imIndexToVect(yb,xb,h);


[Gx,Gy,Gxx,Gxy,Gyy]=getGMat12(w,h);


LB=zeros(imgSize,1); UB=I(:);

A1=[Gx;Gy;Gxx;Gyy;Gxy];
A=[A1;A1];
b=[zeros(size(A1,1),1);A1*I(:)];
f1=[ones(imgSize*2,1)*w1;w2*ones(imgSize*2,1);w2*ones(imgSize,1)*2];
f2=f1;
f1([f_inds,f_inds+imgSize])=0;
f1([b_inds,b_inds+imgSize])=100;
f2([f_inds,f_inds+imgSize])=100;
f2([b_inds,b_inds+imgSize])=0;
f=[f1;f2];

rinds1=find(sum(A~=0,2)==0);
rinds2=find(f==0);
inds=setdiff([1:size(A,1)],[rinds1;rinds2]);
A=A(inds,:);
b=b(inds,:);
f=f(inds);
oA=A;
f=[zeros(imgSize,1);f;f];

of=f;

A=[A,speye(size(A,1)),-speye(size(A,1))];
equConstNum=size(A,1);

LB=zeros(size(f));
UB=[UB;ones(length(LB)-length(UB),1)*3];



for k=1:itrN
[x,lambda,how]=loqo([],f,A,b,LB,UB,[],equConstNum,2);
x=x(1:imgSize);
I1=reshape(x,h,w);


errv=full(abs(oA*x-b));
for j=1:length(sigv)
    emw(:,j)=piw(j)/(2*sigv(j))*exp(-errv/sigv(j));
end

emw=emw./(sum(emw,2)*(sigv));

f(1+length(I(:)):end)=of(1+length(I(:)):end).*[sum(emw,2);sum(emw,2)];


end;

I2=I-I1;
