您的位置: 首页  >  文章  >  Poisson image editing算法实现的Matlab代码解析

Poisson image editing算法实现的Matlab代码解析

分类: 文章 • 2024-11-12 17:10:31
之前发了数篇系列博文来仔细研究Poisson Image Editing算法，每次重新审视和深入，仿佛都能有更为深刻的认识和很大的收获。这应该算是我这个系列的完结篇，会用用Matlab代码一点一点的演示，原文作者到底是如何设计和实现他那个强大且影响深远的算法的。希望你在看本文之前务必参考一下文章来了解算法原理，本文将主要讲解编程实现的问题，对于前面讲过的内容，我不会深究。但我个人总体的感觉是，现在图像处理算法对数学的要求是越来越高了，像泊松融合、泊松抠图这样的算法如果没有偏微分方程（本算法中涉及Poisson Equation）、数值计算（需要高斯-塞德尔迭代法或者共轭梯度法）、场论（涉及散度、梯度和拉普拉斯算子）、矩阵论（涉及大型稀疏矩阵线性方程组）、数学物理方法（涉及带狄利克雷条件的椭圆曲线方程）、最优化理论（涉及到了欧拉-朗格拉日变分法）这些艰深的数学知识，要理解起来实在太困难了（每一步都是一个坎，我感觉有一个地方过不了，那篇经典论文的精华就无法领会，更别说编程实现了Poisson image editing算法实现的Matlab代码解析）。


素材是三张图片，如下


                                                                                   bear.jpg                                                              bear-mask.jpg

Poisson image editing算法实现的Matlab代码解析                         Poisson image editing算法实现的Matlab代码解析


pool-target.jpg

Poisson image editing算法实现的Matlab代码解析


首先，我们清理一下Matlab的环境：

[plain] view plain copy
  1. close all;  
  2. clear;  
  3. clc;  

然后，读入三张图片（注意需要把mask图以二值图的形式读入）
[plain] view plain copy
  1. TargetImg   = imread('pool-target.jpg');  
  2. SourceImg   = imread('bear.jpg');  
  3. SourceMask  = im2bw(imread('bear-mask.jpg'));  

用函数bwboundaries(W,CONN) 来获取二值图中对象的轮廓，其中CONN 为8或4，指示联通性采用4方向邻域点判别还是8方向邻域点判别，默认为8。

[plain] view plain copy
  1. [SrcBoundry, ~] = bwboundaries(SourceMask, 8);  

然后我们把这个轮廓在SourceImg上绘制出来，显示出我们将要剪切的区域。参数 'r' 表示红色。

[plain] view plain copy
  1. figure, imshow(SourceImg), axis image  
  2. hold on  
  3. for k = 1:length(SrcBoundry)  
  4.     boundary = SrcBoundry{k};  
  5.     plot(boundary(:,2), boundary(:,1), 'r', 'LineWidth', 2)  
  6. end  
  7. title('Source image intended area for cutting from');  

所得之结果如下图

Poisson image editing算法实现的Matlab代码解析


设定source将要粘贴在target图中的具体位置，并获取TargetImg的长和宽。

[plain] view plain copy
  1. position_in_target = [10, 225];%xy  
  2. [TargetRows, TargetCols, ~] = size(TargetImg);  


函数find()的作用：b=find(X)，X是一个矩阵，查询非零元素的位置，如果X是一个行向量，则返回一个行向量，否则，返回一个列向量。

[plain] view plain copy
  1. [row, col] = find(SourceMask);  

然后来计算mask框在source图中的大小。
[plain] view plain copy
  1. start_pos = [min(col), min(row)];  
  2. end_pos   = [max(col), max(row)];  
  3. frame_size  = end_pos - start_pos;  

