简单来说,光流法保留特征点,但只计算关键点,不计算描述子。这样可以回避计算和匹配描述子带来的时间,二光流计算本身的时间要少于计算描述子需要的时间。

光流法

  1. 假设条件

    **灰度不变假设:**同一个空间点的像素灰度值,在各个图像中是固定不变的。

  2. 作用

    Untitled

    光流是一种描述像素随着时间,在图像之间运动的方法,如图所示。随着时间的经过,同一个像素会在图像中运动,而我们希望追踪它的运动过程。即得到像素的运动 $(u, v)$。

  3. 分类

    计算部分像素运动的称为稀疏光流,计算所有像素的称为稠密光流。稀疏光流以 Lucas-Kanade 光流为代表,并可以在 SLAM 中用于跟踪特征点位置。

Lucas-Kanade 光流

在 LK 光流中,我们认为来自相机的图像是随时间变化的。图像可以看作时间的函数:$I(t)$。那么,一个在 $t$ 时刻,位于 $(x, y)$ 处的像素,它的灰度 $\boldsymbol { I } ( x , y , t )$。

  1. 推导

    对于 $t$ 时刻位于 $(x, y)$ 处的像素,我们设 $t + \mathrm{d}t$ 时刻,它运动到 $(x + \mathrm{d}x, y + \mathrm{d}y)$ 处。 由于灰度不变,我们有:

    $$ \boldsymbol { I } ( x + \mathrm { d } x , y + \mathrm { d } y , t + \mathrm { d } t ) = \boldsymbol { I } ( x , y , t ) $$

    对左边进行泰勒展开,保留一阶项:

    $$ \boldsymbol { I } ( x + \mathrm { d } x , y + \mathrm { d } y , t + \mathrm { d } t ) \approx \boldsymbol { I } ( x , y , t ) + \frac { \partial \boldsymbol { I } } { \partial x } \mathrm {~d} x + \frac { \partial \boldsymbol { I } } { \partial y } \mathrm {~d} y + \frac { \partial \boldsymbol { I } } { \partial t } \mathrm {~d} t $$

    由于灰度不变假设:

    $$ \frac { \partial \boldsymbol { I } } { \partial x } \mathrm {~d} x + \frac { \partial \boldsymbol { I } } { \partial y } \mathrm {~d} y + \frac { \partial \boldsymbol { I } } { \partial t } \mathrm {~d} t = 0 $$

    两边除以 $\mathrm{d}t$,得:

    $$ \frac { \partial \boldsymbol { I } } { \partial x } \frac { \mathrm { d } x } { \mathrm {~d} t } + \frac { \partial \boldsymbol { I } } { \partial y } \frac { \mathrm { d } y } { \mathrm {~d} t } = - \frac { \partial \boldsymbol { I } } { \partial t } $$

    令:

    $$ u :=\mathrm { d } x / \mathrm { d } t, v:=\mathrm { d } y / \mathrm { d } t \\ \boldsymbol { I } _ { x } := \frac { \partial \boldsymbol { I } } { \partial x } , \boldsymbol { I } _ { y } :=\frac { \partial \boldsymbol { I } } { \partial y }, \boldsymbol { I } _ { t } :=\frac { \partial \boldsymbol { I } } { \partial t } $$

    则:

    $$ \left[ \begin{array} { l l } \boldsymbol { I } _ { x } & \boldsymbol { I } _ { y } \end{array} \right] \left[ \begin{array} { c } u \\ v \end{array} \right] = - \boldsymbol { I } _ { t } $$

    由于该式是带有两个变量的一次方程,仅凭它无法计算出 $u, v$。因此,必须引入额外的约束来计算 $u, v$。在 LK 光流中,我们假设某一个窗口内的像素具有相同的运动。

    考虑一个大小为 $w × w$ 大小的窗口,它含有 $w^2$ 数量的像素。由于该窗口内像素具有 同样的运动,因此我们共有 $w^2$ 个方程:

    $$ \left[ \begin{array} { l l } \boldsymbol { I } _ { x } & \boldsymbol { I } _ { y } \end{array} \right] _ { k } \left[ \begin{array} { c } u \\ v \end{array} \right] = - \boldsymbol { I } _ { t k } , \quad k = 1 , \ldots , w ^ { 2 } $$

    令:

    $$ \boldsymbol { A } = \left[ \begin{array} { c } { \left[ \boldsymbol { I } _ { x } , \boldsymbol { I } _ { y } \right] _ { 1 } } \\ \vdots \\ { \left[ \boldsymbol { I } _ { x } , \boldsymbol { I } _ { y } \right] _ { k } } \end{array} \right] , \boldsymbol { b } = \left[ \begin{array} { c } \boldsymbol { I } _ { t 1 } \\ \vdots \\ \boldsymbol { I } _ { t k } \end{array} \right] $$

    则:

    $$ \boldsymbol { A } \left[ \begin{array} { l } u \\ v \end{array} \right] = - \boldsymbol { b } $$

    由最小二乘法:

    $$ \left[ \begin{array} { l } u \\ v \end{array} \right] ^ { * } = - \left( \boldsymbol { A } ^ { T } \boldsymbol { A } \right) ^ { - 1 } \boldsymbol { A } ^ { T } \boldsymbol { b } $$

