双目图像重叠的视差计算_一文读懂经典双目稠密匹配算法SGM

水深无声 2023-01-01 08:53 113阅读 0赞

![d88830f0aafa760422bd479a8a0f989e.png][]

题图为SGM算法的一个处理结果。

最近来看看一些双目稠密匹配的算法。说来惭愧，SGM在航测领域是很重要的算法（当然也是最好的双目稠密匹配算法之一），自己却没有认真读过，只是大致有些了解。看了论文，再结合网上一些资料，自己做了些论文笔记。想到关于SGM论文网上还没看到比较翔实的博客，就把自己做的笔记再加些解释分享出来了（下文中的引用部分多为我自己的思考）。

其中还有一些是自己不理解的地方，准备后续写代码时研究。也希望大家可以分享自己的见解，或对文章内容进行批评指正。

由于是论文笔记，所以有些中英夹杂，而且有些想想也不好翻译，就不再改了。闲话少说，直接来看论文吧。

## 基础 ##

双目图像稠密匹配的**4个基本步骤为：**

1.  **Matching cost computation**;
2.  **Cost aggregation**: connects the matching cost within a certain neighborhood;
3.  **Disparity computation**: selects the disparity with the lowest matching cost;
4.  **Disparity refinement**: removing peaks, interpolating gaps or increasing the accuracy by sub-pixel interpolation.

这4个都比较好理解，SGM论文也是按照这样来组织的。下面就按照论文的顺序来详细了解下SGM算法。

## 中心思想及要求 ##

**idea**: 使用MI (Mutual Information) 来进行单像素匹配 + 多个一维平滑约束（来拟二维约束）来进行“全局”优化。

**前提**: 已知立体像对间的对极几何关系。

对于普通标定好的双目相机，对极约束为：左图的待匹配像素 p，右图上的对应像素

![Image 1][] ，

![Image 1][] 表示极线方程。

## Pixelwise cost calculation ##

SGM算法不用图像块进行匹配，只考虑当前像素。因为利用图像块进行匹配对应的隐性约束为块内像素的视差是相同的，而这在深度变化（物体边界）的地方是不成立的。

互信息（MI）：defined from the entropy（熵）of two images:

![Image 1][]

其中，

![Image 1][]

> 显然
> 
> ![Image 1][] 就是归一化直方图，
> 
> ![Image 1][]是两张图像的联合分布直方图（左右像素对的灰度值对出现的频率）。
> 
> ![Image 1][]
> 
>   
> 相等时T\[\] = 1。这显然是一个统计，若像素p 及其对应像素的灰度值为(i, k)，则相关量+1。依据此遍历整个左影像进行统计。  
> 论文原文：Calculation of P(I1, I2) is done by counting the number of pixels of **all combinations** of intensities, divided by the number of all correspondences.

上式针对整张图像而非单个像素，不能用作 cost。因此，对于一张图像的联合熵，有其他论文利用泰勒展开计算：

![Image 1][]

其中，公式内第一项为像素 p 的灰度值，第二项为匹配像素的灰度值。

> 这里，联合熵被简化为左图所有像素点（及其对应点）的 h(i1,i2) 之和。

而

![Image 1][] 的定义为：

![Image 1][]

> i, k 为像素灰度值，n 是匹配的左右像素对的数量 *（如何得知？）*。将联合分布直方图做2D高斯平滑后取对数再做平滑（高斯核的大小建议为 7\*7）。  
> 联合分布直方图是 **固定大小**的 256 \* 256 图像，因为只有256个灰度值。整个处理过程如下图所示：

![99153938fd29170ea5b37f7853b74edb.png][]

> 可以看到，由于极线矫正后左右图十分相似，所以得到的联合分布直方图类似于一个对角矩阵。

以上计算并不包含被遮挡的像素*（如何得知哪些像素是被遮挡的？）*。因此，为了避免包含遮挡像素，建议将熵也这么计算：

![Image 1][]

此外，以上计算也只针对两张图像的重叠部分。

