双目立体匹配算法：Patch Match Stereo

一、简介

  我们知道，现有立体匹配算法一般被分类为局部算法、全局算法和半全局算法，其中局部算法和半全局算法是应用最为广泛的。在局部算法中，一个最简单的做法就是采用某种像素相似性度量，比如像素灰度差的绝对值AD，给定左图中的一个点p，在右图中的对应行上（假设输入是已经校正好的图像）搜索与其AD值最小的点q，这样得到的点q就是p在右图中对应的匹配点，p、q的水平坐标差称作视差。然而这种做法所得到的视差图中会包含大量的噪声，即错误的匹配对，原因可能是多方面的，如传感器噪声，左右相机的采集性能差异，图像中存在大面积无纹理、弱纹理或重复纹理，左右相机接收的光照差异（室外环境）等。
  一个更好的做法是不直接匹配单个像素，而是匹配像素点所在的区域，这个区域叫做支撑窗口(support window)，支撑窗口的尺寸可以是固定的，也可以是自适应变化的。如下式所示，由于支撑窗口内的像素提供了更多的信息，因此可以有效降低匹配歧义。

  然而使用支撑窗口的做法是存在问题的，实际上它隐性的遵从了一个假定，即窗口内的所有像素具有相同的视差。然而这个假定在很多情况下并不成立，比如：
  A. 窗口内的像素与中心像素位于不同的表面；
  B. 窗口所捕获的是一个倾斜表面或曲面，即非平行表面（这个平行指的应该是与相机成像平面平行）。
如下图(a)所示，Q点实际上位于一个亚像素的视差平面上，R位于倾斜平面上，S位于一个弧形表面上：

  许多研究人员主要关注的是A中所描述的情况，为了解决这种问题，一个比较有效的做法是对窗口内的像素采用自适应权重进行匹配，如下式所示，这样的方法可以很好的避免edge fattening(边缘平滑)问题：

权重计算所使用的方法可以是类似双边滤波的核函数：

也可以是基于测地线距离的核函数：

  那么如何解决B中的问题呢？我们知道空间中任意一点都可以认为是在一个唯一的平面上，如下式所示，该平面的参数为$($$\pi$${}_1,$$\pi$${}_2,$$\pi$${}_3,$$\pi$${}_4)$ $\in$ ${\mathbb{R}}^4$：

当我们知道相机的基线距离$B$和标定内参($u=fX/Z$, $v=fY/Z$, $d=fB/Z$)，那么上面的方程就可以变换为下式，其中${\pi }^{{}^{\prime }}=({\pi }_1^{{}^{\prime }},{\pi }_2^{{}^{\prime }},{\pi }_3^{{}^{\prime }})$ ∈ ${\mathbb{R}}^3$为视差空间的平面参数。

  基于上面的想法，2011年Michael Bleyer 提出了一个相当新颖的立体匹配算法叫做PMS，该算法的主要思想是对每一个像素计算一个独立的最优3D视差空间平面，如果该像素位于一个曲面上，那么该平面代表曲面在该像素点处的切面，在这个基础上，用于匹配的支撑窗口就是not fronto-parallel support window，简称slanted support window。这样一来问题的挑战就转移到了如何为每一个像素点在所有可能的视差空间平面中挑选出最优的视差平面。
  显然视差平面的数量是无限多的，因此通过遍历所有的视差平面来寻找最优平面是不可能的。Michael Bleyer 找到了一个巧妙的方法解决这个问题，那就是采用Patch Match的思想，Patch Match本身是一个高效求解近似最近邻场(Nearest neighbor filed)的方法，主要包括初始化、空间传播、随机搜索三个步骤。除了Patch Match中的空间传播外，作者还另外提出了视图传播和帧间传播（用于连续的视频帧）能够更好的帮助算法收敛。如下图所示，利用PMS不仅可以直接估计亚像素视差， 还能够精确的重建倾斜表面，甚至能够重建圆形的曲面。

二、算法

1.模型

  算法的目标是计算每一个像素所在的视差平面$f_p$，一旦得到了$f_p$，该像素的视差就可以按照下式来计算，其中$a_{fp}$, $b_{fp}$和$c_{fp}$是平面$f_p$的参数，$p_x$、$p_y$是像素$p$的$x$坐标和$y$坐标。

待求解的视差平面$f_p$是在所有可能的平面中使得聚合匹配代价$m(p,f)$最小的那一个平面：

其中F表示所有视差平面的集合。由于集合内元素的数目是无限的，所以不可能像传统的局部匹配算法那样，通过检查所有的平面来找到聚合代价最小的那一个平面（传统的局部算法是在一个有限的离散视差空间中进行搜索的）。像素$p$在视差平面$f$下的聚合代价是按照下式来计算的：

