Poisson image editing算法实现的Matlab代码解析
来源:互联网 发布:利用广告联盟做淘宝客 编辑:程序博客网 时间:2024/05/16 11:26
之前我发了数篇系列博文来仔细研究Poisson Image Editing算法,每次重新审视和深入,仿佛都能有更为深刻的认识和很大的收获。这应该算是我这个系列的完结篇,会用用Matlab代码一点一点的演示,原文作者到底是如何设计和实现他那个强大且影响深远的算法的。希望你在看本文之前务必参考一下文章来了解算法原理,本文将主要讲解编程实现的问题,对于前面讲过的内容,我不会深究。但我个人总体的感觉是,现在图像处理算法对数学的要求是越来越高了,像泊松融合、泊松抠图这样的算法如果没有偏微分方程(本算法中涉及Poisson Equation)、数值计算(需要高斯-塞德尔迭代法或者共轭梯度法)、场论(涉及散度、梯度和拉普拉斯算子)、矩阵论(涉及大型稀疏矩阵线性方程组)、数学物理方法(涉及带狄利克雷条件的椭圆曲线方程)、最优化理论(涉及到了欧拉-朗格拉日变分法)这些艰深的数学知识,要理解起来实在太困难了(每一步都是一个坎,我感觉有一个地方过不了,那篇经典论文的精华就无法领会,更别说编程实现了)。
图像处理之让手心长出眼睛,其实嘴也可以~
图像的泊松(Poisson)编辑、泊松融合完全详解(3)
图像的泊松(Poisson)编辑、泊松融合完全详解(2)
图像的泊松(Poisson)编辑、泊松融合完全详解(1)
素材是三张图片,如下
bear.jpg bear-mask.jpg
pool-target.jpg
首先,我们清理一下Matlab的环境:
close all;clear;clc;
然后,读入三张图片(注意需要把mask图以二值图的形式读入)
TargetImg = imread('pool-target.jpg');SourceImg = imread('bear.jpg');SourceMask = im2bw(imread('bear-mask.jpg'));
用函数bwboundaries(W,CONN) 来获取二值图中对象的轮廓,其中CONN 为8或4,指示联通性采用4方向邻域点判别还是8方向邻域点判别,默认为8。
[SrcBoundry, ~] = bwboundaries(SourceMask, 8);
然后我们把这个轮廓在SourceImg上绘制出来,显示出我们将要剪切的区域。参数 'r' 表示红色。
figure, imshow(SourceImg), axis imagehold onfor k = 1:length(SrcBoundry) boundary = SrcBoundry{k}; plot(boundary(:,2), boundary(:,1), 'r', 'LineWidth', 2)endtitle('Source image intended area for cutting from');
所得之结果如下图
设定source将要粘贴在target图中的具体位置,并获取TargetImg的长和宽。
position_in_target = [10, 225];%xy[TargetRows, TargetCols, ~] = size(TargetImg);
函数find()的作用:b=find(X),X是一个矩阵,查询非零元素的位置,如果X是一个行向量,则返回一个行向量,否则,返回一个列向量。
[row, col] = find(SourceMask);
然后来计算mask框在source图中的大小。
start_pos = [min(col), min(row)];end_pos = [max(col), max(row)];frame_size = end_pos - start_pos;
如果在position_in_target的位置放置frame将超出Target图的范围,则改变position_in_target,以保证frame不会超出Target图的范围。
if (frame_size(1) + position_in_target(1) > TargetCols) position_in_target(1) = TargetCols - frame_size(1);endif (frame_size(2) + position_in_target(2) > TargetRows) position_in_target(2) = TargetRows - frame_size(2);end
构建一个大小与Target图相等的新的Mask,然后在position_in_target的位置放入SourceMask。
MaskTarget = zeros(TargetRows, TargetCols);MaskTarget(sub2ind([TargetRows, TargetCols], row - start_pos(2) + position_in_target(2), ... col - start_pos(1) + position_in_target(1))) = 1;
此时我们得到的新MaskTarget如下图所示
同前面一样,获取二值图像中对象的轮廓,然后把这个轮廓在TargetImg上绘制出来。
TargBoundry = bwboundaries(MaskTarget, 8);figure, imshow(TargetImg), axis imagehold onfor k = 1:length(TargBoundry) boundary = TargBoundry{k}; plot(boundary(:,2), boundary(:,1), 'r', 'LineWidth', 1)endtitle('Target Image with intended place for pasting Source');
结果如下图所示
根据文章所给出的算法我们知道,对于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三个分量。templt = [0 -1 0; -1 4 -1; 0 -1 0];LaplacianSource = imfilter(double(SourceImg), templt, 'replicate');VR = LaplacianSource(:, :, 1);VG = LaplacianSource(:, :, 2);VB = LaplacianSource(:, :, 3);
然后根据Mask,把上述计算结果贴入TargetImg。
TargetImgR = double(TargetImg(:, :, 1));TargetImgG = double(TargetImg(:, :, 2));TargetImgB = double(TargetImg(:, :, 3));TargetImgR(logical(MaskTarget(:))) = VR(SourceMask(:));TargetImgG(logical(MaskTarget(:))) = VG(SourceMask(:));TargetImgB(logical(MaskTarget(:))) = VB(SourceMask(:));TargetImgNew = cat(3, TargetImgR, TargetImgG, TargetImgB);figure, imagesc(uint8(TargetImgNew)), axis image, title('Target image with laplacian of source inserted');
结果如下图所示
然后,我们要来计算稀疏的临接矩阵。(这个地方需要参考前面给出的博文1)
如果简单从代码来分析,假设我们现在有这样一个mask
mask = 0 0 0 0 0 0 1(1) 1(2) 0 0 0 1(3) 1(4) 1(5) 0 0 0 1(6) 1(7) 0 0 0 0 0 0
那么我们就要建立一个7×7的邻接矩阵,例如(1)是(2)邻接的,所以下面矩阵中[1,2]和[2,1]的位置就1。
>> neighbors = calcAdjancency( mask );>> full(neighbors)ans = 0 1 1 0 0 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 0 0 1 1 0 1 1 0 0 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 0 0 1 1 0
下面我们给出函数calcAdjancency()的实现代码:
function neighbors = calcAdjancency( Mask )[height, width] = size(Mask);[row_mask, col_mask] = find(Mask);neighbors = sparse(length(row_mask), length(row_mask), 0);%下标转索引roi_idxs = sub2ind([height, width], row_mask, col_mask);for k = 1:size(row_mask, 1), %4 邻接点 connected_4 = [row_mask(k), col_mask(k)-1;%left row_mask(k), col_mask(k)+1;%right row_mask(k)-1, col_mask(k);%top row_mask(k)+1, col_mask(k)];%bottom ind_neighbors = sub2ind([height, width], connected_4(:, 1), connected_4(:, 2)); %二分查找函数 i = ismembc2(t, X):返回t在X中的位置,其中X必须为递增的的数值向量 for neighbor_idx = 1: 4, %number of neighbors, adjacent_pixel_idx = ismembc2(ind_neighbors(neighbor_idx), roi_idxs); if (adjacent_pixel_idx ~= 0) neighbors(k, adjacent_pixel_idx) = 1; end end endend
上述代码中第一个应该知道的地方是Matlab的二分查找函数ismembc2(),这一点我的注释已经比较明确,不再赘言。另一个地方是下标转索引的函数sub2ind()。下面这个例子说明了它的用法和意义:
% % e.g. M =% % 11(1) 12(5) 13 14% 21 22(6) 23 24% 31 32(7) 33 34% 41 42 43 44% roi_idxs = sub2ind(size(M), [1, 1, 2, 3], [2, 1, 2, 2])% % roi_idxs =% % 5 1 6 7
回到我们的主干程序,我们调用函数calcAdjancency()来计算MaskTarget 中Ω区域的邻接矩阵。
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.
ResultImgR = PoissonSolver(TargetImgR, MaskTarget, AdjacencyMat, TargBoundry);ResultImgG = PoissonSolver(TargetImgG, MaskTarget, AdjacencyMat, TargBoundry);ResultImgB = PoissonSolver(TargetImgB, MaskTarget, AdjacencyMat, TargBoundry);
最后我们将三个分量合到一起,并显示出最终的融合结果。
ResultImg = cat(3, ResultImgR, ResultImgG, ResultImgB);figure;imshow(uint8(ResultImg));
如下图所示:
更多示例与讨论
在Matlab中建立稀疏矩阵如果简单的用 zeros()或ones()函数的话,当图片稍微大一点(例如我在做下面的手眼融合示例时所用的图),就会产生 out of memory的问题。所以程序中要特别注意,使用sparse()等等专门用于建立稀疏矩阵的函数来避免超内存的问题。
全文完。
===========
如果你是图像处理的同道中人,,欢迎加入图像处理算法交流群(单击链接查看群号)。更多有趣有用的图像处理算法还可以参考我的《数字图像处理原理与实践(Matlab版)》
- Poisson image editing算法实现的Matlab代码解析
- Poisson Image Editing
- 泊松图像编辑(Poisson image editing)
- Poisson Image Editing & Texture Based Terrain Synthesis
- 泊松图像编辑(Poisson image editing)
- 泊松图像编辑(Poisson image editing)
- 泊松图像编辑(poisson Image Editing)
- 图像融合之泊松编辑(Poisson Editing)(2):详解算法和实现
- 图像无缝融合研究(一)Poisson Image Editing
- Poisson Mesh Editing
- opencv3.0中的无缝克隆图像——seamless_cloning(Poisson Image Editing)
- 图像融合之泊松编辑(Poisson Editing)(1):简略语言概述算法
- 图像融合之泊松编辑(Poisson Editing)(1):简略语言概述算法
- LDA算法-matlab代码实现
- KNN算法matlab代码实现
- Apriori算法,MATLAB代码实现
- Kmeans及其改进算法的Matlab解析与实现
- 分水岭算法MATLAB编程代码解析
- Context initialization failed
- 【Java SE】————如何做到跨平台性的?
- C/C++log日志库比较
- [LeetCode] 315
- xcode 下载更新地址
- Poisson image editing算法实现的Matlab代码解析
- Java常见错误汇总
- Oracle与Tomcat端口冲突,报错:Protocol handler initialization failed;Socket bind failed
- USACO Your Ride Is Here 解题日志
- UniCode 下 CString 转 char* 的方法
- Julia: 趣!,13579分别在一本168页书的页码中出现的次数
- hdu(1171)——Big Event in HDU
- Leetcode 270Closest Binary Search Tree Value
- JavaScript优化(三)