SGM(Semi-Global Matching)算法笔记

SGM

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

【论文重点就是基于互信息的代价，代价聚合方式这两部分，后面作者也做了很多视差优化、后处理等很多工作，但是后面部分我看的不是特别仔细，可能表达或者理解有点问题，可以参考下】

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

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

一.逐像素计算代价

1.SGM怎么计算代价？

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

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

1)互信息：

$M{I}_{I1,I2}=H_{I1}+H_{I2}-H_{I1,I2}$
是变量间相关性，可以看成一个随机变量中包含的关于另一个随机变量的信息量，所以互信息越大代表变量之间重复的信息更多，就是匹配的更好；

2)熵：

$H_I=-{\int }_0^1{P}_I(i)\log P_I(i)di$
是随机变量的不确定性，代表信息量，熵越大信息量越大；

3)联合熵：

$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$
是变量之间的不确定性，代表多个随机变量组成的变量系统的信息量，联合熵越大代表这个组合信息量越大；

匹配的结果好，联合熵会比较低，因为一张图可以用另一张预测出来(知道左图每个像素对应右图的位置，直接扭曲过去就行，匹配的好基本位置都能扭曲对)，这样信息量就比较低。

3.概率P是什么意思？

灰度直方图是关于灰度级分布的函数，是对图像中灰度级分布的统计；

1)熵里的P:

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

2)联合熵的P：

联合熵里的 $P_{I1,I2}$ 相当于变成了二维的图像对，范围是(0,0)-(255,255), 概率分布是256*256的，计算的时候一对一对计算像素对出现的频数再除以像素对数。

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

1) $H_{I1,I2}$计算方式：

利用泰勒展开转化成像素和的形式近似联合熵，就是求所有$h_{I1,I2}$的和，自变量$I_{1p}$和$I_{2p}是灰度对$。

$H_{I1,I2}={\sum }_P{h}_{I1,I2}(I_{1p},I_{2p})$

2)$h_{I1,I2}$ 计算方式如下公式：

$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)$

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

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

5.SGM里的熵怎么计算？

1)泰勒展开计算熵:

上面图说的是怎么计算联合熵，除此之外互信息还和两张图的熵有关，可以直接计算按照公式计算，但是有一个问题就是联合熵没考虑遮挡点，所以这里需要 $H_{I1}$ 和 $H_{I2}$ 也不考虑遮挡点，方法是联合熵一样用泰勒展开计算熵。实验发现这么算确实会改善边界效果。

$H_I={\sum }_P{h}_I(I_p)$

2)如何找这些遮挡点？

利用联合概率分布图(联合分布是通过初始视差扭曲右图再计算的，已经把遮挡点过滤掉了，因为只考虑了遮挡点视差)，这里直接对联合概率分布图的行和列求和就行了(很好理解，行和列分别求和对应两张图的0-255像素出现的概率)。

$h_I(i)=-\frac{1}{n}log(P_I(i)\otimes g(i))\otimes g(i)$

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

$M_{I1,I2}={\sum }_{P}m{i}_{I1.I2}(I_{1p},I_{2p})$
$m{i}_{I1,I2}(i,k)=h_{I1}(i)+h_{I2}(k)-h_{I1,I2}(i,k)$

基于互信息的匹配代价如下,代价取互信息的负是因为【互信息越大->匹配效果越好->损失越小】：
$C_{MI}(p,d)=-m{i}_{I_b,f_D(I_m)}(I_{bp},I_{mq}),其中q=e_{bm}(p,d)$

7.分层互信息计算视差

1) 计算视差步骤：

随机生成初始视差，利用初始视差扭曲右图计算联合概率分布，计算 $mi()$，得到新视差图，然后利用新视差图迭代上述步骤；

2) 迭代次数和计算量：

迭代计算几次即可，3次就能得到不错的结果，因为像素很多所以初始视差不好的话也能得到好的概率分布 $P$;
计算的时候采用分层互信息，从小分辨率开始算，一层一层扩大分辨率到最大；
最开始在 $\frac{1}{16}$ 分辨率考虑迭代3次，因为随机视差图不准确多迭代几次；
算法复杂度是 $O(WHD)$ , 可以看下面式子，HMI计算量相当于BT算法的1…14倍，相差不多。

$1+\frac{1}{2^3}+\frac{1}{4^3}+\frac{1}{8^3}+3\frac{1}{16^3}\approx 1.14$

