问题描述
已知n个空间点 P i = [ x i , y i , z i ] T P_i=[x_i,y_i,z_i]^T Pi=[xi,yi,zi]T,其投影的像素坐标 p i = [ u i , v i ] T p_i=[u_i,v_i]^T pi=[ui,vi]T求相机的位姿R,T。
问题分析
根据相机模型,像素点和空间点的位置关系:
s
i
[
u
i
v
i
1
]
=
K
T
[
x
i
y
i
z
i
1
]
s_i\begin{bmatrix}u_i \\v_i \\1\\\end{bmatrix}=KT \begin{bmatrix} x_i \\ y_i \\ z_i\\ 1\\ \end{bmatrix}
si
uivi1
=KT
xiyizi1
即:
s
i
u
i
=
K
T
P
i
s_iu_i=KTP_i
siui=KTPi
由于存在噪声误差,因此以最小化误差平方和为目标构建最小二乘问题:
T
∗
=
a
r
g
min
T
1
2
∑
i
=
1
n
∥
u
i
−
1
s
i
K
T
P
i
∥
2
2
T^{*}=arg \min_{T} \frac{1}{2} \sum_{i=1}^{n} \begin{Vmatrix}u_i-\frac{1}{s_i}KTP_i\end{Vmatrix}_{2}^{2}
T∗=argTmin21i=1∑n
ui−si1KTPi
22
因为这个误差是将3D点的理论投影位置与观测到的实际投影位置之间的误差,因此称为重投影误差:
e
i
=
u
i
−
1
s
i
K
T
P
i
e_i=u_i-\frac{1}{s_i}KTP_i
ei=ui−si1KTPi
按照之前讲过的高斯牛顿法进行求解:
旋转矩阵本身带有约束,即正交且行列式为1。而有约束的优化问题比无约束的优化问题复杂的多。因为李代数的特点,李代数表示的天然满足旋转矩阵的约束,因此通常使用李代数进行表示来求解。
问题求解
步骤概述
- 初始化位姿:使用闭式解法(如EPnP、DLT)获取初始相机位姿 ( R ) 和 ( t )。
- 构建重投影误差:将3D点投影到图像平面,计算与观测值的误差。
- 计算雅可比矩阵:分析误差关于位姿参数的导数,指导优化方向。
- 迭代优化:使用高斯-牛顿或Levenberg-Marquardt算法更新位姿,直至收敛。
详细推导与算法
1. 重投影误差定义
设第
i
i
i 个3D点为
P
i
=
(
X
i
,
Y
i
,
Z
i
)
⊤
\mathbf{P}_i = (X_i, Y_i, Z_i)^\top
Pi=(Xi,Yi,Zi)⊤,对应的图像观测为
p
i
=
(
u
i
,
v
i
)
⊤
\mathbf{p}_i = (u_i, v_i)^\top
pi=(ui,vi)⊤。相机位姿用李代数
ξ
∈
s
e
(
3
)
\boldsymbol{\xi} \in \mathfrak{se}(3)
ξ∈se(3) 表示,对应的变换矩阵为
T
=
exp
(
ξ
∧
)
\mathbf{T} = \exp(\boldsymbol{\xi}^\wedge)
T=exp(ξ∧)。将
P
i
\mathbf{P}_i
Pi 变换到相机坐标系:
P
i
′
=
R
P
i
+
t
=
(
X
i
′
,
Y
i
′
,
Z
i
′
)
⊤
.
\mathbf{P}'_i = \mathbf{R} \mathbf{P}_i + \mathbf{t} = (X'_i, Y'_i, Z'_i)^\top.
Pi′=RPi+t=(Xi′,Yi′,Zi′)⊤.
投影到归一化平面:
Z
i
′
p
^
i
′
=
K
P
i
′
Z_i'\mathbf{\hat{p}}'_i =\mathbf{KP_i'}
Zi′p^i′=KPi′
即:
[
Z
i
′
u
i
Z
i
′
v
i
Z
i
′
]
=
[
f
x
0
c
x
0
f
y
c
y
0
0
1
]
[
X
i
′
Y
i
′
Z
i
′
]
.
\begin{bmatrix} Z'_iu_i \\ Z'_iv_i \\ Z'_i \end{bmatrix} = \begin{bmatrix} f_x & 0 & c_x \\ 0 & f_y & c_y \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} X'_i \\ Y'_i \\ Z'_i \end{bmatrix}.
Zi′uiZi′viZi′
=
fx000fy0cxcy1
Xi′Yi′Zi′
.
解得:
{
u
i
=
f
x
X
i
′
Z
i
′
+
c
x
v
i
=
f
y
Y
i
′
Z
i
′
+
c
y
\begin{cases} u_i=f_x \frac{X_i'}{Z_i'}+c_x\\ v_i=f_y \frac{Y_i'}{Z_i'}+c_y\\ \end{cases}
{ui=fxZi′Xi′+cxvi=fyZi′Yi′+cy
假设相机内参已知,重投影误差为:
e
i
=
p
^
i
′
−
p
i
.
\mathbf{e}_i = \mathbf{\hat{p}}'_i - \mathbf{p}_i.
ei=p^i′−pi.
2. 误差函数与优化目标
最小化所有点的重投影误差平方和:
min
ξ
1
2
∑
i
=
1
n
∥
e
i
∥
2
.
\min_{\boldsymbol{\xi}} \frac{1}{2} \sum_{i=1}^n \|\mathbf{e}_i\|^2.
ξmin21i=1∑n∥ei∥2.
3. 雅可比矩阵计算
误差关于李代数
ξ
\boldsymbol{\xi}
ξ 的雅可比矩阵
J
i
\mathbf{J}_i
Ji 分为两部分:
J
i
=
∂
e
i
∂
ξ
=
∂
e
i
∂
P
i
′
⋅
∂
P
i
′
∂
ξ
.
\mathbf{J}_i = \frac{\partial \mathbf{e}_i}{\partial \boldsymbol{\xi}} = \frac{\partial \mathbf{e}_i}{\partial \mathbf{P}'_i} \cdot \frac{\partial \mathbf{P}'_i}{\partial \boldsymbol{\xi}}.
Ji=∂ξ∂ei=∂Pi′∂ei⋅∂ξ∂Pi′.
- 误差对相机坐标系点的导数(2×3矩阵):
∂ e ∂ P ′ = [ ∂ u ∂ X ′ ∂ u ∂ Y ′ ∂ u ∂ Z ′ ∂ v ∂ X ′ ∂ v ∂ Y ′ ∂ v ∂ Z ′ ] = [ f x Z ′ 0 − f x X ′ Z ′ 2 0 f y Z ′ − f y Y ′ Z ′ 2 ] . \frac{\partial e}{\partial P'} = \left[ \begin{array}{ccc} \frac{\partial u}{\partial X'} & \frac{\partial u}{\partial Y'} & \frac{\partial u}{\partial Z'} \\ \frac{\partial v}{\partial X'} & \frac{\partial v}{\partial Y'} & \frac{\partial v}{\partial Z'} \end{array} \right] = \left[ \begin{array}{ccc} \frac{f_x}{Z'} & 0 & -\frac{f_x X'}{Z'^2} \\ 0 & \frac{f_y}{Z'} & -\frac{f_y Y'}{Z'^2} \end{array} \right]. ∂P′∂e=[∂X′∂u∂X′∂v∂Y′∂u∂Y′∂v∂Z′∂u∂Z′∂v]=[Z′fx00Z′fy−Z′2fxX′−Z′2fyY′]. - 相机坐标系点对位姿的导数(3×6矩阵,考虑李代数扰动):
∂ ( T P ) ∂ δ ξ = [ I − P ′ ∧ 0 T 0 T ] \frac{\partial (TP)}{\partial \delta \xi} = \begin{bmatrix} I & -P'^\wedge \\ 0^T & 0^T \end{bmatrix} ∂δξ∂(TP)=[I0T−P′∧0T]
由于所定义的P只有三维,因此取出结果的前三维,展开后为 2×6 矩阵:
J i = − [ f x Z i ′ 0 − f x X i ′ ( Z i ′ ) 2 − f x X i ′ Y i ′ ( Z i ′ ) 2 f x ( 1 + ( X i ′ ) 2 ( Z i ′ ) 2 ) − f x Y i ′ Z i ′ 0 f y Z i ′ − f y Y i ′ ( Z i ′ ) 2 − f y ( 1 + ( Y i ′ ) 2 ( Z i ′ ) 2 ) f y X i ′ Y i ′ ( Z i ′ ) 2 f y X i ′ Z i ′ ] . \mathbf{J}_i = -\begin{bmatrix} \frac{f_x}{Z'_i} & 0 & -\frac{f_xX'_i}{(Z'_i)^2} & -\frac{f_xX'_i Y'_i}{(Z'_i)^2} & f_x\left(1 + \frac{(X'_i)^2}{(Z'_i)^2}\right) & -\frac{f_xY'_i}{Z'_i} \\ 0 & \frac{f_y}{Z'_i} & -\frac{f_yY'_i}{(Z'_i)^2} & -f_y(1 + \frac{(Y'_i)^2}{(Z'_i)^2}) & \frac{f_yX'_i Y'_i}{(Z'_i)^2} & \frac{f_yX'_i}{Z'_i} \end{bmatrix}. Ji=− Zi′fx00Zi′fy−(Zi′)2fxXi′−(Zi′)2fyYi′−(Zi′)2fxXi′Yi′−fy(1+(Zi′)2(Yi′)2)fx(1+(Zi′)2(Xi′)2)(Zi′)2fyXi′Yi′−Zi′fxYi′Zi′fyXi′ .
4. 迭代优化过程
-
高斯-牛顿法:
- 计算误差和雅可比:对每个点计算 e i \mathbf{e}_i ei 和 J i \mathbf{J}_i Ji。
- 构建线性系统:堆叠所有 J i \mathbf{J}_i Ji 和 e i \mathbf{e}_i ei,得到 J ⊤ J Δ ξ = − J ⊤ e \mathbf{J}^\top \mathbf{J} \Delta \boldsymbol{\xi} = -\mathbf{J}^\top \mathbf{e} J⊤JΔξ=−J⊤e。
- 求解增量:解线性方程 Δ ξ \Delta \boldsymbol{\xi} Δξ。
- 更新位姿: ξ ← ξ + Δ ξ \boldsymbol{\xi} \leftarrow \boldsymbol{\xi} + \Delta \boldsymbol{\xi} ξ←ξ+Δξ。
- 判断收敛:若 Δ ξ \Delta \boldsymbol{\xi} Δξ 足够小或误差不再下降,停止迭代。
-
Levenberg-Marquardt:通过引入阻尼因子 λ \lambda λ 稳定求解:
( J ⊤ J + λ I ) Δ ξ = − J ⊤ e . (\mathbf{J}^\top \mathbf{J} + \lambda \mathbf{I}) \Delta \boldsymbol{\xi} = -\mathbf{J}^\top \mathbf{e}. (J⊤J+λI)Δξ=−J⊤e.
关键点
- 李代数参数化:避免旋转矩阵的约束,简化优化过程。
- 雅可比推导:通过扰动模型或链式法则计算导数,确保优化方向正确。
- 鲁棒性:初始值影响收敛性,建议先用闭式解法初始化。