因此，一维直方图也可以这么计算：

![Image 1][]

> h(i) 即为 P(I1, I2) 第 i 行（列）像素之和。

最终，互信息的**定义**为

![Image 1][]

> 个人觉得，互信息之所以对于光照、噪声等干扰较为鲁棒在于：  
> 1. 极线匹配限制了搜索范围；  
> 2. 像素（对）出现的频率是一个比较鲁棒且全局的特征。论文提到决定 cost 的主要是联合熵（图像的熵近乎常数）。由上图可以看到联合熵基本就是个对角矩阵，亦即，拿到左图像素灰度值p之后，右图对应像素的灰度值基本可以确定就在p的附近，与p差距很大的可能性很小；  
> 3. 额外一点，就是后续对于所得视差的连续性有约束，并且对无纹理区域有后处理（这点后面再说）。

对于像素点 p，若取其视差为 d，则对应的 cost 为：

![Image 1][]

**问题**：想要对匹配图像 Im 进行视差矫正就需要视差图。但我们的目标就是获取视差图。

**解决**：

1.  迭代法：start with a random disparity map for calculating Cmi. And use the cost for matching and calculating a new disparity map. 论文表示大致迭代3次就可以了。
2.  层次法：recursively use the **up-scaled** disparity map (**half resolution**). Start with a random disparity map of 1/16 resolution. 在这个1/16的视差图上重复计算3次（每次迭代3次，然后每放大一倍迭代3次计算出新的视差图。

> 考虑到图像中具有大量像素（图像大小），即使是随机初始化的视差图也可以得到较好的概率分布结果（P）。因此不需要太多次的迭代计算。  
> 既然随机初始化的视差图也可以，那不如直接先对整张影像做一个 **correlation-based matching**.

## Aggregation of costs ##

以上就是论文进行匹配的第一步。接下来就要进行一个“全局”上的优化。

pixelwise matching 还是不够稳定，因此需要加上一些约束来保持同一平面的像素具有相同的视差并惩罚邻域内视差不连续的像素（即平滑）。

所以最小化如下能量方程：

![Image 1][]

上式第一项为所有 MI cost 的和；第二项对 像素p的邻域内出现视差与p的视差相差1 的像素 加上一个惩罚项P1；第二项对 视差相差大于1 的像素加入更大的惩罚P2.

> 后两项的作用相当于机器学习中的正则化项。  
> 由于物方空间中深度不连续的表现之一为像素灰度变化（即有梯度），因此有时
> 
> ![Image 1][]。

**作用**：对于小的视差相差使用较小的惩罚以适应斜面或曲面；对大的视差相差使用大的惩罚来防止深度（视差）不连续。

问题在于全局最小化 E(D) 是一个 NP 完全问题，很难解算。单行（1D）约束可以利用**动态规划**达到多项式时间。

因此很自然地想到优化多个单行的约束条件来拟合2D优化，作者建议至少要选择待优化像素的8个方向。

The aggregated smoothed cost S(p, d) for a pixel *p* and disparity *d* is calculated by summing the costs of all 1D minimum cost paths that end in pixel *p* at disparity *d*.（也就是说S是所有L之和。）

另

![Image 1][] 来表示像素 p，选择视差为d，方向为 r 上的costs:

![Image 1][]

> 意即，当前代价 + 该方向上上个像素的最小代价：视差不变时的代价，视差差1时的代价，视差之差大于1的代价三者之一。  
> These paths through disparity space 投影到左图上是一条直线；投影到右图上则表现为视差变化的折线。

![1e3fb524bcdd0784f5904241444dee57.png][]

考虑到上述代价最后可能非常大，可将之修改为：

![Image 1][]

> 即，减去该方向上，上个像素的最小路径代价。

注意到被减项对于该方向上，像素 p 的**任意视差**都是一个**常数**（上一个像素的视差及其对应的最小cost已经计算好）。

最终，对于像素 p，选择视差 d 的cost 就是：

