SGM(Semi-Global Matching)算法笔记

SGM

论文：《Stereo Processing by Semi-Global Matching and Mutual Information》

半全局双目匹配算法：
逐像素匹配；用互信息来做匹配代价；用多个一维平滑约束近似二维平滑约束进行“全局”优化。
【全局匹配需要考虑所有像素，用一个全局能量函数；局部匹配需要针对局部区域做代价聚合；SGM不是只考虑像素局部区域，也没有考虑所有的像素，只考虑了非遮挡点。】

逐像素匹配由平滑约束支持，该通常表示为全局代价函数。 SGM通过从各个方向的逐路径优化来快速近似全局代价函数。除此之外文章还讨论了遮挡部分检测，子像素优化和多基线匹配。后处理包括去除离群点，缺失点插值等。处理大图像用正交投影融合视差图。

一.逐像素计算代价

1.

输入图像必须知对极几何关系，可以是未矫正的(有的不容易矫正，推扫式图像)；
匹配是用两张图，这里左图是 $b$b (base Image)，右图是 $m$m (match Image)，左图的一个像素位置的p的灰度用 $I_{bp} $代表，它在 $m$m 的可能匹配的点q用 $I_{mq}$Imq​ 表示。$d$d 是极限参数，$q=e_{bm}(p, d)$q=ebm​(p,d) 是p和q之间关于极线的关系。如果两张图矫正了则$e_{bm}(p, d) = [p_x − d, p_y]^T$ebm​(p,d)=[px​−d,py​]T, $d$d 是视差；
匹配区域大小和形状对结果比较重要，一般考虑范围越大鲁棒性会越好，但是在区域内视差连续的假设在不连续的地方其实是不成立的，所以大区域会让边界处和精细结构处变模糊。本论文的方法就只考虑了 $p$p 和 $q$q 两个点来计算匹配代价；
匹配代价用互信息，对光强变化不敏感，互信息用两张图的熵 $H$H 和他们的联合熵定义。

2.什么是熵、联合熵、互信息？

互信息： $MI_{I1,I2} = H_{I1}+H_{I2}-H_{I1,I2}$MII1,I2​=HI1​+HI2​−HI1,I2​
熵： $H_{I}=-\int_{0}^{1} P_{I}(i) \log P_{I}(i) d i$HI​=−∫01​PI​(i)logPI​(i)di

联合熵： $H_{I_{1}, I_{2}}=-\int_{0}^{1} \int_{0}^{1} P_{I_{1}, I_{2}}\left(i_{1}, i_{2}\right) \log P_{I_{1}, I_{2}}\left(i_{1}, i_{2}\right) d i_{1} d i_{2}$HI1​,I2​​=−∫01​∫01​PI1​,I2​​(i1​,i2​)logPI1​,I2​​(i1​,i2​)di1​di2​

熵：是随机变量的不确定性，代表信息量，熵越大信息量越大；
联合熵：是变量之间的不确定性，代表多个随机变量组成的变量系统的信息量，联合熵越大代表这个组合信息量越大；
互信息：是变量间相关性，可以看成一个随机变量中包含的关于另一个随机变量的信息量，所以互信息越大代表变量之间重复的信息更多，就是匹配的更好；
很容易理解匹配的结果好，联合熵会比较低，因为一张图可以用另一张预测出来(知道左图每个像素对应右图的位置，直接扭曲过去就行，匹配的好基本位置都能扭曲对)，这样信息量就比较低。

3.概率P是什么意思？

灰度直方图是关于灰度级分布的函数，是对图像中灰度级分布的统计；
熵里的P也是根据直方图计算的，把图像里所有出现的像素数值都统计出它出现的次数，这个次数除以总像素数就是熵公式里的 $P_{I1}$PI1​ ，图像里像素范围是0-255，概率分布就是一维的256个数的概率；

4.SGM怎么计算联合熵，里面的参数是什么意思？

联合熵里的 $P_{I1,I2}$PI1,I2​ 相当于变成了二维的图像对，范围是(0,0)-(255,255), 概率分布是256*256的，计算的时候一对一对计算像素对出现的频数再除以像素对数。
这里的联合熵利用泰勒展开转化成像素和的形式近似，就是求所有$h_{I1,I2}的和$hI1,I2​的和，自变量$I_{1p}$I1p​和$I_{2p}是灰度对$I2p​是灰度对。

