SLAM入门之视觉里程计

视觉里程计 VO

7 视觉里程计(1)—特征点法

7.1 特征点法

从图像中选取比较有代表性的点，这些点在相机视角发生少量变化后会保持不变，可以在各个图像中找到相同的点。在经典 SLAM 模型中，这些点称为路标（ Landmark）。在视觉 SLAM 中，路标则是指图像特征（ Feature）。特征点由关键点（key-point）和描述子（descriptor）两部分组成。

• 关键点：该特征点在图像里的位置，有些特征点还具有朝向、大小等信息。
• 描述子：通常是一个向量，按照某种人为设计的方式，描述了该关键点周围像素的信息。描述子是按照“外观相似的特征应该有相似的描述子”的原则设计的。因此，只要两个特征点的描述子在向量空间上的距离相近，就可以认为它们是同样的特征点。
• 特征点常提取的是角点，提取算法有Harris 角点、FAST 角点、GFTT 角点等。但有时候从远处看上去是角点的地方，当相机走近之后，可能就不显示为角点了。或者，当旋转相机时，角点的外观会发生变化，我们也就不容易辨认出那是同一个角点。现在已有的稳定的局部图像特征有：SIFT、SURF、ORB等。SIFT（尺度不变特征变换）性能强但计算量太大，很难实时计算。也有的适当降低精度和鲁棒性，以提升计算的速度，如FAST关键点（无描述子），但其不具有方向性。ORB特征是对其改进后的具有代表性的实时图像特征。

ORB 特征

ORB：Oriented FAST and Rotated BRIEF

1. 关键点：Oriented FAST，一种改进的 FAST 角点。FAST 是一种角点，主要检测局部像素灰度变化明显的地方，以速度快著称。针对 FAST 角点不具有方向性和尺度的弱点， ORB 添加了尺度和旋转的描述。尺度不变性由构建图像金字塔，并在金字塔的每一层上检测角点来实现。而特征的旋转是由灰度质心法（ Intensity Centroid）实现的。
2. 描述子：ORB 使用改进的 BRIEF 特征描述，BRIEF 是一种二进制描述子，其描述向量由许多个 0 和 1 组成，这里的 0 和 1 编码了关键点附近两个随机像素（比如 p 和 q）的大小关系：如果 p 比 q 大，则取 1，反之就取 0。如果我们取了 128个这样的 p，q，最后就得到 128 维由 0、 1 组成的向量。BRIEF 使用了随机选点的比较，速度非常快，而且由于使用了二进制表达，存储起来也十分方便，适用于实时的图像匹配。原始的 BRIEF 描述子不具有旋转不变性，因此在图像发生旋转时容易丢失。而 ORB 在FAST 特征点提取阶段计算了关键点的方向，所以可以利用方向信息，计算了旋转之后的“Steer BRIEF”特征使 ORB 的描述子具有较好的旋转不变性。

FAST 检测过程

1. 在图像中选取像素 $p$，假设它的亮度为 $I_p$。
2. 设置一个阈值 $T$（比如， $I_p$ 的 20%）。
3. 以像素 $p$ 为中心，选取半径为 3 的圆上的 16 个像素点。
4. 假如选取的圆上有连续的 $N$ 个点的亮度大于 $I_p+T$ 或小于 $I_p-T$，那么像素 $p$ 可以被认为是特征点（ N 通常取 12，即为 FAST-12)。
5. 循环以上四步，对每一个像素执行相同的操作。

• 金字塔：指对图像进行不同层次的降采样，以获得不同分辨率的图像
  每往上一层，就对图像进行一个固定倍率的缩放，这样就有了不同分辨率的图像。较小的图像可以看成是远处看过来的场景。在特征匹配算法中，可以匹配不同层上的图像，从而实现尺度不变性。例如，如果相机在后退，那么应该能够在上一个图像金字塔的上层和下一个图像的下层中找到匹配。
• 灰度质心：质心是指以图像块灰度值作为权重的中心，在一个小的图像块 B 中，定义图像块的矩，通过矩找到图像块的质心，连接图像块的几何中心 O 与质心 C，得到一个方向向量 $\overrightarrow{OC}$，于是可以定义特征点的方向。

特征匹配

最简单的特征匹配方法就是暴力匹配（ BruteForce Matcher）。即对每一个特征点 $x_t^m$ 与所有的 $x_{t+1}^n$ 测量描述子的距离，然后排序，取最近的一个作为匹配点。描述子距离表示了两个特征之间的相似程度。而对于二进制的描述子（如BRIEF），使用汉明距离作为度量——两个二进制串之间的汉明距离，指的是其不同位数的个数。
快速近似最近邻（ FLANN） 算法更加适合于匹配点数量极多的情况。

匹配好点对后根据点对来估计相机的运动。

• 单目：根据两组 2D 点估计运动，用对极几何来解决
• 双目、 RGB-D：根据两组 3D点估计运动，用 ICP 来解决
• 一组为 3D，一组为 2D：即一些 3D 点和它们在相机的投影位置，通过 PnP 求解。

7.2 对极几何（2D−2D）

7.2.1对极约束

• 几何关系
  

