“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$是输出的结果图像。我们还有一个输入，就是进行图像融合的区域$\Omega$。既然我们把图像看成了函数，那么区域$\Omega$就是整个定义域$S$的一个子集。
$v$是$g$的梯度，$\partial\Omega$是区域$\Omega$的梯度，这两个量都是计算得到的。
我们还要引入一个符号——梯度算子$\triangledown$。$v$是$g$的梯度，可以写成$v=\triangledown g$。
重申一下结果图像求解的思想：在融合区域边缘的颜色和目标图像一样的前提下，让源图像的梯度尽可能不变（最小化差值）。用数学语言就是：
如果理解了求解的思想，那么看懂，或者说推导出这个公式是十分简单的。
这个公式是理想情况：图像在一个连续的二维数集上定义。但实际上，我们图像是离散是，只有在整数的位置处有值。我们要把问题用离散的形式表达出来。这里又要提出一些标记，对于$S$中的每个像素$p$，令 $N_p$ 为其4邻域中仍在 $S$ 的像素集合。也就是说，这个符号表示的某个像素上下左右这四个像素的集合，在这些像素没有跑出整幅图像的前提下。令 $(p,q)$ 为一个像素对，其中 $q \in N_p$ 。
再令 $v_{pq}$ 表示 $pq$ 两点间的像素值之差，即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方程组的代码：
1
2
3
4
5
6
7
| 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替换掉即可。
1
2
3
4
5
6
7
| 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的用法，也没有把整个方法的实现流程完全想好。
18:26
理解了此方法中如何构建最后的求解矩阵。现在学习如何用Eigen求解稀疏矩阵。
18:36
知道如何用Eigen创建矩阵，求解矩阵了。可以考虑写代码了。先想办法把Eigen导入项目，并且测试一下Eigen解方程代码。
18:39
很烦！Eigen还要编译，计划中断。我准备一边编译一边盲写代码（写的代码无法进行测试）。
18:57
一边复习OpenCV，一边写完了拉普拉斯滤波代码。同时Eigen似乎编译完了。
19:09
Eigen编译到一半，出了问题，开始向搜索引擎求教。最烦人的配环境时间开始了。
19:21
气死了，这个库可以不编译直接使用。重新开始Eigen的使用测试。
19:36
写到一半，OpenCV配置又出错了。原来是之前写的OpenCV配置文件少包含了一个库文件。批评一下过去的自己，太坑了吧。
19:47
OK，Eigen测试完了，游戏结束。现在所有不确定的东西都搞定了，前置知识学习完毕。接下来我可以展示王者级别的写代码能力了。
20:26
潇潇洒洒写了几十行，突然发现我不会对OpenCV和Eigen两个库的矩阵进行互相转化，又开始查东西了。时间还很充裕，问题不大。
20:30
原来OpenCV提供了转换函数。好，工作继续。
20:37
虽然我没有花半秒钟在程序的软件工程设计上，但我潜意识里已经把代码的结构、实现顺序给想好了。凭借着强大的功底，我本能地使用了自底向上的思想，先把每个细节的实现函数写好了。现在我准备把整个程序拼起来。
20:46
好，程序拼完了，已经开始编译了。说实话，第一遍能不能过（没有BUG）？肯定过不了。无论多强的高手，都会写BUG，都有大意的时候。BUG本身并不可怕，可怕的是不会调试BUG。编程能力展示结束，现在展示我发现BUG的洞察力和解决BUG的决断力。
20:54
代码还是倒在了OpenCV和Eigen两个库的互相兼容上。三个小时快到了，你会不会觉得我很慌？你会不会觉得我已经写不完了？不好意思，现在我才刚刚登上展示实力的舞台。只有在逆境中绝处逢生，才能体现出一个人的过人之处。
20:58
突然发现代码漏写了！图像有RGB三个颜色通道，该方法要对每个颜色通道单独处理，我忘了拆开颜色通道了！
21:03
OpenCV和Eigen相互转换的库函数不能用，我就自己动手！总算，稍加调试之下，代码可以运行了！只不过代码的运行效果明显不对……这下调试起来就麻烦了。很好，最尴尬的情况出现了。我不会放弃的！
21:17
时间已经过了。但我有大约10分27秒花费在了写博文上，这些时间应该不算。继续计时。
21:45
我发现我对论文方法的理解出现了问题。
22:25
更改对论文理解上的逻辑BUG后，又发现8位颜色通道溢出了。
22:34
历经千辛万苦，程序总算是成功实现了。总耗时约260分钟。多花费了约50%的时间，你肯定说我很菜，说我写代码写得烂。但是这种说法并不对。仔细看一看这份文字实况，你能看出什么？没错，我一直在强调我纯粹的编程实力，并没有强调我有快速把论文转换成代码的能力，也没有强调我有快速学会调用那些坑爹的公共库的能力。我写代码写了多久？19:47~20:26,39分钟；20:30~20:46，16分钟，一共55分钟。第一版本的代码，我花55分钟就写完了。后面又调试了一下，20：54~22：34，100分钟，减去重看论文的21：17~21：45这28分钟，又花了72分钟调试，最后也就改了5%不到的代码。我说这么多数据是为了什么？是为了说明，恰恰是这份实况记录，强调了我编程能力之强，强调了时间其实都耗费在一些编程之外的时间上。虽然没有在规定时间内写完，但我感觉虽败犹荣。
当然，你要说我每个方面的表现都很好吗？也不是这样的。我自己也承认，理解论文的时候急躁了一些，学库函数的时候学的不仔细了一些，导致我这个第一版的程序，出现了一点小小的偏差。不过真要说起来，我今天确实有点急于求成。因为立下了3小时完成的flag，所以我做事火急火燎，忽略了工作质量，结果败给了自己。如果一开始没有这个flag，我能不能3小时内完成呢？你可能会认为我会说：“会”。好，实际情况是，我不立这个flag，我今天会玩一个晚上，不会来做这个项目，哈哈。
参考文献
[1]. Pérez P, Gangnet M, Blake A. Poisson image editing[M]//ACM SIGGRAPH 2003 Papers. 2003: 313-318.
下载很快就会结束。规模这么庞大的编程软件，几分钟就下载完了？其实我们现在只下载了一个下载器，待会才会正式开始下载VS的本体。
刚看到这个界面，第一次用VS的人肯定会有点懵：哎呀！怎么这么多选项啊！该点哪个呀？
所以，我们应该根据自己的需要选择要开发的语言。只进行C++的开发的话，只点击”使用C++的桌面开发“即可。上图中，我多点了几个选项。
