Poisson Image Editing 论文方法实现

发表于 2021-05-10 分类于记录，论文阅读

“Poisson Image Editing”论文方法实现

在上上篇博文中（现在是上上上了），我介绍了一下这篇图像融合的经典论文。今天，我将记录一下这篇论文的C++代码实现。我已经三个多月没碰C++了，手很生。但是，我会凭借着我高超的编程底力和过人的天赋，三小时内完成代码实现。稍微对工作量有一点了解的人都会产生疑问：“你在开玩笑吧？论文实现不比写算法题，你要先去看懂论文，再调一调别人的库函数，写代码之前有很多准备工作要做。你这不可能完成。”肯定有人这样想。但是，我可以说，包括学习前置知识在内，我可以三小时完成实现。这是一个七年编程王者应该有的自信。这篇博文的正文和会按之前的格式对项目的实现做一个比较全面的记录，最后一个章节会写下我完成此项目的实况。

代码仓库：https://github.com/SingleZombie/Gradient-Domain-Image-Processing-Cpp

由于种种原因，在记录了编程的实况后，我就把这篇博文搁置下来了。但马上要写毕业论文了，这篇博文和论文的内容有很多重合之处，我打算先把博文写完。

知识准备

在代码实现前，我们先整理一下实现方法的整体思路，再提取出实现中涉及的知识点，对每个知识点有关的实现方法和具体的实现技术进行介绍，最后对方法的核心——结果图像求解进行一个详细的介绍。

方法思想

整个方法的总流程图如上。和普通的图像复制类似，方法需要输入一幅源图像、一幅目标图像及复制区域（步骤1），输出一幅融合好的图像（步骤4）。为了得到最终的图像，需要先计算出源图像的梯度（步骤2），再根据目标图像在复制区域边缘的像素值和源图像复制区域的梯度值对结果图像求解（步骤3）。实际上，其他步骤都十分简单，方法的关键就在于第三步结果图像求解。

方法中设计的知识点有：

图像输入/输出：这项操作的内容很显然：把操作系统中的图像文件读入到程序的一个数组中，再把一个数组输出成一个图像文件。这个可以通过OpenCV库来轻松实现。图像输入输出是OpenCV最基本的操作，网上随便搜一搜OpenCV的教程就能找到。我自己也有一篇博文讲了这项知识。
图像梯度计算：所谓图像的梯度，就是每个像素与它左边和上边像素的颜色值之差，某处的梯度值表示此处颜色变化速度的快慢。由于后续操作需要源图像的梯度，在实际进行结果图像求解前需要先预处理出图像梯度。实现的时候用一个滤波器对图像做一次滤波即可，OpenCV包括了滤波功能。我的这篇博客简单介绍了OpenCV滤波函数fillter2D的用法。
结果图像求解：图像融合方法的思想是在融合区域边缘的颜色和目标图像一样的前提下，让源图像的梯度尽可能不变。这是一个最优化问题，最后问题可以转换成求解线性方程组。因此，这一步的目的就是通过求解线性方程组来得到我们想要的结果图像。这是本程序最关键、困难的一步，C++的Eigen库提供了高性能的矩阵运算函数。
图像拼接：最后输出前，要把处理过的图像区域和整幅目标图像拼起来。OpenCV提供了方便的图像区域覆盖函数。

本程序设计的图像处理操作都十分基础，用一些简单的OpenCV库函数即可。如前面多次强调的，该程序的关键是第三步。下一节将对第三步进行详细介绍。

结果图像的问题建立与求解

(导向插值标记，来自论文[1])

一般来说，用公式描述一件事是方便描述者，而折磨倾听者的。但是，为了把问题准确无误的描述出来，有时不得不使用严谨的公式标记。

在此问题中，我们把图像看成一个函数。如图所示， $S$ 是一个二维点的集合，即一堆可以用 $(x, y)$ 这样的坐标描述的点。集合 $S$ 就是图像函数的定义域。图像的值域呢？自然是图像的像素值了，这取决于实际情况，比如RGB像素值的范围是 $[0, 2^{24})$ 。

现在，理解了图像其实可以表达为定义域是二维点，值域是颜色空间的值域的一个函数后，就可以继续理解符号标记了。令 $g$ 是源图像， $f^{*}$ 是目标图像， $f$ 是我们把源图像复制到目标区域后，经过图像融合，得到的结果图像。按上一节的话说， $g$ 和 $f^{*}$ 是输入， $f$ 是输出的结果图像。我们还有一个输入，就是进行图像融合的区域 $Ω$ 。既然我们把图像看成了函数，那么区域 $Ω$ 就是整个定义域 $S$ 的一个子集。

$v$ 是 $g$ 的梯度， $\partial Ω$ 是区域 $Ω$ 的梯度，这两个量都是计算得到的。