如果在position_in_target的位置放置frame将超出Target图的范围，则改变position_in_target，以保证frame不会超出Target图的范围。
[plain] view plain copy
  1. if (frame_size(1) + position_in_target(1) > TargetCols)  
  2.     position_in_target(1) = TargetCols - frame_size(1);  
  3. end  
  4.   
  5. if (frame_size(2) + position_in_target(2) > TargetRows)  
  6.     position_in_target(2) = TargetRows - frame_size(2);  
  7. end  

构建一个大小与Target图相等的新的Mask，然后在position_in_target的位置放入SourceMask。

[plain] view plain copy
  1. MaskTarget = zeros(TargetRows, TargetCols);  
  2. MaskTarget(sub2ind([TargetRows, TargetCols], row - start_pos(2) + position_in_target(2), ...  
  3.  col - start_pos(1) + position_in_target(1))) = 1;  

此时我们得到的新MaskTarget如下图所示

Poisson image editing算法实现的Matlab代码解析


同前面一样，获取二值图像中对象的轮廓，然后把这个轮廓在TargetImg上绘制出来。

[plain] view plain copy
  1. TargBoundry = bwboundaries(MaskTarget, 8);  
  2. figure, imshow(TargetImg), axis image  
  3. hold on  
  4. for k = 1:length(TargBoundry)  
  5.     boundary = TargBoundry{k};  
  6.     plot(boundary(:,2), boundary(:,1), 'r', 'LineWidth', 1)  
  7. end  
  8. title('Target Image with intended place for pasting Source');  

结果如下图所示

Poisson image editing算法实现的Matlab代码解析


根据文章所给出的算法我们知道，对于Mask轮廓的内部，我们是不考虑边界项的，此时计算梯度（G）的散度div，就是执行一个拉普拉斯算子，如下

div ( G( Source(x,y) ) )  = -4 f(x,y) + f(x-1,y) + f(x,y-1) + f(x+1,y) + f(x,y+1)

对SourceImg执行拉普拉斯算子，然后提取结果的R、G、B三个分量。

[plain] view plain copy
  1. templt = [0 -1 0; -1 4 -1; 0 -1 0];  
  2. LaplacianSource = imfilter(double(SourceImg), templt, 'replicate');  
  3. VR = LaplacianSource(:, :, 1);  
  4. VG = LaplacianSource(:, :, 2);  
  5. VB = LaplacianSource(:, :, 3);  

然后根据Mask，把上述计算结果贴入TargetImg。

[plain] view plain copy
  1. TargetImgR = double(TargetImg(:, :, 1));  
  2. TargetImgG = double(TargetImg(:, :, 2));  
  3. TargetImgB = double(TargetImg(:, :, 3));  
  4.   
  5. TargetImgR(logical(MaskTarget(:))) = VR(SourceMask(:));  
  6. TargetImgG(logical(MaskTarget(:))) = VG(SourceMask(:));  
  7. TargetImgB(logical(MaskTarget(:))) = VB(SourceMask(:));  
  8.   
  9. TargetImgNew = cat(3, TargetImgR, TargetImgG, TargetImgB);  
  10. figure, imagesc(uint8(TargetImgNew)), axis image, title('Target image with laplacian of source inserted');  

结果如下图所示


Poisson image editing算法实现的Matlab代码解析


然后，我们要来计算稀疏的临接矩阵。（这个地方需要参考前面给出的博文1）

如果简单从代码来分析，假设我们现在有这样一个mask

[plain] view plain copy
  1. mask =  
  2.   
  3.      0     0     0     0     0  
  4.      0     1(1)  1(2)  0     0  
  5.      0     1(3)  1(4)  1(5)  0  
  6.      0     0     1(6)  1(7)  0  
  7.      0     0     0     0     0  

那么我们就要建立一个7×7的邻接矩阵，例如(1)是(2)邻接的，所以下面矩阵中[1,2]和[2,1]的位置就1。


