文章目录
- 1. 运动恢复结构
- 2. 欧式结构恢复
- 2.1. 基本概念
- 2.2. 求解 R 和 T
- 2.3. 欧式结构恢复的歧义性
- 3. 仿射结构恢复
- 3.1. 基本概念
- 3.2. 基于因式分解的仿射结构恢复
- 3.3. 仿射结构恢复的歧义性与局限性
课程视频链接: 计算机视觉之三维重建(深入浅出SfM与SLAM核心算法)—— 6.多视图几何(运动恢复结构)。
1. 运动恢复结构
运动恢复结构(Structure from Motion, SfM)是计算机视觉领域的核心问题,旨在从一组多视角二维图像序列中恢复三维场景的结构信息(3D 点)和相机的运动参数(相机位姿),如下图所示:
运动恢复结构问题的数学模型可表述为::
根据相机模型和已知参数的不同,SfM 可分为三类典型任务:
任务类型 | 相机模型 | 参数 | 求解目标 |
---|---|---|---|
欧式结构恢复 | 透视相机 | 已知相机内参数矩阵 | 相机外参以及 3D 点坐标 |
仿射结构恢复 | 仿射相机 | 相机内外参数均未知 | 投影矩阵和 3D 点坐标 |
透视结构恢复 | 透视相机 | 相机内外参数均未知 | 投影矩阵和 3D 点坐标 |
2. 欧式结构恢复
2.1. 基本概念
欧式结构恢复是最常见的任务,相机内参(焦距、主点等)已通过标定确定,仅需计算相机外参(旋转矩阵 R \mathbf{R} R 和平移向量 T T T)和 3D 场景点坐标。欧式结构恢复问题的数学模型可以表述如下:
参考博客:三维重建 —— 4. 三维重建基础与极几何,已知相机内外参数及两张图像中匹配的像素点对,便可通过三角化方法求解该点的 3D 坐标。三角化问题一般有线性法和非线性法两种求解方法,如下图所示:
在欧式结构恢复中,相机外参是未知的,所以我们无法直接使用三角化求解 3D 点。考虑到基础矩阵 F \mathbf{F} F 与本质矩阵 E \mathbf{E} E 的数学关系( F = K 2 − T E K 1 − 1 \mathbf{F} = \mathbf{K}_2^{-T} \mathbf{E} \mathbf{K}_1^{-1} F=K2−TEK1−1)以及本质矩阵与相机外参的数学关系( E = T × R = [ T × ] R \mathbf{E} = T \times \mathbf{R} = [T_{\times}] \mathbf{R} E=T×R=[T×]R),我们可以使用归一化八点法求出基础矩阵 F \mathbf{F} F,然后求解出 E = K 2 T F K 1 \mathbf{E} = \mathbf{K}_2^T \mathbf{F} \mathbf{K}_1 E=K2TFK1,在从本质矩阵 E \mathbf{E} E 中分解出相机外参 R \mathbf{R} R 和 T,这样我们就可以使用三角化求解场景的 3D 点。
2.2. 求解 R 和 T
根据方程 x 2 T F x 1 = 0 x_2^T \mathbf{F} x_1 = 0 x2TFx1=0 和 E = K 2 T F K 1 \mathbf{E} = \mathbf{K}_2^T \mathbf{F} \mathbf{K}_1 E=K2TFK1,基础矩阵 F \mathbf{F} F 和本质矩阵 E \mathbf{E} E 只能确定到相差一个尺度因子(即对于任意非零常数 k k k, k F k\mathbf{F} kF 同样满足方程)。
为了辅助后续的推导,我们先定义两个矩阵:
W = ( 0 − 1 0 1 0 0 0 0 1 ) , Z = ( 0 1 0 − 1 0 0 0 0 0 ) \mathbf{W} = \begin{pmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{pmatrix}, \mathbf{Z} = \begin{pmatrix} 0 & 1 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} W= 010−100001 ,Z= 0−10100000 则有:
d i a g ( 1 , 1 , 0 ) W = ( 1 0 0 0 1 0 0 0 0 ) ( 0 − 1 0 1 0 0 0 0 1 ) = ( 0 − 1 0 1 0 0 0 0 0 ) = − Z diag(1, 1, 0) \mathbf{W} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{pmatrix} = \begin{pmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} = -\mathbf{Z} diag(1,1,0)W= 100010000 010−100001 = 010−100000 =−Z d i a g ( 1 , 1 , 0 ) W T = ( 1 0 0 0 1 0 0 0 0 ) ( 0 1 0 − 1 0 0 0 0 1 ) = ( 0 1 0 − 1 0 0 0 0 0 ) = Z diag(1, 1, 0) \mathbf{W}^T = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} 0 & 1 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 1 \end{pmatrix} = \begin{pmatrix} 0 & 1 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} = \mathbf{Z} diag(1,1,0)WT= 100010000 0−10100001 = 0−10100000 =Z所以有:
Z = − d i a g ( 1 , 1 , 0 ) W = d i a g ( 1 , 1 , 0 ) W T (1) \mathbf{Z} = -diag(1, 1, 0) \mathbf{W} = diag(1, 1, 0) \mathbf{W}^T \tag{1} Z=−diag(1,1,0)W=diag(1,1,0)WT(1)
先证引理 1:已知 T = ( t x t y t z ) T= \begin{pmatrix} t_x \\ t_y \\ t_z \end{pmatrix} T= txtytz ,则有 [ T × ] = ( 0 − t z t y t z 0 − t x − t y t x 1 ) = k U Z U T [T_{\times}] = \begin{pmatrix} 0 & -t_z & t_y \\ t_z & 0 & -t_x \\ -t_y & t_x & 1 \end{pmatrix} = k \mathbf{U} \mathbf{Z} \mathbf{U}^T [T×]= 0tz−ty−tz0txty−tx1 =kUZUT,其中 U \mathbf{U} U 是单位正交矩阵。
证明如下:
取 u 3 = T ∥ T ∥ = ( t x ∥ T ∥ , t y ∥ T ∥ , t z ∥ T ∥ ) T u_3 = \dfrac{T}{\|T\|} = \begin{pmatrix} \dfrac{t_x}{\|T\|}, \dfrac{t_y}{\|T\|}, \dfrac{t_z}{\|T\|}\end{pmatrix}^T u3=∥T∥T=(∥T∥tx,∥T∥ty,∥T∥tz)T 和 u 1 = ( − t y t x 2 + t y 2 , t x t x 2 + t y 2 , 0 ) T u_1 = \begin{pmatrix} -\dfrac{t_y}{\sqrt{t_x^2 + t_y^2}}, \dfrac{t_x}{\sqrt{t_x^2 + t_y^2}}, 0 \end{pmatrix}^T u1=(−tx2+ty2ty,tx2+ty2tx,0)T,可知 u 3 ⋅ u 1 = 0 u_3 \cdot u_1 = 0 u3⋅u1=0。
令 u 2 = u 3 × u 1 ∥ u 3 × u 1 ∥ u_2 = \dfrac{u_3 \times u_1}{\|u_3 \times u_1\|} u2=∥u3×u1∥u3×u1,取 U = [ u 1 , u 2 , u 3 ] \mathbf{U} = [u_1, u_2, u_3] U=[u1,u2,u3],则有 U T U = I \mathbf{U}^T \mathbf{U} = \mathbf{I} UTU=I,其中 u 1 , u 2 , u 3 u_1, u_2, u_3 u1,u2,u3 为标准正交基。
可知 U T T = ( u 1 T u 2 T u 3 T ) T = ( u 1 T T u 2 T T u 3 T T ) = ( u 1 ⋅ T u 2 ⋅ T u 3 ⋅ T ) = ( u 1 ⋅ T ) ( 1 0 0 ) + ( u 2 ⋅ T ) ( 0 1 0 ) + ( u 3 ⋅ T ) ( 0 0 1 ) \mathbf{U}^T T = \begin{pmatrix} u_1^T \\ u_2^T \\ u_3^T \end{pmatrix} T = \begin{pmatrix} u_1^T T \\ u_2^T T \\ u_3^T T \end{pmatrix} = \begin{pmatrix} u_1 \cdot T \\ u_2 \cdot T \\ u_3 \cdot T \end{pmatrix} = (u_1 \cdot T)\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} + (u_2 \cdot T)\begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix} + (u_3 \cdot T)\begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix} UTT= u1Tu2Tu3T T= u1TTu2TTu3TT = u1⋅Tu2⋅Tu3⋅T =(u1⋅T) 100 +(u2⋅T) 010 +(u3⋅T) 001 ,又因为
{ u 3 ⋅ T = T ∥ T ∥ ⋅ T = ∥ T ∥ u 1 ⋅ T = − t x t y t x 2 + t y 2 + t x t y t x 2 + t y 2 = 0 u 2 ⋅ T = u 2 ⋅ ( ∥ T ∥ u 3 ) = 0 \begin{cases} u_3 \cdot T = \dfrac{T}{\|T\|} \cdot T = \|T\| \\ u_1 \cdot T = -\dfrac{t_x t_y}{\sqrt{t_x^2 + t_y^2}} + \dfrac{t_x t_y}{\sqrt{t_x^2 + t_y^2}} = 0 \\ u_2 \cdot T = u_2 \cdot (\|T\|u_3) = 0 \end{cases} ⎩ ⎨ ⎧u3⋅T=∥T∥T⋅T=∥T∥u1⋅T=−tx2+ty2txty+tx2+ty2txty=0u2⋅T=u2⋅(∥T∥u3)=0所以 U T T = ∥ T ∥ ( 0 0 1 ) \mathbf{U}^T T =\|T\| \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix} UTT=∥T∥ 001 。
对任意向量 ω = ( ω x ω y ω z ) \omega = \begin{pmatrix} \omega_x \\ \omega_y \\ \omega_z \end{pmatrix} ω= ωxωyωz 有:
U T [ T × ] U ω = U T [ T × ( U ω ) ] = ( U T T ) × ( U T U ω ) = ∥ T ∥ ( 0 0 1 ) × ω = ∥ T ∥ ( 0 1 0 − 1 0 0 0 0 0 ) ω \begin{align*} \mathbf{U}^T [T_{\times}] \mathbf{U} \omega &= \mathbf{U}^T [T \times (\mathbf{U} \omega)] = (\mathbf{U}^TT)\times (\mathbf{U}^T \mathbf{U} \omega) = \|T\| \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix} \times \omega \\ &= \|T\| \begin{pmatrix} 0 & 1 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \omega \end{align*} UT[T×]Uω=UT[T×(Uω)]=(UTT)×(UTUω)=∥T∥ 001 ×ω=∥T∥ 0−10100000 ω即有:
[ U T [ T × ] U − ∥ T ∥ ( 0 1 0 − 1 0 0 0 0 0 ) ] ω = 0 \left[\mathbf{U}^T [T_{\times}] \mathbf{U} - \|T\| \begin{pmatrix} 0 & 1 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \right] \omega = 0 UT[T×]U−∥T∥ 0−10100000 ω=0由于 ω \omega ω 是任意向量,所以有:
U T [ T × ] U = ∥ T ∥ ( 0 1 0 − 1 0 0 0 0 0 ) \mathbf{U}^T [T_{\times}] \mathbf{U} = \|T\| \begin{pmatrix} 0 & 1 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} UT[T×]U=∥T∥ 0−10100000 令 σ = ∥ T ∥ \sigma = \|T\| σ=∥T∥,即有:
[ T × ] = σ U Z U T (2) [T_{\times}] = \sigma \mathbf{U} \mathbf{Z} \mathbf{U}^T \tag{2} [T×]=σUZUT(2)证毕。
结合方程 ( 1 ) (1) (1) 和引理 ( 1 ) (1) (1) 有:
[ T × ] = σ U Z U T = − σ U d i a g ( 1 , 1 , 0 ) W U T = σ U d i a g ( 1 , 1 , 0 ) W T U T (3) \begin{align*} [T_{\times}] &= \sigma \mathbf{U} \mathbf{Z} \mathbf{U}^T \\ &= -\sigma \mathbf{U} diag(1, 1, 0) \mathbf{W} \mathbf{U}^T \\ &= \sigma \mathbf{U} diag(1, 1, 0) \mathbf{W}^T \mathbf{U}^T \end{align*} \tag{3} [T×]=σUZUT=−σUdiag(1,1,0)WUT=σUdiag(1,1,0)WTUT(3)取 [ T × ] = − σ U d i a g ( 1 , 1 , 0 ) W U T [T_{\times}] = -\sigma \mathbf{U} diag(1, 1, 0) \mathbf{W} \mathbf{U}^T [T×]=−σUdiag(1,1,0)WUT,则有:
E = [ T × ] R = U d i a g ( − σ , − σ , 0 ) W U T R \mathbf{E} = [T_{\times}] \mathbf{R} = \mathbf{U} diag(-\sigma, -\sigma, 0) \mathbf{W} \mathbf{U}^T \mathbf{R} E=[T×]R=Udiag(−σ,−σ,0)WUTR对 E \mathbf{E} E 做奇异值分解有: E = U Σ V T \mathbf{E} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^T E=UΣVT,对比上述两个式子有:
V T = W U T R \mathbf{V}^T = \mathbf{W} \mathbf{U}^T \mathbf{R} VT=WUTR考虑到 W \mathbf{W} W 和 U \mathbf{U} U 都是单位正交阵,即有:
R = U W T V T (4) \mathbf{R} = \mathbf{U} \mathbf{W}^T \mathbf{V}^T \tag{4} R=UWTVT(4)取 [ T × ] = σ U d i a g ( 1 , 1 , 0 ) W T U T [T_{\times}] = \sigma \mathbf{U} diag(1, 1, 0) \mathbf{W}^T \mathbf{U}^T [T×]=σUdiag(1,1,0)WTUT,同理可得:
R = U W V T (5) \mathbf{R} = \mathbf{U} \mathbf{W} \mathbf{V}^T \tag{5} R=UWVT(5)公式 ( 4 ) (4) (4) 和 ( 5 ) (5) (5) 只保证了矩阵 U W V T \mathbf{U} \mathbf{W} \mathbf{V}^T UWVT 和 U W T V T \mathbf{U} \mathbf{W}^T \mathbf{V}^T UWTVT 是正交的,但是旋转矩阵还要满足 d e t ( R ) = 1 det(\mathbf{R}) = 1 det(R)=1,因此为确保行列式的值为正,将公式 ( 4 ) (4) (4) 和 ( 5 ) (5) (5) 修改成如下:
R = d e t ( U W V T ) U W V T o r R = d e t ( U W T V T ) U W T V T (6) \mathbf{R} = det(\mathbf{U} \mathbf{W} \mathbf{V}^T) \mathbf{U} \mathbf{W} \mathbf{V}^T \quad or \quad \mathbf{R} = det(\mathbf{U} \mathbf{W}^T \mathbf{V}^T) \mathbf{U} \mathbf{W}^T \mathbf{V}^T \tag{6} R=det(UWVT)UWVTorR=det(UWTVT)UWTVT(6)此外,可知:
T × T = [ T × ] T = σ U Z U T T = σ ( u 1 , u 2 , u 3 ) ( 0 1 0 − 1 0 0 0 0 0 ) ( u 1 T u 2 T u 3 T ) T = σ ( − u 2 , u 1 , 0 ) ( u 1 T u 2 T u 3 T ) T = σ ( − u 2 u 1 T + u 1 u 2 T ) T = 0 \begin{align*} T \times T = [T_{\times}] T = \sigma \mathbf{U} \mathbf{Z} \mathbf{U}^T T &= \sigma \left( u_1, u_2, u_3 \right) \begin{pmatrix} 0 & 1 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} u_1^T \\ u_2^T \\ u_3^T \end{pmatrix} T \\ &= \sigma \left( -u_2, u_1, 0 \right) \begin{pmatrix} u_1^T \\ u_2^T \\ u_3^T \end{pmatrix} T \\ &= \sigma \left( -u_2 u_1^T + u_1u_2^T \right) T = 0 \end{align*} T×T=[T×]T=σUZUTT=σ(u1,u2,u3) 0−10100000 u1Tu2Tu3T T=σ(−u2,u1,0) u1Tu2Tu3T T=σ(−u2u1T+u1u2T)T=0即有: − u 2 ( u 1 T T ) + u 1 ( u 2 T T ) = 0 -u_2 (u_1^T T) + u_1 (u_2^T T) = 0 −u2(u1TT)+u1(u2TT)=0,因为 u 1 T u 3 = 0 u_1^T u_3 = 0 u1Tu3=0 且 u 2 T u 3 = 0 u_2^T u_3 = 0 u2Tu3=0,不难得出:
T = ± u 3 (7) T = ±u_3 \tag{7} T=±u3(7)其中, u 3 u_3 u3 为矩阵 U \mathbf{U} U 的第三列。上述推导总结如下图:
根据公式 ( 6 ) (6) (6) 和公式 ( 7 ) (7) (7), R \mathbf{R} R 和 T T T 共有四组可能的解。为了确定正确的解,通常通过三角化一个或多个空间点进行验证。其中,只有确保这些点在各自相机坐标系下深度值 z z z 为正(或者 z z z 坐标均为正的个数最多)的那组解,才是有效的解,如下图所示:
2.3. 欧式结构恢复的歧义性
欧式结构恢复无法恢复场景的绝对尺度,仅能重建与真实场景相似的三维结构(尺度比例一致但大小未知)。若需恢复实际尺度,必须依赖外部先验信息
(例如,场景中包含已知尺寸的物体,如标准高度的行人、车辆长度、建筑层高等,通过比对重建模型与真实尺寸的比例计算尺度因子)或辅助传感器数据
(例如 IMU、深度相机)。这一特性使得单目三维重建在无人驾驶、机器人导航等需精确尺度的场景中必须与其他传感器协同使用。
如下图所示,在欧式结构恢复中,重建的三维场景与真实场景之间存在相似变换歧义
:恢复的结构不仅在尺寸上与真实场景相差一个缩放因子 s s s,还因坐标系定义差异存在旋转矩阵 R \mathbf{R} R 和平移向量 T T T 的偏移。整体而言,重建场景 P r e c P_{rec} Prec 与真实场景 P r e a l P_{real} Preal 的关系可表示为:
P r e c = s ⋅ R ⋅ P r e a l + T P_{rec} = s \cdot \mathbf{R} \cdot P_{real} + T Prec=s⋅R⋅Preal+T这一变换的根源在于算法将第一帧相机的坐标系设为世界坐标系,而真实场景的绝对坐标系可能与相机坐标系存在任意刚体变换(如旋转和平移)。此简化设定虽方便计算,但导致重建结果缺乏绝对空间基准。
3. 仿射结构恢复
3.1. 基本概念
参考博客:三维重建 —— 1. 摄像机几何,仿射相机(如弱透视投影)假设场景深度变化远小于物距,忽略透视效应,导致投影方程线性化,线性模型丢失了深度信息,使解空间扩大。
弱透视投影摄像机的投影矩阵 M = [ m 1 m 2 m 3 ] = [ A 2 × 3 b 2 × 1 0 1 × 3 1 ] = [ m 1 m 2 0 , 0 , 0 , 1 ] M = \begin{bmatrix} m_1 \\ m_2 \\ m_3 \end{bmatrix} = \begin{bmatrix} A_{2\times3} & b_{2\times1} \\ 0_{1\times3} & 1 \end{bmatrix} = \begin{bmatrix} m_1 \\ m_2 \\ 0,0,0,1 \end{bmatrix} M= m1m2m3 =[A2×301×3b2×11]= m1m20,0,0,1 。假设三维点 P P P 的坐标为 X = [ x y z ] X = \begin{bmatrix} x \\ y \\ z \end{bmatrix} X= xyz ,其在图像平面上对应的像素点 p p p 的坐标为 x = [ u v ] x = \begin{bmatrix} u \\ v \end{bmatrix} x=[uv],则有:
x = [ A , b ] [ X 1 ] = A X + b x = [\mathbf{A}, b] \begin{bmatrix} X \\ 1 \end{bmatrix} = \mathbf{A} X + b x=[A,b][X1]=AX+b仿射结构恢复问题的数学模型如下图所示:
3.2. 基于因式分解的仿射结构恢复
基于因式分解的仿射结构恢复分为两个步骤:数据中心化
和因式分解获得运动与结构
。
首先,我们来介绍一下数据中心化,如下图所示:
上图中, x i j x_{ij} xij 表示第 i i i 张图像中的第 j j j 个像素点,该像素点对应于场景中的三维点 X j X_j Xj。
下面我们来介绍使用因式分解来获得运动与结构。假设一共有 m m m 张图像和 n n n 个三维点,根据 x ^ i j = A i X ^ j = A i X j \hat{x}_{ij} = \mathbf{A}_i \hat{X}_j = \mathbf{A}_i X_j x^ij=AiX^j=AiXj,可以列出数据中心化的矩阵方程,如下图所示:
由矩阵 S \mathbf{S} S 和 M \mathbf{M} M 的秩都为 3 可知, D = M S \mathbf{D} = \mathbf{M} \mathbf{S} D=MS 的秩也为 3,因式分解方法如下图所示:
3.3. 仿射结构恢复的歧义性与局限性
仿射结构恢复歧义是指在利用仿射相机模型从多视图图像中重建三维结构时,解的不唯一性。这种歧义表现为恢复的场景结构与真实结构之间存在一个仿射变换的差异,即旋转、平移、缩放的组合。
对于任意可逆矩阵 H ∈ R 3 × 3 \mathbf{H} \in \Bbb{R}^{3 \times 3} H∈R3×3,有 D = M S = ( M H ) ( H − 1 S ) = M ∗ S ∗ \mathbf{D} = \mathbf{M} \mathbf{S} = (\mathbf{M} \mathbf{H})(\mathbf{H}^{-1} \mathbf{S}) = \mathbf{M}^* \mathbf{S}^* D=MS=(MH)(H−1S)=M∗S∗,如下图所示:
歧义对重建的影响包括:
几何属性失真
:仿射变换不保持欧氏空间中的角度、距离和正交性,导致恢复的结构可能发生扭曲(如正方形变为平行四边形)尺度与方向不确定性
:恢复的场景可能被任意缩放或旋转,无法确定绝对尺度或朝向(如建筑物高度的真实值无法直接获取)
仿射结构恢复适合远距离、深度变化小的物体重建(如航拍建筑物),其局限性包括:
- 无法处理大深度变化场景(需透视模型)
- 依赖可见性假设(所有点需在所有视图可见)