79. 三维重建14-立体匹配10，经典视差优化算法Fast ...
一. 前言到目前为止，我的双目立体匹配系列文章已经完成了接近经典视差优化算法的介绍，让我们观察一下现在的进度：
今天的文章，我们先从回顾立体匹配的评价指标开始，我在文章74. 三维重建9-立体匹配5，解析MiddleBurry立体匹配数据集中提到过MiddleBurry立体匹配排行榜的评价指标：
MiddleBurry立体匹配评价指标
这些评价指标大概包括：
[*]在全图上计算视差图和理想视差图之间的均方根误差，及错误像素占比
[*]在无纹理区域，有纹理区域，遮挡区域，非遮挡区域，深度不连续区域共5个区域计算和理想视差图之间的均方根误差，及错误像素占比
[*]在不同视角下进行反向变换，计算变换后的投影误差，即所谓预测误差
在我介绍其他的数据集里，基本上采用了类似的评价指标。我不知道你想过没有，这些评价指标有什么特点？
事实上，这些指标都在强调视差图在视差方向上的准确性，却没有怎么强调视差图与原图之间的贴合性——或者说，它们都在强调Z方向的准确性，却没有重视其在X/Y方向上的准确性。
对于很多三维重建的应用，比如测距、尺寸测量、甚至自动驾驶中的障碍物检测识别，似乎这也没有构成很大的问题。然而有一类应用，对重建后的深度或视差在X/Y方向上的准确性确有很高很高的要求，这就是我在文章手机中的计算摄影1-人像模式(双摄虚化)以及手机中的计算摄影6-旷视技术开放日上展示的手机电影中提到的虚化渲染类应用——事实上我认为绝大多数像素级渲染类的应用，都对场景模型和图像纹理之间的贴合度有很高的要求，而这也意味着我们要求深度图中的内容和原始的彩色图之间也需要紧密贴合。
虚化渲染对边缘贴合有很高要求
如果不能做到这点，我们渲染的结果就会出现非常突兀的瑕疵，比如下面这样的：
边缘不贴合出错
或是这样的：
边缘不贴合出错
很显然，一般的基于MiddleBurry类数据集评价指标开发的立体匹配算法是无法满足这种需求的。我们所需要的算法，必须要从原理上能够保证这种贴边，今天我将要展示的算法就做到了这一点。按照作者的说法，根据这种算法的原理你甚至无法生成不贴边的视差图！
本文同步发表在我的微信公众号和知乎专栏“计算摄影学”，欢迎关注，转载请注明作者和来源，觉得有用的话顺手点个赞哦 看了后还有不懂的，可以点下面卡片发起咨询，或联系我加入我的知识星球进行深入讨论
在正式介绍今天的算法之前，让我们先来看看效果。首先我们来看看一对输入的经过立体校正的图像，及其局部（下方）：
输入原始图
然后我们看看普通的SGM算法得到的视差图，以及用这种视差图进行浅景深虚化渲染的结果。你可以很清楚的看到细节部分出现了严重的不贴边现象：
SGM算法结果和渲染图
下面则是我今天将要介绍的算法的结果。真是不怕不识货，就怕货比货啊，从渲染的结果看上去明显好了很多！
本文算法结果和渲染图
下面是另外一个场景，注意看原图中雕塑的细节部分。由于上采样看细节的缘故，小图显得不那么清晰，可以忽略这一点。
输入原始图
我们先看看一个叫做SPS-StFI的算法结果：
SPS-StFI算法渲染结果
表面上看起来视差图还OK，但我们看看渲染的细节，很明显在雕塑的边缘处出现了奇怪的瑕疵，你可以和上面的原图对比一下，这显然是由于视差图中的边缘错误导致的
SPS-StFI算法渲染结果细节
现在来看看本文将要介绍的算法的结果：
本文算法结果
看细节，很明显主体和虚化背景的边界非常自然，没有出现上面SPS-StFI算法的边缘瑕疵：
本文算法渲染细节
小结一下我们现在得到的信息：
传统的基于类似MiddleBurry这样的数据集的评价指标来设计的双目立体匹配算法，大多数侧重于视差值或者深度值的准确性，却对视差图贴合原图物体边缘的程度不够重视，不适合图形图像渲染类应用。而有一类算法，比如说我今天将要介绍这种算法，特别强调视差图与原图目标边缘的贴合程度，从而比较适合产生渲染类应用所需要的视差图(深度图)。
那么是什么样的算法能够从原理上保证这一点呢？它是如何做到的呢？让我请出今天的主角吧：
对，就是Jon Barron，有的朋友可能已经知道他了。 如果我再提及他参与提出的NeRF(神经辐射场)，就肯定有更多人会哇哦了。
我们今天要谈的是Jon Barron在2015年CVPR上演讲展示的一篇文章Fast Bilateral-Space Stereo for Synthetic Defocus
它首先从原理上保证了视差图的贴边，而且速度还大大超越其他算法，在下图中我们看到它比之前最快的算法还快十几倍！
本文算法速度非常快
那么，它是如何做到的呢？要回答这个问题，我们需要先回顾一下立体匹配的全局优化法的思想。
二. 全局优化算法的思想
在我的文章72. 三维重建7-立体匹配3，立体匹配算法中的视差优化中，我曾经讲过，立体匹配的全局代价优化法的思想是希望寻找到每个像素的最优视差结果，使得全局的、整体的匹配代价最小，这一步被称为视差优化(Disparity Optimization)。于是这个过程就变成了一个最优化某个能量函数的过程，该函数通常写成如下的形式：
等号右边第1项是数据项，它衡量计算出的视差与实际输入图像关系的差异。一般来说，可以用下面的式子来表示，其中C表示代价函数。
这一项用于约束全局代价最小化。 但是代价函数中通常含有噪声和错误，直接最小化求得的结果也会有很多问题，所以还需要第2项平滑项。 这一项一般用于给出某些额外的约束条件，比如通常假设整个图像的视差是平滑变化的。这样视差的大变化只会在场景内视差边缘处产生，一般也和图像内物体边缘高度相关。
很显然，这个全局代价的优化是非常困难和复杂的问题。 比如，一幅1024x1024的图像，如果其有效视差范围为128，那么上面的数据项就是一个有着1024x1024x128=134217728个元素的巨大的立方体。这还没有算上后面的平滑项的复杂度。
其实，虽然如此复杂，但似乎也不是不可以接受，无非就是多花些时间吧。然而，在有些应用场景下，我们需要高速进行深度重建，甚至实时深度重建，比如我提到的手机上的各种渲染应用就是如此。你想象一下，你拿起手机要给女朋友/男朋友拍一张人像照，肯定是希望所见即所得，并且按下快门立马得到结果，而不是等上几分钟都看不到效果吧。
所以，总结下Jon Barron要想解决的问题：用极高的速度解决立体匹配的全局代价优化问题
在计算机科学中，甚至在所有门类的科学中，有一种思想是问题的转换。我们不是说上面的全局代价函数的规模太大，导致优化复杂吗？那么，就先把这种大规模的全局代价函数转换成小规模的全局代价函数吧，这样求解起来不就简单了吗？
我们接着来看看Jon是如何一步步做到这种转换的。
三. 视差图的平滑性与双边滤波
现在重新想想我们的视差图需要满足什么样的特性：视差图总体上是平滑变化的，视差的突然变化只会在场景内物体的边缘处产生。你想到什么了吗？这和什么算法有相通之处？
是的，这和我之前讲过的双边滤波有非常大的相似之处，让我们回忆下4. 数码相机内的图像处理-更多图像滤波中的内容：
高斯滤波只使用了空间距离来衡量像素的权重，而双边滤波则在空间距离的基础上，加入了像素亮度距离。而由于边缘两边的亮度值差异很大，因此越过边缘的像素的滤波权重很小，这就是双边滤波能够保持边缘的秘诀。
双边滤波的图示如下：
双边滤波有这种优异的性质，但其计算量却非常大，这也促使了人们想了各种各样的办法来快速计算双边滤波。其中有一种称为双边网格的方法特别有意思，这是2007年Chen Jiawen等人提出的。
我们看下面的图，这就是Bilateral Grid(双边网格)的思想。一幅2D的灰度图像中的每个像素，按照其x坐标、y坐标、灰度值，被分组放到了一个3D的立方体中，每一小组就是一个格子。这样，原本在2D空间中相邻的像素，到了双边网格空间中，就不相邻了，至少被放到了不同的格子里。
利用这种思想，我们实际上是对原始图像进行了降维，并且使得灰度值差异较大的像素不会相互污染：
利用这种思想，可以实现快速的双边滤波，我们以一个一维图像来举例。首先，我们用刚才的思想把图像投射到
双边网格中。接下来，我们在双边网格空间中进行高斯滤波。最后再把滤波后的双边网格重新投射回原始的图像，如下图所示：
由于双边网格的规模大大小于原始图像，比如一幅八百万像素的图像投射到双边网格中可能只有70*70*10=49000个格子。而高斯滤波又比双边滤波的计算量低很多，并且很容易通过GPU来并行计算。因此通过这样的方式，就可以实现非常快速的双边滤波了。
比如Chen JiaWen等在演讲中提到，一幅两百万像素的图像，如果用原始的CPU上实现的双边滤波需要10分钟，而用他们的双边网格思想在2006年的G80 GPU上实现，只需要9毫秒！ 这就是将难解问题转换成易解问题的威力，让我们再看一遍这个图，即便你研究的领域并非立体匹配，并非图像滤波，今天这篇文章的其他内容你都可以忘记，但我希望你记住这幅图：
在相关文献中，把从原始图像到双边空间的过程称为Splat，在双边空间中进行滤波称为Blur，把滤波后的结果转换到原始图像像素空间中的过程称为Slice。 Splat-Blur-Slice，就是这类方法的三个重要步骤。由于高斯滤波是可以分解的，在双边网格中的高斯滤波可以变换为每个维度上（对于1维图像来说是space和intensity两个维度）分别进行的一维滤波相加的结果（这是一种近似）
那么，这对我们今天将的立体匹配全局代价函数优化有什么启示呢？
Jon Barron在演讲材料中写道：
什么意思？ 前人的工作是将图像信号投射(Splat)到双边空间，在双边空间中滤波(Blur)，再反投射(Slice)滤波后的图像信号到原始像素空间。而Jon Barron的工作则是将所要求解的最小化问题投射(Splat)到双边空间，在双边空间中求解(Solve)该问题，再将求解结果重新投射(Slice)到原始空间。
由于双边空间中的问题规模大大小于原始问题，因此问题求解就变得很快速很容易了。我们在下一节，进一步看看Jon是如何做到的。
四. 利用双边空间优化全局代价函数
我们看看下面Jon构建的代价函数，加号前面一项是平滑项，后面一项是数据项，我们就是要最小化这个式子，
全局匹配代价函数
4.1 平滑项
平滑项
平滑项中的di和dj是一幅视差图中不同位置的像素的视差值。这个时候，我们要把视差图看作展开的一维向量。一个宽高分别为W和H的视差图，展开后就是WHx1的一维向量，于是i和j就是不同像素的编号，或者坐标值。
最神秘的是这里的 \hat{A} 矩阵，为了理解它，我们首先得重新回到双边滤波器的矩阵表达，如果x是一个一维表示的尺寸为WHx1的图像，我们可以用下面的公式来表示双边滤波，等号坐标是滤波后的图像y，尺寸还是WHx1：
双边滤波矩阵表示法
这里的矩阵A是一个对称的矩阵，尺寸是WHxWH。对于彩色图像的双边滤波，A矩阵的每个元素如下：
双边滤波权重
它实际上是表示第j个像素和第i个像素间的在空间和颜色上的相似度。如果两个像素越靠近，并且颜色越相似，其对应的A元素就越大，当然如果是单色图像那么就只是计算空间和亮度上的相似度。
我们看看下面这个3x3的小图像（我们假设其中存储的就是对应像素的视差）
3x3的小图像
按照上面的公式很容易计算得到其双边滤波矩阵如下
权重矩阵A示意
利用这个矩阵，可以把任何我们认为相邻的像素的视差值加权平均到一起，得到最后的滤波结果。而我们前面数据项中的 \hat{A} 矩阵则是双边滤波矩阵A的某种特殊归一化形式，其最大特点是每一行和每一列的和都为1。
4.2 平滑项转换到双边网格空间
我们可以用下面的公式表示双边网格算法，其中S矩阵代表我最开始提到的Splat操作，将像素值投射到双边空间。 \bar{B} 矩阵代表双边空间中进行的模糊操作，而 S^T 矩阵则代表上面提到的Slice操作，将双边空间中的滤波结果重新投射到原始像素空间。在双边网格的构造过程中，每个网格的顶点采用了最近邻方式采样，所以这种方式是对原始双边滤波的一种近似：
双边网格是矩阵分解过程
进行变量替换后，双边滤波就变成了下面这样：
在双边空间中进行滤波
而作者最精妙的一点是，把匹配代价中的平滑项按照类似的方法，投射到了双边空间。为了做到这一点，首先是把原始的平滑项写成矩阵表达形式, 再通过变量替换转换到双边空间
下面对双边网格中的各个变量作出解释，请注意在双边空间中的这些矩阵和向量的规模都是大大低于像素空间中的对应实体的，这也为之和的快速求解打下了基础。
你可能对为什么Jon Barron能够把原始像素空间中的平滑项做上面所说的转换感兴趣。这里面有非常复杂的数学推导过程，都记录在了其论文的附加材料里面。不过很多人直接阅读此份附加材料也是看不懂的（也许你就在这些人里面 ）。如果你对此感兴趣的话，可以私聊我。如果咨询的人足够多，我还打算在我的知识星球中做更详细的讲解，包括针对代价做细致的分析。
4.2 数据项
现在我们来看看数据项。正如前面所说，数据项表示的是像素点的匹配代价
数据项
我们把它也要投射到双边空间中——每个双边网格所对应的代价值是多个原始像素代价值的加权平均：
编号为j的网格对应的匹配代价
比如，这里j代表第j个双边网格中的格子，其中包含了多个以下标i表示的原始像素，w代表像素的权重。这样就把像素空间中的代价fi转换成了双边空间中的代价 g_j
现在我们汇总一下，原始的代价函数为：
像素空间全局代价函数
在双边空间中的代价函数为：
双边空间全局代价函数
为了让计算尽可能简单，而且可优化，这就需要g是一个连续函数，且可导，这也意味着f是一个连续可导的函数，现在我们必须写出g函数。
我在之前的文章中介绍过了很多种不同的匹配代价函数，而作者选择的是下面这个满足所谓Birchfield-Tomasi 标准的函数：
我先在下面给出作者论文中的图示，然后用人话给你解释一下：
单像素匹配代价示意
什么意思？如果两个像素通过某种标准判断是匹配的，那么其匹配代价为0. 否则，其匹配代价正比于代价值偏离匹配范围的大小。 对任何一个像素i，在水平极线上都有可能有多个匹配的像素。其第1个匹配的像素对应的视差我们称为 l_i ， 最后一个匹配的像素对应的视差值我们称为 u_i 。所以只要di在[ l_i , u_i ]之间，我们认为是匹配的，代价为0，否则匹配代价将位于斜率为1或-1的直线上。在实际操作中，我们需要记录每个像素的匹配视差的上下值，其过程我这里略过，感兴趣的私聊我吧。
相比我之前介绍的很多匹配代价函数，这是一种很弱的匹配代价函数。作者选择它并非因为它有什么神奇的地方，仅仅是因为它非常高效，特别是当我们需要最小化之前提到的全局代价函数时这一点尤为重要。
4.3 在双边空间中进行优化
希望你还没有晕掉我们已经得到了双边空间中的数据项和平滑项了，为了进行优化我们需要将数据项表示为连续可导的函数。所以这里采用了线性插值的方法，将在双边空间中离散的g函数变为连续的，原论文公式如下：
双边空间中的全局代价函数
最优化问题有很多求解方法，作者采用的是一种叫做L-BFGS的方法，这是一种在有限的内存空间中进行的优化算法，是Broyden,Fletcher,Goldfarb,Shanno这几位数学家发明的，他们名字的首字母组合就是BFGS
BFGS
我们不需要了解算法的细节，只需要记住许多最优化问题求解的一般流程如下：
优化算法过程
所以我们再把上述代价函数的梯度表达出来，原论文的梯度公式如下
有了上面的(15)和(16)两个式子，我们就可以进行优化求解了。首先用下面的(17)式，我们求解双边空间中的视差值v*，然后把v*重新投射到像素空间中，得到最终的视差值
Solve/Splice
五. 实验和讨论
我已经在文章开头展示了作者演示的效果，我总体上的感觉是这种算法在当目标和背景的颜色不一致时，视差边缘很贴合目标的边缘，我将作者给出的视差图和原始彩色图叠到一起给你看看，很明显能看出来。
不过，由于作者采用的数据项中匹配代价函数过于简单，这也就使得数据项误差较大，导致算法会出现过平滑的现象。表现在结果上，就是会在结果视差图中带入很多原图的纹理。我们来看几张网上的高人实现算法后跑出的结果（参考资料4）：
六. 总结
今天我们学习的这种叫做Fast Bilateral-Space Stereo的算法与之前介绍的其他立体匹配算法有很大的不同。它并不追求深度重建结果在Z方向的绝对精确，反而更强调其在X/Y方向的贴边。因此其结果更适合在图形图像渲染领域中使用，而非在对尺寸、距离等有严格要求的三维重建领域。该算法速度很快，因此一经推出，就在Google的Camera应用中得到了应用，该功能被称为Lens Blur，专用于浅景深虚化，也就是现在大家熟悉的手机上的人像虚化模式。
我们从作者的思路中学到的最关键的洞见是：一个复杂大规模的问题可以通过空间转换，变换为小规模的问题。求解完小规模的问题后，再重新映射回原始的问题空间，得到想要的解。
实际上，关于这篇文章还有很多可以讲的，但限于篇幅原因我不得不略过了。但有两大内容，我非常建议读者进一步学习：
[*]双边网格的原理和具体实现方法，这是本文的基础，也是我接下来还会描述的另外一个算法的基础。而且快速双边滤波在很多领域都有用，因此我觉得掌握双边网格很有用。
[*]本文所略过的Jon Barron的匹配代价平滑项、数据项转换到双边空间的具体推导过程。这是本算法背后的精髓所在，理解了它对理解我下一篇文章也很有帮助。
上面两点，我都会给出公开的参考资料。但坦率说，这些论文和公式你不就是不愿意自己读嘛 所以，如果你真的感兴趣，可以私聊我或者加入我创建的知识星球咨询，咱们通过付费学习的方式把这些重要的问题搞懂，岂不甚好O(∩_∩)O
现在来看看咱们在双目立体匹配上的学习进度：
你可能会想，怎么视差优化算法还没有完？ 是的，因为我们下一次将学习Jon Barron等在2016年ECCV上发表的算法The Fast Bilateral Solver，简称FBS，它是The Fast Bilateral-Space Stereo算法的自然衍生，但里面的数学知识和优化方法都有了很多扩展，得到的效果也很不错。而且，除了立体匹配，还能用于包括图像分割、画笔绘制等在内的很多计算机视觉问题。让我们拭目以待吧
本文同步发表在我的微信公众号和知乎专栏“计算摄影学”，欢迎关注，转载请注明作者和来源，觉得有用的话顺手点个赞哦
七. 参考资料
[*]Jon Barron主页: https://jonbarron.info/
[*]原始论文：Fast Bilateral-Space Stereo for Synthetic Defocus
[*]作者的附加材料，给出了推导过程，很难懂
[*]很好的一份国外大神的解读：Ugo Capeto. On fast bilateral-space stereo for synthetic defocus bybarron et al. https://www.dropbox.com/s/ke46mg8pkyou7s2/fast_bilateral_space_ stereo.pdf?dl=0. Accessed September 28, 2016.
[*]双边网格的论文：Jiawen Chen, Sylvain Paris, and Frédo Durand. Real-time edge-aware image processing with the bilateral grid. In ACM Transactions on Graphics (TOG), volume 26, page 103. ACM, 2007
[*]74. 三维重建9-立体匹配5，解析MiddleBurry立体匹配数据集
[*]手机中的计算摄影1-人像模式(双摄虚化)
[*]手机中的计算摄影6-旷视技术开放日上展示的手机电影
写的很好，很受启发！ 大佬你好， 请问如果只用单目， 我录一个1s的带有motion的video， 是不是也可以用这个方法恢复深度， 相当于我用单目的移动提供了baseline 单目深度重建有很多流派和方法，我之后会讲到。具体到你的问题，通过运动视频恢复深度是有可能的，只不过具体实操时会遇到一些坑。比如，如果依赖于运动来恢复深度，那么当运动幅度太小时，重建误差会很大 确实， 我把这篇源码改成单目video的了， 效果不是很好， baseline很小+生成的depth图含有很多原图纹理， 跟它的data loss 过于疏松有关。 是的
页:
[1]