三维重建 —— 4. 三维重建基础与极几何
文章目录
- 1. 三角化
- 2. 极几何
- 2.1. 基本概念
- 2.2. 极几何特例
- 3. 本质矩阵
- 4. 基础矩阵
- 4.1. 基本概念
- 4.2. 八点法
- 4.3. 归一化八点法
- 5. 单应矩阵
- 6. 总结
课程视频链接:计算机视觉之三维重建(深入浅出SfM与SLAM核心算法)—— 4. 三维重建基础与极几何。
1. 三角化
鉴于单目相机三维重建的困难,引入多目相机。多目相机三维重建需要用到三角化,现在来介绍一下三角化的概念。
假设存在两个相机,不妨记为 O 1 O_1 O1 和 O 2 O_2 O2,且已知从相机 O 1 O_1 O1 到 O 2 O_2 O2 的旋转变换矩阵和平移向量分别为 R \mathbf{R} R 和 T T T,引出下图的问题:
如上图所示,在相机坐标系 O 1 O_1 O1 中存在射线 l l l,在相机坐标系 O 2 O_2 O2 中存在射线 l ′ l' l′,两条射线相交于三维点 P P P,则 P = l × l ′ P = l \times l' P=l×l′。值得注意的是:在计算 l × l ′ l \times l' l×l′ 时,需要先将相机坐标系 O 2 O_2 O2 下的 l ′ l' l′ 先转换到相机坐标系 O 1 O_1 O1 下,利用旋转矩阵 R \mathbf{R} R 和 T T T 即可实现坐标转换。
线性解法如下图所示:
非线性解法如下图所示:
我们将之前的三角化问题更加一般化,引出新的两个问题,如下图所示:
多视图几何的关键问题总结如下:
2. 极几何
2.1. 基本概念
极几何描述了同一场景或者物体的两个视点图像间的几何关系。现在介绍一下极几何的基本概念:
其中, O 1 O 2 O_1O_2 O1O2 为基线, e e e 和 e ′ e' e′ 是极点, l l l 和 l ′ l' l′ 是极线,极平面为点 P P P、 O 1 O_1 O1 和 O 2 O_2 O2 构成的平面。
如下图所示,假设存在另一点 Q Q Q,则会形成新的极线 q e → \overrightarrow{qe} qe 和 q ′ e ′ → \overrightarrow{q'e'} q′e′,我们可以发现 q e → \overrightarrow{qe} qe 和 p e → \overrightarrow{pe} pe 相交于 e e e, q ′ e ′ → \overrightarrow{q'e'} q′e′ 和 p ′ e ′ → \overrightarrow{p'e'} p′e′ 相交于 e ′ e' e′。事实上,相机坐标系 O 1 O_1 O1 对应的图像平面上的所有极线都会相交于点 e e e,相机坐标系 O 2 O_2 O2 对应的图像平面上的所有极线都会相交于点 e ′ e' e′。
2.2. 极几何特例
平行视图的特点是相机 O 2 O_2 O2 相对于相机 O 1 O_1 O1 只有平移(沿图像平面 u u u 轴移动),没有旋转变换。平行视图的性质如下所示:
与平行视图不同的是,前向平移是指相机 O 2 O_2 O2 沿相机 O 1 O_1 O1 的前方进行移动,如下图所示:
现在我们需要根据一个图像中的点 p p p,找到该点在另一个图像上的对应点 p ′ p' p′。我们可以通过极几何约束,将搜索范围缩小到对应的极线上,如下图所示:
3. 本质矩阵
本质矩阵对规范化相机拍摄的两个视点图像间的极几何关系进行代数描述。在规范化相机中,相机内参 α = β = 1 \alpha = \beta = 1 α=β=1, θ = 90 ° \theta = 90° θ=90°, c x = c y = 0 c_x = c_y = 0 cx=cy=0,此时的相机内参矩阵 K = ( 1 0 0 0 1 0 0 0 1 ) \mathbf{K} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} K= 100010001 ,我们规定图像平面位于相机坐标系 Z Z Z 轴的正方向。
如上图所示,规范化相机坐标系通常假设所有点位于 Z = 1 Z=1 Z=1 平面,直接映射像素坐标到该平面。因此像素坐标系上点 p = ( u v ) p = \begin{pmatrix} u \\ v \end{pmatrix} p=(uv) 在规范化相机 O 1 O_1 O1 坐标系下的空间坐标为 ( u v 1 ) \begin{pmatrix} u \\ v \\ 1 \end{pmatrix} uv1 。
因为 p ′ = R p + T p' = \mathbf{R}p + T p′=Rp+T,则 p = R − 1 ( p ′ − T ) = R T ( p ′ − T ) = R T p ′ − R T T p = \mathbf{R}^{-1}(p' - T) = \mathbf{R}^T (p' - T) = \mathbf{R}^T p' - \mathbf{R}^TT p=R−1(p′−T)=RT(p′−T)=RTp′−RTT。所以点 p ′ p' p′ 在相机 O 1 O_1 O1 坐标系下的坐标为 R T p ′ − R T T \mathbf{R}^T p' - \mathbf{R}^TT RTp′−RTT。令 p ′ = ( 0 0 0 ) p' = \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix} p′= 000 ,则可得点 O 2 O_2 O2 在相机 O 1 O_1 O1 坐标系下的坐标为 − R T T - \mathbf{R}^TT −RTT。
现在我们来计算 n ⃗ = O 1 O 2 → × O 1 p ′ → = − R T T × ( R T p ′ − R T ) = − R T T × R T p ′ \vec{n} = \overrightarrow{O_1 O_2} \times \overrightarrow{O_1 p'} = -\mathbf{R}^TT \times (\mathbf{R}^Tp' - \mathbf{R}T) = -\mathbf{R}^TT \times \mathbf{R}^Tp' n=O1O2×O1p′=−RTT×(RTp′−RT)=−RTT×RTp′,可知 n ⃗ \vec{n} n 为平面 O 1 O 2 P O_1O_2P O1O2P 的法向量,即有 [ R T T × R T p ′ ] T ⋅ p = 0 ⇒ [ R T ( T × p ′ ) ] T ⋅ p = 0 ⇒ [ R T [ T x ] p ′ ] T ⋅ p = 0 ⇒ ( p ′ ) T [ T x ] T R ⋅ p = 0 [\mathbf{R}^TT \times \mathbf{R}^Tp']^T \cdot p = 0 \Rightarrow [\mathbf{R}^T (T \times p')]^T \cdot p = 0 \Rightarrow [\mathbf{R}^T[T_{x}]p']^T \cdot p = 0 \Rightarrow (p')^T[T_x]^T\mathbf{R}\cdot p=0 [RTT×RTp′]T⋅p=0⇒[RT(T×p′)]T⋅p=0⇒[RT[Tx]p′]T⋅p=0⇒(p′)T[Tx]TR⋅p=0,注意到 [ T x ] T = − [ T x ] [T_x]^T = -[T_x] [Tx]T=−[Tx],所以有:
( p ′ ) T [ T x ] R p = 0 (p')^T [T_x] \mathbf{R} p = 0 (p′)T[Tx]Rp=0则本质矩阵 E = [ T x ] R \mathbf{E} = [T_x] \mathbf{R} E=[Tx]R,即有:
( p ′ ) T E p = 0 (p')^T \mathbf{E} p = 0 (p′)TEp=0本质矩阵的一些性质如下图所示:
现在我们来证明 l ′ = E p l' = \mathbf{E}p l′=Ep 且 l = E T p ′ l = \mathbf{E}^T p' l=ETp′。
证:
{ ( p ′ ) T l ′ = 0 ( p ′ ) T E p = 0 ⇒ l ′ = E p \begin{cases} (p')^T l' = 0 \\ (p')^T \mathbf{E} p = 0 \end{cases} \Rightarrow l' = \mathbf{E} p {(p′)Tl′=0(p′)TEp=0⇒l′=Ep { l T p = 0 ( p ′ ) T E p = 0 ⇒ l = E T p ′ \begin{cases} l^T p = 0 \\ (p')^T \mathbf{E} p = 0 \end{cases} \Rightarrow l = \mathbf{E}^T p' {lTp=0(p′)TEp=0⇒l=ETp′证毕。
再证: E e = 0 \mathbf{E} e = 0 Ee=0 且 ( e ′ ) T E = 0 (e')^T \mathbf{E} = 0 (e′)TE=0。
证:因为 ( l ′ ) T p ′ = ( E p ) T p ′ = 0 (l')^Tp' = (\mathbf{Ep})^T p' = 0 (l′)Tp′=(Ep)Tp′=0,令 p ′ = e ′ p' = e' p′=e′,则有 ( E p ) T e ′ = 0 (\mathbf{E} p)^T e' = 0 (Ep)Te′=0。当我们取不同的点 p p p,均有 ( E p ) T e ′ = p T E T e ′ = 0 (\mathbf{E} p)^T e' = p^T \mathbf{E}^T e' = 0 (Ep)Te′=pTETe′=0 恒成立,所以有 E T e ′ = 0 \mathbf{E}^T e' = 0 ETe′=0,即 ( e ′ ) T E = 0 (e')^T \mathbf{E} = 0 (e′)TE=0。
同理可得 E e = 0 \mathbf{E} e = 0 Ee=0。
证毕。
最后,我们来证明 E \mathbf{E} E 的秩为 2。
不妨设 T = ( t x t y t z ) T = \begin{pmatrix} t_x \\ t_y \\ t_z \end{pmatrix} T= txtytz 。先证明当 t x , t y , t z t_x, t_y, t_z tx,ty,tz 不全为 0 时, [ T x ] = ( 0 − t z t y t z 0 − t x − t y t x 0 ) \mathbf{[T_x]} = \begin{pmatrix} 0 & -t_z & t_y \\ t_z & 0 & -t_x \\ -t_y & t_x & 0 \\ \end{pmatrix} [Tx]= 0tz−ty−tz0txty−tx0 的秩为 2。
如果 t x , t y . t z t_x, t_y. t_z tx,ty.tz 中有一个或者两个为 0,显然可知 r a n k ( [ T x ] ) = 2 rank(\mathbf{[T_x]}) = 2 rank([Tx])=2。
当 t x , t y . t z t_x, t_y. t_z tx,ty.tz 都不为 0 时,有: ( 0 − t z t y t z 0 − t x − t y t x 0 ) ∼ ( − t y t x 0 0 − t z t y t z 0 − t x ) ∼ ( − t y t x 0 0 − t z t y 0 t x t z t y − t x ) ∼ ( − t y t x 0 0 − t z t y 0 0 0 ) \begin{pmatrix} 0 & -t_z & t_y \\ t_z & 0 & -t_x \\ -t_y & t_x & 0 \\ \end{pmatrix} \sim \begin{pmatrix} -t_y & t_x & 0 \\ 0 & -t_z & t_y \\ t_z & 0 & -t_x \\ \end{pmatrix} \sim \begin{pmatrix} -t_y & t_x & 0 \\ 0 & -t_z & t_y \\ 0 & \dfrac{t_x t_z}{t_y} & -t_x \\ \end{pmatrix} \sim \begin{pmatrix} -t_y & t_x & 0 \\ 0 & -t_z & t_y \\ 0 & 0 & 0 \\ \end{pmatrix} 0tz−ty−tz0txty−tx0 ∼ −ty0tztx−tz00ty−tx ∼ −ty00tx−tztytxtz0ty−tx ∼ −ty00tx−tz00ty0 可知 r a n k ( [ T x ] ) = 2 rank(\mathbf{[T_x]}) = 2 rank([Tx])=2。
因为旋转矩阵 R \mathbf{R} R 是正交矩阵,所以 r a n k ( R ) = 3 rank(\mathbf{R}) = 3 rank(R)=3。
一方面,Sylvester 秩不等式( r a n k ( A B ) ≥ r a n k ( A ) + r a n k ( B ) − n rank(\mathbf{A}\mathbf{B}) ≥ rank(\mathbf{A}) + rank(\mathbf{B}) - n rank(AB)≥rank(A)+rank(B)−n)可知: r a n k ( [ T x ] R ) ≥ r a n k ( [ T x ] ) + r a n k ( R ) − 3 = 2 rank(\mathbf{[T_x]} \mathbf{R}) ≥ rank(\mathbf{[T_x]}) + rank(\mathbf{R}) - 3 =2 rank([Tx]R)≥rank([Tx])+rank(R)−3=2。
另一方面,根据定理 r a n k ( A B ) ≤ min { r a n k ( A ) , r a n k ( B ) } rank(\mathbf{A} \mathbf{B}) ≤ \min \{rank(\mathbf{A}), rank(\mathbf{B})\} rank(AB)≤min{rank(A),rank(B)} 可知, r a n k ( [ T x ] R ) ≤ m i n { 2 , 3 } = 2 rank(\mathbf{[T_x]} \mathbf{R}) ≤ min \{2, 3\} = 2 rank([Tx]R)≤min{2,3}=2。
所以, r a n k ( E ) = r a n k ( [ T x ] R ) = 2 rank(\mathbf{E}) = rank(\mathbf{[T_x]} \mathbf{R}) = 2 rank(E)=rank([Tx]R)=2,即 E \mathbf{E} E 为奇异矩阵,秩为 2。
证毕。
为什么 E \mathbf{E} E 的自由度为 5,而不是 6(3 个旋转自由度 + 3 个平移自由度)呢?
根据对极约束满足 ( p ′ ) T E p = 0 (p')^T \mathbf{E} p = 0 (p′)TEp=0, E \mathbf{E} E 与 k E k \mathbf{E} kE 描述相同的几何关系,这意味着 E \mathbf{E} E 的缩放不影响对极约束,所以 E \mathbf{E} E 需要减去一个自由度,即自由度为 5。例如,我们对矩阵 [ T x ] \mathbf{[T_x]} [Tx] 乘以 1 t z \dfrac{1}{t_z} tz1,再令 t x ′ = t x t z t_x' = \dfrac{t_x}{t_z} tx′=tztx 和 t y ′ = t y t z t_y' = \dfrac{t_y}{t_z} ty′=tzty,则有 [ T x ] ′ = 1 t z [ T x ] = ( 0 − 1 t y ′ 1 0 − t x ′ − t y ′ t x ′ 0 ) \mathbf{[T_x]}' = \dfrac{1}{t_z} \mathbf{[T_x]} = \begin{pmatrix} 0 & -1 & t_y' \\ 1 & 0 & -t_x' \\ -t_y' & t_x' & 0 \\ \end{pmatrix} [Tx]′=tz1[Tx]= 01−ty′−10tx′ty′−tx′0 ,即去除了 t z t_z tz 这个自由度。
4. 基础矩阵
4.1. 基本概念
基础矩阵对一般的透视摄像机拍摄的两个视点的图像间的极几何关系进行代数描述。前面的讨论,我们是基于规范化相机进行的,如果是一般的相机,我们需要想办法变换到规范化相机上。
对于一般相机,我们有 p = K [ I 0 ] P p = \mathbf{K}[\mathbf{I} \quad0]P p=K[I0]P,其中 P P P 为三维空间点, p p p 为像素点。由于 K \mathbf{K} K 是满秩矩阵,所以 K \mathbf{K} K 可逆,则有 K − 1 p = K − 1 K [ I 0 ] P = ( 1 0 0 0 0 1 0 0 0 0 1 0 ) P \mathbf{K}^{-1} p = \mathbf{K}^{-1}\mathbf{K}[\mathbf{I} \quad 0]P = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ \end{pmatrix} P K−1p=K−1K[I0]P= 100010001000 P,如果我们令 p c = K − 1 p p_c = \mathbf{K}^{-1} p pc=K−1p 和 p c ′ = K ′ − 1 p ′ p_c' = \mathbf{K'}^{-1} p' pc′=K′−1p′,则有 p c = ( 1 0 0 0 0 1 0 0 0 0 1 0 ) P p_c = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ \end{pmatrix} P pc= 100010001000 P进一步有:
基础矩阵的性质如下图所示:
基础矩阵的作用如下图所示:
基础矩阵 F \mathbf{F} F 刻画了两幅图像的极几何关系,即相同场景在不同视图中的对应关系,广泛应用于三维重建和多视图匹配。
4.2. 八点法
基础矩阵 F \mathbf{F} F 有 7 个自由度,理论上 7 个点即可求解 F \mathbf{F} F,但计算方法比较复杂。实际中,我们可以选取 8 个点利用最小二乘法进行基础矩阵 F \mathbf{F} F 估计,如下图所示:
值得注意的是,上述最小二乘法估计出来的 F ^ \mathbf{\hat{F}} F^ 是满秩矩阵(秩为 3),而实际的基础矩阵 F \mathbf{F} F 的秩为 2,为此我们需要对 F ^ \mathbf{\hat{F}} F^ 做进一步优化,如下图所示:
八点算法的步骤总结如下:
八点法的精度较低,可能有 10 个像素以上的误差,主要原因有两个:
- 现在的图像的宽高都比较大,这导致 u u , u v uu, uv uu,uv 与 u , v u, v u,v 的数值相差较大,从而导致计算精度下降
- SVD 分解本身存在数值计算问题
4.3. 归一化八点法
5. 单应矩阵
单应矩阵反映了空间平面在两个摄像机下的投影几何,如下图所示:
单应矩阵的推导如下图所示:
单应矩阵估计方法如下图所示:
单应矩阵的性质如下图所示: