单目视觉（6）：SFM之矩阵估计 (五)

SFM之矩阵估计

---

引入

从单幅图像中是不可能获得图像中像素点对应的深度的，因此至少需要多余两幅的相关图像才能通过一定的方法来估计图像中像素点对应的深度信息。由此，专门来几何角度来研究从多幅图像获取相关信息的方法称之为多视角几何学，包括利用几何的方法来获取不同图像在采集过程中，物体投影所满足的一些几何约束条件。

在SFM中，我们也是通过多幅图像来恢复物体的三维信息。因此，该问题其实也就是属于多视角问题，其中的一些参数（矩阵）需要用一些算法来进行估计。

---

双视角模型

双视角模型最理想的模型就是立体视觉了。回想立体视觉模型：

本征矩阵

在相机内参已经得知的情况下(已经标定),假定以左边的图像作为参考，点p$p$在其左边摄像机坐标系统中p$p$,在右边摄像机坐标系统中 p′$p^{\prime }$，这两个坐标系统之间的关系可以使用刚体变换[R,t]$[R,t]$来描述。
以左边摄像机坐标系为世界坐标系，则依据几何关系, O,O’,P三点共面：

OP−→−⋅[OO′−→−×O′P′−→−−]=0p⋅[t×(R⋅p′)]=0$$\begin{gathered} \overrightarrow{OP}\cdot [\overrightarrow{O{O}^{\prime }}\times \overrightarrow{O^{\prime }{P}^{\prime }}]=0 \\ p\cdot [t\times (R\cdot p^{\prime })]=0 \end{gathered}$$

上式也可以改写为：

pT[t×R]p′=pTEp′=0$$p^T[t\times R]p^{\prime }=p^{T}E{p}^{\prime }=0$$

其中E$E$称为本征矩阵。它包含了物理空间中两个摄像机相关的旋转(R)和平移信息(t)。t和R描述了一台摄像机相对于另外一台摄像机在全局坐标系中的相对位置。相对地，也就是将两台摄像机获取的同一个物体的坐标（不同的摄像机坐标系）通过这两个矩阵旋转平移矩阵进行联系起来。
本征矩阵的秩为2（不是满秩），因此理论上上述方程组有无数的可行解。本征矩阵在一个尺度上可以通过5个独立的参数定义，可以使用五点法来进行参数的估计。

基础矩阵

当相机未进行标定时，即存在：p=K−1p^$p=K^{-1}\hat{p}$,p′=K′−1p^′$p^{\prime }=K^{\prime -1}{\hat{p}}^{\prime }$。其中,p^$\hat{p}$和p^′${\hat{p}}^{\prime }$为像素坐标,K$K$和K′$K^{\prime }$为相机3×3$3\times 3$内参矩阵。代入前面的式子可以得到：

p^TK−TEK′−1p^′=p^TFp^′=0$${\hat{p}}^T{K}^{-T}E{K}^{\prime -1}{\hat{p}}^{\prime }={\hat{p}}^{T}F{\hat{p}}^{\prime }=0$$

其中F$F$称为基础矩阵。除了包含E的信息外，还包含了两个摄像机的内参数。由于F包含了这些内参数，因此它可以将同一物体在不同图片中的坐标（不同的像素坐标系）相关联起来。如果图像没有畸变。焦距归一化之后，内参数矩阵变为单位矩阵，此时基础矩阵就是本征矩阵。

基础矩阵也不是满秩的矩阵，这和本征矩阵一样。基础矩阵可以通过7个独立的参数进行定义。假定其中3×3$3\times 3$的F$F$矩阵中最后一个（右下角）一个参数为1，那么可以使用8个匹配好的点对来进行参数估计。

---

参数估计

8点算法求基本矩阵

x↔x'$x\leftrightarrow x\prime$是两幅图像的任意一对匹配点。由于每一组点的匹配提供了计算F$F$的一个线性方程，当给定至少7个点时方程就可以计算出基础矩阵的F$F$。假设匹配点对的像素坐标分别为x=(u,v,1)T，x'=(u',v',1)T$x=(u,v,1{)}^T，x\prime =(u\prime ,v\prime ,1{)}^T$，则对应的方程为

[uv1]⎡⎣⎢f11f21f31f12f22f32f13f23f33⎤⎦⎥⎡⎣⎢u′v′1⎤⎦⎥=0$$\left[\begin{array}{ccc}u& v& 1\end{array}\right]\left[\begin{array}{ccc}{f}_{11}& f_{12}& f_{13}\\ f_{21}& f_{22}& f_{23}\\ f_{31}& f_{32}& f_{33}\end{array}\right]\left[\begin{array}{c}{u}^{\prime }\\ v^{\prime }\\ 1\end{array}\right]=0$$

给定n$n$组点的集合，将每个点得到的方程组改写成向量向量相乘的形式，则有：

