slam 系列—多视图几何

162 阅读3分钟

持续创作,加速成长!这是我参与「掘金日新计划 · 10 月更文挑战」的第7天,点击查看活动详情

运动恢复结构 通过 3 维场景的多张图像,恢复出该场景的 3 维结构信息,以及每张图像对应摄像机参数。

  • 欧式结构恢复
  • 仿射结构恢复
  • 透视结构恢复

这 3 种任务区别在于摄像机理解

欧式结构恢复问题

已知 nn 个 3 维点 Xjj{1,,n}X_j j\in \{1,\cdots,n\}mm 张图像中的对应点像素坐标 xijx_{ij} 也就是 3 维点 XjX_j 在第 ii 个像平面上映射的点,也就是在 m 张图像中的对应点的像素坐标为 xijx_{ij}

  • 这里 m 张图像对应的摄像机的内参数矩阵 Ki(i=1,,m)K_i(i=1,\cdots,m) 是已知的
xij=MiXj=Ki[RiTi]Xjx_{ij} = M_iX_j = K_i\begin{bmatrix}R_i & T_i \end{bmatrix}X_j

其中 i=1,,m;j=1,,ni = 1,\cdots,m; j = 1,\cdots,n

也就是求解

  • nn 个 3 维点 Xj(j=1,,n)X_j(j=1,\cdots,n) 的坐标
  • mm 个摄像机的外参 RiR_iTiT_i (i=1,,m)(i=1,\cdots,m)

三角化

线性法

x1j=M1Xj=K1[I0]Xjx2j=M2Xj=K2[RT]Xjx_{1j} = M_1\blue{X_j} = K_1\begin{bmatrix}I & 0\end{bmatrix}\blue{X_j}\\ x_{2j} = M_2\blue{X_j} = K_2\begin{bmatrix}R & T\end{bmatrix}\blue{X_j}\\

非线性法

pj=arg minp(d(p1j,M1Xj)+d(p2j,M2Xj))p^{*}_j = \argmin_p \left( d(p_{1j},M_1\blue{X_j}) + d(p_{2j},M_2\blue{X_j})\right)
F=KT[T×]RK1=KTEK1\bold{F} = {K^{\prime}}^{-T} [T_{\times}]RK^{-1} = {K^{\prime}}^{-T} \bold{E} K^{-1}
E=T×R=[T×]R\bold{E} = T \times R = [T_{\times}]R
  • 求解基础矩阵 F\bold{F}
  • 求解本质矩阵 E=K2TFK1\bold{E} = K_2^T\bold{F}K_1
  • 分解本质矩阵 E\bold{E} 分解为 RRTT
  • 三角化

本质矩阵 EE 分解为 RRTT

E=[T×]RE = [T_{\times}]R

这里 [T×][T_{\times}] 是一个矩阵,是一个反对称矩阵,这样矩阵具有一定性质,接下来就利用其性质进行分解

a×b=[0azayaz0axayax0][bxbybz]=[a×]ba \times b = \begin{bmatrix} 0 & -a_z & a_y\\ a_z & 0 & -a_x\\ -a_y & a_x & 0\\ \end{bmatrix} \begin{bmatrix} b_x\\ b_y\\ b_z \end{bmatrix} = \begin{bmatrix} a_{\times} \end{bmatrix}b

我们向量叉乘可以写成一个矩阵和一个向量相乘的结果,[a×]\begin{bmatrix} a_{\times} \end{bmatrix} 是一个反对称矩阵

设A为n维方阵,若有A'=-A,则称矩阵A为反对称矩阵。对于反对称矩阵,它的主对角线上的元素全为零,而位于主对角线两侧对称的元素反号。反对称矩阵具有很多良好的性质,如若A为反对称矩阵,则A',λA均为反对称矩阵;若A,B均为反对称矩阵,则A±B也为反对称矩阵;设A为反对称矩阵,B为对称矩阵,则AB-BA为对称矩阵;奇数阶反对称矩阵的行列式必为0。反对称矩阵的特征值是0或纯虚数,并且对应于纯虚数的特征向量的实部和虚部形成的实向量等长且互相正交

