【多视图几何】TUM 课程 第5章 双视图重建：线性方法

课程的 YouTube 地址为：https://www.youtube.com/playlist?list=PLTBdjV_4f-EJn6udZ34tht9EVIW7lbeo4 。视频评论区可以找到课程所使用课件与练习题的下载地址。

课程第5章介绍了如何从两张影像重建同名点的三维空间坐标，其中涉及到本征矩阵、基础矩阵、八点法恢复、四点法恢复。

1. 重建问题

重建问题是指通过影像重建场景中三维点的空间位置，在重建场景中三维点的空间位置之前需要获得影像的空间位置与姿态。故，重建问题涉及两方面的问题：1. 重建各影像的空间位置与姿态；2. 重建场景三维点的空间位置。

1.1 假设

前一章讲述了如何从影像中获得同名点，于是本章的讨论基于以下假设：

1. 已获得各影像之间的特征点对应关系；
2. 场景中任何物体在拍摄过程中不存在移动；
3. 相机已经过标定，影像的内参已知。

1.2 几何模型

在本章中用的几何模型如下图所示（图片来源于 TUM 课程 Slides）：

\(o_1, o_2\) 分别是左右两张影像的成像中心，它们俩的连线 $ (o_1, o_2) $ 被称为基线（Baseline）。

直线 $ (o_1, o_2) $ 分别与左右影像交于 \(e_1, e_2\) 两点，$ e_1, e_2 $ 被称为极点（Epipole），极点可以看做是影像成像中心在其他影像上的投影。

空间中一点 \(X\) 与在左右影像上的投影分别为 \(x_1, x_2\)，\(o_1, o_2, X\) 三点确定的平面被称为核面（Epipolar Plane）。

连接 \(x_1, e_1\) 得到直线 \(l_1\)，\(l_1\) 被称为 \(x_2\) 在左影像上对应的核线（Epipolar Line），同理 \(l_2\) 也是 \(x_1\) 在右影像上对应的核线。

2. 极线约束

像平面坐标 \(\mathbf{x}\) 与三维空间点在世界坐标系下的坐标 \(\mathbf{X_0}\) 的关系如下：

\[\lambda \mathbf{x} = K_f \Pi_0 g \mathbf{X_0} \]

对于左影像，可以将世界坐标系设置为与左影像的相机坐标系重合（\(g\) 为单位阵），于是上式简化为：

\[\lambda \mathbf{x} = K_f \mathbf{X} \]

其中 \(\mathbf{X}\) 是三维空间点在左影像相机坐标系下的非齐次坐标，$ \mathbf{X} = \Pi_0 \mathbf{X_0} $。

方便起见，将像平面坐标 \(\mathbf{x}\) 所使用的单位长度设置为影像的焦距 \(f\)，即每一个坐标值除以 \(f\)，\(K_f\) 是一个对角矩阵，所以这么做相当于 \(\mathbf{x} = K_f^{-1}\mathbf{x}\)，于是上式可以继续化简为：

\[\lambda \mathbf{x} = \mathbf{X} \]

2.1 约束推导

对左右影像建立投影关系：

\[\lambda_1 \mathbf{x_1} = \mathbf{X}, \quad \lambda_2 \mathbf{x_2} = R \mathbf{X} + T \]

两式整合

\[\lambda_2 \mathbf{x_2} = R (\lambda_1\mathbf{x_1}) + T \]

消除式子最右边单独的 \(T\)，因为解线性方程最好整个式子只有一项（齐次）或两项（非齐次），采用的方法是用 \(T\) 的正交补 \(\hat T\) 左乘式子两边：

\[\lambda_2 \hat T \mathbf{x_2} = \lambda_1 \hat T R \mathbf{x_1} \]

现在，消除式子中的 \(\lambda_2\) （\(\lambda_2\) 在 \(\lambda_1, R, T\) 已知的情况下是可以计算出来的，所以现在专心计算 \(\lambda_1, R, T\)，不要理会 \(\lambda_2\)），采用的方法是 \(\mathbf{x_2}^T\) 右乘式子：

\[\mathbf{x_2}^T \hat T R \mathbf{x_1} = 0 \]