这里$W_p$表示中心位置在p处的矩形窗口。如果应用在视频中，这里的矩形窗口就由2D变成了3D，其中第三维代表前一帧或者后一帧。权重函数$w(p,q)$ 采用了自适应权重，目的是为了降低egde-fattenning问题：

式中$\gamma$是一个人为指定的参数，$||I_p-I_q||$表示像素$p$和像素$q$在RGB颜色空间中的$L_1$距离。当然，你也可以使用其他的权重函数，只要它能提升效果就行。函数$\rho (q,q^{\prime })$ 用来度量待匹配像素对之间的相似性，这里采用的方法是：

在上式中，$q^{\prime }$是通过$q$的视差计算出来的，所以$q^{\prime }$值在一个连续的实数区间内，它的颜色和梯度可以通过近邻像素进行线性插值得到。参数$\alpha$用于平衡颜色项和梯度项的比重。截断参数${\tau }_{col}$和${\tau }_{grad}$可以增强匹配代价在遮挡区域内的鲁棒性。

2.基于patch match的视差推理

  有了以上的铺垫，现在最重要的问题就是如何通过最小化匹配代价来为每一个像素找到一个最优的3D视差平面。这个问题可以通过patch match来解决：首先对所有像素的视差平面进行随机初始化或根据某些先验信息进行初始化，然后基于迭代传播(空间传播、视图传播、帧间传播、平面细化)不断地更新所有像素的平面参数。patch match的优点是只要至少有一个像素的初始化视差平面位于或接近最优平面，就可以找到剩余所有像素的最优平面。

2.1随机初始化

  patch match的第一步就是对左右视图的视差平面进行随机初始化，原则上是可以直接对平面$f$的三个参数$a_f$、$b_f$、$c_f$进行随机赋值，但是这样做无法保证在可允许的视差空间中进行均等的采样。因此需要对像素$(x_0,y_0)$设定一个视差搜索区间$[0,d_{max}]$(注意这个区间是连续的)，然后从该区间中随机挑选一个视差值$z_0$作为该像素的初始视差，这样就得到了视差空间中的一个点$P(x_0,y_0,z_0)$，有了点P，只要再得到一个单位法向量$\vec{n}=(n_x,n_y,n_z)$就能计算出平面参数：

如果想要使用fronto-parallel window，可以将法向量设定为$\vec{n}=(0,0,1)$，如果想关闭亚像素估计，可以将$z_0$限制在离散的整数范围内取值。

2.2迭代

  如下图所示，迭代过程中，每一个像素将经历4个阶段：空间传播、视图传播、帧间传播、平面细化，一般先处理左图，然后再处理右图。在偶数次迭代中从左上到右下依次遍历每一个像素，在奇数次迭代中，则顺序相反。

2.2.1空间传播(spatial propogation)

  空间传播的思想基础是在位置上相邻的像素一般具有相近的视差。设当前像素$p$的视差平面为$f_p$，其邻域像素$q$的视差平面为$f_q$，若$m(q,f_q)<m(p,f_p)$，则将像素$q$的视差平面$f_q$赋予像素$p$：$f_p:=f_q$，在偶数次迭代中，$q$为$p$的左侧和上侧的像素，在奇数次迭代中，$q$为$p$的右侧和下侧的像素。

2.2.3视图传播(view propogation)

  视图传播的思想基础是左(右)图中的像素点与右(左)图中的对应匹配点应当具有相同的视差平面。设左(右)图像素$p$的视差平面为$f_p$，其在右(左)图中的对应像素$p^{\prime }$的视差平面为$f_{p^{\prime }}$，若$m(p,f_{p^{\prime }})<m(p,f_p)$，则将像素$p^{\prime }$的视差平面$f_{p^{\prime }}$赋予像素$p$：$f_p:=f_{p^{\prime }}$。

2.2.3帧间传播(temporal propogation)

  帧间传播主要用于视频序列中，其思想基础是当前帧的像素$p$与其相邻帧的像素$q$的坐标相同，那么它们应当拥有相近的视差平面。很明显这个假设在帧率较高的视频中是成立的。设当前帧像素$p$的视差平面为$f_p$，其前一帧或后一帧的像素$q$的视差平面为$f_q$，若$m(p,f_p)<m(q,f_q)$，则将像素$p$的视差平面赋予像素$q$：$f_q:=f_p$。