[plain] view plain copy
  1. >> neighbors = calcAdjancency( mask );  
  2. >> full(neighbors)  
  3.   
  4. ans =  
  5.   
  6.      0     1     1     0     0     0     0  
  7.      1     0     0     1     0     0     0  
  8.      1     0     0     1     0     0     0  
  9.      0     1     1     0     1     1     0  
  10.      0     0     0     1     0     0     1  
  11.      0     0     0     1     0     0     1  
  12.      0     0     0     0     1     1     0  

下面我们给出函数calcAdjancency()的实现代码：


[plain] view plain copy
  1. function neighbors = calcAdjancency( Mask )  
  2.   
  3. [height, width]      = size(Mask);  
  4. [row_mask, col_mask] = find(Mask);  
  5.   
  6. neighbors = sparse(length(row_mask), length(row_mask), 0);  
  7.   
  8. %下标转索引  
  9. roi_idxs = sub2ind([height, width], row_mask, col_mask);  
  10.   
  11. for k = 1:size(row_mask, 1),  
  12.     %4 邻接点  
  13.     connected_4 = [row_mask(k), col_mask(k)-1;%left  
  14.                    row_mask(k), col_mask(k)+1;%right  
  15.                    row_mask(k)-1, col_mask(k);%top  
  16.                    row_mask(k)+1, col_mask(k)];%bottom  
  17.   
  18.     ind_neighbors = sub2ind([height, width], connected_4(:, 1), connected_4(:, 2));  
  19.       
  20.     %二分查找函数 i = ismembc2(t, X)：返回t在X中的位置，其中X必须为递增的的数值向量  
  21.     for neighbor_idx = 1: 4, %number of neighbors,  
  22.         adjacent_pixel_idx =  ismembc2(ind_neighbors(neighbor_idx), roi_idxs);  
  23.         if (adjacent_pixel_idx ~= 0)  
  24.             neighbors(k, adjacent_pixel_idx) = 1;  
  25.         end  
  26.     end   
  27. end  
  28. end  

上述代码中第一个应该知道的地方是Matlab的二分查找函数ismembc2()，这一点我的注释已经比较明确，不再赘言。另一个地方是下标转索引的函数sub2ind()。下面这个例子说明了它的用法和意义：

[plain] view plain copy
  1. % % e.g. M =  
  2. %   
  3. %     11(1) 12(5) 13    14  
  4. %     21    22(6) 23    24  
  5. %     31    32(7) 33    34  
  6. %     41    42    43    44  
  7.   
  8. % roi_idxs = sub2ind(size(M), [1, 1, 2, 3], [2, 1, 2, 2])  
  9. %   
  10. % roi_idxs =  
  11. %   
  12. %      5     1     6     7  

回到我们的主干程序，我们调用函数calcAdjancency()来计算MaskTarget 中Ω区域的邻接矩阵。
[plain] view plain copy
  1. AdjacencyMat = calcAdjancency( MaskTarget );  

然后我们调用PoissonSolver()函数分别对彩色图像的R、G、B三个分量解线性方程组。其核心就是使用迭代法求解型如Ax=b这样的大型稀疏线性方程组。Patrick Perez在原文中使用了高斯-塞德尔迭代法。与此类似，但我们直接调用Matlab中共轭梯度法函数来求解或许更为方便。参考Matlab中给出的帮助信息如下：
X = CGS(A,B) attempts to solve the system of linear equations A*X=B for X.
The N*N coefficient matrix A must be square and the right hand side column vector B must have length N.


[plain] view plain copy
  1. ResultImgR = PoissonSolver(TargetImgR, MaskTarget, AdjacencyMat, TargBoundry);  
  2. ResultImgG = PoissonSolver(TargetImgG, MaskTarget, AdjacencyMat, TargBoundry);  
  3. ResultImgB = PoissonSolver(TargetImgB, MaskTarget, AdjacencyMat, TargBoundry);  

最后我们将三个分量合到一起，并显示出最终的融合结果。

[plain] view plain copy
  1. ResultImg = cat(3, ResultImgR, ResultImgG, ResultImgB);  
  2.   
  3. figure;  
  4. imshow(uint8(ResultImg));  