注意 \(\hat T \mathbf{x_2} = T \times \mathbf{x_2}\) 是两个向量的叉乘，结果垂直于 \(T\) 与 \(\mathbf{x_2}\)，所以 $\mathbf{x_2}^T \hat T \mathbf{x_2} = 0 $。

上面这个式子称作极线约束（Epipolar Constraint）。

2.2 几何意义

极线约束可以写成

\[\mathbf{x_2}^T (T \times (R\mathbf{x_1})) = 0 \]

这是一个三重积（Triple Product）的形式，涉及到三个向量 $ \mathbf{x_2}, T, R\mathbf{x_1} $。三重积的几何意义是体积，即上式表示以 \(\mathbf{x_2}, T, R\mathbf{x_1}\) 为边的平行六面体体积为0，即三个向量共面。

\(R\mathbf{x_1}\) 指左影像上的点旋转到右影像上的坐标。

由 \(\lambda \mathbf{x} = \mathbf{X}\)，整个三个向量共面，这个平面就是核面。

2.3 本征矩阵

方便起见，将核线约束写作

\[\mathbf{x_2}^T E \mathbf{x_1} = 0 \]

其中 $$ E = \hat T R \in \mathbb{R}^3 $$ 是本征矩阵（Essential Matrix）。同时，核线约束（Epipolar Constraint）也被称作 Essential Constraint 或 Bilinear Constraint。

本征矩阵所处的空间可以用以下集合表示：

\[\epsilon \equiv \left\{ \hat T R \left.\right| R \in SO(3), T \in \mathbb{R}^3 \right\} \subset \mathbb{R}^{3\times3} \]

Huang & Faugeras, 1981 表明了作为一个本征矩阵，它的 SVD 分解是如下形式：

\[E = U \Sigma V^T, \quad \Sigma = diag\{\sigma, \sigma, 0\} \]

其中 \(\sigma > 0, U, V \in SO(3)\)。（\(det(U) = +1, det(V) = +1\)）

本征矩阵只有两个特征值，说明它的秩（rank）为2。\(E = \hat T R\)，\(\hat T\) 的秩为2，\(R\) 的秩为3，矩阵乘积的秩小于两个因子矩阵的秩。

3. 八点法

使用本征矩阵进行三维重建的过程是计算本征矩阵--计算影像相对位姿--三角测量获得三维点坐标。

第1章已经说到有下面这个转换关系：

\[u^TAv = (v \otimes u)^T A^s \]

那么极线约束 \(\mathbf{x_2}^TE\mathbf{x_1} = 0\) 就能够转换为

\[(\mathbf{x_1} \otimes \mathbf{x_2})^TE^s = 0 \]

其中

\[E^s = (e_{11}, e_{21}, e_{31}, e_{12,}, e_{22}, e_{32}, e_{13}, e_{23}, e_{33})^T \in \mathbb{R}^9 \]

令 $a \equiv \mathbf{x_1} \otimes \mathbf{x_2} $，极限约束就能转换为一个简洁的线性系统：

\[a^T E^s = 0 \]

\(E^s\) 是 \(a^T\) 的 null space，对于一个点而言，\(a^T\) 是一行向量，秩为1，无法求解。

考虑求解 \(E^s\) 需要多少个点，也就是讨论 \(E^s\) 有多少个自由度，\(E^s \in \mathbb{R}^9\)，去除一个尺度量，\(E^s\) 有八个自由度，\(E\) 的秩为8，所以需要至少8个点求解 \(E^s\)。

寻找8对同名点，得到8个 \(a^T\) 的行向量，将8个行向量上下排列形成 \(8 \times 9\) 的矩阵 \(\chi\)：

\[\chi = (a^1, a^2, a^3, \dots, a^n)^T, n = 8 \]

当 \(n \gt 8\) 时，\(\chi\) 的秩就为 9，列满秩，null space 的维度为0，无解，在这种情况下需要使用最小二乘求解。

当 \(n \lt 8\) 时，\(\chi\) 的秩小于8，null space 的维度大于1，在这种情况下求解得到的 \(E^s\) 有无穷多的解。

