matlab 狄利克雷函数图像,Poisson image editing算法实现的Matlab代码解析
之前我发了数篇系列博文来仔细研究Poisson Image Editing算法,每次重新审视和深入,仿佛都能有更为深刻的认识和很大的收获。这应该算是我这个系列的完结篇,会用用Matlab代码一点一点的演示,原文作者到底是如何设计和实现他那个强大且影响深远的算法的。希望你在看本文之前务必参考一下文章来了解算法原理,本文将主要讲解编程实现的问题,对于前面讲过的内容,我不会深究。但我个人总体的感觉是,现在图像处理算法对数学的要求是越来越高了,像泊松融合、泊松抠图这样的算法如果没有偏微分方程(本算法中涉及Poisson Equation)、数值计算(需要高斯-塞德尔迭代法或者共轭梯度法)、场论(涉及散度、梯度和拉普拉斯算子)、矩阵论(涉及大型稀疏矩阵线性方程组)、数学物理方法(涉及带狄利克雷条件的椭圆曲线方程)、最优化理论(涉及到了欧拉-朗格拉日变分法)这些艰深的数学知识,要理解起来实在太困难了(每一步都是一个坎,我感觉有一个地方过不了,那篇经典论文的精华就无法领会,更别说编程实现了
)。
素材是三张图片,如下
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 image
hold on
for k = 1:length(SrcBoundry)
boundary = SrcBoundry{k};
plot(boundary(:,2), boundary(:,1), 'r', 'LineWidth', 2)
end
title('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);
end
if (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 image
hold on
for k = 1:length(TargBoundry)
boundary = TargBoundry{k};
plot(boundary(:,2), boundary(:,1), 'r', 'LineWidth', 1)
end
title('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
end
end
上述代码中第一个应该知道的地方是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()等等专门用于建立稀疏矩阵的函数来避免超内存的问题。
全文完。
如果你是图像处理的同道中人,欢迎加入图像处理学习群(529549320)。
更多有趣有用的图像处理算法还可以参考我的《数字图像处理原理与实践(Matlab版)》
matlab 狄利克雷函数图像,Poisson image editing算法实现的Matlab代码解析相关推荐
- Poisson image editing算法实现的Matlab代码解析
之前我发了数篇系列博文来仔细研究Poisson Image Editing算法,每次重新审视和深入,仿佛都能有更为深刻的认识和很大的收获.这应该算是我这个系列的完结篇,会用用Matlab代码一点一点的 ...
- Matlab|绘制函数图像
欢迎点击「算法与编程之美」↑关注我们! 本文首发于微信公众号:"算法与编程之美",欢迎关注,及时了解更多此系列文章. 欢迎加入团队圈子!与作者面对面!直接点击! 一.绘制图像的常用 ...
- matlab一般函数的绘制方法,基于MATLAB的函数图像绘制方法
C DOI:10.16707~.cnki.fjpc.2017.01.084 E 晒 亍嚣 基于 MATLAB的函数图像绘制方法 张笑笑 一,童 键 z (1湖南省长沙市第一中学 湖南 长沙 410() ...
- qr-rls算法matlab实现,【预测模型】基于RLS算法进行预测matlab源码
一.简介 1 概述 递归最小二乘(RLS)算法是一种典型的数据处理方法,由著名学者高斯在1795年提出,高斯认为,根据所获得的观测数据来推断未知参数时,未知参数最可能的值是这样一个数据,即它使各项实际 ...
- 视觉SLAM开源算法ORB-SLAM3 原理与代码解析
来源:深蓝学院,文稿整理者:何常鑫,审核&修改:刘国庆 本文总结于上交感知与导航研究所科研助理--刘国庆关于[视觉SLAM开源算法ORB-SLAM3 原理与代码解析]的公开课. ORB-SLA ...
- 基于单层决策树的adaBoost算法思想分析和源代码解析
基于单层决策树的AdaBoost算法思想分析和源代码解析 前言: 上一篇SVM可是废了我好鼻子劲,这一篇咱们来点愉快的东西.我们一定听说过这句俗语:"三个臭皮匠,顶个诸葛亮!" 大 ...
- ICCV2017跟踪算法BACF原理及代码解析
文章和代码下载地址: Galoogahi H K, Fagg A, Lucey S. Learning Background-Aware Correlation Filters for Visual ...
- TPAMI2015跟踪算法KCF原理及代码解析
文章和代码下载地址: http://www.robots.ox.ac.uk/~joao/circulant/ 一.基础公式的由来 相关滤波跟踪器可以表示为最小化的岭回归问题: 表示期望相应,表示正则系 ...
- 利用matlab绘制函数图像
文章目录 一.下载matlab.rar 二.解压缩matlab.rar 三.启动matlab窗口 四.绘制一元函数图像 - 直线或曲线 1.绘制直线 2.绘制曲线 五.绘制二元函数图像 - 平面或曲面 ...
最新文章
- Linux C连接Mysql
- make 编译可执行
- Failed to resolve hostname 192: The name does not resolve for the supplied parameters
- qt 启动画面显示图片_Qt程序起动画面QSplashScreen
- ASP.NET 26个常用性能优化方法
- SQL Server代理(11/12):维护计划作业
- html5离开网页自动暂停,通过html5代码在网页中实现播放和暂停音乐mp3,mav等文
- java 累加 0-9 a-z_JAVA获得包含0-9、a-z、A-Z范围内字符串的的随机数实例
- 摄像头分辨率怎么调整_手机摄像头测试:细数手机摄像头由单摄到多摄有哪些变化...
- Bella Email邮件发送模板
- 研发解决方案介绍#Tracing(鹰眼)
- [置顶] Embedded Server:像写main函数一样写Web Server
- cannot resolve method ‘println(java.lang.String)
- Perhaps JAVA_HOME does not point to the JDK的解决方法
- Keil5下载及安装
- cad图纸服务器共享文件慢,DWG文件打开慢?3个技巧教你实现快速预览!
- JIRA上根据前置任务自动计算到期日之automation实现实例
- 计算机硬件检测与维修理论试题,计算机硬件检测与维修试题10.doc
- 京东到家话费券系统NIO实战
- 圣天诺 加密java_圣天诺Sentinel LDK 7.8壳加密的编译环境是什么?