原理+代码详解 | 稠密重建之SGM/tSGM算法

2021-04-28 16:14 作者:计算机视觉life 0人读过 | 我要投稿

注：本文来自计算机视觉life独家课程视觉三维重建：原理剖析+逐行代码详解中的课件及注释代码。

目录：

立体匹配算法介绍

SGM算法

tSGM算法

SGM/tSGM代码实现

立体匹配算法介绍

全局立体匹配算法

全局立体匹配算法主要是采用了全局的优化理论方法估计视差，建立全局能量函数，通过最小化全局能量函数得到最优视差值；
通过二维相邻像素视差之间的约束（如平滑性约束）而得到更好的匹配效果，但是对内存的占用量大，速度慢不适合实时运行。主要的算法有图割（graph cuts）、信念传播（belief propagation）、动态规划等算法。

局部立体匹配算法

主要是采用局部优化方法进行视差值估计，局部立体匹配算法有 SAD，SSD 等算法，与全局立体匹配算法一样，也是通过能量最小化方法进行视差估计，但是在能量函数中，只有数据项，而没有平滑项；
该算法由于每个像素计算互不干扰可以并行计算，所以可以实时。但由于所基于的局部窗口视差相同的假设在很多情况下并不成立导致匹配效果较差。

半全局立体匹配算法SGM

综合上述局部和全局算法的优缺点，半全局算法依旧采用全局框架，但是在计算能量函数最小化的步骤时使用高效率的一维路径聚合方法来代替全局算法中的二维最小化算法，使用一维最优来近似二维最优，得到的视差图在效果上和全局算法没有太大的差别，但是算法效率却有非常大的提升。

SGM算法

参考文献：

[1]Stereo Processing by Semi-global Matching and Mutual Information
[2]SURE: Photogrammetric Surface Reconstruction from Imagery

SGM算法详细介绍

匹配代价计算

文献[1]中介绍的SGM的代价计算是基于互信息（Mutual Information，MI）的匹配测度计算算法来计算匹配代价，互信息是一种对影像明暗变化不敏感的相关性测度。但由于原理复杂且计算需要迭代效率比较低，在实际应用中，更简单有效的方法如Census变换，故在此不再介绍MI。
Census变换：

使用像素邻域内的局部灰度差异将像素灰度转换为比特串即为census值；
基于Census变换的匹配代价计算方法是计算左右影像对应的两个像素的Census变换值的汉明（Hamming）距离（两个比特串的对应位不相同的数量：先进行异或运算，再统计运算结果中1的个数）

代价聚合

采用全局立体匹配算法，即找到每个像素的最优视差使得整体能量最小。能量方程如下：

$%0Aenergy(X%2CD)%3D%5Csum_%7Bi%7DDataCost(d_i%2Cx_i)%2B%5Csum_%7Bj%3Dneighbor(i)%7DSmoothCost(x_i%2Cx_j)%0A%09%09%09%09%09$

式中， $C$ 为匹配代价，第二项和第三项是平滑项，我们期望的是视差是连续的。所以如果当前像素 $x_b$ 与邻域像素 $x_N$ 视差相差比较小（1个像素）我们会给一个比较小的惩罚 $P_1$ ,如果大于一个像素则给一个大的惩罚 $P_2$ 。但是实际场景中肯定会有一些视差不连续区域相差比较大（比如：前景和背景）如图示：

为了处理这种情况，我们对 $P_2$ 进行调整：

$P_2%3DP%5E%7B'%7D_2*(1%2B%CE%B1*e%5E%7B(-(I_%7Bx_b%7D-I_%7Bx_N%7D)%5E%7B2%7D%2F(2*%5Cbeta%5E%7B2%7D))%7D)$

思考：为什么这样调整第二个惩罚项？

答案：如果像素和它的邻域像素亮度差很大，那么该像素很可能是位于视差非连续区域，则一定程度上允许其和邻域像素的视差差值超过1个像素，对于超过1个像素的惩罚力度就适当减小一点。

具体求解过程中，SGM路径代价聚合的思路：

将像素所有视差下的匹配代价进行像素周围所有路径上的一维聚合得到路径下的路径代价值，然后将所有路径代价值相加得到该像素聚合后的匹配代价值。

一维聚合路径示意图（图中各个箭头方向就是聚合各个路径）：

某一维路径代价计算公式如下：

