基本矩阵 FF

本章讲述如何用数值方法在已知若干对应点的情况下求解基本矩阵 FF

11.1 基本知识

基本矩阵定义如下

x′TFx=0 x^{'T} F x = 0

定义 x=(x,y,1)T\bold{x} = (x,y,1)^{T}，x′=(x′,y′,1)T\bold{x^{'}} = (x^{'},y^{'},1)^{T} 。带入基本矩阵的定义并整理可以得到一个方程组，记为Af=0Af = 0。具体形式参见书 279页

解析解：AA的rank如果是8，那么直接求方程组的特解就可以。但是一般情况下有噪声，所以AA的rank有可能是9，那么就用SDV分解。可以表达成A=UDVTA=UDV^{T}，那么ff 就是VV的最后一列.

11.1.1 奇异性约束

FF是一个奇异矩阵，rank为2. 利用这个性质，我们可以先用解析解求出一个FF, 然后分解为F=UDVTF=UDV^{T},进一步写成F=Udiag(r,s,t)VTF=Udiag(r,s,t)V^{T}, 接着把tt换成0就可以了。利用解析解求FF的时候需要8组对应点，所以这个方法叫8点法

11.1.2 已知7个对应点怎么办

求解Af=0Af=0,然后选基础解系中的两个特解f1,f2f_1,f_2 （书中叫零空间，null space），然后构造一个方程αf1+(1−α)f2\alpha f_1 + (1-\alpha )f_2, 再利用FF是奇异矩阵的性质，得到det(αf1+(1−α)f2)=0det(\alpha f_1 + (1-\alpha )f_2) = 0. 解这个方程就可以得到α\alpha

11.2 正交化的8点法

把8个点的质心平移到原点,把点看成一个向量,然后把向量模长归一化.

11.3 数值最优化

优化目标: 找到一个FF, 使得∣∣Af∣∣||Af||最小,且满足 ∣∣f∣∣=1||f||=1 且 detF=0det F=0 步骤如下:

  1. 通过8点法找初始值F0F_0, 然后找到F0F_0的零向量e0e_0
  2. 记e0=eie_0 = e_i, 通过 eie_i 构造 EiE_i

    E0=[e0000e0000e0]E_0 = \left[ \begin{matrix} e_0 & 0 & 0 \\ 0 & e_0 & 0 \\ 0 & 0 & e_0 \\ \end{matrix} \right]

  3. 设 fi=Eimif_i=E_i m_i 用 595 页的算法求解优化目标
  4. 计算误差 ϵi=Afi\epsilon_i = A f_i
  5. 利用600页Levenberg–Marquardt算法改变eie_i
  6. 重复1到5步直到算法收敛

11.4 几何损失函数

本章主要介绍如何用几何损失函数来计算基本矩阵。

11.4.1 黄金法则 （Gold Standard method）

我们假设图像噪声服从高斯分布，那么我们优化下式：

Σid(xi,x^)2+d(xi′,x^′)2\Sigma_i d(x_i,\hat{x})^2 + d(x_i^{'},\hat{x}^{'})^2

xi↔xi′x_i \leftrightarrow x^{'}_i是已知对应点。x^,x^′\hat{x},\hat{x}^{'}是没有噪声的点。

优化步骤如下：

  1. 根据含有噪声的点xi,xi′x_i,x^{'}_i估计出一个初始的F^\hat{F}
  2. 根据F^\hat{F}求出e′e^{'} (回忆基本矩阵定义 F=[e′]×P′P+F=[e']_{\times}P^{'}P^{+}),然后构造 P=[I∣0],P′=[[e′]×F^∣e′]P=[I|0],P^{'}=[[e']_{\times}\hat{F}|e^{'}]
  3. 根据xi↔xi′x_i \leftrightarrow x^{'}_i 还有F^\hat{F},用三角化求出空间点X^\hat{X}
  4. 根据X^\hat{X}有以下两式成立 xi^=PX^,xi^′=P′X^\hat{x_i} = P\hat{X}, \hat{x_i}^{'} = P^{'}\hat{X}
  5. 把上两式带入优化目标

11.4.2 参数化基本矩阵

非线性优化要求把FF参数化，下面介绍几种方法。

  1. Over-parametrization. 把FF写成F=[t]×MF=[t]_{\times}M M是任意的3×33 \times 3 矩阵
  2. Epipolar parametrization. 把FF 第三列写成前两列的线性组合
  3. Both epipoles as parameters，参见书p286 11.8式

11.4.3 First-order geometric error (Sampson distance)

Sampson distance以 x′Fxi=0x^{'} F x_{i}=0为优化目标，将其写成：

(x′Fxi)2(Fxi)12+(Fxi)22+(Fxi′)12+(Fxi′)12\frac{(x^{'} F x_{i})^2}{(Fx_i)^2_1 + (Fx_i)^2_2 + (Fx^{'}_i)^2_1 + (Fx^{'}_i)^2_1}

(Fxi)j2(Fx_i)^2_j 代表 FxiFx_i的第jj个元素

Symmetric epipolar distance

以 x′Fxi=0x^{'} F x_{i}=0为优化目标，那么我们可以让x′x^{'}尽可能靠近 FxiF x_{i}.所以就有以下损失函数

Σd(x′,Fxi)2+d(x,Fxi′)2\Sigma d(x^{'},Fx_i)^2 + d(x,Fx_i^{'})^2

11.5 对各种损失函数的实验结果

本节选用3中不同的算法来进行对比：

  1. 正则化的8点法
  2. 最小化代数误差
  3. 最小化几何误差

我们直接上结论：

  1. 不要使用没有正则化的8点法，换句话说，用8点法必须要把点的坐标正则化
  2. 8点法式最快最方便的
  3. 如果追求准确，用最小化代数误差。
  4. 用Sampson distance也非常好。

11.6 使用RANSAC的算法

其步骤如下：

  1. 先找出每幅图像中的特征点(SIFT, FAST, ORB)

  2. 计算特征点之间的对应

  3. RANSAC采样：

    1. 任选7个点，用SVD解出一个FF
    2. 用FF和每一个点的坐标，计算一个距离d⊥d_{\bot}. d⊥d_{\bot}可以选择重投影误差，或者Sampson。(多说一句，重投影误差就是指两次投影之间的差别，由于是同一个三维点的两次投影，这两次投影在理论上应该是一样的，所以误差应该为0.)
    3. 找出满足d⊥<td_{\bot} < t的点，其数量记为m，然后利用这些点重新计算FF,重复第二步，直到m足够多。

  4. 利用m对点，再选一个损失函数(比如说重投影和sampson)，用Levenberg–Marquardt algorithm来优化这个损失函数

11.7 一些特殊情况下FF的计算

11.7.1 纯旋转

纯旋转情况下 F=[e′]×F = [e^{'}]_{\times} 这种情况下2对对应点就可以

11.7.2 Planar motion

Planar motion意思是一个刚体沿着与某平面平行的方向(lsl_{s})运动。在这种情况下 FF满足F=[e′]×[ls]×[e]×F=[e^{'}]_{\times} [l_{s}]_{\times} [e]_{\times}

11.7.3 相机已经标定

如果相机已经标定了，那么我们可以在图像坐标系下计算本质矩阵EE, 用$E$d代替FF. 依旧是8点法。在E已知的情况下，把EE分解成UDVTUDV^T, 其中D=diag(a,b,c)D = diag(a,b,c),然后把DD换成D^=diag((a+b)/2,(a+b/2),0)\hat{D} = diag((a+b)/2,(a+b/2),0), 于是答案就是E^=UD^V\hat{E} = U\hat{D}V, 知道了E^\hat{E}和相机内参，就可以得出两个摄像机矩阵

11.8 其他的对应方式

上文讲述的是我们用点对应的方法求解FF,本节我们介绍如何用其他对应方式来求解FF.

线段对应。线段对应不会给基本矩阵FF增加任何约束，因为线段重投影回去是平面，三维空间中平面一直会相交。

曲线对应。假设三维空间中有一条和极平面相切的曲线，该曲线投影在图像上的曲线一定和极线相切。这个性质就是一个约束：如果极线与图像中某曲线相切，则第二幅图像中对应的极线一定与某曲线相切。在这个结论里把曲线换成曲面也成立。

11.9 退化的情况

退化的情况是指若干对对应点无法唯一确定基本矩阵FF. 其具体细节如下 定义 NN 为 Af=0Af=0 中 AA的零空间，也就是该方程的特解。

补充一点基础知识：零空间的维度就是自由度。比如说解方程组Af=0Af=0的时候，方程的解里有3个自由变量，那么零空间的维数就是3（ff本身是有9个变量）

  1. dim(N) = 1，这种情况是对应点n>8n>8,基本矩阵有唯一解。
  2. dim(N) = 2，这种情况是对应点n>7n>7,基本矩阵有1或3个解，主要会发生在对应点都落在某二次曲面上
  3. dim(N) = 3，这种情况是对应点n>6n>6,基本矩阵有2个解系，主要会发生在两种情况：关于两个摄像机光心重合了，或者世界坐标中的点全部在一个平面上。

11.10 求解基本矩阵的几何解释

几何解释就是x′Fx=0x^{'}Fx = 0定义了四维空间中关于点(x,y,x′,y′)(x,y,x^{'},y^{'})的一个二次曲面

11.11 极线搜索

基本矩阵的一个重要作用就是已知第一幅图中的某一点xx，找出第二幅图中该点对应的极线，那么第二幅途中xx的对应点x′x^{'}肯定在极线上。如果我们考虑到噪声问题，那x′x^{'}应该落在极线的附近。所以本章讲述如何用基本矩阵的协方差矩阵来决定x′x^{'}的搜索区域.

我们直接上结论。xx代表一个点，对于基本矩阵FF, 记其协方差矩阵ΣF\Sigma_F, xx对应的极线 l=Fxl=Fx, ll 的均值 mm 记为 lˉ\bar{l}, 我们根据以上几项构建似然函数：

(l−m)TΣl(l−m)=k2(l-m)^T \Sigma_{l} (l-m) = k^2

kk是某一个常数。同时我们用 F2(x)F_2(x)代表 χ22\chi_2^2 分布。我们只需要把k2k^2带入χ22\chi_2^2 分布，就可以得到要搜索的点落在 C 代表的区域内的概率。C如下所示：

C=lˉlˉT−k2Σ1C=\bar{l} {\bar{l}}^T -k^2\Sigma_1

11.11.1 实验验证

本节主要验证一下上一节介绍方法的有效性。先找出一些对应点，当成ground truth. 然后给每个点加上高斯噪声，然后再计算一个FF,通过FF求出极线，并且用上一节方法计算ground truth落在cc中的概率。

11.12 图像校正

校正的意思是说把两个图像旋转到同一个平面，这样就是一对双目立体图像，可以做双目立体匹配。校正分为以下几个部分。

11.12.1 把极点映射到无穷远

如果两幅图的极点都映射到了无穷远，那么他们的极线就平行了。极点本来是 (f,0,1)(f,0,1)，无穷远处的点是 (f,0,1)(f,0,1), 我们有以下矩阵可以做到这一点：

G=[101010−1/f01]G=\left[ \begin{matrix} 1 & 0 & 1\\ 0 & 1 & 0\\ -1/f & 0 & 1 \\ \end{matrix} \right]

如果考虑任意一点x0x_0，那么GG应该变成GRTGRT. TT把x0x_0转换到原点，RR把e′e^{'}旋转到(f,0,1)(f,0,1)

11.12.2 匹配变换

在上一节我们把一个图像已经转选转了某一个角度 （矩阵HH做了这个事），使其极线平行与立体匹配的基线，接下来我们需要把另外一个图像也转一下 (矩阵H′H^{'}做这事)，使其极线也平行与基线。我们不是对第二幅图像重复上一节，而是寻找HH和H′H^{'}的关系。该关系描述如下： 如果$l$he$l^{'}$是两幅图像中对应的极线，那么有下式成立:H−Tl=H′−Tl′H^{-T}l=H'^{-T}l' (HH如果是点变换H−TH^{-T}就是线变换。) 根据该式我们有以下优化函数：

Σid(Hxi,H′xi′)2\Sigma_{i} d(Hx_i,H'x'_{i})^2

同时，HH和H′H^{'}应当满足如下关系：

H=(I+H′E′aT)H′MH=(I+H^{'}E^{'}a^{T})H^{'}M

其证明参见p306

如果我们考虑一个特殊的H′H^{'},它可以把e′e^{'}变换到无穷远点(1,0,0)T(1,0,0)^T，(上文提到的H′H^{'}没有这个性质)，那么HH和H′H^{'} 应该有如下性质：在已知F=[e′]×MF=[e^{'}]_{\times}M的条件下, H=HAH0H=H_A H_0, H0=H′MH_0=H^{'}M, HAH_A是一个仿射变换。 利用这个性质，我们做一点改写：X^′=H′Xi′\hat{X}^{'} = H^{'}X_i^{'}, Xi^=H0Xi\hat{X_i} = H_0X_i, 则优化目标可以变成

Σid(HAXi^,Xi′^)2\Sigma_i d(H_A\hat{X_i},\hat{X^{'}_i})^2

接下来我们参数化，令 Xi^=(xi^,yi^,1)\hat{X_i}=(\hat{x_i},\hat{y_i},1),Xi′^=(xi′^,yi′^,1)\hat{X^{'}_i}=(\hat{x^{'}_i},\hat{y^{'}_i},1), 以上优化目标就变成了

Σi(axi^+byi^+c−xi′^)2+(yi^−y′^)\Sigma_i (a\hat{x_i}+b\hat{y_i}+c-\hat{x^{'}_i})^{2} + (\hat{y_i}-\hat{y^{'}})

11.12.3 总结

本节对算法本身做一个总体的描述：

  1. 找出至少7对对应点。
  2. 计算基本矩阵FF和两个极点e,e′e,e^{'}
  3. 找出H′H^{'},该矩阵把e′e^{'}映射到(1,0,0)T(1,0,0)^T
  4. 优化目标函数Σd(Hxi,H′x′)\Sigma d(Hx_i,H^{'}x^{'})
  5. 利用HH和H′H^{'}分别将两幅图像转换

11.12.5 仿射校正

如果摄像机本身可以本仿射摄像机近似，那么就直接用仿射变换来校正，把上文的基本矩阵换成仿射相机的基本矩阵就可以。