zehua

Zehua

计算机视觉：3D重建理论

零散笔记
23
0
2024-10-22

 我并不擅长计算机视觉方向，因此下面内容只涉及一些基础理论知识。由于我的立体几何、3D空间想象能力非常垃圾，这种课程一直学不明白，因此我对它的兴趣也不大，因此我也没投入过多时间和精力来进一步研究或详细解释。

原本每章内容后都有实例分析，包括代码和结果展示，但由于老师不让，我无法将这些内容直接搬运到这里。

1. 图像投影模型与3D重建理论

1.1 SLAM的概念

SLAM，即同步定位与地图构建（Simultaneous Localization and Mapping），是计算机视觉和机器人学中的一个核心概念。其本质是通过估计每个相机的位置和场景的三维点，实现对场景的重建。换句话说，SLAM让设备能够在未知环境中一边移动一边构建环境地图，同时确定自己在地图中的位置。

1.2 逆向二维图像

逆向投影是3D重建的核心任务，它要从二维图像中提取三维场景的信息。图像是三维场景经过投影后的二维表示，要恢复三维信息，需要逆转这个投影过程。这个过程的关键在于建立一个数学模型，描述三维场景如何投影到二维图像中，然后尝试逆向求解。

1.3 针孔相机模型

1.3.1 模型概述

针孔相机模型是计算机视觉中最基础也是最重要的模型之一。这个模型假设所有的光线都通过一个公共点，即光心（光学中心）。虽然这是一个理想化的假设，但它的优点在于模型简单，易于逆向计算，因此在三维重建中被广泛使用。

1.3.2 坐标系和符号约定

在针孔相机模型中，我们需要建立一个清晰的坐标系统。摄像机坐标系的原点 O_C 是光学中心，坐标为 (0, 0, 0) 。轴方向的设置遵循右手坐标系，其中 X_C 向右， Y_C 向下， Z_C 指向后方（场景深度方向）。这种设置的优势在于 Z 轴指向后方，使得物体深度为正，符合我们的直觉理解。

1.3.3 三维点的投影到归一化焦平面

在空间中，我们用点 U 来表示一个三维点，其坐标为 (U_X, U_Y, U_Z) 。为了理解投影过程，我们引入归一化焦平面的概念——与光心 O_C 距离为1的平面（ Z_C = 1 ），远处的三维点会被投影到此平面上，这个过程构成了整个投影变换的基础。

1.3.4 齐次坐标与非齐次坐标

齐次坐标的定义很简单：在原有坐标后增加一个维度（通常为1），这样做的目的是方便表示投影和变换。对于二维点 m = (m_X, m_Y)^\top ，其齐次坐标为 \bar{m} = (m_X, m_Y, 1)^\top 。齐次坐标的作用在于它能方便矩阵运算，尤其是在投影和变换过程中，让原本非线性的投影变换可以用线性的矩阵运算来表示。相对的，非齐次坐标就是我们熟悉的标准笛卡尔坐标表示法，不包含额外的维度。

1.4 摄像机的线性校准

1.4.1 从归一化焦平面到图像平面

图像平面是我们最终看到的成像平面，它使用像素坐标系，通常以图像左上角为原点，向右为 X 轴（列索引），向下为 Y 轴。我们的目的是将归一化焦平面上的点映射到实际的图像像素坐标上。这个映射通过线性变换实现：

\begin{cases} P _U = f \cdot m _X + U _0 \\ P _V = f \cdot m _Y + V _0 \end{cases}

其中焦距 f 决定了成像的放大倍数， m_X, m_Y 是归一化焦平面上的点，而 (U_0, V_0) 是光学中心在图像平面中的坐标。

1.4.2 摄像机内参矩阵

为了更好的表示上述线性变换，我们引入摄像机内参矩阵：

K = \begin{pmatrix} f & 0 & U _0 \\ 0 & f & V _0 \\ 0 & 0 & 1 \end{pmatrix}

有了这个矩阵，映射关系可以简洁地表示为：

\underline{P} = K \cdot \underline{m}

其中， \underline{P} 是图像平面中的点的齐次坐标。这个简洁的矩阵形式让我们能够方便地进行各种计算和变换。

1.4.3 逆向过程

理解了从归一化焦平面到图像平面的映射后，我们现在想得到其逆向过程。从图像平面到归一化焦平面的变换可以通过内参矩阵的逆来实现：

\underline{M} = K^{-1} \cdot \underline{P}

这个逆向过程在很多应用中都非常重要，比如当我们需要从图像坐标推断三维信息时。

1.4.4 可视锥

可视锥表示相机能够看到的空间范围。通过将图像的四个角点转换到归一化焦平面，然后连接光学中心，我们就能形成视锥。这个概念帮助我们理解相机的视野范围，在场景分析和路径规划中都有重要应用。

1.5 畸变建模与校正

1.5.1 相机畸变的来源

实际相机镜头存在光学缺陷，尤其在广角镜头中，这些缺陷导致图像出现畸变。原本应该是直线的物体在图像中变成了曲线，图像边缘出现拉伸或压缩。这些畸变在精确的视觉测量中必须考虑和校正。

1.5.2 畸变模型

为了数学化地描述畸变，我们引入畸变函数，它将归一化焦平面上的理想点经过畸变函数（从理想图像到畸变图像），得到畸变后的点。这个点位于2D实际畸变聚焦平面上：

\underline{m} _d = d(\underline{m}, k)

其中， k 是畸变参数。一个常用的例子是多项式径向畸变模型：

M _d = \left(1 + k _1 \|m\| _2^2 + k _2 \|m\| _2^4 + \dots \right) m

其中 |m|_2^2 = m_x^2 + m_y^2 表示点到图像中心的距离平方。

1.6 畸变校正的实现

1.6.1 任务描述

畸变校正的目标很明确：将畸变的实际图像校正为理想的无畸变图像。这个过程需要我们理解畸变的数学模型，并能够逆向应用这个模型。

1.6.2 实现步骤

