2D-2D:对极几何(epipolar geometry)
对极约束
现在,假设我们从两张图像中,得到了一对配对好的特征点,像图7-7里显示的那样。如果我们有若干对这样的匹配点,就可以通过这些二维图像点的对应关系,恢复出在两帧之间摄像机的运动。
连线
O
1
p
1
‾
\overline{O_{1} p_{1}}
O1p1 和连线
O
2
p
2
‾
\overline{O_{2} p_{2}}
O2p2 在三维空间中会相交于点
P
P
P。
O
1
O
2
P
O_1O_2P
O1O2P称为极平面(Epipolar plane)。
e
1
,
e
2
,
e_{1}, e_{2},
e1,e2, 称为极点(Epipoles),
O
1
O
2
O_{1} O_{2}
O1O2 被称为基线(Baseline)。极平面与两个像平面
I
1
,
I
2
I_{1}, I_{2}
I1,I2 之间的相交线
l
1
,
l
2
l_{1}, l_{2}
l1,l2 为极线(Epipolar line)。
现在,我们从代数角度来看一下这里出现的几何关系。在第一监的坐标系下,设 P 的空间位置为:
P
=
[
X
,
Y
,
Z
]
T
\boldsymbol{P}=[X, Y, Z]^{T}
P=[X,Y,Z]T
根据针孔相机模型,我们知道两个像素点
p
1
,
p
2
p_{1}, p_{2}
p1,p2 的像素位置为:
s
1
p
1
=
K
P
,
s
2
p
2
=
K
(
R
P
+
t
)
s_{1} \boldsymbol{p}_{1}=\boldsymbol{K} \boldsymbol{P}, \quad s_{2} \boldsymbol{p}_{2}=\boldsymbol{K}(\boldsymbol{R} \boldsymbol{P}+\boldsymbol{t})
s1p1=KP,s2p2=K(RP+t)
这里
K
K
K 为相机内参矩阵,
R
,
t
R, t
R,t 为两个坐标系的相机运动。如果使用齐次坐标,我们也可以把上式写成在乘以非零常数下成立的(up to a scale) 等式
p
1
=
K
P
,
p
2
=
K
(
R
P
+
t
)
\boldsymbol{p}_{1}=\boldsymbol{K} \boldsymbol{P}, \quad \boldsymbol{p}_{2}=\boldsymbol{K}(\boldsymbol{R} \boldsymbol{P}+\boldsymbol{t})
p1=KP,p2=K(RP+t)
现在,取:
x
1
=
K
−
1
p
1
,
x
2
=
K
−
1
p
2
\boldsymbol{x}_{1}=\boldsymbol{K}^{-1} \boldsymbol{p}_{1}, \quad \boldsymbol{x}_{2}=\boldsymbol{K}^{-1} \boldsymbol{p}_{2}
x1=K−1p1,x2=K−1p2
这里的
x
1
,
x
2
\boldsymbol{x}_{1}, \boldsymbol{x}_{2}
x1,x2 是两个像素点的归一化平面上的坐标。代入上式,得:
x
2
=
R
x
1
+
t
\boldsymbol{x}_{2}=\boldsymbol{R} \boldsymbol{x}_{1}+\boldsymbol{t}
x2=Rx1+t
两边同时左乘
t
∧
∘
\boldsymbol{t}^{\wedge} \circ
t∧∘ 回忆
∧
^{\wedge}
∧ 的定义,这相当于两侧同时与
t
\boldsymbol{t}
t 做外积:
t
∧
x
2
=
t
∧
R
x
1
\boldsymbol{t}^{\wedge} \boldsymbol{x}_{2}=\boldsymbol{t}^{\wedge} \boldsymbol{R} \boldsymbol{x}_{1}
t∧x2=t∧Rx1
然后,两侧同时左乘
x
2
T
\boldsymbol{x}_{2}^{T}
x2T
x
2
T
t
∧
x
2
=
x
2
T
t
∧
R
x
1
\boldsymbol{x}_{2}^{T} \boldsymbol{t}^{\wedge} \boldsymbol{x}_{2}=\boldsymbol{x}_{2}^{T} \boldsymbol{t}^{\wedge} \boldsymbol{R} \boldsymbol{x}_{1}
x2Tt∧x2=x2Tt∧Rx1
观察等式左侧,
t
∧
x
2
t^{\wedge} x_{2}
t∧x2 是一个与
t
\boldsymbol{t}
t 和
x
2
\boldsymbol{x}_{2}
x2 都垂直的向量。把它再和
x
2
\boldsymbol{x}_{2}
x2 做内积时,将得
到
0。因此,我们就得到了一个简洁的式子:
x
2
T
t
∧
R
x
1
=
0
\boldsymbol{x}_{2}^{T} \boldsymbol{t}^{\wedge} \boldsymbol{R} \boldsymbol{x}_{1}=0
x2Tt∧Rx1=0
重新代入
p
1
,
p
2
,
\boldsymbol{p}_{1}, \boldsymbol{p}_{2},
p1,p2, 有:
p
2
T
K
−
T
t
∧
R
K
−
1
p
1
=
0
\boldsymbol{p}_{2}^{T} \boldsymbol{K}^{-T} \boldsymbol{t}^{\wedge} \boldsymbol{R} \boldsymbol{K}^{-1} \boldsymbol{p}_{1}=0
p2TK−Tt∧RK−1p1=0
这两个式子都称为对极约束,它以形式简洁著名。它的几何意义是
O
1
,
P
,
O
2
O_{1}, P, O_{2}
O1,P,O2 三者 共面。对极约束中同时包含了平移和旋转。我们把中间部分记作两个矩阵: 基础矩阵 (Fundamental Matrix)
F
F
F 和本质矩阵(Essential Matrix)
E
\boldsymbol{E}
E, 可以进一步简化对极约束:
E
=
t
∧
R
,
F
=
K
−
T
E
K
−
1
,
x
2
T
E
x
1
=
p
2
T
F
p
1
=
0
\boldsymbol{E}=\boldsymbol{t}^{\wedge} \boldsymbol{R}, \quad \boldsymbol{F}=\boldsymbol{K}^{-T} \boldsymbol{E} \boldsymbol{K}^{-1}, \quad \boldsymbol{x}_{2}^{T} \boldsymbol{E} \boldsymbol{x}_{1}=\boldsymbol{p}_{2}^{T} \boldsymbol{F} \boldsymbol{p}_{1}=0
E=t∧R,F=K−TEK−1,x2TEx1=p2TFp1=0
对极约束简洁地给出了两个匹配点的空间位置关系。于是,相机位姿估计问题变为以下两步:
- 根据配对点的像素位置,求出 E \boldsymbol{E} E 或者 F \boldsymbol{F} F
- 根据 E \boldsymbol{E} E 或者 F , \boldsymbol{F}, F, 求出 R , t \boldsymbol{R}, \boldsymbol{t} R,t
本质矩阵
根据定义,本质矩阵
E
=
t
∧
R
E = t ^{\wedge} R
E=t∧R 。它是一个
3
×
3
3 \times 3
3×3 的矩阵,内有 9 个未知数。那么,是
不是任意一个
3
×
3
3 \times 3
3×3 的矩阵都可以被当成本质矩阵呢?从
E
E
E 的构造方式上看,有以下值得
注意的地方:
- 基础矩阵是由对极约束定义的。由于对极约束是等式为零的约束,所以对 E E E 乘以任 意非零常数后,对极约束依然满足。我们把这件事情称为 E E E 在不同尺度下是等价的。
- 根据 E = t ∧ R , E = t ^{\wedge} R , E=t∧R, 可以证明本质矩阵 E E E 的奇异值必定是 [ σ , σ , 0 ] T [\sigma, \sigma, 0]^{T} [σ,σ,0]T 的形式。这称为本质矩阵的内在性质。
- 另一方面,由于平移和旋转各有三个自由度,故 t ∧ R t ^{\wedge} R t∧R 共有六个自由度。但由于尺度等价性, E E E 实际上有五个自由度。
E
E
E 具有五个自由度的事实,表明我们最少可以用五对点来求解
E
E
E。但是,E 的内在性质是一种非线性性质,在求解线性方程时会带来麻烦,因此,也可以只考虑它的尺度等价性,使用八对点来估计
E
E
E,这就是经典的八点法
(Eight-point-algorithm)。
考虑一对匹配点,它们的归一化坐标为:
x
1
=
[
u
1
,
v
1
,
1
]
T
,
x
2
=
[
u
2
,
v
2
,
1
]
T
x_{1}=\left[u_{1}, v_{1}, 1\right]^{T}, \boldsymbol{x}_{2}=\left[u_{2}, v_{2}, 1\right]^{T}
x1=[u1,v1,1]T,x2=[u2,v2,1]T 。根据对
极约束,有:
(
u
1
,
v
1
,
1
)
(
e
1
e
2
e
3
e
4
e
5
e
6
e
7
e
8
e
9
)
(
u
2
v
2
1
)
=
0
\left(\begin{array}{ll} u_{1}, v_{1}, 1 \end{array}\right)\left(\begin{array}{lll} e_{1} & e_{2} & e_{3} \\ e_{4} & e_{5} & e_{6} \\ e_{7} & e_{8} & e_{9} \end{array}\right)\left(\begin{array}{l} u_{2} \\ v_{2} \\ 1 \end{array}\right)=0
(u1,v1,1)⎝⎛e1e4e7e2e5e8e3e6e9⎠⎞⎝⎛u2v21⎠⎞=0
我们把矩阵
E
\boldsymbol{E}
E 展开,写成向量的形式:
e
=
[
e
1
,
e
2
,
e
3
,
e
4
,
e
5
,
e
6
,
e
7
,
e
8
,
e
9
]
T
\boldsymbol{e}=\left[e_{1}, e_{2}, e_{3}, e_{4}, e_{5}, e_{6}, e_{7}, e_{8}, e_{9}\right]^{T}
e=[e1,e2,e3,e4,e5,e6,e7,e8,e9]T
那么对极约束可以写成与 e 有关的线性形式:
[
u
1
u
2
,
u
1
v
2
,
u
1
,
v
1
u
2
,
v
1
v
2
,
u
2
,
v
2
,
1
]
⋅
e
=
0
\left[u_{1} u_{2}, u_{1} v_{2}, u_{1}, v_{1} u_{2}, v_{1} v_{2}, u_{2}, v_{2}, 1\right] \cdot \boldsymbol{e}=0
[u1u2,u1v2,u1,v1u2,v1v2,u2,v2,1]⋅e=0
这八个方程构成了一个线性方程组。它的系数矩阵由特征点位置构成,大小为
8
×
9
8 \times 9
8×9。
接下来的问题是如何根据已经估得的本质矩阵
E
,
\boldsymbol{E},
E, 恢复出相机的运动
R
,
t
\boldsymbol{R}, \boldsymbol{t}
R,t 。这个过 程是由奇异值分解 (SVD) 得到的。设
E
E
E 的 SVD 分解为
E
=
U
Σ
V
T
E = U \Sigma V ^{T}
E=UΣVT
其中
U
,
V
U , V
U,V 为正交阵,
Σ
\Sigma
Σ 为奇异值矩阵。根据
E
E
E 的内在性质,我们知道
Σ
=
diag
(
σ
,
σ
,
0
)
\Sigma =\operatorname{diag}(\sigma, \sigma, 0)
Σ=diag(σ,σ,0)在SVD分解中,对于任意一个
E
E
E, 存在两个可能的
t
,
R
t,R
t,R与它对应:
t
1
∧
=
U
R
Z
(
π
2
)
Σ
U
T
,
R
1
=
U
R
Z
T
(
π
2
)
V
T
t
2
∧
=
U
R
Z
(
−
π
2
)
Σ
U
T
,
R
2
=
U
R
Z
T
(
−
π
2
)
V
T
\begin{array}{l} t _{1}^{\wedge}= U R _{Z}\left(\frac{\pi}{2}\right) \Sigma U ^{T}, \quad R _{1}= U R _{Z}^{T}\left(\frac{\pi}{2}\right) V ^{T} \\ t _{2}^{\wedge}= U R _{Z}\left(-\frac{\pi}{2}\right) \Sigma U ^{T}, \quad R _{2}= U R _{Z}^{T}\left(-\frac{\pi}{2}\right) V ^{T} \end{array}
t1∧=URZ(2π)ΣUT,R1=URZT(2π)VTt2∧=URZ(−2π)ΣUT,R2=URZT(−2π)VT
其中
R
Z
(
π
2
)
R _{Z}\left(\frac{\pi}{2}\right)
RZ(2π) 表示沿
Z
Z
Z 轴旋转 90 度得到的旋转矩阵。同时,由于
−
E
- E
−E 和
E
E
E 等价,所 以对任意一个
t
t
t 取负号,也会得到同样的结果。因此,从
E
E
E 分解到
t
,
R
t, R
t,R 时,一共存在四个可能的解。幸运的是,只有一种解中,
P
P
P 在两个相机中都具有正的深度。因此,只要把任意一点代入四种解中,检测该点在两个相机下的深度,就可以确定哪个解是正确的了。
剩下的问题还有一个: 根据线性方程解出的
E
,
E,
E, 可能不满足
E
E
E 的内在性质一它 的奇异值不一定为
σ
,
σ
,
0
\sigma, \sigma, 0
σ,σ,0 的形式。这时,在做 SVD 时,我们会刻意地把 \Sigma 矩阵调整成上面的样子。通常的做法是,对八点法求得的 E 进行 SVD 分解后,会得到奇异值矩阵
Σ
=
diag
(
σ
1
,
σ
2
,
σ
3
)
,
\Sigma =\operatorname{diag}\left(\sigma_{1}, \sigma_{2}, \sigma_{3}\right),
Σ=diag(σ1,σ2,σ3), 不妨设
σ
1
≥
σ
2
≥
σ
3
\sigma_{1} \geq \sigma_{2} \geq \sigma_{3}
σ1≥σ2≥σ3。 取:
E
=
U
diag
(
σ
1
+
σ
2
2
,
σ
1
+
σ
2
2
,
0
)
V
T
E = U \operatorname{diag}\left(\frac{\sigma_{1}+\sigma_{2}}{2}, \frac{\sigma_{1}+\sigma_{2}}{2}, 0\right) V ^{T}
E=Udiag(2σ1+σ2,2σ1+σ2,0)VT
当给定的点数多于八对时,我们可以计算一个最小二乘解。
min
e
∥
A
e
∥
2
2
=
min
e
e
T
A
T
A
e
\min _{e}\|A e\|_{2}^{2}=\min _{e} e^{T} A^{T} A e
emin∥Ae∥22=emineTATAe
不过,当可能存在误匹配的情况时,倾向于使用随机采样一致性
(Random Sample Concensus, RANSAC)来求,而不是最小二乘。
单应矩阵
除了基本矩阵和本质矩阵,我们还有一种称为单应矩阵(Homography) H H H 的东西, 它描述了两个平面之间的映射关系。若场景中的特征点都落在同一平面上(比如墙,地面等),则可以通过单应性来进行运动估计。这种情况在无人机携带的俯视相机,或扫地机携带的顶视相机中比较常见。
单应性在SLAM中具重要意义。当特征点共面,或者相机发生纯旋转的时候,基础矩阵的自由度下降,这就出现了所谓的退化(degenerate)。现实中的数据总包含一些噪声,这时候如果我们继续使用八点法求解基础矩阵,基础矩阵多余出来的自由度将会主要由噪声决定。为了能够避免退化现象造成的影响,通常我们会同时估计基础矩阵 F F F 和单应矩阵 H H H , 选择重投影误差比较小的那个作为最终的运动估计矩阵。