function [I1,I2]=sepRefItrLS(I,xf,yf,xb,yb,iI)
%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]=sepRefItrLS(I,xf,yf,xb,yb);
%    I - input superposition image.
%    xf,yf- coordinates of labeled foreground gradients.
%    xb,yb- coordinates of labeled background gradients.
%    iI- (optional) initial guess. If no guess is provided, 
%        I is used as initial guess. 
%    I1,I2- output layers.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


if (~exist('iI'))
  iI=I;
end



%Number of EM iterations
EMitrN=15;

%Number of inner M step iterations=Number if iterative reweighted least
%squares iterations to solve each L1 optimization problem.
%Note that in case MitrN is set to one, the resulting rewighting is
%equivalent to iterative reweighted least squares directly on the sparse
%mixture cost function 
MitrN=2;

%Mixture parameters
sigv=[0.0100    0.0520];
piw=[0.4000    0.6000];

w1=2;
w2=1;


[h,w]=size(I);
imgSize=h*w;
LB=zeros(imgSize,1); UB=I(:);

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

%keyboard

%building constraints matrix
[Gx,Gy,Gxx,Gxy,Gyy]=getGMat12(w,h);
A1=[Gx;Gy;Gxx;Gyy;Gxy];
%A1=[Gx;Gy];
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)];
%f1=[ones(imgSize*2,1)*w1];

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;

f1([b_inds+imgSize*2,b_inds+imgSize*3,b_inds+imgSize*4])=4;
f2([f_inds+imgSize*2,f_inds+imgSize*3,f_inds+imgSize*4])=4;

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;
ob=b;
of=f;

df=spdiags(f,0,length(f),length(f));
A=df*oA; b=df*ob;

if (max(abs(iI(:)-I(:)))<0.05)|(max(abs(iI(:)))<0.05 )
   x=(A'*A)\(A'*b);
else
   x=iI(:);
end

for j=1:EMitrN
   %keyboard 
   %M step: solving L1 optimization problem using iterative reweighted
   %least squares
   for i=1:MitrN
      e=abs(oA*x-ob);
      e=max(e,0.00001);
      e=1./e;
      E=spdiags(e,0,length(f),length(f));
      x=(A'*E*A)\(A'*E*b);
   end

  
   %E step- reweighting equations. 
   e=abs(A*x-b);
   terrv=full(abs(oA*x-ob));

   for j=1:length(sigv)
      emw(:,j)=piw(j)/(2*sigv(j))*exp(-terrv/sigv(j));
   end

   emw=emw./(sum(emw,2)*(sigv));
   f=of.*sum(emw,2);
   df=spdiags(f,0,length(f),length(f));
   A=df*oA; b=df*ob;
   
end


I1=reshape(x,h,w);
I2=I-I1;