畸变校正的实现需要仔细的步骤规划。首先要定义参数：理想的摄像机内参矩阵 K_{\text{ideal}} 、畸变的摄像机内参矩阵 K_{\text{real}} ，以及失真参数 k 。

对于每个理想图像的像素坐标，我们执行以下步骤。首先将像素坐标转换到归一化焦平面：

\underline{m} _{\text{ideal}} = K _{\text{ideal}}^{-1} \cdot \underline{P} _{\text{ideal}}

然后应用畸变函数：

\underline{m} _d = d(\underline{m} _{\text{ideal}}, k)

接着映射回实际图像坐标系：

\underline{P} _{\text{real}} = K _{\text{real}} \cdot \underline{m} _d

最后进行插值处理。由于 \underline{P}_{\text{real}} 的坐标可能为非整数，我们需要使用双线性插值等方法来生成校正后的图像。

2. 二维刚性变换和单应性

2.1 二维刚性变换

二维刚性变换是计算机视觉中的基础变换，它包括平移和旋转两种基本操作，这些变换保持物体的形状和大小不变。

2.1.1 旋转

旋转变换的核心是理解不同坐标系之间的关系。我们从一个参考系中选取一个向量，然后将其转换到另一个坐标系中：

\mathbf{U}^c = \overrightarrow{O _c U}^c \quad \mathbf{U}^w = \overrightarrow{O _w U}^w
\mathbf{R} _{wc} \underline{\mathbf{U}}^c = \mathbf{R} _{wc} \cdot \overrightarrow{O _c U}^c = \overrightarrow{O _w U}^w

这里 \mathbf{R}_{wc} 是一个正交矩阵，它描述了从坐标系 c 到坐标系 w 的旋转关系。

2.1.2 平移

平移变换相对简单，它描述了两个坐标系原点之间的位移：

\mathbf{T} _{wc} = \overrightarrow{O _w O _c}^{w}

 

2.1.3 刚性变换公式

将旋转和平移结合起来，我们得到完整的刚性变换公式：

\mathbf{U}^w = \mathbf{R} _{wc} \cdot \mathbf{U}^c + \mathbf{T} _{wc}

这个公式的正确性可以通过向量运算来证明：

\mathbf{R} _{wc} \cdot \mathbf{U}^c + \mathbf{T} _{wc} = \mathbf{R} _{wc} \cdot \overrightarrow{O _c U}^c + \overrightarrow{O _w O _c}^w = \overrightarrow{O _c U}^w + \overrightarrow{O _w O _c}^w = \overrightarrow{O _w U}^w = \mathbf{U}^w

 

2.1.4 齐次坐标

为了将旋转和平移统一在一个矩阵运算中，我们再次使用齐次坐标：

\underline{\mathbf{U}}^w = \begin{bmatrix} \mathbf{U}^w \\ 1 \end{bmatrix}

 

\mathbf{M} _{wc} = \begin{bmatrix} \mathbf{R} _{wc} & \mathbf{T} _{wc} \\ 0 & 1 \end{bmatrix} = \begin{bmatrix} r _{11} & r _{12} & r _{13} & t _{x} \\ r _{21} & r _{22} & r _{23} & t _{y} \\ r _{31} & r _{32} & r _{33} & t _{z} \\ 0 & 0 & 0 & 1 \end{bmatrix}

 

2.1.5 反变换

任何变换都有其逆变换，对于刚性变换：

\mathbf{M} _{cw} = \mathbf{M} _{wc}^{-1}

2.1.6 变换的组合性

刚性变换具有良好的组合性质，多个变换可以通过矩阵乘法组合：

\mathbf{M} _{ab} \cdot \mathbf{M} _{bc} = \mathbf{M} _{ac}

2.2 单应性

2.2.1 平面场景假设

在处理平面场景时，我们需要一个关键的假设和表达式：

\mathbf{U} _i^A = z _i^A \cdot \underline{\mathbf{m}} _{Ai}

这个方程的含义很直观：点 \mathbf{U}_i^A 可以由 \underline{\mathbf{m}}_{Ai} 来表示，怎么表示呢？只需要乘以它的深度 z_i^A 即可。因为 \underline{\mathbf{m}}_{Ai} 是单位深度的点，所以通过深度缩放就能得到实际的三维点。

2.2.2 寻找 \underline{\mathbf{m}}_{Ai} 和 \underline{\mathbf{m}}_{Bi} 之间的对应关系

光有上面的方程还不够，我们需要找到 \underline{\mathbf{m}}_{Ai} 和 \underline{\mathbf{m}}_{Bi} 之间的对应关系。通俗来讲，就是要进行坐标对应变换。

法线关键公式

为了建立这种对应关系，我们需要先回顾一个关于平面和法线的重要性质。在参考系A中，平面P的方程为： ax + by + cz + d = 0 ，其中 a, b, c 是平面法向量分量， d 是常数项，代表平面P和原点 O_A 的相对距离。在向量形式中，平面方程可以化简为：

\mathbf{n} _A^\top \mathbf{U} _i^A + d = 0

这里 \mathbf{n}_A^\top 代表向量P在参考系A中的法向量。通过这个平面方程的向量形式，我们得到了一个带有法向量的很重要的公式。

利用变量代换得到深度表达式

现在我们将 \mathbf{U} i^A = z i^A \cdot \underline{\mathbf{m}} _{A,i} 带入上式中：

\mathbf{n} _A^\top \cdot z _i^A \cdot \underline{\mathbf{m}} _{A,i} + d = 0 \quad \Rightarrow \quad z _i^A = -\dfrac{d}{\mathbf{n} _A^\top \cdot \underline{\mathbf{m}} _{A,i}}

这样，我们就把 \underline{\mathbf{m}}_{Ai} 给引进来了，其中 z_i^A = -\dfrac{d}{\mathbf{n}_A^\top \cdot \underline{\mathbf{m}}_{A,i}} 代表了深度。换句话说，我们利用 \mathbf{U}_i^A 的两个方程，将 \mathbf{U}_i^A 替换掉了，得到了 z_i^A 深度的表达式。但是光有关于 \underline{\mathbf{m}}_{Ai} 的方程还不够，我们还需要从 \underline{\mathbf{m}}_{Bi} 入手。

 

 

 