如下图所示：

Poisson image editing算法实现的Matlab代码解析


更多示例与讨论

在Matlab中建立稀疏矩阵如果简单的用 zeros()或ones()函数的话，当图片稍微大一点（例如我在做下面的手眼融合示例时所用的图），就会产生 out of memory的问题。所以程序中要特别注意，使用sparse()等等专门用于建立稀疏矩阵的函数来避免超内存的问题。

Poisson image editing算法实现的Matlab代码解析


全文完。

===========

更多有趣有用的图像处理算法还可以参考我的《数字图像处理原理与实践（Matlab版）》


Poisson image editing算法实现的Matlab代码解析


素材是三张图片，如下


                                                                                   bear.jpg                                                              bear-mask.jpg

Poisson image editing算法实现的Matlab代码解析                         Poisson image editing算法实现的Matlab代码解析


pool-target.jpg

Poisson image editing算法实现的Matlab代码解析


首先，我们清理一下Matlab的环境：

[plain] view plain copy
  1. close all;  
  2. clear;  
  3. clc;  

然后，读入三张图片（注意需要把mask图以二值图的形式读入）
[plain] view plain copy
  1. TargetImg   = imread('pool-target.jpg');  
  2. SourceImg   = imread('bear.jpg');  
  3. SourceMask  = im2bw(imread('bear-mask.jpg'));  

用函数bwboundaries(W,CONN) 来获取二值图中对象的轮廓，其中CONN 为8或4，指示联通性采用4方向邻域点判别还是8方向邻域点判别，默认为8。

[plain] view plain copy
  1. [SrcBoundry, ~] = bwboundaries(SourceMask, 8);  

然后我们把这个轮廓在SourceImg上绘制出来，显示出我们将要剪切的区域。参数 'r' 表示红色。

[plain] view plain copy
  1. figure, imshow(SourceImg), axis image  
  2. hold on  
  3. for k = 1:length(SrcBoundry)  
  4.     boundary = SrcBoundry{k};  
  5.     plot(boundary(:,2), boundary(:,1), 'r', 'LineWidth', 2)  
  6. end  
  7. title('Source image intended area for cutting from');  

所得之结果如下图

Poisson image editing算法实现的Matlab代码解析


设定source将要粘贴在target图中的具体位置，并获取TargetImg的长和宽。

[plain] view plain copy
  1. position_in_target = [10, 225];%xy  
  2. [TargetRows, TargetCols, ~] = size(TargetImg);  


函数find()的作用：b=find(X)，X是一个矩阵，查询非零元素的位置，如果X是一个行向量，则返回一个行向量，否则，返回一个列向量。

[plain] view plain copy
  1. [row, col] = find(SourceMask);  

然后来计算mask框在source图中的大小。
[plain] view plain copy
  1. start_pos = [min(col), min(row)];  
  2. end_pos   = [max(col), max(row)];  
  3. frame_size  = end_pos - start_pos;  

如果在position_in_target的位置放置frame将超出Target图的范围，则改变position_in_target，以保证frame不会超出Target图的范围。
[plain] view plain copy
  1. if (frame_size(1) + position_in_target(1) > TargetCols)  
  2.     position_in_target(1) = TargetCols - frame_size(1);  
  3. end  
  4.   
  5. if (frame_size(2) + position_in_target(2) > TargetRows)  
  6.     position_in_target(2) = TargetRows - frame_size(2);  
  7. end  

构建一个大小与Target图相等的新的Mask，然后在position_in_target的位置放入SourceMask。

[plain] view plain copy
  1. MaskTarget = zeros(TargetRows, TargetCols);  
  2. MaskTarget(sub2ind([TargetRows, TargetCols], row - start_pos(2) + position_in_target(2), ...  
  3.  col - start_pos(1) + position_in_target(1))) = 1;  

此时我们得到的新MaskTarget如下图所示

