The Fast Bilateral Solver 快速双边求解器 双边网格可视化展示

The Fast Bilateral Solver 快速双边求解器

一. 前言

双目相机拍到的原始图后计算出来的深度图a，右图则是通过FBS算法得到的图像b。

使用这个算法优化计算后的深度图我们还需要两个图像信息，一个是置信图c，c中的元素越接近1，a中相应元素的正确概率就越高。还有是原始图像的RGB像素图d。

二. 算法的原理与优化步骤

2.1 创建双边空间优化全局代价函数

下面文章作者构建的代价函数，加号前面一项是平滑项，后面一项是数据项，接下来我们就是要最小化这个式子：
$$\underset{x}{\mathrm{minimize}}\frac{\lambda }{2}\sum _{i,j}{W}_{i,j}(x_i-x_j{)}^2+\sum _i{c}_i(x_i-t_i{)}^2\text{\hspace{1em}(1)}$$

• $x$ 像素空间的输出图
• $y$ 双边空间的输出图
• $c$ 输入的置信图
• $t$ 输入的待优化图像
• $W$ 双边平滑权重

平滑项$\frac{\lambda }{2}{\sum }_{i,j}^{}{W}_{i,j}(x_i-x_j{)}^2$中的 $x_i$ 和 $x_j$ 是一幅视差图中不同位置的像素的视差值。这个时候，我们要把视差图看作展开的一维向量。一个宽高分别为 w 和 h 的视差图，展开后就是$(w\ast h,1)$的一维向量，于是 $i$ 和 $j$ 就是不同像素的编号，或者坐标值。

$W$矩阵,是在YUV域的双边滤波矩阵,如果x是一个一维表示的尺寸为$(w\ast h,1)$的图像，我们可以用下面的公式来表示双边滤波：
$$y=\frac{Wx}{W\vec{1}}\text{\hspace{1em}(2)}$$
$y$即是滤波后的图像，为$(w\ast h,1)$的一维向量。
这里的矩阵$W$是一个对称的矩阵，尺寸是$(w\ast h,w\ast h)$。对于彩色图像的双边滤波，$W$矩阵的每个元素如下:
$$W_{i,j}=exp(-\frac{||(p_i^x,p_j^y)-(p_i^x,p_j^y)|{|}^2}{2{\sigma }_{xy}^2}-\frac{(p_l^x-p_l^y{)}^2}{2{\sigma }_{xy}^2}-\frac{||(p_i^u,p_j^v)-(p_i^u,p_j^v)|{|}^2}{2{\sigma }_{xy}^2})\text{\hspace{1em}(3)}$$
这么说还是太抽象了，让我们可视化一下，下面图像是一张 16x16 的黑白图像，每个像素用黑白方块表示。所以双边滤波矩阵的大小为$(256,256)$：

如上图所示，首先一张 16X16 的灰度图变为$(256,1)$的向量，左乘一个滤波矩阵，得到一个新的$(256,1)$的向量，reshape一下就得到了滤波之后的图像。大家看懂了吗？如果能回答下面的问题，就可以说是看懂了

1. 图中的滤波矩阵分别都是什么类型的滤波矩阵？
2. 两个滤波矩阵有什么不同？(颜色的亮度代表亲和力的强弱，越亮的代码亲和力越强)

2.2 由双边滤波到双边网格

好了，现在我们知道什么是滤波矩阵了，那么问题又来了，如果我们的图片是 1280X720 的图像，那将是一个巨大的计算灾难(1280x720,1280x720)。这时就用到了我们一开始提到的思想：将大规模问题转换为小规模问题。其中有一种称为双边网格的方法，这是2007年Chen Jiawen等人提出的。

依然将上面的 16x16 的黑白灰度图像中的每个像素，按照其x坐标、y坐标、灰度值I，被分组放到了一个小盒子中，这时候坐标和灰度值就相当于小盒子的门牌号。实际在我们操作时，不可能让小盒子中只放一个像素，将一幅2D的灰度图像中的每个像素。所以可能会让$(x,y,I)\in ((x_1,x_2),(y_1,y_2),(I_1,I_2))$的像素被分到同1个小盒子中，这时候小盒子中可能会存在一个或多个像素。还是很抽象我们依然使用上面$(16,16)$图像来说明：

