泊松抠图

泊松抠图MATLAB代码        (高手代码)

图像处理 2010-09-27 09:54:25 阅读93 评论0   字号：大中小 订阅 

poison matting matlab code 

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

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

A=imread('ppmm2.bmp');

imshow(A);title(‘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); title(‘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); title(‘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); title(‘NewAlpha’)

figure,imshow(uint8(NewImage)); title(‘uint8(NewImage)’)

figure,imshow(uint8(RawFore)); title(‘uint8(RawFore’)

figure,imshow(uint8(RawBack)); title(‘uint8(RawBack)’)