$O_1，O_2，P$ 三点确定一个平面称为极平面 $e_1，e_2$ 称为极点（ Epipoles），$O_1{O}_2$ 被称为基线，极平面与两个像平面 $I_1,I_2$ 之间的相交线 $l_1,l_2$ 为极线（ Epipolar line）。
特征点$p_1$是$O_{1}P$所在直线上所有点的投影，特征点$p_2$是$O_{2}P$所在直线上所有点的投影，又应为$p_1$和$p_2$匹配，所以就能推断出 P 的空间位置，以及相机的运动

• 代数关系

  • ${x\text{x}x}_2^{T}t\text{t}t$^${Rx\text{Rx}Rx}_1=0$

  • ${p\text{p}p}_2^T{K\text{K}K}^{-T}t\text{t}t$^$R\text{R}R{K\text{K}K}^{-1}{p\text{p}p}_1=0$

  • 令$E\text{E}E=t\text{t}t$^$R\text{R}R$ ：本质矩阵$E\text{E}E$(Essential Matrix)

  • 令$F\text{F}F={K\text{K}K}^{-T}{EK\text{EK}EK}^{-1}$ ：基础矩阵$F\text{F}F$( Fundamental Matrix)

  • 简化得 ${x\text{x}x}_2^{T}E\text{E}E{x\text{x}x}_1={p\text{p}p}_2^{T}F\text{F}F{p\text{p}p}_1=0$

对极约束简洁地给出了两个匹配点的空间位置关系。相机位姿估计问题变为以下两步

1. 根据配对点的像素位置求出 $E\text{E}E$ 或者 $F\text{F}F$。（实践中常用$E\text{E}E$）
2. 根据 $E\text{E}E$ 或者 $F\text{F}F$ 求出 $R，t\text{R，t}R，t$。

7.2.2 本质矩阵

$E\text{E}E=t\text{t}t$^$R\text{R}R$

• 这是一个是一个 3 × 3 的矩阵，内有 9 个未知数。采用八点法估计 $E\text{E}E$，即使用8对点。
• 估得的本质矩阵 $E\text{E}E$ 后 奇异值分解（ SVD） 得到相机的运动 $R，t\text{R，t}R，t$。

${t\text{t}t}_1$^ $={UR\text{UR}UR}_z(\frac{\pi }{2}){\Sigma U\text{ΣU}\Sigma U}^T$	${R\text{R}R}_1={UR\text{UR}UR}_z^T(\frac{\pi }{2}){V\text{V}V}^T$
${t\text{t}t}_2$^ $={UR\text{UR}UR}_z(-\frac{\pi }{2}){\Sigma U\text{ΣU}\Sigma U}^T$	${R\text{R}R}_2={UR\text{UR}UR}_z^T(-\frac{\pi }{2}){V\text{V}V}^T$

从 $E\text{E}E$ 分解到 $R，t\text{R，t}R，t$ 时，一共存在 4 个可能的解，只有第一种解中 P 在两个相机中都具有正的深度。因此，只要把任意一点代入 4 种解中，检测该点在两个相机下的深度，就可以确定哪个解是正确的了。

7.2.3 单应矩阵

单应矩阵(Homography) $H\text{H}H$，它描述了两个平面之间的映射关系。若场景中的特征点都落在同一平面上（比如墙、地面等），则可以通过单应性来进行运动估计

• 单应性在 SLAM 中具有重要意义。当特征点共面或者相机发生纯旋转时，基础矩阵的自由度下降，这就出现了所谓的退化（ degenerate）。现实中的数据总包含一些噪声，这时候如果继续使用八点法求解基础矩阵，基础矩阵多余出来的自由度将会主要由噪声决定。为了能够避免退化现象造成的影响，通常我们会同时估计基础矩阵 $F\text{F}F$ 和单应矩阵 $H\text{H}H$，选择重投影误差比较小的那个作为最终的运动估计矩阵。

尺度不确定性

$E\text{E}E$ 具有尺度等价性，所以分解的 $R,t\text{R,t}R,t$ 也具有尺度等价性，但因为$R\text{R}R\in SO(3)$ 自身具有约束，所以 $t\text{t}t$ 具有一个尺度，分解过程中对 $t\text{t}t$ 乘以任意非零常数分解都成立，因此常对t进行归一化，使其长度为1。这直接导致了单目视觉的尺度不确定性。例如，程序中输出的 $t\text{t}t$ 第一维约为 0.822。这个 0.822 究竟是指 0.822 米还是 0.822 厘米，我们是没法确定的。

• 初始化 ： 对两张图像的 $t\text{t}t$ 归一化相当于固定了尺度。虽然不知道它的实际长度是多少，但以这时的 $t\text{t}t$ 为单位 1，计算相机运动和特征点的 3D 位置。这被称为单目 SLAM 的初始化
• 另一种方法是令初始化时所有的特征点平均深度为 1，也可以固定一个尺度。相比于令 t 长度为 1 的做法，把特征点深度归一化可以控制场景的规模大小，使计算在数值上更稳定些。

初始化的纯旋转问题

• 单目初始化不能只有纯旋转，必须要有一定程度的平移。从 $E\text{E}E$ 分解到 $R,t\text{R,t}R,t$ 的过程中，如果相机发生的是纯旋转，导致 $t\text{t}t$ 为零，那么，得到的$E\text{E}E$ 也将为零，这将导致我们无从求解 $R\text{R}R$。让相机进行左右平移以顺利地进行初始化。