$H_{I1,I2}=\sum_Ph_{I1,I2}(I_{1p},I_{2p})$HI1,I2​=∑P​hI1,I2​(I1p​,I2p​)

$h_{I1,I2}$hI1,I2​ 计算方式如下，里面n就是像素个数(或者说灰度对个数)，$P_{I1,I2}(i,j)$PI1,I2​(i,j) 是256*256的灰度对概率分布，里面的 i 和 j 不是坐标而是灰度值。$g(i,k)$g(i,k) 是二维高斯卷积核，这里需要对二维的概率
分布图做高斯平滑去噪声。
$h_{I_{1}, I_{2}}(i, k)=-\frac{1}{n} \log \left(P_{I_{1}, I_{2}}(i, k) \otimes g(i, k)\right) \otimes g(i, k)$hI1​,I2​​(i,k)=−n1​log(PI1​,I2​​(i,k)⊗g(i,k))⊗g(i,k)
联合分布里的P用下面的方式计算，其中T[]函数判断[]里面等式是否成立，返回1或者0。
$P_{I1,I2}(i,k)=\frac{1}{n}\sum_{P}T[(i,k)=(I_{1p},I_{2p})]$PI1,I2​(i,k)=n1​∑P​T[(i,k)=(I1p​,I2p​)]

注意：
联合熵里的P计算的时候需要初始视差，根据初始视差可以把每个左图位置 $I_b$Ib​ 在对应右图找到唯一匹配位置 $I_m$Im​ 组成一对，而不是把右图的所有像素都匹配成一对。所以如果图是10 * 10的就除以100;
可以从Fig.1看出是利用视差来对右图做扭曲操作使之对应像素和左图位置尽量一致，然后计算联合分布就直接计算每个像素位置就可以，P就是二维的联合分布概率图，以概率图每个位置做高斯平滑，取对数，取负数，再高斯平滑；
可以看出分布图位置对角线的元素较多，是因为左右图对应像素匹配，灰度值基本相同；
实验发现小高斯核效果和大的基本相同，但是效率高，这里选择 7x7；
计算对数的时候，如果是0就替换成一个很小的数。

5.SGM里的熵怎么计算？

上面图说的是怎么计算联合熵，除此之外互信息还和两张图的熵有关，可以直接计算，但是有一个问题就是联合熵没考虑遮挡点，所以这里需要 $H_{I1}$HI1​ 和 $H_{I2}$HI2​ 也不考虑遮挡点，可以像联合熵一样用泰勒展开计算熵。实验发现这么算确实会改善边界效果；
问题是如何找这些遮挡点？**这里的计算方法是利用联合概率分布图(联合分布是通过初试视差扭曲右图在计算的，已经把遮挡点过滤掉了)，这里直接对联合概率分布图的行和列求和就行了(很好理解，行和列分别求和就是两张图的0-255像素出现的概率)。
$H_I=\sum_Ph_I(I_p)$HI​=∑P​hI​(Ip​)
$h_I(i)=-\frac{1}{n}log(P_I(i){\otimes}g(i)){\otimes}g(i)$hI​(i)=−n1​log(PI​(i)⊗g(i))⊗g(i)

6.基于互信息的代价怎么计算？

$M_{I1,I2}=\sum_P mi_{I1.I2}(I_{1p},I_{2p})$MI1,I2​=∑P​miI1.I2​(I1p​,I2p​)
$mi_{I1,I2}(i,k)=h_{I1}(i)+h_{I2}(k)-h_{I1,I2}(i,k)$miI1,I2​(i,k)=hI1​(i)+hI2​(k)−hI1,I2​(i,k)

基于互信息的匹配代价如下,代价取互信息的负是因为互信息越大->匹配效果越好->损失越小：
$C_{MI}(p,d)=-mi_{I_b,f_D(I_m)}(I_{bp},I_{mq}),其中q=e_{bm}(p,d)$CMI​(p,d)=−miIb​,fD​(Im​)​(Ibp​,Imq​),其中q=ebm​(p,d)

7.分层互信息计算视差