我们还要引入一个符号——梯度算子 $▽$ 。 $v$ 是 $g$ 的梯度，可以写成 $v = ▽ g$ 。

重申一下结果图像求解的思想：在融合区域边缘的颜色和目标图像一样的前提下，让源图像的梯度尽可能不变（最小化差值）。用数学语言就是：

如果理解了求解的思想，那么看懂，或者说推导出这个公式是十分简单的。

这个公式是理想情况：图像在一个连续的二维数集上定义。但实际上，我们图像是离散是，只有在整数的位置处有值。我们要把问题用离散的形式表达出来。这里又要提出一些标记，对于 $S$ 中的每个像素 $p$ ，令 $N_{p}$ 为其4邻域中仍在 $S$ 的像素集合。也就是说，这个符号表示的某个像素上下左右这四个像素的集合，在这些像素没有跑出整幅图像的前提下。令 $(p, q)$ 为一个像素对，其中 $q \in N_{p}$ 。

再令 $v_{p q}$ 表示 $p q$ 两点间的像素值之差，即p到q方向的梯度值。则上面那个离散化的式子可以被转换成：

能看懂之前的式子，看懂新的式子也不难。

经过一系列数学推导，上面这个离散域的最优化问题可以被转换成一个线性方程组：

这个方程组就比较难懂了。不过，可以用一个直观的方法描述出来：待求解的结果图像的每个像素与周围四个像素的像素值之差（梯度）的求和等于源图像与周围四个点的像素值之差（梯度）的求和。公式中结果图像的每个像素与周围的梯度的求和分了两部分，是因为如果这个像素来自图像边界，则这个值是指定的（来自约束条件，边界像素值等于目标图像的像素值），在等式右边；如果像素不来自边界，则这个值是未知的，在等式左边。

有了联立的线性方程组，问题就变成了纯粹的数学问题，有很多的工具和方法来求解。

程序实现

程序结构图

在理想情况下，程序的结构图如下：（目前的代码很乱，没有严格按照这个图来）

整体分成输入、输出、和图像处理模块，非常简明。

如方法思想中所介绍的，图像处理可以分成图像梯度计算和目标图像求解。而求解又可以分成方程组建立和方程组求解。

这里仅介绍图像处理模块的实现细节。该模块整体的伪代码如下：

这份伪代码只能算是把整个流程用英文表达出来了而已，前文的梯度计算、方程组建立、方程组求解等内容均有体现。接下来几节介绍这三个重要步骤的细节。其他诸如获取图像的一个区域、填充图形一个区域的操作十分简单，可以用OpenCV轻松实现。OpenCV用Rect表示一个矩形区域，假设其有实例rect。令mat为任意一个图像，则mat[rect]就是图像的一个区域，可以执行读写操作。

图像梯度计算

梯度计算是学数字图像处理时的一个基本作业，由于有库函数的帮助，实现起来很方便，其伪代码如下：

拉普拉斯滤波器就是

1
2
3

[[0, -1, 0],
[-1, 4, -1],
[0, -1, 0]]

这样一个3X3的滤波器。手动建立一个滤波器后，调用OpenCV的filter2D就可以完成滤波了。

方程组建立

问题的方程组就是“结果图像的问题建立与求解”一节中最后一个公式所表示的方程组。只有彻底理解了那个公式中每一项的由来，才能把方程组建好。方程组求解时用到了矩阵，这其实就是把方程组的所有有关参数塞入了一个二维的存储结构中，没有太深奥的东西。

方程组主要有3项：系数矩阵项，边缘像素值项和梯度项。系数矩阵是方程组的左端项，后两者之和就是方程组的右端项。这一步的伪代码如下：

系数矩阵项lhs只由融合区域的大小决定。设矩阵的第i行表示第i个像素有关的方程，则第i列就是表示i像素自己的系数，由公式可知系数值是 $| N_{p} |$ ；而其他至多4个在这个像素上下左右，且在融合区域内部的像素的系数是 $- 1$ 。

边缘像素值来自原问题的限制条件：融合区域的边缘像素值等于目标图像在此处的像素值。因此，这一部分的值要根据目标图像及融合区域得到。由于边缘项的出现没有什么规律，实际实现时可以在计算系数矩阵的同时计算边缘项：对每个像素四周判断，如果这一项在边缘上，则更新返回的边缘项；否则更新系数矩阵。

梯度项则完全来自源图像区域的梯度，甚至每个元素的位置都没有变，只要把它的形状改成列数为1的矩阵。

方程组求解

有了所有参数，调库求解方程组就是几行的事情了。下面是Eigen求解Ax=B方程组的代码：