Poisson image editing算法实现的Matlab代码解析


同前面一样，获取二值图像中对象的轮廓，然后把这个轮廓在TargetImg上绘制出来。

[plain] view plain copy
  1. TargBoundry = bwboundaries(MaskTarget, 8);  
  2. figure, imshow(TargetImg), axis image  
  3. hold on  
  4. for k = 1:length(TargBoundry)  
  5.     boundary = TargBoundry{k};  
  6.     plot(boundary(:,2), boundary(:,1), 'r', 'LineWidth', 1)  
  7. end  
  8. title('Target Image with intended place for pasting Source');  

结果如下图所示

Poisson image editing算法实现的Matlab代码解析


根据文章所给出的算法我们知道，对于Mask轮廓的内部，我们是不考虑边界项的，此时计算梯度（G）的散度div，就是执行一个拉普拉斯算子，如下

div ( G( Source(x,y) ) )  = -4 f(x,y) + f(x-1,y) + f(x,y-1) + f(x+1,y) + f(x,y+1)

对SourceImg执行拉普拉斯算子，然后提取结果的R、G、B三个分量。

[plain] view plain copy
  1. templt = [0 -1 0; -1 4 -1; 0 -1 0];  
  2. LaplacianSource = imfilter(double(SourceImg), templt, 'replicate');  
  3. VR = LaplacianSource(:, :, 1);  
  4. VG = LaplacianSource(:, :, 2);  
  5. VB = LaplacianSource(:, :, 3);  

然后根据Mask，把上述计算结果贴入TargetImg。

[plain] view plain copy
  1. TargetImgR = double(TargetImg(:, :, 1));  
  2. TargetImgG = double(TargetImg(:, :, 2));  
  3. TargetImgB = double(TargetImg(:, :, 3));  
  4.   
  5. TargetImgR(logical(MaskTarget(:))) = VR(SourceMask(:));  
  6. TargetImgG(logical(MaskTarget(:))) = VG(SourceMask(:));  
  7. TargetImgB(logical(MaskTarget(:))) = VB(SourceMask(:));  
  8.   
  9. TargetImgNew = cat(3, TargetImgR, TargetImgG, TargetImgB);  
  10. figure, imagesc(uint8(TargetImgNew)), axis image, title('Target image with laplacian of source inserted');  

结果如下图所示


Poisson image editing算法实现的Matlab代码解析


然后，我们要来计算稀疏的临接矩阵。（这个地方需要参考前面给出的博文1）

如果简单从代码来分析，假设我们现在有这样一个mask

[plain] view plain copy
  1. mask =  
  2.   
  3.      0     0     0     0     0  
  4.      0     1(1)  1(2)  0     0  
  5.      0     1(3)  1(4)  1(5)  0  
  6.      0     0     1(6)  1(7)  0  
  7.      0     0     0     0     0  

那么我们就要建立一个7×7的邻接矩阵，例如(1)是(2)邻接的，所以下面矩阵中[1,2]和[2,1]的位置就1。


[plain] view plain copy
  1. >> neighbors = calcAdjancency( mask );  
  2. >> full(neighbors)  
  3.   
  4. ans =  
  5.   
  6.      0     1     1     0     0     0     0  
  7.      1     0     0     1     0     0     0  
  8.      1     0     0     1     0     0     0  
  9.      0     1     1     0     1     1     0  
  10.      0     0     0     1     0     0     1  
  11.      0     0     0     1     0     0     1  
  12.      0     0     0     0     1     1     0  

下面我们给出函数calcAdjancency()的实现代码：