这里计算 $mi()$mi() 之前需要知道视差图来扭曲右图计算联合概率分布，但是我们的目的就是求视差图，所以开始就随机生成一个，用这个随机生成的视差图计算一个新视差图，在利用新的计算一个更新的；
迭代计算几次即可，3次就能得到不错的结果，因为像素很多所以初始视差不好的话也能得到好的概率分布 $P$P;
计算的时候采用分层互信息，从小分辨率开始算，一层一层扩大分辨率到最大；
最开始在 $\frac{1}{16}$161​ 分辨率考虑迭代3次，因为最开始是随机生成的视差图；
算法复杂度是 $O(WHD)$O(WHD) , 可以看下面约等式子，HMI计算量相当于BT算法的1…14倍，相差不多。

$1+\frac{1}{2^3}+\frac{1}{4^3}+\frac{1}{8^3}+3\frac{1}{16^3}\approx1.14$1+231​+431​+831​+31631​≈1.14

二.代价聚合

逐像素匹配不是很精确，为了防止噪声干扰，惩罚每个位置的邻域视差变化来加入了平滑性约束。

1.能量函数：

$\begin{aligned} E(D)=& \sum_{\mathbf{p}}\left(C\left(\mathbf{p}, D_{\mathbf{p}}\right)+\sum_{\mathbf{q} \in N_{\mathbf{p}}} P_{1} \mathbf{T}\left[\left|D_{\mathbf{p}}-D_{\mathbf{q}}\right|=1\right]\right.\\ &\left.+\sum_{\mathbf{q} \in N_{\mathbf{p}}} P_{2} \mathrm{T}\left[\left|D_{\mathbf{p}}-D_{\mathbf{q}}\right|>1\right]\right) \end{aligned}$E(D)=​p∑​⎝⎛​C(p,Dp​)+q∈Np​∑​P1​T[∣Dp​−Dq​∣=1]+q∈Np​∑​P2​T[∣Dp​−Dq​∣>1]⎠⎞​​

1.1各个符号的含义：

D：disparity map，代表视差图；
p，q：表示某个像素位置；
Np：p的邻域像素，8连通；
C(p,Dp):视差为Dp时候，该像素位置的代价cost；
P1、P2：惩罚项，分别用于像素p和p的邻域像素q视差之差等于1和大于1的像素；
T[]：这个中括号返回0和1，里面式子成立返回1，否则为0.

1.2各项含义：