二.代价聚合

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

1.能量函数：

$\begin{array}{rl}E(D)=& \sum _{\mathbf{p}}(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]\\ & +\sum _{\mathbf{q}\in N_{\mathbf{p}}}{P}_2\mathrm{T}\left[\left|D_{\mathbf{p}}-D_{\mathbf{q}}\right|>1\right])\end{array}$

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.

2)各项含义：

第一项是所有像素位置匹配代价和；
第二项对每个像素和邻域视差大小相差等于1的做惩罚，是个较小的常量；
第三项对每个像素和邻域视差大小相差大于1的做惩罚，是个较大的常量。

3)P1、P2两个惩罚项有什么用？

P1针对较小的视差变化，让结果适应倾斜和弯曲表面，让结果更平滑；
P2针对较大的视差变化，保留不连续性，让结果能适应遮挡或者边界这种不连续的地方从而保存边界。
不连续性常出现在灰度变化的地方, 调整$P_2$ 为灰度的梯度$P_2=\frac{P_2^{{}^{\prime }}}{|I_{bp}-I_{bq}|}$, 需要保证$P1<P2$。

2.怎么做代价聚合

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

上面左图横轴是 $(x，y)$ 表示整张图所有像素位置，纵坐标是视差搜索范围 $d$，实线表示从$L_r$路径起始点到点 $p$的路径上经过动态规划聚合后的每个点的视差 $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,mi{n}_i{L}_r(p-r,i)+P_2)-mi{n}_k{L}_r(p-r,k)$

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

三. 视差计算

1.最小代价作为视差：

和其他局部方法一样用代价聚合后每个像素选择最小代价的D作为视差(winner take all)；

2.亚像素估计:

之前的到的结果都是整数精度要求高的场景可能不满足要求，通过相邻视差处的相邻代价来拟合二次曲线重新计算视差；

3.一致性检查：

右图作为基础图用左图来匹配，计算得到的视差图$D_m$，和之前计算的视差图 $D_b$可以用来找遮挡区域和错误匹配区域，对应位置差距不满足要求的视差设置为无效；

