closed-form 抠图背景和前景重建实现方法

在www.alphamatting.com这个网站中，可以下载到closed-from抠图的代码，在代码中有关于背景和前景的重建代码，matlab代码如下：

[F,B]=solveFB(I,alpha);%I是输入的原始图像，alpha是整个图像的抠图结果，F和B分别是返回的前景 和别经值

function [F,B]=solveFB(I,alpha)
  [h,w,c]=size(I);

%得到图像的大小高度为h，宽度为w，通道为c，RGB图像中c=3
  mask=(alpha>=0.02).*(alpha<=0.98);

mask取两个极端情况，alpha的值很小，一个就是alpha的值很大接近到1，那么就可以认为分别是背景和前景

mask中保存的是不确定区域，在三分图中有背景区域、前景区域、和不确定区域

solveFB函数就要从不确定区域中求的F和B的值
  [Gx,Gy,Gd1,Gd2]=getGMatByMask(w,h,mask);

getGMatByMask这个函数很复杂，

  G=[Gx;Gy;Gd1;Gd2];梯度计算Gx、GyGd1、Gd2都是相同的大小，imagesize*imagesize
  Ga=G*alpha(:);Ga是个向量，大小为iamgesize*4的行向量
  wgf=abs(Ga).^0.5+0.003*repmat((1-alpha(:)),4,1);
  wgb=abs(Ga).^0.5+0.003*repmat(alpha(:),4,1);


  
  wf=(alpha(:)>0.98)*100+0.03*alpha(:).*(alpha(:)>0.7)+0.01*(alpha(:)<0.02);
  wb=(alpha(:)<0.02)*100+0.03*(1-alpha(:)).*(alpha(:)<0.3)+0.01*(alpha(:)>0.98);
  
  
  wgf=spdiags(wgf(:),0,length(wgf),length(wgf));
  wgb=spdiags(wgb(:),0,length(wgb),length(wgb)); 
  wf=spdiags(wf(:),0,length(wf),length(wf));
  wb=spdiags(wb(:),0,length(wb),length(wb)); 
  
  
  for t=1:c
    tI=I(:,:,t);
    Ag=[wgf*G,sparse(size(G,1),size(G,2));sparse(size(G,1),size(G,2)),wgb* ...
        G];
    bg=zeros(size(Ag,1),1);
    Ai=[wf,sparse(w*h,w*h);sparse(w*h,w*h),wb];
    bi=[wf*tI(:).*(alpha(:)>0.02);wb*tI(:).*(alpha(:)<0.98)];
    As=[spdiags(alpha(:),0,w*h,w*h),spdiags(1-alpha(:),0,w*h,w*h)];
    bs=tI(:);
    A=[Ai;As;Ag];
    b=[bi;bs;bg];
    x=(A'*A)\(A'*b);
    F(:,:,t)=reshape(x(1:w*h),h,w);
    B(:,:,t)=reshape(x(w*h+1:end),h,w);


  end  

下面这个代码是微分代码，

function [Gx,Gy,G3,G4]=getGMatByMask(w,h,mask)
imgSize=w*h;


dS=[1,-1];
filtSizeS=1;




%indsGx1=[]; indsGx2=[]; valsGx=[];
%indsGy1=[]; indsGy2=[]; valsGy=[];
indsGx1=zeros(imgSize*2,1);
indsGx2=zeros(imgSize*2,1);
valsGx=zeros(imgSize*2,1);
indsGy1=zeros(imgSize*2,1);
indsGy2=zeros(imgSize*2,1);
valsGy=zeros(imgSize*2,1);
indsG31=zeros(imgSize*2,1);
indsG32=zeros(imgSize*2,1);
valsG3=zeros(imgSize*2,1);
indsG41=zeros(imgSize*2,1);
indsG42=zeros(imgSize*2,1);
valsG4=zeros(imgSize*2,1);%没一个都是图像大小的2倍


indy=0; indx=0; ind3=0; ind4=0;


for x=1:w-1,
  for y=1:h,
    if ((~mask(y,x))&(~mask(y,x+1)))
      continue
    end  
    for disp=0:filtSizeS,
      indx=indx+1;
      indsGx1(indx)=imIndexToVect(y,x,h);
      indsGx2(indx)=imIndexToVect(y,x+disp,h);
      valsGx(indx)=dS(disp+1);
    end
  end
end
for x=1:w,
  for y=1:h-1,
    if ((~mask(y,x))&(~mask(y+1,x)))
      continue
    end
    for disp=0:filtSizeS,
      indy=indy+1;
      indsGy1(indy)=imIndexToVect(y,x,h);
      indsGy2(indy)=imIndexToVect(y+disp,x,h);
      valsGy(indy)=dS(disp+1);
    end;
  end;
end


for x=1:w-1,
  for y=1:h-1,
    if ((~mask(y,x))&(~mask(y+1,x+1)))
      continue
    end
    for disp=0:filtSizeS,
      ind3=ind3+1;
      indsG31(ind3)=imIndexToVect(y,x,h);
      indsG32(ind3)=imIndexToVect(y+disp,x+disp,h);
      valsG3(ind3)=dS(disp+1);
    end
  end
end


for x=1:w-1,
  for y=2:h,
    if ((~mask(y,x))&(~mask(y-1,x+1)))
      continue
    end
    for disp=0:filtSizeS,    
      ind4=ind4+1;
      indsG41(ind4)=imIndexToVect(y,x,h);
      indsG42(ind4)=imIndexToVect(y-disp,x+disp,h);
      valsG4(ind4)=dS(disp+1);
    end;
  end;
end;
%'done inds'
indsGx1=indsGx1(1:indx);
indsGx2=indsGx2(1:indx);
valsGx=valsGx(1:indx);
indsGy1=indsGy1(1:indy);
indsGy2=indsGy2(1:indy);
valsGy=valsGy(1:indy);
indsG31=indsG31(1:ind3);
indsG32=indsG32(1:ind3);
valsG3=valsG3(1:ind3);
indsG41=indsG41(1:ind4);
indsG42=indsG42(1:ind4);
valsG4=valsG4(1:ind4);






Gx=sparse(indsGx1,indsGx2,valsGx,imgSize,imgSize);
Gy=sparse(indsGy1,indsGy2,valsGy,imgSize,imgSize);
G3=sparse(indsG31,indsG32,valsG3,imgSize,imgSize);
G4=sparse(indsG41,indsG42,valsG4,imgSize,imgSize);