红色框是$(x,y,I)\in ((x_1,x_2),(y_1,y_2),(I_1,I_2))$；蓝色框是小盒子
Splat矩阵记录了像素到小盒子的变化，所以它的转置矩阵就是小盒子返回到红色框。
之后我们称红色框为网格宽度，小盒子为格点，格点组成了新的向量，特别像数据聚类，从而实现了数据量的压缩。

现在我们仅仅是得到一个小规模的数据的变化矩阵，依然没有生成滤波矩阵，所以我们需要对每个格点生成滤波矩阵，每个格点中其实隐藏了3个纬度的信息，所以需要对三个纬度的信息分别做滤波，为了方便计算使用$[1,2,1]$滤波核

每个$B$矩阵每行不超过3个非零元素，在大多数文献中，双边网格是用双线性插值实现的，并且是“密集的”即使没有填充顶点，也会创建顶点，模糊矩阵近似为每个维度的模糊矩阵的外积。相比之下，本文的简化网格是“稀疏的”，没有分配像素的顶点永远不会被创建。这使得网格可以扩展到比密集的双边网格更高的维度。

请注意，因为我们的模糊操作是模糊的和而不是一个结果，每个模糊可以有效地并行计算。
为了更好地理解这些作为矩阵因式分解的双边表示，请参见图：

从上往下分别为网格宽度是8，4，2的双边网格矩阵。
由于双边网格的规模远远小于原始图像，比如一幅100万像素的图像投射到双边网格中可能只有10万个格子。而高斯滤波又比双边滤波的计算量低很多。因此通过这样的方式，就可以实现非常快速的双边滤波了。
总结一下：
把从原始图像到双边空间的过程称为Splat，在双边空间中进行滤波称为Blur，把滤波后的结果转换到原始图像像素空间中的过程称为Slice。 Splat-Blur-Slice，就是这类方法的三个重要步骤。由于高斯滤波是可以分解的，在双边网格中的高斯滤波可以变换为每个维度上（对于1维图像来说是space和intensity两个维度，对于RGB图像是x,y,r,g,b五个纬度）分别进行的一维滤波相加的结果（这是一种近似）:

$$W\approx \hat{W}=S^{T}BS\text{\hspace{1em}(4)}$$

• $S$ splat
• $B$ blur
• $S^T$ slice

现在我们终于可以回到主线任务中了，最小化公式$(1)$。

2.3 双随机矩阵

文中使用了一种迭代算法，该算法将矩阵的每一行和每列的$\infty$范数渐近地缩放为1(双随机矩阵)。这种缩放算法保持了原始矩阵的对称性，并以1/2的渐近速率显示出快速的线性收敛。对于1范数的情况，收敛是线性的，其速率依赖于缩放矩阵的谱。使用缩放算法通过减少旋转的需要来帮助直接求解。特别是对于对称矩阵。下面是双随机对称矩阵的求解步骤：

然后公式$(1)$整理为：
$$\underset{y}{\mathrm{minimize}}=\frac{1}{2}{y}^T(\lambda (D_m-D_n\overline{B}{B}_n)+diag(S_c))y-(S(diag(c)t){)}^T+\frac{1}{2}(diag(c)t{)}^{T}t，其中y=Sx\text{\hspace{1em}(5)}$$
然后我们用标准二次型重写他：
$$\underset{y}{\mathrm{minimize}}\frac{1}{2}{y}^{T}Ay-b^{T}y+c\text{\hspace{1em}(6)}$$
$$A=\lambda (D_m-D_n\overline{B}{B}_n)+diag(S_c)b=S(diag(c)t)c=\frac{1}{2}(diag(c)t{)}^{T}t$$
通过对该损失函数求导并将其设为零，我们可以看到最小化该二次形式等价于求解稀疏线性系统。

2.4 共轭梯度法求解