7.3 三角测量 (三角化 Triangulation)

• 对极几何约束估计了相机运动后需要用相机的运动估计特征点的空间位置

三角测量是指，通过在两处观察同一个点的夹角，从而确定该点的距离。理论上直线 $O_1{p}_1$ 与 $O_2{p}_2$ 在场景中会相交于一点 $P$，该点即两个特征点所对应的地图点在三维场景中的位置。然而由于噪声的影响，这两条直线往往无法相交。因此，可以通过最二小乘法求解。

• 设 $x_1，x_2$ 为两个特征点的归一化坐标，已知 $R,t\text{R,t}R,t$ ，求解两特征点深度$s_1，s_2$
  • $s_1{x}_1\text{x1}{x}_1=s_2{Rx\text{Rx}Rx}_2+t\text{t}t$
  • $s_1{x}_1\text{x1}{x}_1$^$x_1\text{x1}{x}_1=0=s_2{x}_1\text{x1}{x}_1$^${Rx\text{Rx}Rx}_2+x_1\text{x1}{x}_1$ ^ $t\text{t}t$
• 右侧可看成 $s_2$ 的一个方程，可以根据它直接求得 $s_2$，然后求出 $s_1$,得到了两帧下的点的深度，确定了它们的空间坐标。当然，由于噪声的存在，我们估得的 $R,t\text{R,t}R,t$ 不一定精确使上式为零，所以更常见的做法是求最小二乘解而不是零解。

三角测量的矛盾：增大平移，会导致匹配失效；而平移太小，则三角化精度不够

7.4 PnP（3D-2D)

• 描述了当我们知道n 个 3D 空间点以及它们的投影位置时，如何估计相机所在的位姿。如果两张图像中，其中一张特征点的 3D 位置已知，那么最少只需三个点对（需要至少一个额外点验证结果）就可以估计相机运动。
• PnP 问题解法用三对点估计位姿的 P3P，直接线性变换（DLT），EPnP（Efficient PnP）等。还能用非线性优化的方式，构建最小二乘问题并迭代求解，也就是万金油式的 Bundle Adjustment 。

1. 直接线性变换

   • 设某个空间点 $P$ 的齐次坐标为 $P=(X，Y，Z，1{)}^T$。在图像 $I_1$ 中，投影到特征点 $x_1=(u_1,v_1,1{)}^T$（以归一化平面齐次坐标表示）。此时相机的位姿 $R,t\text{R,t}R,t$ 是未知的。与单应矩阵的求解类似，定义增广矩阵 $[R|t\text{R∣t}R|t]$ 为一个 3 × 4 的矩阵，包含了旋转与平移信息。(和 SE(3) 中的变换矩阵 $T\text{T}T$ 是不同的)
   • 定义这个矩阵的行向量分别为$t_1=(t_1,t_2,t_3,t_4{)}^T,t_2=(t_5,t_6,t_7,t_8{)}^T,t_3=(t_9,t_{10},t_{11},t_{12}{)}^T$ $t\text{t}t$ 是待求的变量，一共12维，最少六对匹配点可求解出该矩阵。
     $$\left(\begin{array}{ccc}{P\text{P}P}_1^T& 0& -u_1{P\text{P}P}_1^T\\ 0& {P\text{P}P}_1^T& -v_1{P\text{P}P}_1^T\\ .& .& .\\ .& .& .\\ {P\text{P}P}_N^T& 0& -u_N{P\text{P}P}_N^T\\ 0& {P\text{P}P}_N^T& -v_N{P\text{P}P}_N^T\end{array}\right)\left(\begin{array}{c}{t\text{t}t}_1\\ {t\text{t}t}_2\\ {t\text{t}t}_3\end{array}\right)=0$$
   • 在 DLT 求解中，我们直接将 $T\text{T}T$ 矩阵看成了 12 个未知数，忽略了它们之间的联系。因为旋转矩阵 $R\text{R}R\in SO(3)$，用 DLT 求出的解不一定满足该约束，它是一个一般矩阵。平移向量比较好办，它属于向量空间。对于旋转矩阵 $R\text{R}R$，我们必须针对 DLT 估计的 $T\text{T}T$ 的左边3 × 3 的矩阵块，寻找一个最好的旋转矩阵对它进行近似。这可以由 QR 分解完成 [3, 48]，相当于把结果从矩阵空间重新投影到 $SE(3)$ 流形上，转换成旋转和平移两部分。