![Image 1][]

选择使 S(p, d) 最小的视差d. 若要到亚像素精度，可以利用 neighboring costs 来拟合一个二次曲线。

> 问题：  
> 每个 d 都试一遍？ 边缘像素怎么办？ 路径多长？（有多长算多长？）  
> （论文）可以先将图片在一个方向上计算完，保存起来。算下一个方向时再加上去。

## Multi-baseline matching ##

对左图计算完视差图后，还可以以同样的方法来对右图进行视差图的计算，在将二者融合起来。也可以利用多个立体像对来计算多张稠密视差图。

**问题**：areas which are not seen by **all images** will become invalid.

**解决**：取多视差图的**并集**而非交集。

> 对于左图来说，它的视差图是恒定值（参考依据视差计算深度的公式），因此可以通过多图匹配来优化。

对左图计算一个视差图后，也可对**右图**计算一个视差图。然后进行**一致性校验**：

![Image 1][]

如此，确保 one to one mapping (**uniqueness constrain**).

## Disparity Refinement ##

之后，是对所得的视差图进行优化处理。个人觉得这是双目稠密匹配很重要的一步，因为对于很多算法，其得到的视差图之间的差别并不会特别大。重要的是对所得的视差图进行后处理，提高深度计算的精度。

## 移除极值 ##

采用**联通域分割法**，将视差相差在1个像素以内的像素4邻域归为一个联通域来分割所得的视差图。

根据场景事先确定一个大小阈值，去除像素数量过少的联通域。

## Intensity Consistent Disparity Check ##

这步是SGM后处理比较核心的一步。

室内场景的前后景之间会有视差不连续。但前述的能量方程对于发生视差不连续的位置没有偏好，这种不连续可能在错误的地方出现，因此会造成错误的前后景物体边界，或者一个斜面。

先前的将**像素梯度**与惩罚值P2关联可以部分解决该问题。但在一些视差不连续的**不对称**的情况（如前景物体没有完全出现，只在出现的一边产生视差不连续）下可能失败*（原因待解）*。

室内环境中，后景常常为无纹理区域，如墙等。为解决这个问题，先提出以下三个**假设**：

1.  无纹理区域内部没有视差不连续现象（即视差变化伴会随梯度变化）；
2.  无纹理区域包含一些弱纹理（否则也无法进行像素匹配）；
3.  无纹理区域的表面可以用一个平面表示（参考1，如此则区域内视差恒定。否则为两个可以分割开的区域）。

**解决步骤：**先对图像进行分割，得到前后景区域。问题就变为如何正确选择后景（无纹理区域）的边界和视差，较好的分割出前景物体。

1.  利用 fixed bandwidth Mean Shift Segmentation 方法分割原始图像，得到各个无纹理区域（假设1），忽略过小区域，所得区域表示为
    
    ![Image 1][]。
2.  利用联通域分割法分割步骤1得到的无纹理区域（假设2），忽略过小区域，所得分割区域表示为
    
    ![Image 1][]；
3.  针对2得到的区域拟合平面
    
    ![Image 1][] （假设3）。将得到的各个可能平面推广到所在的整个无纹理区域，取能量方程最小的平面。
    
    1.  遮挡判断子步骤：对于像素p，若其匹配像素q有了另一个匹配，且新的匹配像素的视差大于像素p，则像素p认为被遮挡；
4.  该无纹理区域内的视差都替换为最佳拟合平面所对应的视差；

![Image 1][]

在这种情况下，只有大小足够的无纹理区域会被修正。