得到由八个点组成的矩阵 \(\chi\) 之后就能求解 \(E^s\) 了：

1. 将 \(\chi\) 进行 SVD 分解，\(\chi = U_{\chi} \Sigma_{\chi} V_{\chi}^T\)，\(V_{\chi}\) 的第9列对应着将 \(\|\chi E^s\|\) 最小化的 \(E^s\)，将 \(E^s\) 变换回 \(E\)。
2. 拿到 \(E\) 之后，继续进行 SVD 分解，因为“现实世界是有误差的”，所以分解结果是 \(E = U diag\{\sigma_1, \sigma_2, \sigma_3\}V^T\) 的形式。现在需要把这个“现实世界”的 \(E\) 投影到本征矩阵的空间中去，做法是设置 \(\sigma_1 = 1, \sigma_2 = 1, \sigma_3 = 0\)，强制降秩，然后 \(E\) 就是一个秩为2的本征矩阵了。这个时候再检查一下 \(U, V\) 的行列式值，按照 2.3 中说的，他们都应该为 \(+1\)，不行就给 \(U, V\) 加一个负号。
3. 继续上面的分解结果，计算相对位姿 \(T, R\)，$$ R = U R_Z^T(\pm{\pi \over 2}) V^T, \quad \hat T = U R_Z(\pm{\pi \over 2})\Sigma U^T $$，得到了两个 \(R\)，两个 \(T\)，组合一下，有四种情况，再从这四种情况中挑选正确的，这里就需要一些先验知识了。

4. 空间点三维重建

最初是从公式 $$ \lambda_2 \mathbf{x_2} = R (\lambda_1\mathbf{x_1}) + T $$ 推导出核线约束的，现在计算得到了相对位姿 \(R, T\) 之后用这个公式继续计算景深 \(\lambda_1, \lambda_2\) 以得到空间点的三维坐标（其实只需要计算其中之一就行了）。

整理一下，假设我们现在在两张影像上存在 \(n\) 对同名点，需要求解的也就是这 \(n\) 个空间点的三维坐标，对于每一个点都能列出下面的公式：

\[\lambda_2^j \mathbf{x_2^j} = R (\lambda_1^j\mathbf{x_1^j}) + \gamma T, \quad j = 1, \dots, n \]

其中 \(\gamma \in \mathbb{R}^{+}\) 是一个比例系数，因为在八点法求解的时候存在一个比例系数的放松。

对于一个点只需要计算 \(\lambda_1, \lambda_2\) 其中之一就足够了，现在想办法去掉 \(\lambda_2\)，当然是用 \(\mathbf{\hat {x}_2^j}\) 左乘式子的左右两侧：

\[\lambda_1^j\mathbf{\hat {x}_2^j} R \mathbf{x_1^j} + \gamma \mathbf{\hat {x}_2^j} T = 0, \quad j = 1, \dots, n \]

将这 \(n\) 个方程写成 \(n\) 个线性系统的形式：

\[\begin{bmatrix}\mathbf{\hat {x}_2^j} R \mathbf{x_1^j}, \quad \mathbf{\hat {x}_2^j} T \end{bmatrix} \begin{bmatrix} \lambda_1^j \\ \gamma \end{bmatrix} = 0, \quad j = 1, \dots, n \]

然后将这 \(n\) 个线性系统合为一个：

\[M \vec{\lambda} = 0 \]

其中

\[\vec{\lambda} \equiv (\lambda_1^1, \lambda_1^2, \dots, \lambda_1^j, \gamma)^T \in \mathbb{R}^{n+1} \]

\[M \equiv \begin{bmatrix} \mathbf{\hat {x}_2^1} R \mathbf{x_1^1} & 0 & 0 & 0 & 0 & \mathbf{\hat {x}_2^1}T \\ 0 & \mathbf{\hat {x}_2^2} R \mathbf{x_1^2} & 0 & 0 & 0 & \mathbf{\hat {x}_2^2}T \\ 0 & 0 & \ddots & 0 & 0 & \vdots \\ 0 & 0 & 0 & \mathbf{\hat {x}_2^{n-1}} R \mathbf{x_1^{n-1}} & 0 & \mathbf{\hat {x}_2^{n-1}}T \\ 0 & 0 & 0 & 0 & \mathbf{\hat {x}_2^n} R \mathbf{x_1^n} & \mathbf{\hat {x}_2^n}T \end{bmatrix} \]