2. P3P

   1. 仅使用三对匹配点，对数据要求较少。它的输入数据为三对 3D-2D 匹配点。记 3D点为 A，B，C， 2D 点为 a， b， c，其中小写字母代表的点为大写字母在相机成像平面上的投影。此外， P3P 还需要使用一对验证点，以从可能的解出选出正确的那一个（类似于对极几何情形）。记验证点对为 D − d，相机光心为 O。求解思路是利用三角形相似性质，求解投影点 a，b，c 在相机坐标系下的 3D 坐标，最后把问题转换成一个 3D 到 3D 的位姿估计问题。
   2. $(1-u)y^2-u{x}^2-cos⟨b,c⟩y+2uxycos⟨a,b⟩+1=0$
      $(1-w)x^2-w{y}^2-cos⟨a,c⟩x+2wxycos⟨a,b⟩+1=0$
      由于我们知道 2D 点的图像位置，三个余弦角cos ⟨a,b⟩ ，cos ⟨b,c⟩ ，cos ⟨a,c⟩ 是已知的。同时， $u=\frac{B{C}^2}{A{B}^2}$，$w=\frac{A{C}^2}{A{B}^2}$ 可以通过A，B，C 在世界坐标系下的坐标算出，变换到相机坐标系下之后，并不改变这个比值。该式中的 x， y 是未知的，随着相机移动会发生变化。因此，该方程组是关于 x，y 的一个二元二次方程（多项式方程）。求解这个解析解比较麻烦，需要用吴消元法。该方程最多得到4个解，需要用验证点来验证最有可能的解。
   3. P3P 只利用三个点的信息。当给定的配对点多于 3 组时，难以利用更多的信息。如果 3D 点或 2D 点受噪声影响，或者存在误匹配，则算法失效。
   4. 在 SLAM 当中，通常的做法是先使用 P3P/EPnP 等方法估计相机位姿，然后构建最小二乘优化问题对估计值进行调整（Bundle Adjustment）。

3. Bundle Adjustment
   我们可以把 PnP 问题构建成一个定义于李代数上的非线性最小二乘问题,是一个最小化重投影误差（Reprojection error）
   

• 考虑 n 个三维空间点 $P$ 和它们的投影 $p$，我们希望计算相机的位姿 $R,t\text{R,t}R,t$，它的李代
  数表示为 $\xi \text{ξ}\xi$。假设某空间点坐标为 $P_i=[X_i{Y}_i,Z_i{]}^T$ ，其投影的像素坐标为 $u_i=[u_i,v_i{]}^T$ 。根据像素点和空间点投影位置关系可得：
  $s_i\left[\begin{array}{c}{u}_i\\ v_i\\ 1\end{array}\right]=K\text{K}Kexp(\xi \text{ξ}\xi$^$)\left[\begin{array}{c}{X}_i\\ Y_i\\ Z_i\\ 1\end{array}\right]$ ——>$s_i{u\text{u}u}_i=K\text{K}Kexp(\xi \text{ξ}\xi$^$){P\text{P}P}_i$
  中间隐含着的齐次坐标到非齐次的转换。由于相机位姿未知以及观测点的噪声，该等式存在一个误差,是将像素坐标与 3D 点按照当前估计的位姿进行投影得到的位置相比较得到的误差就是重投影误差
  $${\xi \text{ξ}\xi }^{\ast }=\underset{\xi }{\mathrm{argmin}}\frac{1}{2}\sum _{i=1}^n||{u\text{u}u}_i-\frac{1}{s_i}K\text{K}Kexp({\xi \text{ξ}\xi }^{>}){P\text{P}P}_i|{|}^2$$
  可以构建无约束的优化问题，通过 G-N, L-M 等优化算法进行求解

• 优化位姿：根据链式法则
  $$\begin{array}{rl}\frac{\partial e\text{e}e}{\partial \delta \xi \text{ξ}\xi }& =\underset{\delta \xi \to 0}{\lim }\frac{e\text{e}e(\delta \xi \text{ξ}\xi ⨁\xi \text{ξ}\xi )}{\delta \xi \text{ξ}\xi }=\frac{\partial e\text{e}e}{\partial {P\text{P}P}^{\prime }}\frac{\partial {P\text{P}P}^{\prime }}{\partial \delta \xi \text{ξ}\xi }\\ & \\ & =-\left[\begin{array}{cccccc}\frac{f_x}{Z^{\prime }}& 0& \frac{f_x{X}^{\prime }}{Z^{\prime 2}}& -\frac{f_x{X}^{\prime }{Y}^{\prime }}{Z^{\prime 2}}& f_x+\frac{f_x{X}^2}{Z^{\prime 2}}& -\frac{f_x{Y}^{\prime }}{Z^{\prime }}\\ 0& \frac{f_y}{Z^{\prime }}& -\frac{f_y{Y}^{\prime }}{Z^{\prime 2}}& -f_y-\frac{f_y{Y}^{\prime 2}}{Z^{\prime 2}}& \frac{f_y{X}^{\prime }{Y}^{\prime }}{Z^{\prime 2}}& \frac{f_y{X}^{\prime }}{Z^{\prime }}\end{array}\right]\end{array}$$

• 优化特征点：根据链式法则
  $$\begin{array}{rl}\frac{\partial e\text{e}e}{\partial P\text{P}P}& =\frac{\partial e\text{e}e}{\partial {P\text{P}P}^{\prime }}\frac{\partial {P\text{P}P}^{\prime }}{\partial P\text{P}P}\\ & \\ & =-\left[\begin{array}{ccc}\frac{f_x}{Z^{\prime }}& 0& -\frac{f_x{X}^{\prime }}{Z^{\prime 2}}\\ 0& \frac{f_y}{Z^{\prime }}& \frac{f_y{Y}^{\prime }}{Z^{\prime 2}}\end{array}\right]R\text{R}R\end{array}$$
  BA优化
  节点：第二个相机的位姿节点 $\xi \text{ξ}\xi \in SE(3)$，以及所有特征点的空间位置 $\mathbf{P}\in {\mathbb{R}}^3$
  边：每个 3D 点在第二个相机中的投影，以观测方程来描述：$z_j=h(\xi \text{ξ}\xi ,{P\text{P}P}_j)$
  

