function layer(im,alpha,id)
	run('vlfeat-0.9.16/toolbox/vl_setup');%add the VL Feat library
	level=2;%weight of spatial coordinate in the feature vector
	if(size(alpha,3)>1)%change alpha matte to grayscale in case it is not
		alpha=rgb2gray(alpha);
	end
	alpha=reshape(alpha,[],1);
	[m n d]=size(im);
	nn=[10];%k=10 in KNN search
	row=zeros(n*m*sum(nn),2,'uint32');
	[a b]=ind2sub([m n],1:m*n);
	feature=[reshape(im,m*n,d)';[a;b]/sqrt(m*m+n*n)*level+rand(2,m*n)*1e-6];%feature vector is (r,g,b,x,y)
	now=0;
	for i=1:size(nn,1)%connect each pixel to its k nearest neighbors in the feature space
		ind=vl_kdtreequery(vl_kdtreebuild(feature),feature,feature,'NUMNEIGHBORS',nn(i),'MAXNUMCOMPARISONS',nn(i)*2);
		a=reshape(repmat(uint32(1:m*n),nn(i),1),[],1);
		b=reshape(ind,[],1);
		row(now+1:now+m*n*nn(i),:)=[min(a,b) max(a,b)];
		feature(d+1:d+2,:)=feature(d+1:d+2,:)/100;
		now=now+m*n*nn(i);
	end
	row=unique(row,'rows');

	weight=(1-2*abs(reshape(alpha,m,n)-0.5));%compute the weight based on the known alpha matte
	
	[a b]=ndgrid(1:m,1:n-1);
	row2(1:m*(n-1),1)=sub2ind([m n],reshape(a,1,[]),reshape(b,1,[]));
	row2(1:m*(n-1),2)=row2(1:m*(n-1),1)+m;
	[a b]=ndgrid(1:m-1,1:n);
	row2((1:(m-1)*n)+m*(n-1),1)=sub2ind([m n],reshape(a,1,[]),reshape(b,1,[]));
	row2((1:(m-1)*n)+m*(n-1),2)=row2((1:(m-1)*n)+m*(n-1),1)+1;
	row2=[row2;double(row)];

	ll=m*(n-1)+(m-1)*n;
	value=min(weight(row2(:,1)),weight(row2(:,2)))+0.001;%+0.001 to avoid zero weight
	A2=sparse(double([row2(:,1);row2(:,1)+m*n]),double([row2(:,2);row2(:,2)+m*n]),[value;value],m*n*2,m*n*2);
	A2=A2+A2';
	D2=spdiags(sum(A2,2),0,n*m*2,n*m*2);
	L2=D2-A2;
	lambda2=1;
	alpha=min(max(alpha,0),1);
	row2=zeros(m*n*4,2);
	value2=zeros(m*n*4,1);
	row2(1:m*n*2,1)=[1:m*n 1:m*n];
	row2((1:m*n*2)+m*n*2,1)=[1:m*n 1:m*n]+m*n;
	row2(1:m*n,2)=1:m*n;
	row2((1:m*n)+m*n,2)=(1:m*n)+m*n;
	row2((1:m*n)+m*n*2,2)=1:m*n;
	row2((1:m*n)+m*n*3,2)=(1:m*n)+m*n;
	value2(1:m*n)=lambda2*alpha.*alpha;
	value2((1:m*n)+m*n)=lambda2*alpha.*(1-alpha);
	value2((1:m*n)+m*n*2)=lambda2*(1-alpha).*alpha;
	value2((1:m*n)+m*n*3)=lambda2*(1-alpha).*(1-alpha);
	data=sparse(row2(:,1),row2(:,2),value2,m*n*2,m*n*2);
	B=zeros(m*n*2,1);
	Fg=zeros(m,n,3);
	Bg=zeros(m,n,3);
	tmp=L2+data;
	LL=ichol(tmp);
	for j=1:3%for rgb channels
		j
		B(1:m*n)=lambda2*alpha.*reshape(im(:,:,j),m*n,1);
		B((1:m*n)+m*n)=lambda2*(1-alpha).*reshape(im(:,:,j),m*n,1);
		out=pcg(tmp,B,1e-10,2000,LL,LL');
		Fg(:,:,j)=reshape(out(1:m*n),m,n);
		Bg(:,:,j)=reshape(out((1:m*n)+m*n),m,n);
	end
	Fg=min(max(Fg,0),1);
	Bg=min(max(Bg,0),1);
	imwrite(Fg,sprintf('output/%02d_foreground.png',id),'png');
	imwrite(Bg,sprintf('output/%02d_background.png',id),'png');
	FgAlpha=Fg.*repmat(reshape(alpha,m,n),[1 1 3])+repmat(reshape(1-alpha,m,n),[1 1 3]);
	BgAlpha=Bg.*repmat(reshape(1-alpha,m,n),[1 1 3]);
	figure;
	imshow(FgAlpha);
	imwrite(FgAlpha,sprintf('output/%02d_foreground_times_alpha.png',id),'png');
	figure;
	imshow(BgAlpha);
	imwrite(BgAlpha,sprintf('output/%02d_background_times_alpha.png',id),'png');
	imwrite(im.*repmat(reshape(alpha,m,n),[1 1 3])+repmat(reshape(1-alpha,m,n),[1 1 3]),sprintf('output/%02d_image_times_alpha.png',id),'png');
end