$D_p=\{\begin{array}{ll}{D}_{bp}& if|D_{bp}-D_{mq}|\le 1\\ D_{inv}& otherwise.\end{array}$

四.多基线匹配

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

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

$D_{\mathbf{p}}=\frac{{\sum }_{k\in V_{\mathrm{p}}}{D}_{k_{\mathrm{p}}}}{{\sum }_{k\in V_{\mathrm{p}}}{t}_k}$

【论文里加权求和是写的上面的公式，好像不太对，应该把求和符号提出来吧1/2+2/3!=(1+2)/(2+3)】

V p = { k ∣ ∣ D t k t k − m e d i D i p t i ∣ ≤ 1 t k } 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$ 来提高鲁棒性，不符合规则的都算作无效。而且如果算法是分层计算的互信息，这里也是分层融合作为初始视差传给下一层，复杂度是 $O(KWHD)$.

五.视差优化

1. Removal of Peaks,去除峰值

用联通域分割的方法，以视差相差不到1个像素的4联通邻域划分，把整个视差图分成若干块，像素总数少于预先设置的阈值的联通域就设置为无效；

2. Intensity Consistent Disparity Selection,强度一致性视差选择

2.1 要解决什么问题？

解决视差不连续放错地方导致边界错误的问题，主要是找无纹理区域的边界视差。

室内环境前景经常在低纹理或者无纹理背景前。能量函数对视差不连续的位置没有偏好，可能把视差的不连续放在错误的地方。
代价聚合步骤中根据灰度梯度调整惩罚项 $P_2$ 可以正确的把视差变化放到【前景目标旁边】，因为灰度梯度在这些位置比无纹理区域更吻合。

2.2 SGM可以处理好、不能处理好的情况：

可以处理好路径上的无纹理区域两端有纹理的情况；
不能处理好路径上遇到前或背景区域(一边有一边无)，或者一直无纹理的情况。

SGM的能量函数沿着1D路径单独计算，只有在无纹理区域的两侧都是纹理区域时，才能判断是视差不连续。
无纹理区域有不同的形状和大小并且会超出图像边界(室内墙壁比较常见)。
根据1D路径的位置和方向可能在无纹理部分周围遇到带纹理的前后景目标，这种情况可以确定正确的视差不连续。
还可能在路径上遇到前景或者背景纹理，或者在无纹理的区域离开图像，这种情况无法确定视差的不连续性。
合并所有不连续的路径可能会让无纹理背景前的前景对象周围具有模糊的不连续性。

2.3 针对不能处理好的情况,对无纹理区域做出一些假设:

a.无纹理的区域内不会出现视差不连续性；
b.无纹理化区域里也可以看到一些纹理；
c.无纹理区域的表面可用一个平面近似表示。

第一个假设通常都是正确的，因为深度不连续会导致一些亮度变化；
第二个假设是必需的，因为完全没纹理的背景表面不能匹配出视差；
第三个假设最弱，视差恒定就作为一个平面，不然就划分割成两个平面。

2.4 处理流程

为了区分前后景边界可以根据三个假设找到后景。
a.均值漂移分割原图：
分割原图得到多个无纹理区域$S_i$，因为前后景外观差异导致灰度梯度不同容易分开。

识别无纹理区域，在 $I_b$ 上用固定带宽的(Mean Shift)均值漂移分割算法(根据假设1)。
辐射带宽 ${\sigma }_r$设置为 $P_1$，经常设置为4。把低于平滑度惩罚的灰度变化作为噪声。为了提高速度，空间带宽 ${\sigma }_s$ 设置为一个很低的值。
把所有小于特定阈值(100像素)的分割忽略掉，因为SGM可以很好的处理小的无纹理区域。

b.联通域分割视差：
根据视差变化对$S_i$ 做联通域分割为多个 $S_{ik}$，因为无纹理区域也有一些纹理可以再分。

需要解决的问题是无纹理区域内模糊的不连续性。因此，只要背景表面包含一点纹理就能预测无纹理区域包含的不正确的前景目标视差和正确的背景视差(根据假设2)。所以在每个分割片段 $S_i$ 里的视差应该是正确的(不正确就不能很好利用无纹理区这点仅有的纹理)。可以通过继续分割每个$S_i$的视差图来判断正确的 $S_i$视差的几种可能性。
可以通过允许相邻视差差距在1个像素值以内的联通域分割算法发完成，这样把每个分割 $S_i$ 分为几个小的 $S_{ik}$ 联通域。

c.得到拟合平面：
找最佳拟合平面。

通过 $S_{ik}$ 的视差计算可能的拟合平面$F_{ik}$ (基于假设3选择平面)。小于等于12个像素的分割被忽略掉，因为这么小的一般不对。评估每一个可能的平面，用可能的拟合平面替换掉 $S_i$ 所有像素并计算 $S_i$ 所有未遮挡像素的 $E_{ik}$。在 $S_i$ 内，如果另一个像素和像素q匹配有更高的视差值，则像素p被遮挡。
首先把p映射用$q=e_{bm}(p,D_p^{{}^{\prime }})$ 映射到匹配图里, 然后找q在基础图的对极线 $e_{bm}(q,d)$，$d>D_p^{{}^{\prime }}$ 。如果对极线通过了一个视差大于d的像素则认为像素p是遮挡的。

d.替换视差：
用最佳拟合平面的视差替换无纹理区域视差，因为假设3提出无纹理区域是1个平面，需要视差恒定，这里为了把错误的边界替换成恒定视差(从而保证视差的连续性)。

$$F_i=F_{ik}^{{}^{\prime }}with{k}^{{}^{\prime }}=argmi{n}_k(E_{ik})$$
$$D_p^{{}^{\prime }}=\{\begin{array}{ll}{F}_i(p)& ifp\in S_i\\ D_P& otherwise.\end{array}$$对每一个固定的分割 $S_i$ 选择代价$E_{ik}$最小的拟合平面$F_{ik}$。所有$S_i$内的视差都被所选表面的值代替，以使视差选择和基础图像的灰度一致。

注：
本论文方法使用图像分割和平面拟合优化初始视差；相对其他方法，SGM已经得到比较精确的结果，【只修改了特定大小(比较大的)的无纹理区域的值】。因此，只处理关键区域，不会破坏可能匹配好的区域。另一个区别是，通过考虑初始视差图像中固有的少量假设来选择我们所考虑的区域的视差。没有增加复杂度。计算最佳拟合平面需要访问所有分割的像素。测试所有假设都需要访问所有假设的所有段的所有像素（即最大N）。遮挡测试要求每个像素最多经历D个视差。因此，复杂度的上限是O（WHDN）。
分割后的像素通常仅占整个图像的一小部分，并且一个段的假设N的最大数量通常很小，通常仅为1

3. Discontinuity Preserving Interpolation,不连续处插值

一致性检查、融合多张视差图、峰值过滤这几步会导致一些视差无效，视差图出现空洞，需要用插值的方式获得稠密视差图；

3.1 无效视差包括:

遮挡: 遮挡不能用遮挡物插值(会太光滑或者不连续)，要用背景进行插值；
误匹配: 可以直接用相邻像素直接插值；

因为要插值遮挡区域，所以和跟遮挡区域相邻的误匹配区域也用背景插值；

3.2 怎么判断是遮挡还是误匹配

遮挡和误匹配可以当成连续性检查的一部分来区分，图中穿过遮挡像素 $P1$ 的对极线跟 $D_m$ 的视差函数不相交；误匹配造成的不连续点$P_2$则相交；所以【可以根据无效像素对应极线和$D_m$交点判断是遮挡还是误匹配】；

从有效区域逐步传播到无效区域的方式来插值，类似SGM的8个方向，需要保存8个 $v_{pi}$ 值，最后计算的视差图如下:

$$D_p^{{}^{\prime }}=\{\begin{array}{ll}seclo{w}_i{v}_{pi}& ifpisoccluded,\\ me{d}_i{v}_{pi}& ifpismismatched,\\ D_p& otherwise.\end{array}$$
上面式子第一项是选择8个值第二小的视差来确保从较低深度的背景插值；
第二项是选择中值【不选均值而选择中值的好处是当误匹配处是遮挡处边界时可以保留视差的不连续性】；
最后用中值滤波去除剩下的不规则区域并且平滑视差图。

六.处理超大图像

带重叠的切分大图分别处理，重叠越大越好(图太大就切分完再切分)。

SGM的方法需要很多临时的内存来保存代价矩阵 $C[]$，代价聚合后的矩阵 $S[]$，融合前的视差图，$1000\ast 1000$ 的图像在视差搜索范围100的时候内存可能就不够了；
解决方法是划分成若干块，多基线匹配的视差图融合之前先合并成一张大视差图；
带重叠的划分图像块【因为SGM只使用一侧像素来靠近边界，边界不重叠会导致精度低甚至出现误匹配，尤其是位于图像中间的边界和边界的低纹理区域】；
重合的边界加权平均，分块尽可能大；
如果是特别大的类似航空图像，需要先划分成大块单独计算，每大块划分成小块进行上述操作，可以多机同时进行。

七.融合多个视差图

多相机的图映射到一个坐标系，坐标系划分网络，每个网络里找中位数融合投影，再通过插值得到最终结果。

多个相机拍摄的图像相匹配重建场景，在公共平面上进行2.5D正投影，选择与所有相机光学中心平行的公共平面，建立坐标系$R_o,T_o$，z轴与平面垂直，x和y轴在平面上划分为等间距单元；
通过重建所有像素，用$R_o,T_o$转换并把z值保存在所属的像素位置把每章视差图分别转换为正交投影【这里说建立坐标系 $R_o,T_o$ 感觉写成找到多个旋转平移参数转换到一个公共坐标系更合适】；
正投影变化会导致出现遮挡，每个摄像机单独映射完之后选择落入每个单元的中位数融合正投影；
融合之前不在单个视差图中进行插值来补充缺失的视差有利于消除异常值，因为有可能在其他视差图里填补缺失会比插值更准；
正投影会导致新的缺失，所以插值其实不能避免，正投影之后用于路径上插值的关于遮挡和误匹配的信息会丢失；
因此需要用其他方法插值2.5D的投影数据，对高度数据分段(每段里的数据可以在一定范围浮动)，每个高度代表一个物理表面；
分割的段内和段间都可能缺失，段内的用相邻有效像素进行反距离插值(InverseDistanceWeighted)，段间的缺失根据缺失位置周围的某些像素来计算(所有靠近缺失处的分割段，它们靠近缺失处的边界中有效的像素计算均值，均值低的用来计算插值)；
这么做插值会更平滑，但是又通过推断背景保持了深度的不连续性，反距离插值复杂度比路径插值要高但是只在融合结果上计算可以接受。

参考：
1
2
3
4