$%5Cbegin%7Baligned%7D%0AL_%7B%5Cmathbf%7Br%7D_%7Bi%7D%7D%5Cleft(%5Cmathbf%7Bx%7D_%7Bb%7D%2C%20d%5Cright)%3D%26%20C%5Cleft(%5Cmathbf%7Bx%7D_%7Bb%7D%2C%20d%5Cright)%2B%5Cmin%20%5Cleft(L_%7Br%7D%5Cleft(%5Cmathbf%7Bx%7D_%7Bb%7D-%5Cmathbf%7Br%7D_%7Bi%7D%2C%20d%5Cright)%2C%5Cright.%5C%5C%0A%26%20L_%7B%5Cmathbf%7Br%7D_%7Bi%7D%7D%5Cleft(%5Cmathbf%7Bx%7D_%7Bb%7D-%5Cmathbf%7Br%7D_%7Bi%7D%2C%20d-1%5Cright)%2BP_%7B1%7D%2C%20%5C%5C%0A%26%20L_%7B%5Cmathbf%7Br%7D_%7Bi%7D%7D%5Cleft(%5Cmathbf%7Bx%7D_%7Bb%7D-%5Cmathbf%7Br%7D_%7Bi%7D%2C%20d%2B1%5Cright)%2BP_%7B1%7D%2C%20%5C%5C%0A%26%5Cleft.L_%7B%5Cmathbf%7Br%7D_%7Bi%7D%7D%5Cleft(%5Cmathbf%7Bx%7D_%7Bb%7D-%5Cmathbf%7Br%7D_%7Bi%7D%2C%20i%5Cright)%2BP_%7B2%7D%5Cright)%20%5C%5C%0A-%26%20%5Cmin%20_%7Bk%7D%20L_%7Br%7D%5Cleft(%5Cmathbf%7Bx%7D_%7Bb%7D-%5Cmathbf%7Br%7D_%7Bi%7D%2C%20k%5Cright)%0A%5Cend%7Baligned%7D$

式子中 $L%7B%5Cmathbf%7Br%7D%7Bi%7D%7D%5Cleft(%5Cmathbf%7Bx%7D%7Bb%7D%2C%20d%5Cright)%3CC%7Bmax%7D(x_b%2Cd)%2BP_2$

则最终代价值（所有路径之和）：

$S(x_b%2Cd)%3D%5Csum_%7Br_i%7DL_%7B%5Cmathbf%7Br%7D_%7Bi%7D%7D%5Cleft(%5Cmathbf%7Bx%7D_%7Bb%7D%2C%20d%5Cright)$

视差计算
代价聚合后，每个像素直接找聚合代价最小对应的视差值就是我们所要求的视差值。
视差优化

视差一致性检查
将左右图互换，得到R-L视差图。与L-R视差图对比，根据左图的视差找到在右图中的对应视差，如果两者小于阈值则认为是准确的，反之是错的把该值剔除。如下图：

亚像素计算
上述计算的视差图是像素级别的，在实际使用中精度是无法满足需求的。SGM通过在上述计算的最小代价的视差层附近进行插值找到亚像素级的精度。如下图：

tSGM算法

与SGM基本相同，区别主要是在代价聚合的时候：

使用金字塔图像计算视差（由粗糙到精细即从低分辨率到高分辨率计算匹配代价）
每个像素的视差范围都不同，只在真值附近搜索大大减少搜索范围和内存占用，如下图：
每个像素的视差范围计算方法：
如果当前像素的视差值$d(x)$无效，在搜索窗口 $win$ =31, 反之 $win$ =7
以当前像素为中心，窗口大小 $win$ 搜索所有有效视差存储在 $max$ 中， $disps$ 的大小是 $numdisp$ ,中值是 $disps$ ,最大值,最小值 $min%0A$ :
$%5Bd_%7Bmax%7D%2Cd_%7Bmin%7D%5D%3D%5Cbegin%7Bcases%7D%0A%5Bmin(width*2%2F3%2C32)%2C%20-min(width*2%2F3%2C32)%5D%2C%26numdisp%3C3%5C%5C%0A%5Bdisp%2B(numDisp%2B1)%2F2%2C%20disp-numDisp%2F2%5D%2C%263%3D%3Cnumdisp%3CmaxNumDisp%20%5C%5C%0A%5Bdisp%2B(maxNumDisp*(max*2%2B1-disp)%2B1)%2FnumDisp%EF%BC%8C%5C%5Cdisp-(maxNumDisp*(disp-min*2)%2B1)%2FnumDisp%5D%20%26numdisp%3E%3DmaxNumDisp%5C%5C%0A%5Cend%7Bcases%7D$
代价聚合时，SGM由于每个像素视差范围都一样，所以各个搜索路径都有对应的视差，tSGM由于每个像素视差范围不一样有的可能都没有重叠范围，所以之前的代价聚合计算方法需要调整，把没有重合的情况考虑进去。
调整如下：
$if%20%5Cquad%20d%3Ed_%7Bmax%7D(x_b-r_i)$
$L_%7Br_i%7D(x_b%2Cd)%3DC_%7Br_i%7D(x_b%2Cd)%2BP_2$
$if%5Cquad%20d%3Cd_%7Bmin%7D(x_b-r_i)$
$L_%7Br_i%7D(x_b%2Cd)%3DC_%7Br_i%7D(x_b%2Cd)%2BP_2$
$else%3A$
$L_%7Br_i%7D(x_b%2Cd)%3DSGM$