78. 三维重建12-立体匹配9，经典算法PatchMatchStereo

一. 前言

我在77. 三维重建12-立体匹配8，经典算法ADCensus中画了一个学习路线图：

那么，今天咱们就进入经典视差优化算法的学习。

本文同步发表在我的微信公众号和知乎专栏“计算摄影学”，欢迎扫码关注，转载请注明作者和来源

我们先来看看，很多立体匹配算法没有解决好的一个问题，即亚像素匹配问题。我以前说过，立体匹配是一个让计算机做连连看游戏的过程——给定R图上一点，我们在T图的同一行上搜索和R图点匹配的同名点。这种搜索方式里，我们会将范围局限在一个有限的范围内，例如0到dmax，并且只是以整数视差来进行。

很多立体匹配算法会假设在某个像素的固定大小形状的支持窗范围内，视差都是一致的。当然在实践中，这种假设是不靠谱的——很多时候方形支持窗内有视差完全不同的像素。于是，就有了很多算法采用自适应的支持窗来解决这个问题，比如我们上一篇文章讲的ADCesus就是这样做的：

然而，自适应窗口依然是假设其内部的视差是一致的，而另外一种情况则会打破这种假设，即如果目标像素位于一个弯曲面或者斜面的情况，比如下图中的R点位于斜面上，S点位于弯曲面上。我们用普通的支持窗（图中红色线段）时，很显然窗内的像素都不在同一个视差上。

另外，由于搜索是在整数视差上进行，因此如果目标的视差刚好位于两个相邻整数之间，搜索出的结果也不够准确——这是自适应窗口无法解决的问题，比如上图中的Q点，它的视差我们看到大概是在1.6左右，不管是普通的正对相机的支持窗，还是自适应支持窗，都无法得到Q点的小数视差值。我在文章73. 三维重建8-立体匹配4，利用视差后处理完善结果中描述的亚像素插值能够部分解决此问题，但它建立的抛物线插值模型只是对真实目标表面的近似。

那么，有没有更好的方式呢？有的，那即是对每个像素所在的平面进行建模，这样我们的匹配过程就变成了求该像素的最佳平面参数的过程。

我们来看看，采用传统的“正面朝向相机”的支持窗，及采用斜面支持窗时，效果的区别，很明显在这种场景下斜面支持窗获得了更好的效果。

那么，该如何在匹配过程中对计算每个像素所在的平面呢？这里面的关键问题是：

1. 传统匹配搜索得到的是整数型的视差值，我们最多搜dmax次即可完整匹配
2. 而一个像素所在的平面可能有无穷多个, 平面的参数是3个浮点数，此时我们该如何确定这个平面呢？

这就引出了我们今天介绍的方法，这是由Michael Bleyer等在2011年发表的方法，看名字我们知道，它采用了基于PatchMatch的思想来解决斜面支持窗时的立体匹配问题。我们在下一节来仔细学习它。

二. PatchMatch Stereo的核心思想

2.1 问题的描述

首先让我们看清楚空间投影的模型，话不多说，如下图所示：

上图中，任意一点P所在的平面参数fp是未知的。我们需要做的就是在所有可能的平面中搜索，直到确定平面参数，一旦确定了平面参数，就可以利用图中视差和平面参数关系，求出所关心的点p的视差值。

这个过程，用论文上的公式，可以写作下式。其中F 代表所有可能的平面

其中m表示一个代价函数，表示如下：

对于某个像素来说，代价是在一个方形窗口内进行的，下图中q是我们关注的像素p点的邻域像素：

其中，自适应权重w(p,q)表示为：

而q点对应的代价则同时考虑了颜色和梯度，并用一个参数α来平衡

看起来，该算法走的是我在文章70. 三维重建5-立体匹配1，立体匹配算法总体理解中提到的全局法的路子，来求取视差图，并且不用做代价聚合。

