北邮鲁鹏老师的课程《计算机视觉之三维重建(深入浅出sfm和SLAM核心算法)》 笔记
6,多视几何
6.1,运动恢复结构问题
问题:
n
n
n个点
X
j
X_j
Xj在
m
m
m张图像中对应点
x
i
j
=
M
i
X
j
x_{ij}=M_iX_j
xij=MiXj,已知
x
i
j
x_{ij}
xij,求
M
i
,
X
j
M_i,X_j
Mi,Xj
三种任务:
- 欧式结构恢复:内参数 K i K_i Ki已知 (实际常用)
- 仿射结构恢复:仿射相机,内外参数未知
- 透视结构恢复:透视相机,内外参数未知
6.2,欧式结构
问题:
x
i
j
,
K
i
⇒
X
j
,
R
i
,
T
i
x_{ij}, K_i\Rightarrow X_j,R_i,T_i
xij,Ki⇒Xj,Ri,Ti
考虑两视图的情形,求解算法:
- 归一化八点法求解基础矩阵 F F F
- 求解本质矩阵 E = K 2 T F K 1 E=K_2^TFK_1 E=K2TFK1
- 分解本质矩阵 E → R 、 T E\rightarrow R、T E→R、T
- 三角化得到 X i j X_{ij} Xij
本质矩阵分解:
定义:
W
=
[
0
−
1
0
1
0
0
0
0
1
]
Z
=
[
0
1
0
−
1
0
0
0
0
0
]
W=\left[ \begin{array}{c} 0&-1&0 \\ 1&0&0 \\ 0&0&1 \end{array}\right]\quad Z=\left[ \begin{array}{c} 0&1&0 \\ -1&0&0 \\ 0&0&0 \end{array}\right]
W=⎣⎡010−100001⎦⎤Z=⎣⎡0−10100000⎦⎤
则
Z
=
−
d
i
a
g
(
1
,
1
,
0
)
W
=
d
i
a
g
(
1
,
1
,
0
)
W
T
Z=-diag(1,1,0)W=diag(1,1,0)W^T
Z=−diag(1,1,0)W=diag(1,1,0)WT
[
T
×
]
=
k
U
Z
U
T
[T_\times]=kUZU^T
[T×]=kUZUT,
U
U
U单位正交阵
证明1:
在垂直与向量 T T T的平面上,取两个互相垂直的单位向量 u 1 , u 2 u_1,u_2 u1,u2,跟 T T T的单位向量 u 3 u_3 u3一起构成单位正交矩阵 U = [ u 1 , u 2 , u 3 ] U=[u_1,u_2,u_3] U=[u1,u2,u3]
则: [ T × ] U = T × [ u 1 , u 2 , u 3 ] = [ T × u 1 , T × u 2 , T × u 3 ] = k [ u 2 , − u 1 , 0 ] = k U Z [T_\times]U=T\times[u_1,u_2,u_3]=[T\times u_1,T\times u_2,T\times u_3]=k[u_2,-u_1,0]=kUZ [T×]U=T×[u1,u2,u3]=[T×u1,T×u2,T×u3]=k[u2,−u1,0]=kUZ
证明2:
反对称矩阵的性质,参考https://geek-docs.com/linear-algebra/matrix/fanduichen-matrix.html
不考虑符号尺度,
[
T
×
]
=
U
d
i
a
g
(
1
,
1
,
0
)
W
U
T
=
U
d
i
a
g
(
1
,
1
,
0
)
W
T
U
T
[T_\times]=Udiag(1,1,0)WU^T=Udiag(1,1,0)W^TU^T
[T×]=Udiag(1,1,0)WUT=Udiag(1,1,0)WTUT
E
=
[
T
×
]
R
=
U
d
i
a
g
(
1
,
1
,
0
)
W
U
T
R
=
U
d
i
a
g
(
1
,
1
,
0
)
W
T
U
T
R
E=[T_\times]R=Udiag(1,1,0)WU^TR=Udiag(1,1,0)W^TU^TR
E=[T×]R=Udiag(1,1,0)WUTR=Udiag(1,1,0)WTUTR
SVD分解:
E
=
U
d
i
a
g
(
1
,
1
,
0
)
V
T
E=Udiag(1,1,0)V^T
E=Udiag(1,1,0)VT
可求得:
R
=
d
e
t
(
U
W
V
T
)
U
W
V
T
R=det(UWV^T)UWV^T
R=det(UWVT)UWVT或
R
=
d
e
t
(
U
W
T
V
T
)
U
W
T
V
T
R=det(UW^TV^T)UW^TV^T
R=det(UWTVT)UWTVT
T
=
±
u
3
T=\pm u_3
T=±u3
得到四组解,选择重建的点的z轴为正的那组解。
欧式结构的恢复歧义:恢复的结构与真实场景存在相似变换
6.3,仿射结构
投影(透视):
M
=
[
A
b
v
1
]
M=\left[ \begin{array}{ll} A&b \\ v&1 \end{array}\right]
M=[Avb1]
仿射(弱透视):
M
=
[
A
b
0
1
]
M=\left[ \begin{array}{ll} A&b \\ 0&1 \end{array}\right]
M=[A0b1]
x
E
=
(
m
1
,
m
2
)
T
=
[
A
b
]
X
=
A
X
E
+
b
x^E=(m_1,m_2)^T=[A\ b]X=AX^E+b
xE=(m1,m2)T=[A b]X=AXE+b
问题:
x
i
j
⇒
X
j
,
A
i
,
b
i
x_{ij}\Rightarrow X_j,A_i,b_i
xij⇒Xj,Ai,bi
求解:
- 数据中心化
数据的中心 x ˉ i = 1 n ∑ k x i k , X ˉ = 1 n ∑ i X i \bar{x}_i=\frac{1}{n}\sum_k{x_{ik}},\bar{X}=\frac{1}{n}\sum_i{X_i} xˉi=n1∑kxik,Xˉ=n1∑iXi,将数据减去其中心 x ^ i j = x i j − x ˉ i , X ^ i = X i − X ˉ \hat{x}_{ij}=x_{ij}-\bar{x}_i,\hat{X}_i=X_i-\bar{X} x^ij=xij−xˉi,X^i=Xi−Xˉ,可得: x ^ i j = A i X ^ j \hat{x}_{ij}=A_i\hat{X}_j x^ij=AiX^j
写成矩阵形式:
D = [ x ^ 11 x ^ 12 ⋯ x ^ 1 n x ^ 21 x ^ 22 ⋯ x ^ 2 n ⋱ x ^ m 1 x ^ m 2 ⋯ x ^ m n ] = [ A 1 A 2 ⋮ A m ] [ X ^ 1 X ^ 2 ⋯ X ^ n ] = M S D=\left[ \begin{array}{ll} \hat{x}_{11}&\hat{x}_{12}&\cdots&\hat{x}_{1n} \\ \hat{x}_{21}&\hat{x}_{22}&\cdots&\hat{x}_{2n} \\ &&\ddots&\\ \hat{x}_{m1}&\hat{x}_{m2}&\cdots&\hat{x}_{mn} \\ \end{array}\right]= \left[ \begin{array}{c}A_1\\A_2\\\vdots\\A_m\end{array}\right] \left[ \begin{array}{c}\hat{X}_1&\hat{X}_2&\cdots&\hat{X}_n\end{array}\right] =MS D=⎣⎢⎢⎡x^11x^21x^m1x^12x^22x^m2⋯⋯⋱⋯x^1nx^2nx^mn⎦⎥⎥⎤=⎣⎢⎢⎢⎡A1A2⋮Am⎦⎥⎥⎥⎤[X^1X^2⋯X^n]=MS - 矩阵分解
对 D D D进行奇异值分解,取不为0的前三个特征值(理论上只有三个特征值不为0), D = U 3 W 3 V 3 T D=U_3W_3V_3^T D=U3W3V3T,得 M = U 3 , S = W 3 V 3 T M=U_3,S=W_3V_3^T M=U3,S=W3V3T,从而得到 A i , X j A_i,X_j Ai,Xj。
仿射结构的恢复歧义:差一个仿射变换
D
=
M
∗
S
∗
=
M
H
⋅
H
−
1
S
D=M^*S^*=MH\cdot H^{-1}S
D=M∗S∗=MH⋅H−1S
4,透视结构
相机模型:
x
i
j
=
M
i
X
j
x_{ij}=M_iX_j
xij=MiXj
x
=
M
X
=
M
H
−
1
H
X
=
M
∗
X
∗
x=MX=MH^{-1}HX=M^*X^*
x=MX=MH−1HX=M∗X∗相差一个H变换
代数方法
两视图:
x
1
j
=
M
1
X
j
=
K
1
[
I
,
0
]
X
j
x
2
j
=
M
2
X
j
=
K
2
[
R
,
T
]
X
j
x_{1j}=M_1X_j=K_1[I,0]X_j \\ x_{2j}=M_2X_j=K_2[R,T]X_j
x1j=M1Xj=K1[I,0]Xjx2j=M2Xj=K2[R,T]Xj
- 求解 F F F矩阵
- 分解 F → M 1 , M 2 F\rightarrow M_1,M_2 F→M1,M2
令 M 1 ∗ = M 1 H − 1 = [ I , 0 ] , M 2 ∗ = M 2 H − 1 = [ A , b ] M_1^*=M_1H^{-1}=[I,0],M_2^*=M_2H^{-1}=[A,b] M1∗=M1H−1=[I,0],M2∗=M2H−1=[A,b]
x ′ = [ A b ] X = A [ I 0 ] X + b = A x + b x'=[A\ b]X=A[I\ 0]X+b=Ax+b x′=[A b]X=A[I 0]X+b=Ax+b
x ′ × b = A x × b x'\times b=Ax\times b x′×b=Ax×b
0 = − x ′ T ( x ′ × b ) = − x ′ T ( A x × b ) = x ′ T ( b × A x ) = x ′ T [ b × ] A x = x ′ T F x 0=-x'^T(x'\times b)=-x'^T(Ax\times b)=x'^T(b\times Ax)=x'^T[b_\times]Ax=x'^TFx 0=−x′T(x′×b)=−x′T(Ax×b)=x′T(b×Ax)=x′T[b×]Ax=x′TFx
F = [ b × ] A F=[b_\times]A F=[b×]A
F T b = 0 F^Tb=0 FTb=0,得到 b b b为 F T F^T FT最小奇异值的右奇异向量( b b b为极点)
A = − [ b × ] F A=-[b_\times]F A=−[b×]F
- 三角化
捆绑调整BA
因式分解:都可见重建点少
代数法:两两算,有累积误差
BA,最小化重投影误差:
E
(
M
,
X
)
=
∑
i
=
1
m
∑
j
=
1
n
D
(
x
i
j
,
M
i
X
j
)
2
E(M,X)=\sum_{i=1}^m\sum_{j=1}^nD(x_{ij}, M_iX_j)^2
E(M,X)=∑i=1m∑j=1nD(xij,MiXj)2
实际用于sfm最后一步,用分解或代数法的解作为其初始解
补充
PnP
Perspective-n-Points,已知
K
,
X
j
,
x
j
K,X_j,x_j
K,Xj,xj,求相机位姿。
P3P:
已知相机内参数
K
K
K,空间三个点
A
,
B
,
C
A,B,C
A,B,C及其对应图像点
a
,
b
,
c
⇒
R
,
T
a,b,c\Rightarrow R,T
a,b,c⇒R,T
求解:
在相机坐标系下: p = K P p=KP p=KP,所以图像点对应的相机坐标系下的三维坐标 P = K − 1 p P=K^{-1}p P=K−1p
可得: K − 1 a = o a → , K − 1 b = o b → , K − 1 a = o c → K^{-1}a=\overrightarrow{oa},K^{-1}b=\overrightarrow{ob},K^{-1}a=\overrightarrow{oc} K−1a=oa,K−1b=ob,K−1a=oc(直线的方向,相差一个尺度)
从而可求得: c o s < o a → , o b → > , c o s < o a → , o c → > , c o s < o b → , o c → > cos<\overrightarrow{oa}, \overrightarrow{ob}>,cos<\overrightarrow{oa}, \overrightarrow{oc}>,cos<\overrightarrow{ob}, \overrightarrow{oc}> cos<oa,ob>,cos<oa,oc>,cos<ob,oc>
求解方程:
{ O A 2 + O B 2 − 2 O A ⋅ O B c o s < o a → , o b → > = A B 2 O B 2 + O C 2 − 2 O B ⋅ O C c o s < o b → , o c → > = B C 2 O A 2 + O C 2 − 2 O A ⋅ O C c o s < o a → , o c → > = A C 2 \left\{\begin{array}{c} OA^2+OB^2-2OA\cdot OBcos<\overrightarrow{oa}, \overrightarrow{ob}>=AB^2\\ OB^2+OC^2-2OB\cdot OCcos<\overrightarrow{ob}, \overrightarrow{oc}>=BC^2\\ OA^2+OC^2-2OA\cdot OCcos<\overrightarrow{oa}, \overrightarrow{oc}>=AC^2\\ \end{array}\right. ⎩⎪⎨⎪⎧OA2+OB2−2OA⋅OBcos<oa,ob>=AB2OB2+OC2−2OB⋅OCcos<ob,oc>=BC2OA2+OC2−2OA⋅OCcos<oa,oc>=AC2
得到四组 O A , O B , O C OA,OB,OC OA,OB,OC的解,结合另一对点对应D来确定唯一解。
根据上面求得的方向和长度,可以确定 A , B , C A,B,C A,B,C在相机坐标系下的坐标,结合其世界坐标,可求得 R , T R,T R,T
RANSAC
RANdom SAmple Consensus
- 随机均匀采样获取模型求解所需最小子集,个数为 s s s;
- 估计模型参数;
- 计算剩余样本跟当前模型的一致性,统计内点数作为模型分;
- 重复1-3步 N N N次,输出得分最高的模型;
样本的外点率为
e
e
e,
N
N
N次采样中,至少有一次采样的点全部是内点的概率为
p
p
p,则有:
(
1
−
(
1
−
e
)
s
)
N
=
1
−
p
(1-(1-e)^s)^N=1-p
(1−(1−e)s)N=1−p
所以,给定
p
p
p,至少需要
N
N
N次迭代:
N
=
log
(
1
−
p
)
/
log
(
1
−
(
1
−
e
)
s
)
N=\log(1-p)/\log(1-(1-e)^s)
N=log(1−p)/log(1−(1−e)s)