Af=⎡⎣⎢⎢x′1x1⋮x′nxnx′1y1⋮x′nynx′1⋮x′ny′1x1⋮y′nxny′1y1⋮y′nyny′1⋮y′nx1⋮xny1⋮yn1⋮1⎤⎦⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢f11f12f13f21f22f23f31f32f33⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥=0$$A\mathbf{f}=\left[\begin{array}{ccccccccc}{x}_1^{\prime }{x}_1& x_1^{\prime }{y}_1& x_1^{\prime }& y_1^{\prime }{x}_1& y_1^{\prime }{y}_1& y_1^{\prime }& x_1& y_1& 1\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ x_n^{\prime }{x}_n& x_n^{\prime }{y}_n& x_n^{\prime }& y_n^{\prime }{x}_n& y_n^{\prime }{y}_n& y_n^{\prime }& x_n& y_n& 1\end{array}\right]\left[\begin{array}{c}{f}_{11}\\ f_{12}\\ f_{13}\\ f_{21}\\ f_{22}\\ f_{23}\\ f_{31}\\ f_{32}\\ f_{33}\end{array}\right]=\mathbf{0}$$

根据线性代数中的知识，对于线性方程组，对于不是方阵的系数矩阵，如果矩阵中存在冗余项，最后经过化简后得到的秩不多于未知数的个数，则方程组有非零解。基础矩阵的行秩不会超过8，未知数为个数为7（9个约束中除去尺度和秩为2的约束）， 也就是说只要找到8个不相关的点对就可以可以直接用线性算法解得我们需要的基础矩阵。
如果由于点坐标存在噪声则矩阵A的秩可能大于8（此时为列秩）。这时候就需要求最小二乘解。

min∑j∥∥XjF−X′j∥∥22$$min\sum _j{\Vert \begin{array}{c}{X}_{j}F-X_j^{\prime }\end{array}\Vert }_2^2$$

这里也可以用SVD$SVD$来求解，f$\mathbf{f}$的解就是系数矩阵A$A$最小奇异值对应的奇异向量，也就是A$A$奇异值分解后A=UDVT$A=UD{V}^T$中矩阵V$V$的最后一列矢量，这是在解矢量f在约束∥f∥$\Vert \begin{array}{c}\mathbf{f}\end{array}\Vert$下取∥Af∥$\Vert \begin{array}{c}A\mathbf{f}\end{array}\Vert$最小的解。由于基本矩阵有一个重要的特点就是奇异性，F$F$矩阵的秩是2。如果基本矩阵是非奇异的，那么所计算的对极线将不重合。所以在上述算法解得基本矩阵后，会增加一个奇异性约束。最简便的方法就是修正上述算法中求得的矩阵F$F$。设最终的解为F'$F\prime$，令detF′=0下求得Frobenius范数（二范数）∥F−F′∥$\Vert \begin{array}{c}F-F^{\prime }\end{array}\Vert$最小的F'$F\prime$。这种方法的实现还是使用了SVD$SVD$分解，若F=UDVT$F=UD{V}^T$，此时的对角矩阵D=diag(r,s,t)$D=diag(r,s,t)$，满足r≥s≥t$r\ge s\ge t$，则F'=Udiag(r,s,0)VT$F\prime =Udiag(r,s,0)V^T$最小化范数∥F−F′∥$\Vert \begin{array}{c}F-F^{\prime }\end{array}\Vert$，也就是最终的解。

以上算法是解基本矩阵的基本方法，称为8点算法。

错误的匹配点对矩阵估计的影响很大，可以使用RANSAC算法来剔除误匹配，以找到基础矩阵：
（1）随机选择8个点；
（2）用这8个点估计基础矩阵F$F$；
（3）用这个F$F$进行验算，计算用F$F$验算成功的点对数n；
重复多次，找到使n$n$最大的F$F$作为基础矩阵。

5点法求本征矩阵

5点算法的大概思路就是把本征方程展开，再加上本质矩阵的充要条件。化简后得到一个高次多项式方程，该多项式实数解中的一个确定出的本质矩阵是系统的真实本质矩阵，对伪解进行剔除。该方法的具体原理可见论文David Nistér. An efficient solution to the five-point relative pose problem. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 26(6):756–770, 2004.和Henrik Stewenius. Calibrated fivepoint solver.

note: 实际中使用这些方法估计出的本质矩阵的结果并不太令人满意，尤其当匹配点对趋于共面时和平移量较小时，这些方法获得的结果完全就不能用。如果用于标定的匹配点对，是在一个平面上的话，这些方法将完全失效，因为平面是一种退化结构或者叫临界曲线（Critical surfaces）。在实际使用中，建议尽量避免匹配点对来自于同一个平面，且必须有平移量，这样估计出的本质矩阵会好很多！

旋转平移矩阵求解

对获得的本质矩阵进行使用SVD$SVD$分解而得到旋转矩阵R$R$和单位平移量t$t$，通常有四对旋转平移组合，进行cheirality check（三维点在相机前方约束），留下一组物理可实现的R,t$R,t$。即完成双目的相互位置标定。

四点法单应性矩阵

单应性矩阵描述的是两个平面之间的关系：

p=H∗p′$$p=H\ast p^{\prime }$$

note: 在实际使用中，建议尽量避免匹配点对来自于同一个平面，且必须有平移量，这样估计出的本质矩阵会好很多！

---

reference：
[1] 计算机视觉：一种现代方法, 福赛思.
[2] (学习笔记)投影与三维视觉——本征矩阵和基础矩阵
[3] 基本矩阵的基本解法之8点算法
[4] 三维重建（一）外极几何，基础矩阵及求解
[5] 基于opencv3的本质矩阵（essential matrix）估计算例
[6] 四点求解单应性矩阵