这里有点值得说明就是 pTFp=0{p^{\prime}}^TFp = 0 这里是无法确定 FF 的符号和尺度的,所以 F-FkF-kF 都是满足上面式子,从而得出 E=KTFKE = {K^{\prime}}^TFK 中的 EE 对于我们也是无法确定尺度和符号的。

定义两个矩阵

W=[010100001]Z=[010100001]W = \begin{bmatrix} 0 & -1 & 0\\ 1 & 0 & 0\\ 0 & 0 & 1\\ \end{bmatrix} Z = \begin{bmatrix} 0 & -1 & 0\\ 1 & 0 & 0\\ 0 & 0 & 1\\ \end{bmatrix}
Z=diag(1,1,0)W=diag(1,1,0)WTZ = diag(1,1,0)W = diag(1,1,0)W^T

[T×]\begin{bmatrix} T_{\times} \end{bmatrix} 可以写成 [T×]=kUZUT\begin{bmatrix} T_{\times} \end{bmatrix} = kUZU^T 其中 UU 是单位正交矩阵

[T×]=UZUT=Udiag(1,1,0)WUT=Udiag(1,1,0)WTUT\begin{bmatrix} T_{\times} \end{bmatrix} = UZU^T = Udiag(1,1,0)WU^T = Udiag(1,1,0)W^TU^T
E=[T×]RE=Udiag(1,1,0)WUTRE=Udiag(1,1,0)(WUTR)E = \begin{bmatrix} T_{\times} \end{bmatrix} R\\ E = Udiag(1,1,0)WU^TR\\ E = Udiag(1,1,0)(WU^TR)\\

对于 EE 进行 SVD 分解 E=Udiag(1,1,0)VTE = Udiag(1,1,0)V^T

反对称矩阵的逆是他的转置

VT=WUTRR=UWTVTV^T = WU^TR\\ R = UW^TV^T

同时我们也可以得到 R=UWVTR = UWV^T

通过上面推导我们得到 RRUWVTUWV^TUWTVTUW^TV^T 两种解中的一种解。

注意: E 的这个因式分解只保证了矩阵 UWVTUWV^TUWTVTUW^TV^T 是正交的。其为旋转矩阵还需要保证行列式的值为正 R=det(UWVT)UWVTR = det(UWV^T)UWV^Tdet(UWTVT)UWTVTdet(UW^TV^T)UW^TV^T

T×T=[T×]T=UZUTT=0T \times T = \begin{bmatrix} T_{\times} \end{bmatrix}T = UZU^TT = 0

从而得到 T=±u3T = \pm u_3 也就是 UU 的第三列

R=(det(UWVT))UWVT  T=u3R=(det(UWVT))UWVT  T=u3R=(det(UWTVT))UWTVT  T=u3R=(det(UWTVT))UWTVT  T=u3R = (det(UWV^T))UWV^T \; T=u_3\\ R = (det(UWV^T))UWV^T \; T=-u_3\\ R = (det(UW^TV^T))UW^TV^T \; T=u_3\\ R = (det(UW^TV^T))UW^TV^T \; T=-u_3\\

到现在为止有 4 种情况,暂时还无法判断哪种情况是正确的,可以通过实际计算点然后观察从以上 4 种情况中选择正确的一组解。

找到一对解来 3 为重建,对多个点进行三角化,选择在两个摄像机坐标系下 z 坐标均在正的个数最多的那组 RRRR

从图像不能估计出绝对尺度,需要一些对照的先验信息,和真实解差异还旋转、平移和缩放,也就是差一个相似变换。这种重构称为度量重构,恢复场景与真实场景之间差相似变换

仿射结构摄像机

  • 内、外参数未知 弱透视投影摄像机
  • 仿射结构摄像机

弱透视投影

x=fzxx^{\prime} = \frac{f^{\prime}}{z}x
y=fzyy^{\prime} = \frac{f^{\prime}}{z}y
x=fz0xy=fz0yx^{\prime} = \frac{f^{\prime}}{z^0}x\\ y^{\prime} = \frac{f^{\prime}}{z^0}y

就是等比放大,x 和 y 不再与深度有关系,就是放大关系,这就是仿射摄像机,