Eigen::SparseLU<Eigen::SparseMatrix<float>> solver;
solver.analyzePattern(A);
solver.factorize(A);
Eigen::VectorXf tmpRes = solver.solve(b);Eigen::SparseLU<Eigen::SparseMatrix<float>> solver;
solver.analyzePattern(A);
solver.factorize(A);
Eigen::VectorXf tmpRes = solver.solve(b);

我只有在这一部分贴了代码，因为其他步骤都是比较灵活的，每个人都有自己的实现方式。只有这一步调库的写法是固定的。

求解完方程组，程序其实就基本写完了。

UPD：SparseLU并不是求解此问题的最佳方法。SimplicialLDLT适用于正定对称矩阵，在此问题中有更高的性能（我并不知道其中的原因，但是官网上明确写出这个方法适用于2D泊松问题，我就把结论拿来用了，哈哈）。只要把上述代码的SparseLU替换掉即可。

Eigen::SimplicialLDLT<Eigen::SparseMatrix<float>> solver;
solver.analyzePattern(A);
solver.factorize(A);
Eigen::VectorXf tmpRes = solver.solve(b);Eigen::SparseLU<Eigen::SparseMatrix<float>> solver;
solver.analyzePattern(A);
solver.factorize(A);
Eigen::VectorXf tmpRes = solver.solve(b);

结果展示

如图，(a)是源图像，(b)是目标图像，(c)是直接复制的结果，(d)是图像融合的结果。从视觉上来看，图像融合的效果还是不错的。如果能够用套索工具代替方框工具选择区域。用套索工具的话，难点在于实现套索工具本身，这已经脱离了本程序的主要任务了。有了套索的选定区域后，就是获取边缘项烦了一点，大部分的步骤还是一样的。

UPD:

程序里发现了一个BUG！！！方程右端的边缘项应该是相加，我写成了直接赋值。上面的结果图片是错的，正确的图片如下：

后续工作

在下面的”实况记录“一节中，你会发现，我5个小时不到就完成了没有bug的程序，基本完成了整个毕业设计。得知了这一事实，你能得到哪些结论？

首先，你会说：”哇，你好强啊，这么快就把一个看上去那么复杂的项目写完了。“确实，这说得很正确。但这对于我来说是基本操作，甚至我还对自己有些不够满意，写代码的时候急躁了一点，不然可以更快写完。比起学东西，写论文和博客等文档的时间，实际编程的时间确实占比很小。

层次高一点的人，会评论道：”你这项目太水了吧！你干了什么？不久求个线性方程组？高中生都会解这种问题，叫一个大一新生学一学也能写出求解线性方程组的程序。“这说得不错。确实这个项目实现起来很简单。但是，你要注意到，这是一个科研性质的毕业设计。什么是科研？写代码是科研吗？不是。科研是要针对一个问题，提出某些解法，或者说改进现有的解法。重点是思路。我刚刚也说了，我的大多数时间其实花在了看文章、学习上面。这十分合理。

对科研有一点经验的人，又能找到我的漏洞了：”你刚刚说科研的定义，说得很好。但你做了什么呢？你复现了别人的方法。这是科研吗？这有任何创新吗？有任何改进吗？“能说出这种话的人，一定是高手。确实，本文的内容是不太够的。我的创新点在其他地方。这篇文章，只是我毕业设计的核心部分，还有一些内容没有展示出来。敬请期待以后的博客。

真正善于见微知著的人，还能得到其他的结论：”这就是本科教育吗？本科论文就这么水吗？“我只能说：“是的。”本科教育基本没有给学生们提供接触科研的机会，而学校为了方便大家毕业，不会在毕业设计上为难大家。没有了GPA，奖学金这些功利的东西，大家也没有认真做这个项目的理由了。本科毕业设计，就是一个从学校、导师到学生，基本没有人会认真对待的东西，从动机的角度来说，十分合理。

要是你再想得透彻一点，会发现这一切都是不是事。还是我最喜欢强调的，重点是你学到什么，你能得到什么。对于我来说，本科毕设强化了我的科研能力，让我找回了一点写代码的熟练度，锻炼了逻辑表达能力。虽然我花费的很多心思没有任何功利的作用，没有人会注意到我做的那些东西，但我知道我能得到什么就够了，不需要得到别人的认可。本科教育也是这样。本科究竟是为了什么呢？说来可笑，一部分去就业，一部分人继续科研。无论是工业界的知识，还是科研能力，本科一概不教。对于大学来说，能通过高考把优秀的人才招进来，本科设立的意义就已经达到了。既然如此，毕业设计怎么搞不都一样了吗。

实况记录

3.22 18:13

现在，游戏开始。先稍微讲一下我已经做好了哪些准备：我已经建立好了VS项目，导入了之前了OpenCV环境，并成功实现了普通的图像复制。我下好了别人实现论文的代码，以及代码中将用到的Eigen数学库。我还没有学过Eigen的用法，也没有把整个方法的实现流程完全想好。