http://blog.sina.com.cn/s/blog_5dd2e9270100xdcp.html

代码不能运行，不知道咋回事

poison matting matlab code

注：显示图像后，单击鼠标左键连出完全背景区域，双击左键结束；

图像改变后，再单击鼠标左键连出完全前景区域，双击左键结束。即可得到结果。

A=imread(‘ppmm2.bmp‘);

imshow(A);

hold on;

Threshold = 1;

UniBack=[0

0

255];

% choose background object, and B results into the Foreground Object and the

% Boundary Region

% choose the outer line of the boundary region. Single left click the mouse to specify vertice.Double left click or single right click to finish;

OuterBW = roipoly(A);

B(:,:,1) = immultiply(A(:,:,1),OuterBW);

B(:,:,2) = immultiply(A(:,:,2),OuterBW);

B(:,:,3) = immultiply(A(:,:,3),OuterBW);

Back(:,:,1) = immultiply(A(:,:,1),~OuterBW);

Back(:,:,2) = immultiply(A(:,:,2),~OuterBW);

Back(:,:,3) = immultiply(A(:,:,3),~OuterBW);

imshow(B);

% choose foreground object, and C results into a Circle of Interest

% choose the inner line of the boundary region. Single left click the mouse to specify vertice.Double left click or single right click to finish;

InnerBW = roipoly(B);

C(:,:,1) = immultiply(B(:,:,1),~InnerBW);

C(:,:,2) = immultiply(B(:,:,2),~InnerBW);

C(:,:,3) = immultiply(B(:,:,3),~InnerBW);

Fore(:,:,1) = immultiply(A(:,:,1),InnerBW);

Fore(:,:,2) = immultiply(A(:,:,2),InnerBW);

Fore(:,:,3) = immultiply(A(:,:,3),InnerBW);

imshow(C);

hold off;

RawAlpha = (double(OuterBW)+double(InnerBW))/2;

I = double(A);

% Method of "averaging" to get raw fore and background colors

for i = 1:size(RawAlpha,1)

for j = 1:size(RawAlpha,2)

RawFore(i,j,1)=0;

RawFore(i,j,2)=0;

RawFore(i,j,3)=0;

RawBack(i,j,1)=0;

RawBack(i,j,2)=0;

RawBack(i,j,3)=0;

if RawAlpha(i,j) == 1

RawFore(i,j,:)=Fore(i,j,:);

end;

if RawAlpha(i,j) == 0

RawBack(i,j,:)=Back(i,j,:);

end;

if RawAlpha(i,j) == 0.5

r=1;

while 1

TestFore = InnerBW(max(i-r,1):min(i+r,size(A,1)),max(j-r,1):min(j+r,size(A,2)));

if size(find(TestFore),1)

[i1, j1] = find(TestFore);

i2 = i1 + max(i-r,1) -1;

j2 = j1 + max(j-r,1) -1;

Rs = double(Fore(i2,j2,:));

Ds(:,1)=diag(Rs(:,:,1));

Ds(:,2)=diag(Rs(:,:,2));

Ds(:,3)=diag(Rs(:,:,3));

if size(Ds,1) == 1

RawFore(i,j,:) =Ds;

else

RawFore(i,j,:) = sum(Ds) / size(Ds,1);

end;

clear Ds;

break;

else

r=r+1;

end;

end;

r=1;

while 1

TestBack = ~OuterBW(max(i-r,1):min(i+r,size(A,1)),max(j-r,1):min(j+r,size(A,2)));

if size(find(TestBack),1)

[i1, j1] = find(TestBack);

i2 = i1 + max(i-r,1) -1;

j2 = j1 + max(j-r,1) -1;

Rs = double(Back(i2,j2,:));

Ds(:,1)=diag(Rs(:,:,1));

Ds(:,2)=diag(Rs(:,:,2));

Ds(:,3)=diag(Rs(:,:,3));

if size(Ds,1) == 1

RawBack(i,j,:) =Ds;

else

RawBack(i,j,:) = sum(Ds) / size(Ds,1);

end;

clear Ds;

break;

else

r=r+1;

end;

end;

end;

end;

end;

%此处没有用高斯滤波

Denorm = RawFore - RawBack;

%red channel

I1 = I(:,:,1);

Denorm1 = Denorm(:,:,1);

for i=1:size(Denorm1,1)

for j= 1: size(Denorm1,2)

if Denorm1(i,j)==0

Denorm1(i,j)=1;

end;

end;

end;

OldAlpha = RawAlpha;

NewAlpha = RawAlpha;

h1=0;

while 1

for i=1:size(OldAlpha,1)

for j=1:size(OldAlpha,2)

NewAlpha(i,j) = OldAlpha(i,j);

if RawAlpha(i,j) == 0.5

Roui = ((I1(i+1,j) + I1(i-1,j) - 2 * I1(i,j)) * Denorm1(i,j) - (I1(i+1,j) - I1(i,j)) * (Denorm1(i+1,j) - Denorm1(i,j)))/(Denorm1(i,j) * Denorm1(i,j));

Rouj = ((I1(i,j+1) + I1(i,j-1) - 2 * I1(i,j)) * Denorm1(i,j) - (I1(i,j+1) - I1(i,j)) * (Denorm1(i,j+1) - Denorm1(i,j)))/(Denorm1(i,j) * Denorm1(i,j));

Rou = Roui + Rouj;

NewAlpha(i,j) = (OldAlpha(i+1,j) + NewAlpha(i-1,j) + OldAlpha(i,j+1) + NewAlpha(i,j-1) - Rou) / 4;

if NewAlpha(i,j)<0

NewAlpha(i,j)=0;

end;

if NewAlpha(i,j)>1

NewAlpha(i,j)=1;

end;

end;

end;

end;

% imshow(uint8(NewAlpha*255));

DifferenceAlpha = abs(NewAlpha - OldAlpha);

OldAlpha = NewAlpha;

if sum(sum(DifferenceAlpha)) < Threshold

break;

end;

h1=h1+1;

end;

for i=1:size(A,1)

for j=1:size(A,2)

if OldAlpha(i,j)==0

NewImage(i,j,:)=UniBack‘;

else

NewImage(i,j,1)=UniBack(1)*(1-OldAlpha(i,j))+RawFore(i,j,1)*OldAlpha(i,j);

NewImage(i,j,2)=UniBack(2)*(1-OldAlpha(i,j))+RawFore(i,j,2)*OldAlpha(i,j);

NewImage(i,j,3)=UniBack(3)*(1-OldAlpha(i,j))+RawFore(i,j,3)*OldAlpha(i,j);

end;

end;

end;

figure,imshow(NewAlpha);

figure,imshow(uint8(NewImage));

figure,imshow(uint8(RawFore));

figure,imshow(uint8(RawBack));

原文：http://www.cnblogs.com/natalie/p/5631548.html