[plain] view plain copy
  1. function neighbors = calcAdjancency( Mask )  
  2.   
  3. [height, width]      = size(Mask);  
  4. [row_mask, col_mask] = find(Mask);  
  5.   
  6. neighbors = sparse(length(row_mask), length(row_mask), 0);  
  7.   
  8. %下标转索引  
  9. roi_idxs = sub2ind([height, width], row_mask, col_mask);  
  10.   
  11. for k = 1:size(row_mask, 1),  
  12.     %4 邻接点  
  13.     connected_4 = [row_mask(k), col_mask(k)-1;%left  
  14.                    row_mask(k), col_mask(k)+1;%right  
  15.                    row_mask(k)-1, col_mask(k);%top  
  16.                    row_mask(k)+1, col_mask(k)];%bottom  
  17.   
  18.     ind_neighbors = sub2ind([height, width], connected_4(:, 1), connected_4(:, 2));  
  19.       
  20.     %二分查找函数 i = ismembc2(t, X)：返回t在X中的位置，其中X必须为递增的的数值向量  
  21.     for neighbor_idx = 1: 4, %number of neighbors,  
  22.         adjacent_pixel_idx =  ismembc2(ind_neighbors(neighbor_idx), roi_idxs);  
  23.         if (adjacent_pixel_idx ~= 0)  
  24.             neighbors(k, adjacent_pixel_idx) = 1;  
  25.         end  
  26.     end   
  27. end  
  28. end  

上述代码中第一个应该知道的地方是Matlab的二分查找函数ismembc2()，这一点我的注释已经比较明确，不再赘言。另一个地方是下标转索引的函数sub2ind()。下面这个例子说明了它的用法和意义：

[plain] view plain copy
  1. % % e.g. M =  
  2. %   
  3. %     11(1) 12(5) 13    14  
  4. %     21    22(6) 23    24  
  5. %     31    32(7) 33    34  
  6. %     41    42    43    44  
  7.   
  8. % roi_idxs = sub2ind(size(M), [1, 1, 2, 3], [2, 1, 2, 2])  
  9. %   
  10. % roi_idxs =  
  11. %   
  12. %      5     1     6     7  

回到我们的主干程序，我们调用函数calcAdjancency()来计算MaskTarget 中Ω区域的邻接矩阵。
[plain] view plain copy
  1. AdjacencyMat = calcAdjancency( MaskTarget );  

然后我们调用PoissonSolver()函数分别对彩色图像的R、G、B三个分量解线性方程组。其核心就是使用迭代法求解型如Ax=b这样的大型稀疏线性方程组。Patrick Perez在原文中使用了高斯-塞德尔迭代法。与此类似，但我们直接调用Matlab中共轭梯度法函数来求解或许更为方便。参考Matlab中给出的帮助信息如下：
X = CGS(A,B) attempts to solve the system of linear equations A*X=B for X.
The N*N coefficient matrix A must be square and the right hand side column vector B must have length N.


[plain] view plain copy
  1. ResultImgR = PoissonSolver(TargetImgR, MaskTarget, AdjacencyMat, TargBoundry);  
  2. ResultImgG = PoissonSolver(TargetImgG, MaskTarget, AdjacencyMat, TargBoundry);  
  3. ResultImgB = PoissonSolver(TargetImgB, MaskTarget, AdjacencyMat, TargBoundry);  

最后我们将三个分量合到一起，并显示出最终的融合结果。

[plain] view plain copy
  1. ResultImg = cat(3, ResultImgR, ResultImgG, ResultImgB);  
  2.   
  3. figure;  
  4. imshow(uint8(ResultImg));  

如下图所示：

Poisson image editing算法实现的Matlab代码解析


更多示例与讨论

在Matlab中建立稀疏矩阵如果简单的用 zeros()或ones()函数的话，当图片稍微大一点（例如我在做下面的手眼融合示例时所用的图），就会产生 out of memory的问题。所以程序中要特别注意，使用sparse()等等专门用于建立稀疏矩阵的函数来避免超内存的问题。

Poisson image editing算法实现的Matlab代码解析


全文完。

===========

更多有趣有用的图像处理算法还可以参考我的《数字图像处理原理与实践（Matlab版）》


Poisson image editing算法实现的Matlab代码解析

相关推荐