7.5 ICP(3D-3D)

假设有一组配对好的 3D 点 ：$P\text{P}P=${$p_1,...,p_n$}，${P\text{P}P}^{\prime }=${$p_1^{\prime },...,p_n^{\prime }$}
现在，想要找一个欧氏变换 $R,t\text{R,t}R,t$，使得：$\forall i,{p\text{p}p}_i={Rp\text{Rp}Rp}_i^{\prime }+t\text{t}t$
这个问题可以用迭代最近点（Iterative Closest Point, ICP）求解。

1. SVD方法（线性代数的求解）

   1. 首先，定义两组点的质心：$p\text{p}p=\frac{1}{n}{\sum }_{i=1}^n({p\text{p}p}_i)$ ， ${p\text{p}p}^{\prime }=\frac{1}{n}{\sum }_{i=1}^n({p\text{p}p}_i^{\prime })$
      然后计算每个点的去质心坐标：${q\text{q}q}_i={p\text{p}p}_i-p\text{p}p$，${q\text{q}q}_i^{\prime }={p\text{p}p}_i^{\prime }-{p\text{p}p}^{\prime }$
   2. 根据优化问题求解旋转矩阵
      $${R\text{R}R}^{\ast }=\underset{\mathbf{R}}{\mathrm{argmin}}\frac{1}{2}\sum _{i=1}^n||{q\text{q}q}_i-{Rq\text{Rq}Rq}_i^{\prime }|{|}^2$$
      先定义3×3矩阵 $\mathbf{W}$： $W\text{W}W={\sum }_{i=1}^n{q\text{q}q}_i{q\text{q}q}_i^{\prime T}$
      对 $\mathbf{W}$ 进行SVD分解:$W\text{W}W={U\sum V\text{U∑V}U\sum V}^T$
      其中， $\sum \text{∑}\sum$为奇异值组成的对角矩阵，对角线元素从大到小排列，而 $U\text{U}U$ 和 $V\text{V}V$ 为正交矩阵。当 $W\text{W}W$ 满秩时，$R\text{R}R={UV\text{UV}UV}^T$
   3. 根据$R\text{R}R$计算$t\text{t}t$：${t\text{t}t}^{\ast }=p\text{p}p-{Rp\text{Rp}Rp}^{\prime }$

2. 非线性优化方法
   目标函数为： $\underset{\xi \text{ξ}\xi }{\min }=\frac{1}{2}{\sum }_{i=1}^n||{p\text{p}p}_i-\exp (\xi$^ $){p\text{p}p}_i^{\prime }|{|}_2^2$
   单个误差项关于位姿导数：$\frac{\partial e\text{e}e}{\partial \delta \xi \text{ξ}\xi }=-(\exp (\xi$^$){p\text{p}p}_i^{\prime }{)}^{⨀}$

• 在RGB-D SLAM 中，由于一个像素的深度数据可能测量不到，所以我们可以混合着使用 PnP和 ICP 优化：对于深度已知的特征点，用建模它们的 3D-3D 误差；对于深度未知的特征点，则建模 3D-2D 的重投影误差。于是，可以将所有的误差放在同一个问题中考虑，使得求解更加方便。

8 视觉里程计(2)—直接法

• 特征点法的缺点

1. 关键点的提取与描述子的计算非常耗时
2. 使用特征点时，忽略了除特征点以外的所有信息
3. 相机有时会运动到特征缺失的地方，往往这些地方没有明显的纹理信息

• 解决以上问题的思路：

1. 保留特征点，但只计算关键点，不计算描述子。同时，使用光流法来跟踪特征点的运动。这样可以回避计算和匹配描述子带来的时间，但光流本身的计算需要一定时间；
2. 只计算关键点，不计算描述子。同时，使用直接法来计算特征点在下一时刻图像的位置。这同样可以跳过描述子的计算过程，而且直接法的计算更加简单。
3. 既不计算关键点、也不计算描述子，而是根据像素灰度的差异， 直接计算相机运动。

• 特征点法：最小化重投影误差——需要精确地知道空间点在两个相机中投影后的像素位置
• 直接法：最小化光度误差——根据像素的亮度信息，估计相机的运动，可以完全不用计算关键点和描述子，既避免了特征的计算时间，也避免了特征缺失的情况。只要场景中存在明暗变化（可以是渐变，不形成局部的图像梯度）。
• 根据使用像素的数量，直接法分为稀疏、稠密和半稠密三种。相比于特征点法只能重构稀疏特征点（稀疏地图），直接法还具有恢复稠密或半稠密结构的能力。

8.1 光流 ( Optical Flow )

光流法可以加速基于特征点的视觉里程计算法，避免计算和匹配描述子的过程，但要求相机运动较慢（或采集频率较高）

• 分类：

稀疏光流	计算部分像素运动	Lucas-Kanade
稠密光流	计算所有像素	Horn-Schunck

本节主要讲LK光流