> 解释一下，对于后景，作者做了两个假设：一是其为一个较大的均匀区域，内部梯度变化很小；二是这是一个平面。  
> 对于深度变化（物体边界），作者的假设为前后景的外观不一致，因此二者之间存在较大的梯度。  
> 因此，利用前面的假设1，对原始图像进行分割时，就可以比较好的分割出前后景，而且认为 **这个分割结果是正确的**。而且，如果前景物体较小，或者前景物体纹理复杂被分割为多个图像块时，会被忽略。只有较大的后景区域会被处理。  
> 以这个分割结果去找出视差图对应的无纹理区域
> 
> ![Image 1][] ，考虑到假设2，后景区域内部可能也会有微小的视差（因为有微小的梯度），因此，容忍其有1个像素的视差相差，再在视差图上对该无纹理区域进行视差上的分割。
> 
>   
> 由于假定由原始图分割得到的无纹理区域
> 
> ![Image 1][] 的边界是正确的，再考虑到前述的问题：该区域内包含了错误的前后景边界（即较大的视差变化，该变化应该出现在
> 
> ![Image 1][] 的边界处），因此我们要对这个区域内的视差进行校正。结合假设3，该区域内的视差要被校正为一个恒定值，最好的办法就是根据视差分割的结果拟合多个平面，取让能量方程最小的那个视差。

## Discontinuity Preserving Interpolating ##

后处理完的视差图可能因为遮挡或误匹配而存在缺失区域。对这些区域的视差进行插值前要对其性质进行分类，以使用不同的插值方法：

1.  **遮挡**：用背景的视差进行外插值，而不使用前面的遮挡物的视差；
2.  **误匹配**：可以利用其周围像素进行内插值。

**注意**：对与遮挡相连的误匹配区域，同样利用背景像素进行外插值。

**区分**：由于遮挡会照成视差不连续。被遮挡位置的像素的极线不与disparity function（即视差随位置变化的方程）相交；误匹配像素则相反。

![472e9db08b43e45b7f6b084b90dfcd72.png][]

如上图，点p1处就是因为遮挡而造成点缺失；而点p2则是因为误匹配。

插值仍然是从像素周围的8个方向计算出8个可能的视差值，然后

![Image 1][]

该式的第一项是选择待选视差中第二小 (second lowest) 的视差。第二项为选择待选视差的中值，这样**有利于**在物体边界处保持视差不连续。

> 该插值方法独立于立体相对匹配算法。

最后可以再用**中值滤波**对得到的视差图进行去噪处理。

最后，作者还介绍了对于超大幅的航空影像进行分块处理的方法和加速技巧，以及解释了文章中不先对右影像整张影像进行极线矫正的原因（对于推扫式的航空/航天影像，其极线为曲线，因此无法进行矫正，只能实时计算）。但对于其他普通立体像对，这些内容都是不必要的，就不再解释了。

如果觉得有帮助，还请点个赞 : )

## 参考 ##

Stereo Processing by Semi-Global Matching and Mutual Informationcore.ac.uk Accurate and Efficient Stereo Processing by Semi-Global Matching and Mutual Informationwww.robotic.dlr.de Stereo Matching文献笔记之（九）：经典算法Semi-Global Matching（SGM）之神奇的HMI代价计算~blog.csdn.net

![f97f7969cce6d3a096fcace69a912a14.png][]

Stereo Matching文献笔记之（十）：经典算法Semi-Global Matching（SGM）之碉堡的动态规划~blog.csdn.net

![1447e19ecc54900bc0c7377cac1cc549.png][]

[d88830f0aafa760422bd479a8a0f989e.png]: /images/20221120/3d34f2103b854c919ccac6587937fefc.png
[Image 1]: 
[99153938fd29170ea5b37f7853b74edb.png]: /images/20221120/b63beaa7272246fe9a3ed9b34eba4db3.png
[1e3fb524bcdd0784f5904241444dee57.png]: /images/20221120/d769bd20649d4575bb5d2356341b4bed.png
[472e9db08b43e45b7f6b084b90dfcd72.png]: /images/20221120/7afcd9e87b97411b8c887a79733e0cfe.png
[f97f7969cce6d3a096fcace69a912a14.png]: /images/20221120/75802282927443359b565ef7944ec0c7.png
[1447e19ecc54900bc0c7377cac1cc549.png]: /images/20221120/3e03610bf1a74d0c8f3eac2ef157b207.png