有趣的是，虽然我们认为这个过程是全局优化的路子，但上面的操作却是对一个一个像素单独进行的，所以作者将这个方法分类为了局部法。这里由于潜在的平面数量几乎是无限多的，所以我们很难直接去优化论文中的(2)式，这也正是本文的贡献点——作者巧妙的采用了PatchMatch的思想，来优化该代价。我们看看它是如何做的。

2.2 用PatchMatch思想来快速获取平面参数，再获取视差

2.2.1 随机初始化

首先，我们对所有像素随机初始化其所在的平面。为了让初始化平面反映出有效的视差值，这里面有一点点讲究，并不是随便分配3个平面参数即可。具体的过程如下图所示：

让我们用MiddleBurry的MotorCycle数据做一下实验，可以看到这一步结束后的左右两幅视差图只是一些随机的噪点，这符合我们的预期。

然而，我们前面指出当前这个算法用到了PatchMatch的思想——到底是什么思想呢？作者观察到了一个重要的情况：虽然每个像素的平面参数都是随机分配的，但是由于大数定律，在成千上万像素中，总会有少数的点它们的平面参数特别接近真实情况。

因此，接下来我们要做的，就是把这些特殊点的平面参数传播给它的邻域像素。

2.2.2 信息迭代传播

接下来我们进入到信息的迭代传播过程，这里还是需要先了解背后的直觉，如下图所示。

正是基于这种直觉，1次传播迭代分为三步：

1. 空间传播
2. 平面优化
3. 视角传播

而作者做了3次迭代，每次迭代对每个像素执行上述三步。对偶数次迭代，我们按扫描顺序先扫描图像左上角的点，一直扫描执行到图像右下角的点。对奇数次迭代，则相反，从图像右下角的点一直执行到图像左上角的点。

我在下面这个视频中可视化了这个过程，戴上耳机欣赏哦：

http://mpvideo.qpic.cn/0bc3xqaaiaaag4aocfizi5rvbpgdas6aabaa.f10002.mp4?dis_k=443fdfbd13bc72eeb1a4c8128ff5c973&dis_t=1669795502&vid=wxv_2567740766828085251&format_id=10002&support_redirect=0&mmversion=false

现在我们来看看具体过程：

A. 先来看第一步，空间传播

如下图所示，我们比较p点和其邻域点的平面参数，看看哪一个代入到p点的坐标时上面公式(2)得到的代价值更小。如果邻域点的平面参数代入时，得到的代价更小，那么就用邻域点的平面参数替代p点的平面参数。如果当前的迭代是偶数次时，我们比较p点和其上、左邻域点，而如果当前迭代是奇数次时，我们比较p点和其右、下邻域点。

让我们只做1次空间传播，看看效果吧：

非常奇妙，简单的做了空间传播，我们就能大概从视差图中看出形状了，别忘了刚才它们还是一堆随机数哦！

我们再来看看做两次和做三次空间传播的效果，可以看到一次比一次好：

看看细节对比就更明显了：

B. 接下来进行平面优化

记住我们并不是做完3次空间传播，再继续进行平面优化。而是对每个像素，每做一次空间传播，紧接着就做平面优化和视角传播。

注意，论文里面是先做空间传播，然后做视角传播，再做平面优化。但我这边找到的两份代码都是先做空间传播，再做平面优化，然后做视角传播。所以我也就按照代码的顺序来讲了。

前面的空间传播，不可能直接算出像素点的最优平面参数。我们只能说经过空间传播，像素点的平面参数更加接近正确值了。那么，对任意一个像素，我们还可以在一定范围内随机改变它所在的平面参数，看看代价是否降低。如果是，那么就接受新的平面参数——这就是所谓的平面优化。具体过程如下图所示：

平面优化的过程的一个关键技巧是随机干扰的范围。对于点的视差值，一开始随机干扰的范围为 [−dmax/2,dmax/2] ，而对于法向量的干扰则为[-1, 1]。接下来就是循环执行上图中的3步，并在每一个循环过程中减半随机干扰的范围。直到视差干扰范围值 |Δd|<0.1

这种干扰范围指数下降的方式能够使我们一开始进行较大范围的调整，而后当平面参数接近正确时又快速降低调整范围，不至于跳过最优平面。