接下来我们找B坐标系下的点 \underline{\mathbf{m}}_{Bi}

我们从刚性变换公式 \mathbf{U}^w = \mathbf{R}_{wc} \cdot \mathbf{U}^c + \mathbf{T}_{wc} 入手。可见从c投影到w只需要对 \mathbf{U}^c 进行变换即可。也就是说，为了得到 \underline{\mathbf{m}}_{Bi} ，只需要对 \underline{\mathbf{m}}_{Ai} 进行刚性变换即可：

\underline{\mathbf{m}}_{B,i} = \Pi \left( \mathbf{R}_{BA} \mathbf{U}_i^A + \mathbf{t}_{BA} \right)

其中 \Pi(\cdot) 是投影函数。将前面得到的深度表达式代入：

\underline{\mathbf{m}}_{B,i}= \Pi \left( \mathbf{R}_{BA} \left( -\dfrac{d}{\mathbf{n}_A^\top \cdot \underline{\mathbf{m}}_{A,i}} \right) \cdot \underline{\mathbf{m}}_{A,i} + \mathbf{t}_{BA} \right)

将上公式左右两边都乘 -\dfrac{\mathbf{n}_A^\top \cdot \underline{\mathbf{m}}_{A,i}}{d} ：

\underline{\mathbf{m}}_{B,i}= \Pi \left( \mathbf{R}_{BA} \cdot \underline{\mathbf{m}}_{A,i} - \dfrac{\mathbf{n}_A^\top \cdot \underline{\mathbf{m}}_{A,i}}{d} \cdot \mathbf{t}_{BA} \right)
\underline{\mathbf{m}}_{B,i} = \Pi \left( \left( \mathbf{R}_{BA} - \dfrac{\mathbf{t}_{BA} \cdot \mathbf{n}_A^\top}{d} \right) \cdot \underline{\mathbf{m}}_{A,i} \right)

这样就得到了各自归一化平面上A点到B点的对应关系。

问题：上述公式中左右两边都乘了 -\dfrac{\mathbf{n}_A^\top \cdot \underline{\mathbf{m}}_{A,i}}{d} ，为什么等式还能保持不变？这是因为投影函数 \Pi(\cdot) 的特点是它是一个比例不变的操作——它只看方向和相对位置，不看绝对尺度。因此，即使我们在右边乘上这个系数，也不会影响等式成立的条件，因为投影结果相同。

2.2.3 寻找 \underline{\mathbf{P}}_{A,i} 和 \underline{\mathbf{P}}_{B,i} 之间的对应关系

现在我们已经建立了归一化平面上点的对应关系，接下来要将其推广到图像平面。我们已知：

