一. 前言
到目前为止，我的双目立体匹配系列文章已经完成了接近经典视差优化算法的介绍，让我们观察一下现在的进度：
今天的文章，我们先从回顾立体匹配的评价指标开始，我在文章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就是不同像素的编号，或者坐标值。
最神秘的是这里的矩阵，为了理解它，我们首先得重新回到双边滤波器的矩阵表达，如果x是一个一维表示的尺寸为WHx1的图像，我们可以用下面的公式来表示双边滤波，等号坐标是滤波后的图像y，尺寸还是WHx1：
双边滤波矩阵表示法
这里的矩阵A是一个对称的矩阵，尺寸是WHxWH。对于彩色图像的双边滤波，A矩阵的每个元素如下：
双边滤波权重
它实际上是表示第j个像素和第i个像素间的在空间和颜色上的相似度。如果两个像素越靠近，并且颜色越相似，其对应的A元素就越大，当然如果是单色图像那么就只是计算空间和亮度上的相似度。
我们看看下面这个3x3的小图像（我们假设其中存储的就是对应像素的视差）
3x3的小图像
按照上面的公式很容易计算得到其双边滤波矩阵如下
权重矩阵A示意
利用这个矩阵，可以把任何我们认为相邻的像素的视差值加权平均到一起，得到最后的滤波结果。而我们前面数据项中的 矩阵则是双边滤波矩阵A的某种特殊归一化形式，其最大特点是每一行和每一列的和都为1。
4.2 平滑项转换到双边网格空间
我们可以用下面的公式表示双边网格算法，其中S矩阵代表我最开始提到的Splat操作，将像素值投射到双边空间。B矩阵代表双边空间中进行的模糊操作，而ST矩阵则代表上面提到的Slice操作，将双边空间中的滤波结果重新投射到原始像素空间。在双边网格的构造过程中，每个网格的顶点采用了最近邻方式采样，所以这种方式是对原始双边滤波的一种近似：
双边网格是矩阵分解过程
进行变量替换后，双边滤波就变成了下面这样：
在双边空间中进行滤波
而作者最精妙的一点是，把匹配代价中的平滑项按照类似的方法，投射到了双边空间。为了做到这一点，首先是把原始的平滑项写成矩阵表达形式, 再通过变量替换转换到双边空间
下面对双边网格中的各个变量作出解释，请注意在双边空间中的这些矩阵和向量的规模都是大大低于像素空间中的对应实体的，这也为之后的快速求解打下了基础。
你可能对为什么Jon Barron能够把原始像素空间中的平滑项做上面所说的转换感兴趣。这里面有非常复杂的数学推导过程，都记录在了其论文的附加材料里面。不过很多人直接阅读此份附加材料也是看不懂的（也许你就在这些人里面