• 一个在 t 时刻，位于 (x， y) 处的像素的灰度可以写成：$I\text{I}I(x,y,t)$
• 灰度不变假设：同一个空间点的像素灰度值，在各个图像中是固定不变的。
• 光流推导：
  • $I\text{I}I(x+dx,y+dy,t+dt)=I\text{I}I(x,y,t)$
  • 泰勒展开保留一阶：$I\text{I}I(x+dx,y+dy,t+dt)\approx I\text{I}I(x,y,t)+\frac{\partial I\text{I}I}{\partial x}dx+\frac{\partial I\text{I}I}{\partial y}dy+\frac{\partial I\text{I}I}{\partial t}dt$
  • 因为灰度不变，所以：$\frac{\partial I\text{I}I}{\partial x}dx+\frac{\partial I\text{I}I}{\partial y}dy+\frac{\partial I\text{I}I}{\partial t}dt=0$
  • 两边同除以$dt$,得：$\frac{\partial I\text{I}I}{\partial x}\frac{dx}{dt}+\frac{\partial I\text{I}I}{\partial y}\frac{dx}{dt}=-\frac{\partial I\text{I}I}{\partial t}$
  • 记像素在 x 轴上运动速度 $u=\frac{dx}{dt}$，像素在 y 轴上运动速度 $v=\frac{dy}{dt}$ ,${I\text{I}I}_x=\frac{\partial I\text{I}I}{\partial x}$ 为图像在该点处 x 方向的梯度，${I\text{I}I}_y=\frac{\partial I\text{I}I}{\partial y}$ 为图像在该点处 y 方向的梯度，${I\text{I}I}_t$为图像灰度对时间的变化量
    $\left[\begin{array}{cc}{I\text{I}I}_x& {I\text{I}I}_y\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right]=-{I\text{I}I}_t$ 该式是有两个变量的一次方程，无法解出$u,v$
  • 假设某个$w\times w$窗口内的像素具有相同的运动，含有$w^2$个像素，得到$w^2$个方程
    $\left[\begin{array}{c}[{I\text{I}I}_x,{I\text{I}I}_y{]}_1\\ .\\ .\\ [{I\text{I}I}_x,{I\text{I}I}_y{]}_k\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right]=-\left[\begin{array}{c}{I\text{I}I}_{t1}\\ .\\ .\\ {I\text{I}I}_{tk}\end{array}\right]$ ——>令$A\text{A}A=\left[\begin{array}{c}[{I\text{I}I}_x,{I\text{I}I}_y{]}_1\\ .\\ .\\ [{I\text{I}I}_x,{I\text{I}I}_y{]}_k\end{array}\right]$，$b\text{b}b=\left[\begin{array}{c}{I\text{I}I}_{t1}\\ .\\ .\\ {I\text{I}I}_{tk}\end{array}\right]$，得 $A\text{A}A\left[\begin{array}{c}u\\ v\end{array}\right]=-b\text{b}b$
  • 这是一个关于 u，v 的超定线性方程，传统解法是求最小二乘解。得到了像素在图像间的运动速度 u, v后当 t 取离散的时刻，就可以估计某块像素在若干个图像中出现的位置。
    ${\left[\begin{array}{c}u\\ v\end{array}\right]}^{\ast }=-({A\text{A}A}^{T}A\text{A}A{)}^{-1}{A\text{A}A}^{T}b\text{b}b$
• 角点具有更好的辨识度，边缘次之，区块最少
• LK 光流跟踪能够直接得到特征点的对应关系。这个对应关系就像是描述子的匹配，但实际上大多数时候只会碰到特征点跟丢的情况，而不太会遇到误匹配，这应该是光流相对于描述子的一点优势。但是，匹配描述子的方法在相机运动较大时仍能成功，而光流必须要求相机运动是微小的。从这方面来说，光流的鲁棒性比描述子差一些。
• 通过光流跟踪的特征点，用 PnP、 ICP 或对极几何来估计相机运动。

8.2 直接法（Direct Methods）

已知某个空间点 P 的世界坐标为 [X，Y，Z]和它在两个相机上成像的非齐次像素坐标${p\text{p}p}_1，{p\text{p}p}_2$。我们的目标是求第一个相机到第二个相机的相对位姿变换。

特征点法中，通过匹配描述子，知道${p\text{p}p}_1，{p\text{p}p}_2$的像素位置，就可以定位空间点然后计算重投影的位置。
直接法的思路是根据当前相机的位姿估计值，来寻找 ${p\text{p}p}_2$ 的位置。若相机位姿不够好， ${p\text{p}p}_2$ 的外观和 ${p\text{p}p}_1$ 会有明显差别（没有用投影，而是用的亮度）。为了减小这个差别，我们优化相机的位姿，来寻找与 ${p\text{p}p}_1$ 更相似的 ${p\text{p}p}_2$ 。