2.2.4平面细化(plane refinement)

  平面细化的目的是通过更新像素$p=(x_0,y_0)$的视差平面$f_p$进一步降低匹配代价。由于平面$f_p$可以由一个点$(x_0,y_0,z_0)$和一个法向量$\vec{n}$表示，所以平面的更新可以通过点坐标和法向量的更新来得到。设$z_0$的最大允许变化量为${\Delta }_{z_0}^{max}$，$\vec{n}$的最大允许变化量为${\Delta }_n^{max}$，然后从$[-{\Delta }_{z_0}^{max},{\Delta }_{z_0}^{max}]$中随机选择一项数值$\Delta z_0$来计算$\Delta z^{\prime }:=z_0+\Delta z_0$ ,，从$[-{\Delta }_n^{max},{\Delta }_n^{max}]$随机选择三项数值来计算法向量$\overrightarrow{n^{\prime }}:=u(\vec{n}+\Delta n)$，式中$u()$用于将向量归一化为单位向量，这样就得到了一个新的平面$f_{p^{\prime }}$，若$m(p,f_{p^{\prime }})<m(p,f_p)$，则将视差平面$f_{p^{\prime }}$赋予像素$p$，即$f_p:=f_{p^{\prime }}$。
  这一步需要经过多次迭代才能完成，初始时设定${\Delta }_{z_0}^{max}:=d_{max}/2$ ，${\Delta }_n^{max}:=1$，然后每次迭代完后，将${\Delta }_{z_0}^{max}:={\Delta }_{z_0}^{max}/2$ ，${\Delta }_n^{max}:={\Delta }_n^{max}/2$，直到${\Delta }_{z_0}^{max}<0.1$迭代停止。假如像素当前的视差平面是完全错误的，初次迭代使用较大的搜索范围就有可能将其转移到正确平面附近，之后的迭代中不断的缩小搜索范围，使得视差平面在当前平面附近捕获视差细节，能够适用于弧形表面。

3.后处理

  后处理部分首先对左右视差图执行左右一致性检验，在小于设定阈值的情况下该视差点就会被设定为无效点，然后进行填充。填充的方法是搜索无效点左侧最近邻和右侧最近邻的有效像素点，设它们的视差平面分别为$f^l$和$f^r$，选择其中视差较小的那一个平面赋给当前无效点。之所以选择小的那一个是因为无效点多为遮挡区域的点，而遮挡区域往往处于背景之中。这样的缺点就是容易导致视差图中产生水平条纹，为了减轻水平条纹，可以对视差做加权中值滤波，滤波核的参数$\gamma$与匹配阶段所使用的参数相同。

4.为全局算法构建数据项

  全局算法对于遮挡区域和无纹理区域的处理要优于局部算法，所以本文提出的匹配代价还可以应用到全局算法当中，只不过此时将无法在使用连续的视差平面。下图显示了局部算法对于纹理匮乏的图像匹配失败，而全局算法却很好的处理了这种情况。

三、代码实现

由于该算法代码量较大，所以这里就不展示代码了，大家可以去百度或谷歌自行查找代码。

四、实验

  这里选择了MiddleBurry数据集中图像进行实验，效果如下图，可以看到该算法生成的视差图是稠密的，而且边缘保持的非常好，整体视差精度较高。

五、总结

1.PMS算法不同于传统的局部算法，它不是直接估计视差，而是估计视差平面，而且利用patch match思想在无限多的视差平面中来推理出最优视差平面。因此它的两个最大的优点是：(1)能够直接计算亚像素视差；(2)能够处理倾斜表面和曲面，这是其他局部算法所不具备的。
2.PMS算法的性能比较有限，虽然在middleburry上表现的确惊艳，但是在处理室外场景时鲁棒性不是特别好，特别是对于较高的图像噪声、大面积的弱纹理或重复纹理等(当然其他的算法也不能很好的处理)，其效果一般弱于SGM和ELAS等算法。
3.PMS算法的速度比较慢，因为它的所有处理流程都是顺序性的，不能并行处理，目前已经有一些算法改进了PMS中的传播方式使其能够并行处理，并应用于GPU，这里给出论文名称《Massively Parallel Multiview Stereopsis by Surface Normal Diffusion》。

参考资料：
1.PatchMatch Stereo - Stereo Matching with Slanted Support Windows，Michael Bleyer, Christoph Rhemann, Carsten Rother.
2.Stereo Matching—State-of-the-Art and Research Challeng，Michael Bleyer, Christian Breiteneder.
3.Massively Parallel Multiview Stereopsis by Surface Normal Diffusion, Silvano Galliani, Katrin Lasinger, Konrad Schindler.