立体匹配经典算法：PatchMatchStereo

点击上方“3D视觉工坊”，选择“星标”

干货第一时间送达

作者丨HawkWang

来源丨计算摄影学 

一. 前言

学习路线图：

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

我们先来看看，很多立体匹配算法没有解决好的一个问题，即亚像素匹配问题。我以前说过，立体匹配是一个让计算机做连连看游戏的过程——给定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次迭代，每次迭代对每个像素执行上述三步。对偶数次迭代，我们按扫描顺序先扫描图像左上角的点，一直扫描执行到图像右下角的点。对奇数次迭代，则相反，从图像右下角的点一直执行到图像左上角的点。

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

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

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

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

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

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

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

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

B. 接下来进行平面优化

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

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

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

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