让我们看看，经过一次空间传播和平面优化后，效果如何吧。从下图你可以明显看到，当做了平面优化后，整个视差图好了很多，特别是地面已经变平了：

细节对比：

如果经过3次迭代呢？从下图可看出3次迭代后视差图更加准确，凸显出了更多的细节，比如后面的椅背：

C. 迭代的最后一步是视角传播

正如前面所说，两个视角的匹配点可能位于同一个平面，因此我们可以在视角之间传播平面信息。具体来说，假如我们在左视图有一个点p，我们在右视图找到所有和点p匹配的点p'(可能多于1个)。如果某一个p'投影到左视图后的的平面是fp'，并且这个平面信息相比点p有更低的代价值（公式(3)），那么我们就用fp'来代替原有的fp。

让我们来比较下加入这个步骤的效果：

放大图，我们看到加入视角传播后摩托车把手更加有细节了：

2.2.3 视差后处理

正如一开始的流程图所示，视差优化后就进入到了视差后处理的部分

首先，还是经典的左右一致性检测，确定错误匹配点。

这个操作的结果如下图所示：

接下来，对于错误匹配的像素p，我们确定其左右两侧最近的正确匹配的像素，设它们为a和b，那么这两个像素所在的平面我们称为fa和fb。现在利用fa和fb，以及p的x/y坐标，我们可以计算出新的视差值。接下来我们取这两者之间更小的那个，作为p点的新视差值。由于我这边选择的摩托车的图集，基距大，最近视差大，因此有很多遮挡区域。可以明显看出在右视图的视差图上出现了很多错误，即便是填充也没有挽救回来。而左视图虽然看起来还不错，但也出现了横向的拉丝效应。

所以作者最终又采用了一个加权中值滤波来消除水平的条纹。下面是最终的结果。不过我对这个滤波器的效果存疑，虽然我们看到摩托车的细节变好了，但本来平滑的地面视差似乎被这个滤波器破坏了，变得略显凹凸不平，有点遗憾。

三. 总结和参考

PatchMatch与ADCensus等算法不一样的地方的地方在于，视差的计算不是直接进行的，而是通过平面参数计算而来的——立体匹配的过程不是在一维的水平极线上进行搜索得到视差值，而是在3D空间中搜索最佳平面，再通过平面参数反算出视差。这意味着在匹配过程中就可以得到亚像素精度的视差值和正确的平面。比如下图就能明显看出这种算法的优势。

在同一篇论文中，作者还描述了一种基于上面的思想进行全局代价优化的方式，类似于我之前将的思想：先建立全局代价函数，然后对它进行优化得到视差图：

其中数据项的计算就会用到上面PatchMatch的思想：

这种方式也取得了非常好的效果，比如下面这个效果，想必你很眼熟这幅图了：

然而，PatchMatch立体匹配算法在未经优化时，运行速度非常慢，因为很多操作都是对一个个像素单独进行的。所以至少要将其并行化才能快速跑起来。

不管如何，本文踏出了将PatchMatch思想用于立体匹配的第一步。后面又有很多人在本文的基础上研究了各种各样用PatchMatch思想进行立体匹配的方案，我就在此不表了。

现在看看我们的进度，如下图所示。你可以看到视差优化这部分我们只完成了一半，这是因为我在下一篇文章中还会介绍一种更震撼的优化算法，是什么呢？让我买个关子吧，你可以留言猜一猜：

本文同步发表在我的微信公众号和知乎专栏“计算摄影学”，欢迎扫码关注，转载请注明作者和来源

四. 参考资料

• Michael Bleyer , Christoph Rhemann , Carsten Rother, PatchMatch Stereo – Stereo Matching with Slanted Support Windows, 2011
• Connelly Barnes, PatchMatch: a randomized correspondence algorithm for structural image editing, 2009
• 李迎松解读PatchMatchStereo: https://ethanli.blog.csdn.net/article/details/107192399
• 一个比较好的实现: https://github.com/ibergonzani/patch-match-stereo