第一项是所有像素位置匹配代价和；
第二项对每个像素和邻域视差大小相差等于1的做惩罚，是个较小的常量；
第三项对每个像素和邻域视差大小相差大于1的做惩罚，是个较大的常量。
P1、P2两个惩罚项有什么用？
P1针对较小的视差变化，让结果适应倾斜和弯曲表面，让结果更平滑；
P2针对较大的视差变化，保留不连续性，让结果能适应遮挡或者边界这种不连续的地方从而保存边界。
不连续性常出现在强度变化的地方, 调整$P_2$P2​ 为强度的梯度$P_2=\frac{P_2^{'}}{|I_{bp}-I_{bq}|}$P2​=∣Ibp​−Ibq​∣P2′​​, 需要保证$P1<P2$P1<P2。

2.怎么做代价聚合

有了上面的能量函数，就需要找一个视差图可以使之最小，2D的全局最小化是NP问题，但是1D的最小化可以用动态规划(DP)解决；但是有一些问题，2D图像每行单独用1D动态规划方法容易受条纹影响，因为一个方向的强约束很难和其他方向的弱约束或者没有约束结合在一起。所以有了一个新思路，聚合各个方向的1D匹配代价；
计算方法就是当计算p点的代价的时候沿着它的邻域(一般选8或者16邻域)方向以它为终点来做1D的动态规划，如下面的右图。

上面左图横轴是 $(x，y)$(x，y) 表示整张图所有像素位置，纵坐标是视差搜索范围 $d$d，实线表示从${L_r}$Lr​路径起始点到点 $p$p的路径上经过动态规划聚合后的每个点的视差 $d$d.

3.代价聚合公式

$L_r(p,d)=C(p,d)+min(L_r(p-r,d),L_r(p-r,d-1)+P1,L_r(p-r,d+1)+P1,min_iL_r(p-r,i)+P_2)-min_kL_r(p-r,k)$Lr​(p,d)=C(p,d)+min(Lr​(p−r,d),Lr​(p−r,d−1)+P1,Lr​(p−r,d+1)+P1,mini​Lr​(p−r,i)+P2​)−mink​Lr​(p−r,k)
公式里$L_r$Lr​表示的是当前优化的路径，$L_r(p-r,d)$Lr​(p−r,d)代表的是这条路径上相邻点视差为 $d$d 时候的代价聚合值(其实就是前一个点的)；
公式第一项是该像素原始代价；
第二项是这条路径相邻点视差为d、视差为d-1的聚合值+P1、视差为d+1的聚合值+P1、代价最小的聚合值+P2这四个数的最小值；
第三项为这条路径相邻点的最小视差值，这里减去是因为防止代价聚合值越来越大，因为会一直加。
$S(p,d)=\sum_rL_r（p,d）$S(p,d)=∑r​Lr​（p,d）
最后代价聚合值为所有路径的聚合值之和。

3. 视差计算

最小代价作为视差：
和其他局部方法一样用代价聚合后每个像素选择最小代价的D作为视差(winner take all)；
亚像素估计:
之前的到的结果都是整数精度要求高的场景可能不满足要求，通过相邻视差处的相邻代价来拟合二次曲线重新计算视差；

一致性检查：
右图作为基础图用左图来匹配，计算得到的视差图$D_m$Dm​，和之前计算的视差图 $D_b$Db​可以用来找遮挡区域和错误匹配区域，对应位置差距不满足要求的视差设置为无效；
$D_p=\begin{cases}D_{bp}& if|D_{bp}-D_{mq}|\leq1 \\ D_{inv}& otherwise.\end{cases}$Dp​={Dbp​Dinv​​if∣Dbp​−Dmq​∣≤1otherwise.​

4.多基线匹配

通过计算并组合基础图像和所有匹配图像之间对应像素的匹配代价，该算法可以扩展到多基线匹配。 但是这种方式匹配像素的时候需要考虑遮挡问题(代价聚合之前)是不稳定的。【这里的意思是先组合多个匹配结果后聚合不稳定】
因此，先用基础图像和所有匹配图像分别匹配并代价聚合出结果，再使用【一致性检查】，以消除遮挡时的错误匹配和其他不匹 配。 最后，每对图像都用单独的缩放比例将得到的视差图像融合在一起。

基础图$I_b$Ib​ 和匹配图 $I_{mk}$Imk​ 匹配出的视差为 $D_k$Dk​，如果所有图像之间都互相矫正(Rectify)过(就是所有图都和所有光学平面平行)，每个视差图 $D_{mk}$Dmk​都有一个和基线线性相关的因子 $t_k$tk​, 这些视差图需要用因子归一化。
融合视差图的时候用 $t_k$tk​ 来加权求和，但是因为会有一些离群点，所以多张图一个像素位置只考虑融合和中值差距不超过1的图。

$D_{\mathbf{p}}=\frac{\sum_{k \in V_{\mathrm{p}}} D_{k_{\mathrm{p}}}}{\sum_{k \in V_{\mathrm{p}}} t_{k}}$Dp​=∑k∈Vp​​tk​∑k∈Vp​​Dkp​​​
【论文里加权求和是写的上面的公式，好像不太对，应该把求和符号提出来吧1/2+2/3!=(1+2)/(2+3)】

$V_p= \{k||\frac{D_{tk}}{t_k}-med_i\frac{D_{ip}}{t_i}|\leq\frac{1}{t_k}\}$Vp​={k∣∣tk​Dtk​​−medi​ti​Dip​​∣≤tk​1​}
通过使用中值和加权的方式会鲁棒，如果图像足够多的话还可以进一步筛选每个位置的有效像素得到 $V_p$Vp​ 来提高鲁棒性，不符合规则的都算作无效。而且如果算法是分层计算的互信息，这里也是分层融合作为初始视差传给下一层，复杂度是 $O(KWHD)$O(KWHD).

5.视差优化

1)Removal of Peaks,去除峰值
2)Intensity Consistent Disparity Selection,强度一致性视差选择
3)Discontinuity Preserving Interpolation,不连续处插值

6.处理超大图像

7.融合多个视差图

参考：
1
2
3
4