多层光流

我们把光流写成了优化问题,就必须假设优化的初始值靠近最优值,才能在一定程度上保障算法的收敛。如果相机运动较快,两张图像差异较明显,那么单层图像光流法容易达到一个局部极小值。这种情况可以通过引入图像金字塔来改善。

图像金字塔是指对同一个图像进行缩放,得到不同分辨率下的图像,如图82所示。以原始图像作为金字塔底层,每往上一层,就对下层图像进行一定倍率的缩放,就得到了一个金字塔。

然后,在计算光流时,先从顶层的图像开始计算,然后把上一层的追踪结果,作为下一层光流的初始值。由于上层的图像相对粗糙,所以这个过程也称为由粗至精(Coarse-to-fine)的光流,也是实用光流法的通常流程。

直接法

单层光流法

考虑某个空间点 $P$ 和两个时刻的相机,$P$ 的世界坐标为 $[X,Y,Z]$,它在两个相机平面上成像的像素坐标为 $p_1, p_2$。

Untitled

我们的目标是求第一个相机到第二个相机的相对位姿变换。我们以第一个相机为参照系,设第二个相机的旋转和平移为 $R,t$(对应李群为$T$)。同时,两相机的内参相同,记为 $K$。为清楚起见,我们列写完整的投影方程:

$$ \begin{aligned} \boldsymbol { p } _ { 1 } = & { \left[ \begin{array} { c } u \\ v \\ 1 \end{array} \right] _ { 1 } = \frac { 1 } { Z _ { 1 } } \boldsymbol { K } \boldsymbol { P } } \\ \boldsymbol { p } _ { 2 } = & { \left[ \begin{array} { l } u \\ v \\ 1 \end{array} \right] _ { 2 } = \frac { 1 } { Z _ { 2 } } \boldsymbol { K } ( \boldsymbol { R } \boldsymbol { P } + \boldsymbol { t } ) = \frac { 1 } { Z _ { 2 } } \boldsymbol { K } ( \boldsymbol { T } \boldsymbol { P } ) _ { 1 : 3 } } \end{aligned} $$

其中 $Z_1$ 是P的深度,$Z_2$ 是P在第二个相机坐标系下的深度。由于 $T$ 需要与其次坐标相乘,因此乘完之后需要取出前 $3$ 个元素。

此时,我们将光度误差 $e = \boldsymbol { I } _ { 1 } \left( \boldsymbol { p } _ { 1 } \right) - \boldsymbol { I } _ { 2 } \left( \boldsymbol { p } _ { 2 } \right)$ 作为最小化的目标,因此,目标函数可以写为:

$$ \min _ { T } J ( \boldsymbol { T } ) = \| e \| ^ { 2 } $$

当我们有 $N$ 个空间点 $P_i$,则相机位姿估计问题可以写为:

$$ \min _ { \boldsymbol { T } } J ( \boldsymbol { T } ) = \sum _ { i = 1 } ^ { N } e _ { i } ^ { \mathrm { T } } e _ { i } , \quad e _ { i } = \boldsymbol { I } _ { 1 } \left( \boldsymbol { p } _ { 1 , i } \right) - \boldsymbol { I } _ { 2 } \left( \boldsymbol { p } _ { 2 , i } \right) $$

我们定义两个中间变量:

$$ \begin{aligned} \boldsymbol { q } & = \boldsymbol { T } \boldsymbol { P } \\ \boldsymbol { u } & = \frac { 1 } { Z _ { 2 } } \boldsymbol { K } \boldsymbol { q } \end{aligned} $$

这里的 $q$ 为 $P$ 在第二个相机坐标系下的坐标,$u$ 为其在相机平面上的像素坐标。显然,$q$ 是 $T$ 的函数,$u$ 是 $q$ 的函数,从而也是 $T$ 的函数。考虑李代数的左扰动模型,利用一阶泰勒展开:

$$ e ( \boldsymbol { T } ) = \boldsymbol { I } _ { 1 } \left( \boldsymbol { p } _ { 1 } \right) - \boldsymbol { I } _ { 2 } ( \boldsymbol { u } ) $$

