论文地址:Stereo Processing by Semiglobal Matching and Mutual Information
问题提出
局部算法基于局部窗口内视差一致的假设太强,不符合实际情况;而全局算法在全局能量函数的约束下计算量大,算法效率不高。为此提出一个半全局匹配算法,利用图像中非遮挡点进行立体匹配,在全局算法的构架下采用一维路径聚合(动态规划)方法来求解全局能量函数最小化问题,用一维最优来近似二维最优解,提高了算法计算效率。
算法框架
1、 基于互信息的代价计算
互信息是一种对图像光照变化比较鲁棒的相关性测度,通过两张图像各自的熵和联合熵来定义(有关于互信息和熵的已经在另外一篇文章有介绍—“信息论(熵、条件熵、联合熵、互信息)”)。在图像中,熵是基于灰度概率分布得到的统计量,图像的熵越大说明图像的纹理越丰富。两幅经过配准的图像的联合熵很小,因为一幅图像可以通过另外一幅图像进行预测,所以两幅校准图像的互信息很大。
计算两幅图的互信息需要将两幅图像的视差图,在文章中使用了分层迭代的方法,对图像进行降采样得到图像金字塔,为最上层的图像对随机生成一张视差图,然后计算得到的代价体作为初始值提供给下一层作为初始值,依次计算每一层的代价体。
基于互信息计算代价体的方法较为复杂,在实际应用中,也被证明效果并不佳。一般采用Census变换来替代互信息的方式计算代价体。
2、基于Census变换的代价计算
Census算法是基于窗口计算局部特征的方法,分别计算左右窗口中邻域像素和中心像素的大小关系,展开为1维后使用汉明距离计算两个窗口之间的距离。窗口特征计算如式1所示:
C
s
(
x
,
y
)
=
⨂
i
=
−
n
′
n
′
⨂
j
=
−
m
′
m
′
ξ
(
I
(
x
,
y
)
,
I
(
x
+
i
,
y
+
j
)
)
(1)
C_{s}(x, y)=\bigotimes_{i=-n^{\prime}}^{n^{\prime}} \bigotimes_{j=-m^{\prime}}^{m^{\prime}} \xi(I(x, y), I(x+i, y+j))\tag1
Cs(x,y)=i=−n′⨂n′j=−m′⨂m′ξ(I(x,y),I(x+i,y+j))(1)
上式表示计算
(
x
,
y
)
(x,y)
(x,y)处,窗口大小为
n
′
∗
m
′
n^{\prime}*m^{\prime}
n′∗m′的特征,其中,
ξ
(
x
,
y
)
\xi(x, y)
ξ(x,y)如式2所示:
ξ
(
x
,
y
)
=
{
0
if
x
≤
y
1
if
x
>
y
(2)
\xi(x, y)=\left\{\begin{array}{ll} 0 & \text { if } x \leq y \\ 1 & \text { if } x>y \end{array}\right.\tag2
ξ(x,y)={01 if x≤y if x>y(2)
census算法基于汉明距离计算点
(
x
,
y
)
(x,y)
(x,y)和视差为d的右视图点
(
x
−
d
,
y
)
(x-d, y)
(x−d,y)的匹配代价如式3所示:
C
(
x
,
y
,
d
)
=
Hamming
(
C
s
l
(
x
,
y
)
,
C
s
r
(
x
−
d
,
y
)
)
(3)
\mathrm{C}(x, y, d)=\operatorname{Hamming}\left(C_{s l}(x, y), C_{s r}(x-d, y)\right)\tag3
C(x,y,d)=Hamming(Csl(x,y),Csr(x−d,y))(3)
Census算法基于窗口内的相对亮度差计算匹配代价,对左右视图光照不一致和噪声具有较强的鲁棒性,对弱纹理区域也具有较好的效果,但Census算法对于重复的纹理区域会产生歧义,而AD算法是基于单像素邻域亮度来计算匹配代价,可以在一定程度上缓解重复纹理的歧义性问题。Census算法是基于局部窗口的运算,可以进行多线程并行计算。经过代价计算,可以得到一个初始的代价体,经过代价聚合才能得到得到更准确的代价体,聚合后的代价体可以通过WTA与亚像素插值得到视差图。代价体如下图所示:
3、多路径代价聚合*
基于Census变换计算得到初始的代价体,Census变换计算代价时只考虑了局部的相关性,计算出的结果对噪声很敏感,无法用于计算最优视差。为此,SGM基于全局算法的架构下,设计了全局能量函数,通过一维多路径代价聚合求取全局能量函数的近似最优解。全局能量函数定义如式4所示:
E
(
d
)
=
E
data
(
d
)
+
E
smooth
(
d
)
(4)
E(d)=E_{\text {data }}(d)+E_{\text {smooth }}(d)\tag4
E(d)=Edata (d)+Esmooth (d)(4)
其中,
E
d
a
t
a
E_{data}
Edata为数据项,是反应视差图对应的总体匹配代价的测度,
E
s
m
o
o
t
h
E_{smooth}
Esmooth为平滑项,是为了让视差图满足某些条件假设的约束(场景表面连续性约束),平滑项会对相邻像素视差变化过大的像素点进行惩罚(加大能量值)。全局能量最小化问题是个NP完全问题,常用的能量最优策略有图割法、置信度传播法、合作优化法等,SGM中使用的是一维动态规划方法,用一维路径聚合的方式来近似二维最优,效率更高。SGM中的能量函数具体表达式如式5所示:
E
(
D
)
=
∑
p
C
(
p
,
D
p
)
+
∑
q
∈
N
p
P
1
T
[
∣
D
p
−
D
q
∣
=
1
]
+
∑
q
∈
N
p
P
2
T
[
∣
D
p
−
D
q
∣
>
1
]
]
(5)
E(\mathrm{D})=\sum_{\mathrm{p}} \mathrm{C}\left(\mathrm{p}, \mathrm{D}_{\mathrm{p}}\right)+\sum_{\mathrm{q} \in N_{\mathrm{p}}} P_{1} T[|\mathrm{D}_{\mathrm{p}}-\mathrm{D}_{\mathrm{q}}|=1]+\sum_{\mathrm{q} \in N_{\mathrm{p}}} P_{2} T\left[\left|\mathrm{D}_{\mathrm{p}}-\mathrm{D}_{\mathrm{q}}\right|>1\right]]\tag5
E(D)=p∑C(p,Dp)+q∈Np∑P1T[∣Dp−Dq∣=1]+q∈Np∑P2T[∣Dp−Dq∣>1]](5)
其中,
∑
p
C
(
p
,
D
p
)
\sum_{\mathrm{p}} \mathrm{C}\left(\mathrm{p}, \mathrm{D}_{\mathrm{p}}\right)
∑pC(p,Dp)对应式 4 中的
E
d
a
t
a
E_{data}
Edata,其余两项对应式 4 中的
E
s
m
o
o
t
h
E_{smooth}
Esmooth,
P
1
P_{1}
P1是对p点邻域内像素最优视差和p点最优视差相差为为1的点进行的小惩罚(如q为p邻域内的点,p点在视差为
d
p
d_{p}
dp时代价最小,即p点的最优视差为
d
p
d_{p}
dp,q点的最优视差为
d
q
d_{q}
dq,若
d
p
−
d
p
=
1
d_{p}-d_{p}=1
dp−dp=1,则对p点进行
P
1
P_{1}
P1惩罚),
P
2
P_{2}
P2是对邻域内像素与p点视差超过1的点进行
P
2
P_{2}
P2 大惩罚,一般而言,
P
2
>
P
1
P_{2}>P_{1}
P2>P1,
P
2
P_{2}
P2一般是按式 6 方式来自适应更新:
P
2
=
P
2
′
∣
I
b
p
−
I
b
q
∣
,
P
2
>
P
1
(6)
P_{2}=\frac{P_{2}^{\prime}}{\left|I_{b \mathrm{p}}-I_{b \mathrm{q}}\right|}, \quad P_{2}>P_{1}\tag6
P2=∣Ibp−Ibq∣P2′,P2>P1(6)
其中
P
2
′
P_{2}^{\prime}
P2′ 为
P
2
P_{2}
P2的初始值,一般设置为远大于
P
1
P_{1}
P1 的数值,结合
P
2
P_{2}
P2表达式,
P
2
T
[
∣
D
p
−
D
q
∣
>
1
]
]
P_{2} T\left[\left|\mathrm{D}_{\mathrm{p}}-\mathrm{D}_{\mathrm{q}}\right|>1\right]]
P2T[∣Dp−Dq∣>1]]表示的含义为:如果这两个像素的颜色值很相近,理论两点应该处于同一个视差平面,但是两个像素的视差却相差大于1,很反常,因此给予一个很大的惩罚。实际上,式5表示的能量函数还是NP完全问题,为了高效计算聚合代价,SGM使用动态规划的方式来近似求解该问题,如式7所示:
L
r
(
p
,
d
)
=
C
(
p
,
d
)
+
min
{
L
r
(
p
−
r
,
d
)
L
r
(
p
−
r
,
d
−
1
)
+
P
1
L
r
(
p
−
r
,
d
+
1
)
+
P
1
min
i
L
r
(
p
−
r
,
i
)
+
P
2
}
−
min
i
L
r
(
p
−
r
,
i
)
(7)
L_{r}(\mathrm{p}, d)=\mathrm{C}(\mathrm{p}, d)+\min \left\{\begin{array}{l} L_{r}(\mathrm{p}-\mathrm{r}, d) \\ L_{r}(\mathrm{p}-\mathrm{r}, d-1)+P_{1} \\ L_{r}(\mathrm{p}-\mathrm{r}, d+1)+P_{1} \\ \underset{i}{\min} L_{r}(\mathrm{p}-\mathrm{r}, i)+P_{2} \end{array}\right\}-\min _{i} L_{r}(\mathrm{p}-\mathrm{r}, i)\tag7
Lr(p,d)=C(p,d)+min⎩
⎨
⎧Lr(p−r,d)Lr(p−r,d−1)+P1Lr(p−r,d+1)+P1iminLr(p−r,i)+P2⎭
⎬
⎫−iminLr(p−r,i)(7)
其中r表示路径聚合方向上单位步长,
p
−
r
p-r
p−r 表示路径上前一个像素。第一项为匹配代价,属于数据项;第二项是平滑项,
L
r
(
p
−
r
,
d
)
L_{r}(\mathrm{p}-\mathrm{r}, d)
Lr(p−r,d) :表示如果路径上前一个像素的最优视差和当前像素的最优视差相同,则不加惩罚;
L
r
(
p
−
r
,
d
−
1
)
+
P
1
和
L
r
(
p
−
r
,
d
+
1
)
+
P
1
L_{r}(\mathrm{p}-\mathrm{r}, d-1)+P_{1}和L_{r}(\mathrm{p}-\mathrm{r}, d+1)+P_{1}
Lr(p−r,d−1)+P1和Lr(p−r,d+1)+P1 :表示如果路径上前一个像素的最优视差和当前像素的最优视差相差1,则做
P
1
P_{1}
P1惩罚、
min
i
L
r
(
p
−
r
,
i
)
+
P
2
\underset{i}{\min} L_{r}(\mathrm{p}-\mathrm{r}, i)+P_{2}
iminLr(p−r,i)+P2:表示表示如果路径上前一个像素的最优视差和当前像素的最优视差相差大于1,则做
P
2
P_{2}
P2惩罚,取三者中的最小值;第三项是为了保证新的路径代价不超过数值上限,即:
L
≤
C
max
+
P
2
(8)
L \leq C_{\max }+P_{2}\tag8
L≤Cmax+P2(8)
所有路径聚合代价计算的公式如式9所示:
S
(
p
,
d
)
=
∑
r
L
r
(
p
,
d
)
(9)
\mathrm{S}(\mathrm{p}, d)=\sum_{r} L_{r}(\mathrm{p}, d)\tag9
S(p,d)=r∑Lr(p,d)(9)
经过多路径代价聚合后的代价体(h,w,d)能够更好地表示左右图像的立体匹配代价,SGM将全图(二维)的能量最小化问题转换成了N个方向(一维)的子问题的最优解问题,求得的是全局能量函数的近似解而非真实最优解,提升了算法效率。
4、视差计算与视差优化
(1)视差值计算与子像素拟合
在SGM算法中采用WTA算法从聚合代价体中计算每个像素最小的代价对应的最优视差值。得到最优视差后一般通过二次曲线优化的方式来得到更精确的最优视差值的。
如在代价体中,沿着d方向将p点的代价体抽取出来,p点的不同视差值的代价值数组为
(
c
m
i
n
.
.
.
c
d
−
1
、
c
d
、
c
d
+
1
.
.
.
c
d
m
a
x
)
(c_{min}...c_{d-1}、c_{d}、c_{d+1}...c_{d_{max}})
(cmin...cd−1、cd、cd+1...cdmax), d为整数,但最优的视差
d
∗
d^{*}
d∗可能是一个小数,因此用d点和d点前后两点来的代价值拟合二次曲线,求得二次曲线的最优值
d
∗
d^{*}
d∗做为p点的最优视差。
(2)左右一致性检查
因为遮挡和重复区域出现,会导致算法中会出现一些错误点。左右一致性检查具体步骤为将右视图作为参考视图,计算视差图
I
d
i
s
p
−
m
I_{disp-m}
Idisp−m,使用左视图计算得到的视差图
I
d
i
s
p
−
b
I_{disp-b}
Idisp−b计算左视图中p点的视差值
b
p
b_{p}
bp与在右视图中的对应点
m
q
m_{q}
mq的位置,利用右视图计算得到的视差图
I
d
i
s
p
−
m
I_{disp-m}
Idisp−m计算
D
m
q
D_{mq}
Dmq,若两者视差差值超过阈值,则将p点设置为错误点,如式10所示
D
p
=
{
D
b
p
if
∣
D
b
p
−
D
m
q
∣
≤
T
D
invalid
otherwise
(10)
D_{\mathrm{p}}=\left\{\begin{array}{l} D_{b \mathrm{p}} \text { if }\left|D_{b \mathrm{p}}-D_{m \mathrm{q}}\right| \leq T \\ D_{\text {invalid }} \text { otherwise } \end{array}\right. \tag{10}
Dp={Dbp if ∣Dbp−Dmq∣≤TDinvalid otherwise (10)
在计算右视图的视差图时,可以通过左视图代价体与右视图代价体的关系直接得到右视图的代价体,再通过子像素拟合与WTA得到右视图的视差图,而不用再从头开始计算。通过左视图计算右视图代价体的公式如式11所示:
c
o
s
t
R
(
x
,
y
,
d
)
=
c
o
s
t
L
(
x
+
d
,
y
,
d
)
(11)
cost_{R}(x,y,d) = cost_{L}(x+d,y,d)\tag{11}
costR(x,y,d)=costL(x+d,y,d)(11)
注:不管拿哪个视图当做参考视图,假设d的定义就是
d
=
x
l
−
x
r
d= x_{l}-x_{r}
d=xl−xr
例如:由
d
=
x
l
−
x
r
d= x_{l}-x_{r}
d=xl−xr知,右视图
(
x
,
y
)
(x,y)
(x,y)点对应左视图(x+d,y)点。
c
o
s
t
R
(
x
,
y
,
d
)
cost_{R}(x,y,d)
costR(x,y,d)表示右视图中
(
x
,
y
)
(x, y)
(x,y)点与左视图
(
x
+
d
,
y
)
(x+d, y)
(x+d,y)点的匹配代价。而
c
o
s
t
L
(
x
+
d
,
y
,
d
)
cost_{L}(x+d,y,d)
costL(x+d,y,d)表示的是左视图
(
x
+
d
,
y
)
(x+d, y)
(x+d,y)与右视图
(
x
,
y
)
(x, y)
(x,y)点的匹配代价,两个点都是一样的,因此匹配代价也是相同的。
(3)剔除小连通域
剔除视差图中连通域极小的区域,同一个连通域内的像素视差之差小于预设的阈值。
(4)唯一性检测
每个像素的最优视差和次优视差之差小于阈值,则直接剔除。
(5)抑制噪声
利用中值滤波和双边滤波来保持视差图边缘。
小结
SGM算法是基于平行窗口来计算匹配代价,然后多路径聚合来进行代价聚合,在此基础上通过WTA来计算得到视差图,最后通过左右一致性检查等算法来优化视差图。SGM算法适合cuda加速,但对倾斜面的重建效果欠优。