\left\{ \begin{aligned} \underline{\mathbf{m}}_{A,i} = K_A^{-1} \cdot \underline{\mathbf{P}}_{A,i}\\ \underline{\mathbf{m}}_{B,i} = K_B^{-1} \cdot \underline{\mathbf{P}}_{B,i} \end{aligned} \right.

从 \underline{\mathbf{P}}_{B,i} = K_B \cdot \underline{\mathbf{m}}_{B,i} 出发，将前面得到的 \underline{\mathbf{m}}_{B,i} 带入：

\underline{\mathbf{P}}_{B,i} = K_B \cdot \Pi \left( \left( \mathbf{R}_{BA} - \dfrac{\mathbf{t}_{BA} \cdot \mathbf{n}_A^\top}{d} \right) \cdot \underline{\mathbf{m}}_{A,i} \right)
\underline{\mathbf{P}}_{B,i} = K_B \cdot \Pi \left( \left( \mathbf{R}_{BA} - \dfrac{\mathbf{t}_{BA} \cdot \mathbf{n}_A^\top}{d} \right) \cdot K_A^{-1} \cdot \underline{\mathbf{P}}_{A,i} \right)

这里我们需要回顾一个重要性质：

K \cdot \Pi \left( \begin{bmatrix} a \\ b \\ c \end{bmatrix} \right) = \Pi \left( K \cdot \begin{bmatrix} a \\ b \\ c \end{bmatrix} \right)

利用此性质，可得：

\underline{\mathbf{P}}_{B,i} = \Pi \left( K_B \cdot \left( \mathbf{R}_{BA} - \dfrac{\mathbf{t}_{BA} \cdot \mathbf{n}_A^\top}{d} \right) \cdot K_A^{-1} \cdot \underline{\mathbf{P}}_{A,i}\right)

2.2.4 得到单应性矩阵 \mathbf{H}_{AB}

现在我们可以定义单应性矩阵：

\mathbf{H}_{AB} = K_B \cdot \left( \mathbf{R}_{BA} - \dfrac{\mathbf{t}_{BA} \cdot \mathbf{n}_A^\top}{d} \right) \cdot K_A^{-1}

因此得到了点映射关系：

\left\{ \begin{aligned} &\underline{\mathbf{P}}_{B,i} = \Pi \left( \mathbf{H}_{BA} \cdot \underline{\mathbf{P}}_{A,i} \right) \quad \quad A \Rightarrow B\\ &\underline{\mathbf{P}}_{A,i} = \Pi \left( \mathbf{H}_{BA}^{-1} \cdot \underline{\mathbf{P}}_{B,i} \right) = \Pi \left( \mathbf{H}_{AB} \cdot \underline{\mathbf{P}}_{B,i} \right) \quad \quad B \Rightarrow A \\ \end{aligned} \right.

通过单应性矩阵，我们可以将某点从一个相机图片坐标系变换到另一个相机图片坐标系，实现了点映射关系。

2.2.5 单应性矩阵估计求解

单应性矩阵的形式为：

\mathbf{H}_{AB} = \begin{bmatrix} h_1 & h_4 & h_7 \\ h_2 & h_5 & h_8 \\ h_3 & h_6 & h_9 \end{bmatrix}

这是一个齐次矩阵，它有9个参数 h_1 到 h_9 。但齐次矩阵在尺度上具有冗余性，这会导致自由度的丢失。

简单的解法——参数化

为了解决这个问题，我们可以采用参数化的方法，将要估计的参数设为自由参数：

\mathbf{H}_{AB} = \begin{bmatrix} h_1 & h_4 & h_7 \\ h_2 & h_5 & h_8 \\ h_3 & h_6 & 1 \end{bmatrix}
\mathbf{h} = \begin{bmatrix} h_1 \\ \vdots \\ h_8 \end{bmatrix}

如何估计 \mathbf{h} ？

在这种参数化情况下，只要我们了解对应点，就可以求得 h_1 到 h_8 。从基本关系式出发：

\underline{\mathbf{P}}_{A,i} = \Pi \left( \begin{bmatrix} h_1 & h_4 & h_7 \\ h_2 & h_5 & h_8 \\ h_3 & h_6 & 1 \end{bmatrix} \cdot \underline{\mathbf{P}}_{B,i} \right)

由于 \underline{\mathbf{P}}_{A,i} 和 \underline{\mathbf{P}}_{B,i} 是齐次坐标，我们将其展开：

\begin{bmatrix} P_{A,i,x} \\ P_{A,i,y} \\ 1 \end{bmatrix} = \Pi \left( \begin{bmatrix} h_1 & h_4 & h_7 \\ h_2 & h_5 & h_8 \\ h_3 & h_6 & 1 \end{bmatrix} \cdot \begin{bmatrix} P_{B,i,x} \\ P_{B,i,y} \\ 1 \end{bmatrix} \right)

展开后得到：

\left\{ \begin{aligned} P_{A,i,x} = \dfrac{h_1 \cdot P_{B,i,x} + h_4 \cdot P_{B,i,y} + h_7}{h_3 \cdot P_{B,i,x} + h_6 \cdot P_{B,i,y} + 1} \\ P_{A,i,y} = \dfrac{h_2 \cdot P_{B,i,x} + h_5 \cdot P_{B,i,y} + h_8}{h_3 \cdot P_{B,i,x} + h_6 \cdot P_{B,i,y} + 1} \end{aligned} \right.

移项整理后：

\left\{ \begin{aligned} P_{A,i,x} \cdot \left( h_3 \cdot P_{B,i,x} + h_6 \cdot P_{B,i,y} + 1 \right) = h_1 \cdot P_{B,i,x} + h_4 \cdot P_{B,i,y} + h_7 \\ P_{A,i,y} \cdot \left( h_3 \cdot P_{B,i,x} + h_6 \cdot P_{B,i,y} + 1 \right) = h_2 \cdot P_{B,i,x} + h_5 \cdot P_{B,i,y} + h_8 \end{aligned} \right.

最终可以写成矩阵形式：

\begin{bmatrix} P_{B,i,x} & 0 & -P_{A,i,x} \cdot P_{B,i,x} & P_{B,i,y} & 0 & -P_{A,i,x} \cdot P_{B,i,y} & 1 & 0 \\ 0 & P_{B,i,x} & -P_{A,i,y} \cdot P_{B,i,x} & 0 & P_{B,i,y} & -P_{A,i,y} \cdot P_{B,i,y} & 0 & 1 \end{bmatrix} \begin{bmatrix} h_1 \\ h_2 \\ h_3 \\ h_4 \\ h_5 \\ h_6 \\ h_7 \\ h_8 \end{bmatrix} = \begin{bmatrix} P_{A,i,x} \\ P_{A,i,y} \end{bmatrix}

因为有8个未知数，需要八个独立的线性方程，而每对对应点可以提供两个对应方程，因此需要至少四对对应点。即需要四个匹配 (P_{A,i}, P_{B,i}) \quad i = 1, 2, 3, 4 ，然后通过线性最小二乘法求解： \mathbf{h}^* = \arg\min_{\mathbf{h}} \sum_{i=1}^{4} \left\lVert M_i \mathbf{h} - P_{A,i} \right\rVert_2^2 。

 

3. 使用RANSAC算法进行稳健的单应性估计

3.1 目标

图像对齐与拼接是计算机视觉中的基础问题。通过估计两幅图像之间的单应性（Homography），我们可以实现图像的自动拼接，这在全景图生成、图像配准等应用中非常重要。

3.2 自动建立对应关系——SIFT算法

在实际应用中，我们需要自动建立图像间的对应关系。使用SIFT等算法在两幅图像中检测特征点（由老师提供的代码），无需手动标记对应点，利用算法自动建立图像间的对应关系。通过兴趣点检测，我们可以找到两幅图像中最相似的点对。

但需要注意的是，这些点对并不一定正确对应，可能会出现错误匹配（离群点）。在这种情况下不可以直接用对应关系，我们将使用另一种算法RANSAC来自动评估对应点之间的正确性，并得到最理想的H矩阵并输出。

3.3 RANSAC算法进行稳健估计

RANSAC（Random Sample Consensus，随机抽样一致性）是一种在存在离群点（错误点）的情况下估计模型参数（H）的稳健算法。它的核心思想是通过反复随机抽样，寻找最符合的模型。

RANSAC流程

整个算法需要重复N次（迭代次数根据经验或计算确定）。每次迭代包含以下步骤：

①随机选取4对匹配点。4是估计单应性矩阵所需的最小匹配点数。从所有的匹配点中随机选四个，由于不确定哪个对应关系正确，所以后续中有一个估计评判标准（欧几里得距离）。

②估计单应性矩阵 H^k 。使用选取的4对匹配点，通过DLT算法（上个实验做过，其目的与作用是在已知对应点的情况下，将一个相机视角转换到另一个相机视角）估计单应性矩阵。

③计算误差并评估模型。对于所有匹配点（包括未选取的），将第二幅图像的点 P_{Bi} 通过估计的 H^k 转换，即 H^k P_{Bi} 。计算变换后的点（估计点）与第一幅图像实际点 P_{A_i} 之间的欧氏距离。定义代价函数时，我们使用二值核函数 \phi_c(d) ：当距离 d < \tau 时，认为匹配正确，代价为0；当距离 d \geq \tau 时，认为匹配错误，代价为1。总代价为 L^k = \sum_{i} \phi_c(|P_{A_i} - H^k P_{B_i}|) 。

④更新最佳模型。如果当前代价 L^k 小于之前的最小代价 L ，则更新 L 和对应的 H 。

最终输出具有最小代价的单应性矩阵 H 。

3.4 阈值的选择

阈值 \tau 是判断匹配是否为内点的距离阈值，它的选择对算法性能有重要影响。通常根据图像分辨率和匹配精度选择，一般在0.5到3个像素之间。选择过大会增加错误匹配，过小会忽略正确匹配。

3.5 为什么不用传统的二次代价函数

这涉及到敏感性问题。二次代价函数（如最小二乘法）对离群点非常敏感，如果某个点的误差很大，会导致代价函数值过大，这时即使其他点的误差很小也没有用。

相比之下，二值核函数具有更好的稳健性。它对那些特别大、离谱的点不敏感（都等于1），能够有效抑制离群点的影响，使得估计结果更稳健。当然，除了二值核函数，还存在其他稳健核函数，如Huber核、Lorentzian核等，它们可以在一定程度上兼顾误差大小和稳健性。

3.6 RANSAC算法的局限性

RANSAC算法虽然强大，但也有其局限性。当模型参数数量增加时，所需的随机采样次数会指数增长，计算成本显著提高。因此，RANSAC适用于参数数量较少的情况，如直线拟合、基础矩阵和单应性估计等。

 

 

4. 立体视觉中的对极几何

到目前为止，我们已经研究了平面场景的情况，使用了单应性来描述两个视图之间的关系。然而，对于一般的三维场景，平面假设不再成立。为此，我们需要引入对极几何这个更一般的框架。

4.1 对极几何

对极几何可以通过一个示意图很好地解释。考虑两个相机，分别位于参考系1和参考系2，相机1的光心为 O_1 ，相机2的光心为 O_2 ，空间中的一点 U 投影到两个相机的图像平面上，得到点 \underline{m}_1 和 \underline{m}_2 。

在这个设置下，我们面临的问题是：在一般情况下，我们无法对点 U 做出任何假设（与之前的平面场景不同）。我们需要找到一种方法，在不知道 U 的情况下，建立 \underline{m}_1 和 \underline{m}_2 之间的关系。这就是对极几何要解决的核心问题。

4.2 对极平面和对极线

对极平面

点 U 和两个光心 O_1 、 O_2 定义了一个平面，这意味着点 m_1 、 m_2 、 O_1 、 O_2 共面，我们称这个平面为对极平面 "Contrainte épipolaire = coplanarité"，即 \underline{m}_1 、 \underline{m}_2 、 O_1 、 O_2 共面。这个共面性条件是整个对极几何的基础，在立体视觉中，基础矩阵 F 和本质矩阵 E 都依赖于这个共面性条件来计算。

对极线

对极平面与两个相机的图像平面相交，分别得到对极线 l_1 和 l_2 。这里有一个重要的几何约束： m_2 是三维点 U 在第二个图像平面的投影，但根据对极几何的约束， m_2 必须位于对极线 l_2 上。

这给我们提供了一个强大的工具：给定点 m_1 的位置，我们可以通过基础矩阵 F 确定对应的对极线 l_2 。基础矩阵 F 捕捉了两个相机之间的相对姿态和内在参数信息。这个关系表明，给定点 m_1 ，我们可以计算出 m_2 必须位于的对极线 l_2 。

4.3 对极约束

现在让我们形式化对极约束，建立 m_1 和 m_2 之间的数学关系。首先定义必要的向量：

\left\{ \begin{aligned} &\mathbf{\underline{m}_1} \text{ is the vector from the optical center } O_1 \text{ to the image point } \underline{m}_1 \quad \overrightarrow{O_1 m_1}^{1} \\ &\mathbf{\underline{m}_2} \text{ is the vector from the optical center } O_2 \text{ to the image point } \underline{m}_2 \quad \overrightarrow{O_2 m_2}^{2} \\ &\mathbf{t_{12}} = \overrightarrow{O_1 O_2}^{1} \text{ is the translation vector between the two camera optical centers} \end{aligned} \right.

接下来定义对极平面的法向量：

\left\{ \begin{aligned} &\text{In reference frame } 1, \quad \overrightarrow{\mathbf{n}_1}^{1} = \underline{\mathbf{m}}_1 \times \mathbf{t}_{12} \\ &\text{In reference frame } 2, \quad \overrightarrow{\mathbf{n}_2}^{2} = \mathbf{R}_{21} \overrightarrow{\mathbf{n}_1}^{1}, \text{ where } \mathbf{R} \text{ is the rotation matrix between the cameras} \end{aligned} \right.

这里需要注意两个重要点。首先， \times 表示两个向量之间的叉积运算。叉积的结果是一个向量，它垂直于运算的两个向量，方向由右手定则决定，大小为这两个向量构成的平行四边形的面积。

其次，法向量的坐标系变换不用考虑平移部分，因为单位法向量并不是坐标位置，方向向量在旋转过程中大小不变，不受平移的影响。总而言之，法向量只考虑旋转矩阵，而点则需要考虑旋转加平移。

继续推导法向量在参考系2中的表达：

\overrightarrow{\mathbf{n}_2}^{2} = \mathbf{R}_{21} \cdot \overrightarrow{\mathbf{n}_1}^{1} = \mathbf{R}_{21} \cdot \left( \underline{\mathbf{m}}_1 \times \mathbf{t}_{12} \right) = \mathbf{R}_{21} \cdot \underline{\mathbf{m}}_1 \times \mathbf{R}_{21} \cdot \mathbf{t}_{12}

由于之前我们已知 \mathbf{t}_{21} = \mathbf{R}_{21} \cdot \mathbf{t}_{12} ，所以上式变为：

\overrightarrow{\mathbf{n}_2}^{2} = \mathbf{t}_{21} \times \left( \mathbf{R}_{21} \cdot \underline{\mathbf{m}}_1 \right)

现在让我们回顾叉积运算的一个重要性质。叉积可以表示为矩阵运算：

\mathbf{a} \times \mathbf{b} = \begin{bmatrix} a_x \\ a_y \\ a_z \end{bmatrix} \times \begin{bmatrix} b_x \\ b_y \\ b_z \end{bmatrix} = \begin{bmatrix} a_y b_z - a_z b_y \\ a_z b_x - a_x b_z \\ a_x b_y - a_y b_x \end{bmatrix}_{3 \times 1} \Rightarrow \left[\mathbf{a}\right]_{\times} = \begin{bmatrix} 0 & -a_z & a_y \\ a_z & 0 & -a_x \\ -a_y & a_x & 0 \end{bmatrix}
\mathbf{a} \times \mathbf{b} = \left[\mathbf{a}\right]_{\times} \mathbf{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}

利用上述性质，我们可以看出叉积运算可以变成矩阵运算，因此：

\overrightarrow{\mathbf{n}_2}^{2} = \mathbf{t}_{21} \times \left( \mathbf{R}_{21} \cdot \underline{\mathbf{m}}_1 \right) = \left[ \mathbf{t}_{21} \right]_{\times} \cdot \mathbf{R}_{21} \cdot \underline{\mathbf{m}}_1

因为 \overrightarrow{\mathbf{n}_2}^{2} 是 \mathbf{m}_2 所在平面的法线，所以有 \mathbf{m}_2^\top \cdot \overrightarrow{\mathbf{n}_2}^{2} = 0 ，可得出：

\mathbf{m}_2^\top \cdot \left( \left[ \mathbf{t}_{21} \right]_{\times} \cdot \mathbf{R}_{21} \right) \cdot \underline{\mathbf{m}}_1 = 0

4.4 本质矩阵（matrice essentielle）

基于上面的推导，我们可以定义本质矩阵：

\mathbf{E}_{21} = \left[ \mathbf{t}_{21} \right]_{\times} \cdot \mathbf{R}_{21} \quad \Rightarrow \quad \text{matrice essentielle}

它包含了两个相机之间的相对旋转 \mathbf{R} 和平移 \mathbf{t} 的信息。有了本质矩阵，对极约束可以简洁地表示为：

\underline{\mathbf{m}}_2^\top \cdot \mathbf{E}_{21} \cdot \underline{\mathbf{m}}_1 = 0

自由度分析

本质矩阵的自由度：

5 \text{ degre de liberte} \\ \downarrow\\ 5 \text{ DDL} \left( \begin{array}{c} 3 \, \mathbf{R} _{21} \quad \text{rotation} \\ \quad 2 \, \mathbf{t} _{21} \quad \text{translation} \end{array} \right)\\ \downarrow\\ \quad \quad \| \mathbf{t} _{21} \|_2 \quad \text{ inconnue}
自由度：\left\{ \begin{aligned} &\text{The rotation matrix } \mathbf{R} \text{ has 3 degrees of freedom} \\ &\text{The translation vector } \mathbf{t} \text{ has 2 degrees of freedom (since the scale is unknown)} \\ &\text{Therefore, } \mathbf{E} \text{ has 5 degrees of freedom} \end{aligned} \right.

自由度（degree of freedom, DoF）是指描述本质矩阵所需的独立参数数量。在几何和线性代数中，自由度反映了系统在不受限制的情况下可以独立变化的方向或方式。旋转矩阵具有3个自由度，描述了三维空间中的旋转。平移向量理论上在三维空间中有3个自由度，但是本质矩阵中的平移向量一般只关注方向，对长度忽略（未知），所以平移向量只剩下2个有效的自由度，描述了平移的方向。

4.5 基础矩阵（Matrix Fundamental）

在很多实际情况下，相机的内参是未知的。这时我们需要引入基础矩阵。继续对前面的公式进行变换：

\underline{\mathbf{m}}_2^\top \cdot \mathbf{E}_{21} \cdot \underline{\mathbf{m}}_1 = 0

我们已知归一化坐标和像素坐标的关系：

\left\{ \begin{aligned} \underline{\mathbf{m}}_2 = K^{-1} \cdot \underline{\mathbf{P}}_2 \\ \underline{\mathbf{m}}_1 = K^{-1} \cdot \underline{\mathbf{P}}_1 \end{aligned} \right.

代入后得到：

\underline{\mathbf{P}}_2^\top \cdot (K^{-1})^\top \cdot \mathbf{E}_{21} \cdot K^{-1} \cdot \underline{\mathbf{P}}_1 = 0

当相机内参未知或未被考虑时，我们引入基础矩阵 \mathbf{F} 来覆盖 K ：

\mathbf{F}_{21}= (K^{-1})^\top \cdot \mathbf{E}_{21} \cdot K^{-1}

基础矩阵的性质：

\mathbf{F}_{21} : \text{ matrice fondamentale} \quad \Rightarrow \quad 7 \text{ DDL}\quad \left\{ \begin{aligned} & \text{- matrice homogène} \\ & \text{- rang}(\mathbf{F}_{21}) = 2 \quad \Rightarrow \quad \det(\mathbf{F}_{21}) = 0 \end{aligned} \right.

基础矩阵 \mathbf{F} 是齐次矩阵，可以乘以任意非零标量而不改变其性质。另一个重要性质是 \mathbf{F} 的秩为2，这意味着 \det(\mathbf{F}_{21}) = 0 。

有了基础矩阵，对极约束变为：

\underline{\mathbf{P}}_2^\top \cdot \mathbf{F}_{21} \cdot \underline{\mathbf{P}}_1 = 0

设 \mathbf{L}_2 = \mathbf{F}_{21} \cdot \underline{\mathbf{P}}_1 = \begin{bmatrix} a \\ b \\ c \end{bmatrix} ，则有：

\underline{\mathbf{P}}_2^\top \cdot \mathbf{L}_2 = 0 \quad \Leftrightarrow \quad a P_{2,x} + b P_{2,y} + c = 0

这就是相机2的图像平面中的直线方程，即对极线。

4.6 本质和基础矩阵的估计

根据相机是否已经校准，我们有不同的估计策略。如果相机已经校准，我们需要估计本质矩阵 \mathbf{E} （5个自由度），可以使用5点对应算法。如果相机未校准，我们需要估计基础矩阵 \mathbf{F} （7个自由度），可以使用7点对应算法。在实际应用中，常用8点对应算法，这个算法故意忽略约束条件 \det(\mathbf{F}) = 0 ，简化了求解过程。

4.7 算法（8点对应算法）

8点对应算法的步骤如下：

首先收集匹配点对。虽然 \mathbf{F} 有7个自由度，但在算法中忽略了秩为2的约束，因此需要至少8对匹配点来估计 \mathbf{F} 。

然后构建线性方程组。对于每一对匹配点 (m_1, m_2) ，构建方程：

\underline{\mathbf{m}}_2^\top \cdot \mathbf{E}_{21} \cdot \underline{\mathbf{m}}_1 = 0
\underline{\mathbf{P}}_2^\top \cdot \mathbf{F}_{21} \cdot \underline{\mathbf{P}}_1 = 0

最后求解线性系统，将方程组表示为 \underline{\mathbf{P}}_2^\top \cdot \mathbf{L}_2 = 0 。

RANSAC算法步骤

为了处理匹配点对中的离群点（错误匹配），稳健地估计 \mathbf{F} ，我们再次使用RANSAC算法。首先随机采样，利用8点对应算法来估计 \mathbf{F} 。然后评估模型，利用估计得到的 \mathbf{F} 来计算所有匹配点对的对极约束误差，即点到对应对极线的距离。根据设定的距离阈值，判断哪些匹配点是内点。重复上述过程，直到找到内点数量最多的模型。

5. 束调整

束调整是一种同时优化摄像机参数（包括位置、姿态和内参）和场景中三维点位置的技术。其核心思想是通过最小化三维点在图像上的重投影误差，使得优化后的模型与实际观测更加吻合。最关键五个字：最小化投影误差。

5.1 两个摄像机的情况

5.1.1 数据

在两个摄像机的情况下，我们的输入数据是N个对应点对：

(P_{A,i}, P_{B,i})_{i=1,\dots,N} \implies N \text{ correspondences}

5.1.2 要估计的参数

我们需要估计的参数包括摄像机的姿态以及三维点云数据集：

\mathbf{R}_{W1} \quad \mathbf{t}_{W1} \quad \mathbf{R}_{W2} \quad \mathbf{t}_{W2} \quad \left\{\mathbf{U}^w_i \right\}_{i=1,\dots,N}

5.1.3 损失函数

束调整的损失函数定义为所有重投影误差的总和：

\mathcal{L} \left( \mathbf{R}_{w1}, \mathbf{t}_{w1}, \mathbf{R}_{w2}, \mathbf{t}_{w2}, \left\{ \mathbf{U}^w_i \right\}_{i=1,\dots,N} \right)= \sum_{i=1}^{N} \left( \left\lVert P_{1,i} - K_1 \Pi \left( \mathbf{R}_{w1}^\top \mathbf{U}_i^{w} - \mathbf{R}_{w1}^\top \mathbf{t}_{w1} \right) \right\rVert_2^2 +\left\lVert P_{2,i} - K_2 \Pi \left( \mathbf{R}_{w2}^\top \mathbf{U}_i^{w} - \mathbf{R}_{w2}^\top \mathbf{t}_{w2} \right) \right\rVert_2^2 \right)

K_A 和 K_B 是摄像机A和B的内参矩阵， \Pi(\cdot) 是投影函数，将三维点投影到二维平面上。

\mathbf{R}_{w1}^\top 和 \mathbf{R}_{w2}^\top 等价于 \mathbf{R}_{1w} 和 \mathbf{R}_{2w} ，即将点从世界坐标系转换到摄像机坐标系。

\mathbf{R}_{w1}^\top \mathbf{t}_{w1} 等价于 \mathbf{t}_{1w} ，表示平移向量。

变换过程可以表示为：

\mathbf{U}_i^{1} = \mathbf{R}_{w1}^\top \cdot \mathbf{U}_i^{w} - \mathbf{R}_{w1}^\top \cdot \mathbf{t}_{w1}

也就是将 \mathbf{U}_i^{w} 变换到 \mathbf{U}_i^{1} ，即从世界坐标系变换到相机坐标系。损失函数中的每一项都是做差：相机A或B中的图像坐标（实际）减去三维空间旋转变换得来的估计图像坐标，这就是重投影误差。

5.2 多个摄像机的情况

5.2.1 数据

在多摄像机的情况下，每张图像中检测到的点为：

\left\{ \left\{ P_{m,i} \right\}_{i=1,\dots,N_m} \right\}_{m=1,\dots,M}

这些点在不同视角下的图像中可以形成轨迹（tracks）。第 m 个摄像机检测到的点，其中 N_m 是第 m 个摄像机检测到的点的数量。

我们还需要维护点的对应关系：

\left\{ \text{p2d-id}_m, \ \text{p3d-id}_m \right\}_{m=1,\dots,M}

其中 \text{p2d-id}_m 是二维点在图像中的索引， \text{p3d-id}_m 是对应的三维点在点云中的索引。它们的大小尺寸都是 C_m \times 1 。

5.2.2 要估计的参数

在多摄像机情况下，要估计的参数包括所有相机的外参：

\left\{ \left( \mathbf{R}_{wm}, \mathbf{t}_{wm} \right) \right\}_{m=1,\dots,M}

以及所有三维点的位置：

\left\{ \mathbf{U}_i^{w} \right\}_{i=1,\dots,N}

5.2.3 损失函数

代价函数扩展为对所有摄像机和所有检测到的点进行误差计算，将投影点与实际观测点之间的距离最小化：

\mathcal{L}(x) = \sum_{m=1}^{M} \sum_{c=1}^{C_m} \left\| P_{m,\ \text{p2d-id}_m(c)} - K_m \Pi \left( \mathbf{R}_{wm}^\top \mathbf{U}_{\text{p3d-id}_m(c)}^{w} - \mathbf{R}_{wm}^\top \mathbf{t}_{wm} \right) \right\|_2^2

其中 C_m 是第 m 台摄像机的观测数量， \mathbf{U}_{\text{p3d-id}_m(c)}^{w} 是与观测对应的三维点。

我们可以简单地将上述代价函数简化成更一般的形式：

\mathcal{L}(x) = \sum_{i=1}^{N} \left\| f_i(x) \right\|_2^2 \quad \left\{ \begin{array}{l} x \in \mathbb{R}^D \\ f_i : \mathbb{R}^D \rightarrow \mathbb{R}^B \end{array} \right.

这里 x 是所有待优化的参数（摄像机参数和三维点坐标）， f_i(x) 是第 i 个残差函数，表示第 i 个观测的重投影误差。我们的目标是找到 x ，使得 \mathcal{L}(x) 最小化。这是一个非线性最小二乘问题，通常使用迭代的方法求解。

5.3 高斯牛顿算法

高斯牛顿算法是用于非线性最小二乘问题的一种迭代优化算法，其基本形式为： \delta_{k+1} = \delta_k + d_k 。

算法的核心思想是对 f_i 进行线性化：

\text{Linearisation de } f_i: \quad f_i(x_k + d_k) \approx f_i(x_k) + \mathbf{J}_i(x_k)\cdot d_k

这里 \delta x 是参数的增量，需要求解。对于每次迭代，我们在当前估计 x_k 附近对 f_i(x) 进行泰勒展开，并忽略高阶项。

各个变量的维度关系为：

f_i(\delta_k + d_k) \in \mathbb{R}^B
f_i(\delta_k) \in \mathbb{R}^B \quad \mathbf{J}_i(\delta_k) \in \mathbb{R}^{B \times D} \quad d_k \in \mathbb{R}^D

其中雅可比矩阵为：

\mathbf{J}_i(x_k) = \frac{\partial f_i(x_k + d_k)}{\partial d_k} \bigg|_{d_k=0}

这代表了在点 x_k 处函数 f_i 对于 d_k 的偏导数，并且此偏导数是在 d_k = 0 的条件下计算的。雅可比矩阵描述了在点 x_k 处函数 f_i 的线性变化率。

线性最小二乘法

将所有残差的线性化组合起来，我们得到线性最小二乘问题：

L_k(d_k) = \sum_{i=1}^{N} \left\| f_i(x_k) + \mathbf{J}_i(x_k) \cdot d_k \right\|_2^2

定义组合的雅可比矩阵和残差向量：

\mathbf{J}_k = \begin{bmatrix} J_1(x_k) \\ J_2(x_k) \\ J_3(x_k) \\ \vdots \\ J_N(x_k) \end{bmatrix} \quad \mathbf{b}_k = \begin{bmatrix} f_1(x_k) \\ f_2(x_k) \\ \vdots \\ f_N(x_k) \end{bmatrix}

其中 \mathbf{b}_k 是所有残差的组合。线性最小二乘问题变为：

L_k(d_k) = \lVert b_k + J_k \cdot d_k \rVert_2^2

通过最小化 L_k(\delta x) ，我们可以得到正规方程：

J_k^T \cdot J_k \cdot d_k = -J_k^T \cdot b_k

其中 b_k \propto \text{ gradient} 。左边的矩阵 \mathbf{J}_k^T \mathbf{J}_k 是海森矩阵的近似，右边的向量 -\mathbf{J}_k^T \mathbf{b}_k 是梯度的负值。求解这个线性系统，得到参数更新 d_x 。

Levenberg-Marquardt算法

Levenberg-Marquardt算法在高斯-牛顿算法的基础上引入阻尼因子 \lambda ，使得优化过程在接近解时具有高斯-牛顿的快速收敛特性，而在远离解时具有梯度下降的稳定性。这是一种常用于非线性最小二乘问题的迭代优化算法。

目标函数变为：

L_k(d_k) = \lVert b_k + J_k d_k \rVert_2^2 + \lambda \lVert d_k \rVert_2^2

对应的正规方程为：

(J_k^T J_k + \lambda I_k)d_k = -J_k^T b_k

这里 \lambda 是阻尼因子，它的作用可以这样理解：

\begin{cases} \text{Si } \lambda = 0 & \Rightarrow \text{Gauss-Newton} \\ \text{Si } \lambda \rightarrow +\infty & \Rightarrow \lambda d_k \rightarrow -J_k^T b_k \quad \text{descente de gradient} \end{cases}

如果新的代价函数值降低了（说明更新有效），则减小 \lambda ，使算法更接近高斯-牛顿法，加快收敛。如果代价函数值没有降低，则增大 \lambda ，使算法更接近梯度下降法，保证稳定性。

5.4 算法步骤总结

在实际应用中，Levenberg-Marquardt算法的步骤如下：

初始化：

  • 设定初始参数x和阻尼因子 \lambda

  • 计算初始代价函数 L_{\min}

迭代：

  • 计算雅可比矩阵 \mathbf{J} 和残差 \mathbf{b}

  • 求解线性系统：

    (J^T J + \lambda I_d) d = -J^T b

  • 更新参数：

    x′=x+d

  • 计算新的代价函数 L'

判断更新效果：

  • 如果 L' < L_{\min} （代价函数降低）：

    • 接受更新： x = x' ， L_{\min} = L'

    • 减小 \lambda：\lambda = \lambda / 2

    • 继续迭代

  • 否则（代价函数未降低）：

    • 拒绝更新，不改变 x

    • 增大 \lambda：\lambda = 2\lambda

    • 检查 \lambda 是否超过最大值，若超过则停止迭代

终止条件：

  • 当 \lambda 超过预设的最大值，或者参数更新的幅度小于阈值时，停止迭代

#图像处理 8
反问题：Wiener-Hunt 方法
应用随机梯度算法解决信道均衡问题