最后就是求解$Ax=B$。本文的求解算法是预条件共轭梯度法，我们首先来介绍一下共轭梯度法共轭梯度法（Conjugate Gradient），它是介于最速下降法与牛顿法之间的一个方法，它仅需利用一阶导数信息，但克服了最速下降法收敛慢的缺点，又避免了牛顿法需要存储和计算Hesse矩阵并求逆的缺点。多用于解维度高的稀疏线性方程组，其优点是所需存储量小，具有步收敛性，稳定性高，而且不需要任何外来参数。
共轭梯度法怎么理解？共轭梯度法对最优梯度法进行了修正，搜索方向为共轭方向，将负梯度方向旋转了一个角度，每次往最优方向需要在负梯度方向进行修正。
什么是共轭方向那？设$Q$是正定矩阵，若对于两不同的非零向量$d_i^T,d_j$满足：
$$d_i^{T}Q{d}_j=0\text{\hspace{1em}(6)}$$
则称$d_i^T,d_j$关于$Q$共轭，$d_i^T$与$d_j$是共轭方向。
再说说什么是共轭共轭向量组：如果一组方向 $d_0,d_1,...d_{n-1}$中任意两个满足向量$d_i,d_j$满足公式$d_i^{T}Q{d}_j=0$,即为共轭向量组，且共轭向量组中的向量一定线性无关。

下面分别为梯度下降法，共轭梯度法，预处理共轭梯度法

对于n阶线性方程组，通常最多n+1次迭代可以获得收敛。但是在极端大的线性方程组，此时，采用共轭梯度法，有时可能需要上万次的迭代才能收敛，虽然每次收敛的计算量并不大，但是整体求解也会花费较多时间。此时共轭梯度法能否快速收敛与系数矩阵A密切相关。
对于共轭梯度法存在以下定理：
如果$A=I+B$为 nxn 的对称正定矩阵，且$rank(B)=r$，则共轭梯度法至多$r+1$步收敛。其中I指单位矩阵，rank是矩阵的秩。
从该定理可知，共轭梯度法能否快速收敛，主要取决于$A$是否“接近于”单位矩阵。特别地，当$A$就是单位矩阵时，$rank(B)=0$，此时一步即可收敛，当然这是很显然的。基于上述定理，多种基于共轭梯度法法衍生的预处理共轭梯度法（PCG）得以广泛地应用
采用不同的M，最终影响共轭梯度法的收敛性。很显然，对于整个原始方程的PCG求解效率来说，首先这个预处理方程需要比较容易求解，否则求解该预处理方程就已经花费了较多时间，则整个原始方程求解效率自然不会高。其次，M的选取需要确实能够减少迭代步数。这两方面在实际中往往是冲突的。
对角线预处理：对角线预处理指的就是M矩阵为原始系数矩阵A的对角线矩阵。一般情况下，这种方法求解预处理方程十分简单，但是通常对于系数矩阵严格对角占优（即对角线元素大于该行其他元素之和）的情况下才比较有效。

三. 应用

双边求解器是一种灵活和快速诱导边缘感知平滑的技术。FBS的不仅可以用于优化初始的视差图，还能做视差图或深度图图的超分辨率重建。这里有一个低质量的输入视差图，以及一个参考图像，通过FBS可以得到高质量的输出图像。

另外一个应用是用户引导的灰度图像上色。这里也有低质量的输入图像y，即用户指定的最终色彩图像的骨架，原始的灰度图像则作为参考图像。再指定一个置信度图c，我们一样通过FBS求解输出彩色图像x，并且要求x自身相对参考图像来说是双边平滑的，同时还在高置信度的区域和用户输入y尽可能一致。

参考文档

[1] Barron, Jonathan T，B Poole ，The Fast Bilateral Solver， DOI：10.48550/arXiv.1511.03296

[2]JT Barron，A Adams，YC Shih，C Hernandez Fast bilateral-space stereo for synthetic defocus DOI : 10.1109/CVPR.2015.7299076

[3]J Chen，S Paris，F Durand，Real-time edge-aware image processing with the bilateral grid DOI：10.1145/1276377.1276506

[4]BORA，UAR，DANIEL，RUIZ，PHILIP，A.，KNIGHT，A Symmetry Preserving Algorithm for Matrix Scaling https://doi.org/10.1137/110825753

[5]JL Nazareth，Conjugate gradient method DOI：10.1002/wics.13