$$ \frac { \partial e } { \partial \boldsymbol { T } } = \frac { \partial \boldsymbol { I } _ { 2 } } { \partial \boldsymbol { u } } \frac { \partial \boldsymbol { u } } { \partial \boldsymbol { q } } \frac { \partial \boldsymbol { q } } { \partial \delta \boldsymbol { \xi } } \delta \boldsymbol { \xi } $$

其中 $\xi$ 为 $T$ 的左扰动。

  1. $\partial I_2 / \partial u$ 为 $u$ 处的像素梯度。

  2. $\partial u / \partial q$ 为投影方程关于相机坐标系下的三维点的导数,记 $q = [X,Y,Z]^T$,则导数为

    $$ \frac { \partial \boldsymbol { u } } { \partial \boldsymbol { q } } = \left[ \begin{array} { c c c } \cfrac { \partial u } { \partial X } & \cfrac { \partial u } { \partial Y } & \cfrac { \partial u } { \partial Z } \\ \cfrac { \partial v } { \partial X } & \cfrac { \partial v } { \partial Y } & \cfrac { \partial v } { \partial Z } \end{array} \right] = \left[ \begin{array} { c c c } \cfrac { f _ { x } } { Z } & 0 & - \cfrac { f _ { x } X } { Z ^ { 2 } } \\ 0 & \cfrac { f _ { y } } { Z } & - \cfrac { f _ { y } Y } { Z ^ { 2 } } \end{array} \right] $$

  3. $\partial \boldsymbol { q } / \partial \delta \boldsymbol { \xi }$ 为变换后的三位点对变换的导数

    $$ \frac { \partial \boldsymbol { q } } { \partial \delta \boldsymbol { \xi } } = \left[ \boldsymbol { I } , - \boldsymbol { q } ^ { \wedge } \right] $$

在实践中,由于后两项只与三维点 $q$ 有关,而与图像无关,我们经常把它合并在一起:

$$ \frac { \partial \boldsymbol { u } } { \partial \delta \boldsymbol { \xi } } = \left[ \begin{array} { c c c c c c } \frac { f _ { x } } { Z } & 0 & - \frac { f _ { x } X } { Z ^ { 2 } } & - \frac { f _ { x } X Y } { Z ^ { 2 } } & f _ { x } + \frac { f _ { x } X ^ { 2 } } { Z ^ { 2 } } & - \frac { f _ { x } Y } { Z } \\ 0 & \frac { f _ { y } } { Z } & - \frac { f _ { y } Y } { Z ^ { 2 } } & - f _ { y } - \frac { f _ { y } Y ^ { 2 } } { Z ^ { 2 } } & \frac { f _ { y } X Y } { Z ^ { 2 } } & \frac { f _ { y } X } { Z } \end{array} \right] $$

于是,我们推导出误差相对于李代数的雅可比矩阵:

$$ \boldsymbol { J } = - \frac { \partial \boldsymbol { I } _ { 2 } } { \partial \boldsymbol { u } } \frac { \partial \boldsymbol { u } } { \partial \delta \boldsymbol { \xi } } $$

对于N个点的问题,我们可以用这种方法计算优化问题的雅可比矩阵,然后使用高斯牛顿法或列文伯格一马夸尔特方法计算增量,迭代求解。

直接法的特点

最后,我们总结直接法的优缺点。大体上,它的优点如下:

  1. 可以省去计算特征点、描述子的时间。
  2. 只要求有像素梯度即可,不需要特征点。因此,直接法可以在特征缺失的场合下使用。比较极端的例子是只有渐变的一幅图像。它可能无法提取角点类特征,但可以用直接法估计它的运动。在演示实验中,我们看到直接法对随机选取的点亦能正常工作。这一点在实用中非常关键,因为实用场景很有可能没有很多角点可供使用。
  3. 可以构建半稠密乃至稠密的地图,这是特征点法无法做到的。

它的缺点也很明显:

  1. 非凸性。直接法完全依靠梯度搜索,降低目标函数来计算相机位姿。其目标函数中需要取像素点的灰度值,而图像是强烈非凸的函数。这使得优化算法容易进入极小,只在运动很小时直接法才能成功。针对于此,金字塔的引入可以在一定程度上减小非凸性的影响。
  2. 单个像素没有区分度。和它像的实在太多了!于是我们要么计算图像块,要么计算复杂的相关性。由于每个像素对改变相机运动的“意见”不一致,只能少数服从多数,以数量代替质量。所以,直接法在选点较少时的表现下降明显,我们通常建议用500个点以上。
  3. 灰度值不变是很强的假设。如果相机是自动曝光的,当它调整曝光参数时,会使得图像整体变亮或变暗。光照变化时也会出现这种情况。特征点法对光照具有一定的容忍性,而直接法由于计算灰度间的差异,整体灰度变化会破坏灰度不变假设,使算法失败。针对这一点,实用的直接法会同时估计相机的曝光参数0,以便在曝光时间变化时也能工作。