问题描述
已知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.
关键点
- 李代数参数化:避免旋转矩阵的约束,简化优化过程。
- 雅可比推导:通过扰动模型或链式法则计算导数,确保优化方向正确。
- 鲁棒性:初始值影响收敛性,建议先用闭式解法初始化。