这个线性系统的求法与八点法相同，将 \(M\) 进行 SVD 分解，\(\vec{\lambda}\) 取特征值最小的对应的 \(V\) 的列向量。

5. 四点法

四点法只需要使用4个点就能进行三维重建，当然这是需要建立在某些先验知识（或限制条件）之上。

四点法是假设有4对同名点在同一平面上。四点共面，这一个面可以用法向量 \(N\) 进行描述，该平面上一点在左影像相机坐标系中的坐标为 \(\mathbf{X_1}\)，左影像成像中心距离平面的距离为 \(d\)：

\[N^T \mathbf{X_1} = d \quad \Leftrightarrow \quad {1 \over d}N^T\mathbf{X_1} = 1 \]

左右影像相机坐标系对于一空间点的坐标有如下的关系：

\[\mathbf{X_2} = R \mathbf{X_1} + T = R \mathbf{X_1} + T {1 \over d}N^T\mathbf{X_1} = (R + {1 \over d}TN^T)\mathbf{X_1} \equiv H\mathbf{X_1} \]

\[H = R + {1 \over d}TN^T \in \mathbb{R}^{3\times 3} \]

将像平面坐标代入（\(\lambda \mathbf{x} = \mathbf{X}\)）

\[\lambda_2 \mathbf{x_2} = H \lambda_1 \mathbf{x_1} \]

同样的方式，去掉不需要的 \(\lambda_2\)

\[\mathbf{\hat {x_2}} H \mathbf{x_1} = 0 \]

转换为线性系统

\[a^TH^s = 0 \]

\[H^s \in \mathbb{R}^9, \quad a \equiv \mathbf{x_1} \otimes \mathbf{\hat {x_2}} \in \mathbb{R}^{9 \times 3} \]

随后可以将对应多对同名点的多个 \(a^T\) 组成 \(\chi\) 以求解 \(H^s\) ：

\[\chi = (a^1, \dots, a^n)^T \in \mathbb{R}^{3n \times 9} \]

这里需要注意 \(\mathbf{\hat {x_2}}\) 的秩为2，也就是一个 \(a^T\) 只能提供2个秩，对于 \(\chi \in \mathbb{R}^{3n \times 9}\)，它的秩 \(rank(\chi) = \min\{2n, 9\}\)。

为了达到 \(H^s\) 八个自由度的要求，需要有4个点，以使得秩满足要求。（多余4个点就需要用最小二乘了）

求解 \(H^s\) 的方法就不赘述了，与求解 \(E^s\) 的方法一致。

将 \(H^s\) 分解为 \(R, T, N\) 的方法，现在还不会，留一个坑在这里。等我来填。

6. 基础矩阵

在相机未标定的情况下，没有办法获得空间点的相机成像平面坐标系坐标 \(\mathbf{x}\)，这个时候只能获得空间点的影像像素坐标 $ \mathbf{x'} $。

可以利用极线约束推导出这种情况下的约束条件：

\[\mathbf{x_2}^T E \mathbf{x_1} = 0 \]

\[\mathbf{x_2'} = K_s \mathbf{x_2}, \quad \mathbf{x_1'} = K_s \mathbf{x_1} \]

\[K_s \equiv \begin{bmatrix} s_x & s_{\theta} & o_x \\ 0 & s_y & o_y \\ 0 & 0 & 1 \end{bmatrix} \]

极线约束可以写作

\[(K_s^{-1}\mathbf{x_2'})^TE(K_s^{-1}\mathbf{x_1'}) = 0 \]

\[\Rightarrow \mathbf{x_2'}^T (K_s^{-1})^TEK_s^{-1}\mathbf{x_1'} = 0 \]

\[\Rightarrow \mathbf{x_2'}^TF\mathbf{x_1'} = 0 \]

\[F \equiv (K_s^{-1})^TEK_s^{-1} \]

\(F\) 就是所谓的基础矩阵（Fundamental Matrix）。