• 推导

  • 以第一个相机为参照系，设第二个相机旋转和平移为 $R,t\text{R,t}R,t$（对应李代数为 $\xi \text{ξ}\xi$）。相机的内参相同为 $K\text{K}K$，$Z_1$ 是 P 的深度， $Z_2$ 是 P 在第二个相机坐标系下的深度。投影方程：
    ${p\text{p}p}_1={\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]}_1={\frac{1}{Z}}_{1}KP\text{KP}KP$
    ${p\text{p}p}_2={\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]}_2={\frac{1}{Z}}_{2}K\text{K}K(RP+t\text{RP+t}RP+t)={\frac{1}{Z}}_{2}K\text{K}K(\exp (\xi \text{ξ}\xi$^$)P\text{P}P{)}_{1:3}$

  • 解一个优化问题———最小化光度误差，也就是 P 的两个像的亮度误差：$e={I\text{I}I}_1({p\text{p}p}_1)-{I\text{I}I}_2({p\text{p}p}_2)$——>$\underset{\xi \text{ξ}\xi }{\min }J\text{J}J(\xi \text{ξ}\xi )=||e|{|}^2$

  • 假设一个空间点在各个视角下，成像的灰度是不变的。N 个空间点 $P_i$的条件下整个相机位姿估计问题变为：
    $$\underset{\xi \text{ξ}\xi }{\min }J\text{J}J(\xi \text{ξ}\xi )=\sum _{i=1}^N{e}_i^T{e}_i，e_i={I\text{I}I}_1({p\text{p}p}_1,i)-{I\text{I}I}_2({p\text{p}p}_2,i)$$

  • 优化变量是相机位姿 $\xi \text{ξ}\xi$

  • $\frac{\partial e\text{e}e}{\partial \delta \xi \text{ξ}\xi }={\lim }_{\delta \xi \to 0}\frac{e\text{e}e(\delta \xi \text{ξ}\xi ⨁\xi \text{ξ}\xi )}{\delta \xi \text{ξ}\xi }$

  • $$\begin{array}{rl}e\text{e}e(\delta \xi \text{ξ}\xi ⨁\xi \text{ξ}\xi )& ={I\text{I}I}_1(\frac{1}{Z_1}KP\text{KP}KP)-{I\text{I}I}_2(\frac{1}{Z_2}K\text{K}K\exp (\delta {\xi }^{\wedge })\exp ({\xi }^{\wedge })P\text{P}P)\\ & \approx {I\text{I}I}_1(\frac{1}{Z_1}KP\text{KP}KP)-{I\text{I}I}_2(\frac{1}{Z_2}K\text{K}K(1+\delta {\xi }^{\wedge })\exp ({\xi }^{\wedge })P\text{P}P\\ & ={I\text{I}I}_1(\frac{1}{Z_1}KP\text{KP}KP)-{I\text{I}I}_2(\frac{1}{Z_2}K\text{K}K\exp ({\xi }^{\wedge })P\text{P}P+\frac{1}{Z_2}K\text{K}K\delta {\xi }^{\wedge }\exp ({\xi }^{\wedge })P\text{P}P)\end{array}$$

  • 设$q\text{q}q$ 为 $P\text{P}P$ 在扰动之后，位于第二个相机坐标系下的坐标，而 $u\text{u}u$ 为它的像素坐标
    $q\text{q}q=\delta \xi \text{ξ}\xi$^$\exp (\xi \text{ξ}\xi$ ^ $)P\text{P}P$，$u\text{u}u=\frac{1}{Z_2}Kq\text{Kq}Kq$

  • 一阶泰勒展开得
    $$\begin{array}{rl}e\text{e}e(\delta \xi \text{ξ}\xi ⨁\xi \text{ξ}\xi )& ={I\text{I}I}_1(\frac{1}{Z_1}KP\text{KP}KP)-{I\text{I}I}_2(\frac{1}{Z_2}K\text{K}K\exp ({\xi }^{\wedge })P\text{P}P+u\text{u}u)\\ & \approx {I\text{I}I}_1(\frac{1}{Z_1}KP\text{KP}KP)-{I\text{I}I}_2(\frac{1}{Z_2}K\text{K}K\exp ({\xi }^{\wedge })P\text{P}P)-\frac{\partial {I\text{I}I}_2}{\partial u\text{u}u}\frac{\partial u\text{u}u}{\partial q\text{q}q}\frac{\partial q\text{q}q}{\partial \delta \xi \text{ξ}\xi }\delta \xi \text{ξ}\xi \\ & =e(\xi \text{ξ}\xi )-\frac{\partial {I\text{I}I}_2}{\partial u\text{u}u}\frac{\partial u\text{u}u}{\partial q\text{q}q}\frac{\partial q\text{q}q}{\partial \delta \xi \text{ξ}\xi }\delta \xi \text{ξ}\xi \end{array}$$

  • $\frac{\partial {I\text{I}I}_2}{\partial u\text{u}u}$为 $u\text{u}u$ 处的像素梯度；

  • $\frac{\partial u\text{u}u}{\partial q\text{q}q}$为投影方程关于相机坐标系下的三维点的导数;

  • $\frac{\partial q\text{q}q}{\partial \delta \xi \text{ξ}\xi }$为变换后的三维点对变换的导数;

  • 后两项只与三维点 $\pm q$ 有关，而与图像无关，把它合并在一起：
    $\frac{\partial u\text{u}u}{\partial \delta \xi \text{ξ}\xi }=\left[\begin{array}{cccccc}\frac{f_x}{Z^{\prime }}& 0& \frac{f_x{X}^{\prime }}{Z^{\prime 2}}& -\frac{f_x{X}^{\prime }{Y}^{\prime }}{Z^{\prime 2}}& f_x+\frac{f_x{X}^2}{Z^{\prime 2}}& -\frac{f_x{Y}^{\prime }}{Z^{\prime }}\\ 0& \frac{f_y}{Z^{\prime }}& -\frac{f_y{Y}^{\prime }}{Z^{\prime 2}}& -f_y-\frac{f_y{Y}^{\prime 2}}{Z^{\prime 2}}& \frac{f_y{X}^{\prime }{Y}^{\prime }}{Z^{\prime 2}}& \frac{f_y{X}^{\prime }}{Z^{\prime }}\end{array}\right]$
    -所以误差相对于李代数的雅可比矩阵：
    $J\text{J}J=-\frac{\partial {I\text{I}I}_2}{\partial u\text{u}u}\frac{\partial u\text{u}u}{\partial \delta \xi \text{ξ}\xi }$

  • 对于 N 个点的问题，可以用这种方法计算优化问题的雅可比，然后使用 G-N 或L-M 计算增量，迭代求解

• 分类

  P 是一个已知位置的空间点，它是怎么来的呢？

1. P 来自于稀疏关键点，我们称之为稀疏直接法。通常我们使用数百个至上千个关键点，并且像 L-K 光流那样，假设它周围像素也是不变的。这种稀疏直接法不必计算描述子，并且只使用数百个像素，因此速度最快，但只能计算稀疏的重构。
2. P 来自部分像素。如果像素梯度为零，整一项雅可比就为零，不会对计算运动增量有任何贡献。因此，可以考虑只使用带有梯度的像素点，舍弃像素梯度不明显的地方。这称之为半稠密（Semi-Dense）的直接法，可以重构一个半稠密结构。
3. P 为所有像素，称为稠密直接法。稠密重构需要计算所有像素（一般几十万至几百万个），因此多数不能在现有的 CPU 上实时计算，需要 GPU 的加速。但是，如前面所讨论的，梯度不明显的点，在运动估计中不会有太大贡献，在重构时也会难以估计位置

8.3 直接法的讨论

• 相比于特征点法，直接法完全依靠优化来求解相机位姿。像素梯度引导着优化的方向。如果我们想要得到正确的优化结果，就必须保证大部分像素梯度能够把优化引导到正确的方向。
  

  • 在$I_1$中测的已知深度的$P$灰度值为229，又得到了一张新的图像，需要估计它的相机位姿。这个位姿是由一个初值不断地优化迭代得到的。假设我们的初值比较差，在这个初值下，空间点 P 投影后的像素灰度值是 126。于是，这个像素的误差为 229 − 126 = 103。为了减小这个误差，我们希望微调相机的位姿，使像素更亮一些。
  • 沿 $u$ 轴往前走一步，该处的灰度值变成了 123，即减去了 3。同样地，沿 $v$ 轴往前走一步，灰度值减 18，变成 108。在这个像素周围，我们看到梯度是 [−3，−18]，为了提高亮度，我们会建议优化算法微调相机，使 P 的像往左上方移动。这个过程中，我们用像素的局部梯度近似了它附近的灰度分布，不过请注意真实图像并不是光滑的，所以这个梯度在远处就不成立了。
  • 但是，优化算法不能只听这个像素的一面之词，还需要听取其他像素的建议。综合听取了许多像素的意见之后，优化算法选择了一个和我们建议的方向偏离不远的地方，计算出一个更新量 exp(ξ^)。加上更新量后，图像从 $I_2$ 移动到了 $I_2^{{}^{\prime }}$，像素的投影位置也变到了一个更亮的地方。通过这次更新，误差变小了。
  • 直接法的梯度是直接由图像梯度确定的，因此我们必须保证沿着图像梯度走时，灰度误差会不断下降。然而，图像通常是一个很强烈的非凸函数，很容易由于图像本身的非凸性（或噪声）落进一个局部极小值中，无法继续优化。只有当相机运动很小，图像中的梯度不会有很强的非凸性时，直接法才能成立。
  • 在例程中，我们只计算了单个像素的差异，并且这个差异是由灰度直接相减得到的。然而，单个像素没有什么区分性，周围很可能有好多像素和它的亮度差不多。所以，我们有时会使用小的图像块（patch），并且使用更复杂的差异度量方式，例如归一化相关性（ NCC）等。

  8.4 直接法的优缺点

  • 优点：
    • 可以省去计算特征点、描述子的时间
    • 只要求有像素梯度即可，无须特征点。可以在特征缺失的场合下使用，比较极端的例子是只有渐变的一张图像。它可能无法提取角点类特征，但可以用直接法估计它的运动
    • 可以构建半稠密乃至稠密的地图
  • 缺点：
    • 非凸性——直接法完全依靠梯度搜索，降低目标函数来计算相机位姿。其目标函数中需要取像素点的灰度值，而图像是强烈非凸的函数。这使得优化算法容易进入极小，只在运动很小时直接法才能成功。
    • 单个像素没有区分度. 找一个和他像的实在太多了！——于是我们要么计算图像块，要么计算复杂的相关性。由于每个像素对改变相机运动的“意见”不一致。只能少数服从多数，以数量代替质量
    • 灰度值不变是很强的假设. 如果相机是自动曝光的，当它调整曝光参数时，会使得图像整体变亮或变暗。光照变化时亦会出现这种情况。特征点法对光照具有一定的容忍性，而直接法由于计算灰度间的差异，整体灰度变化会破坏灰度不变假设，使算法失败。针对这一点，目前的直接法开始使用更细致的光度模型标定相机，以便在曝光时间变化时也能让直接法工作