在SLAM的视觉里程计中，比较常用的就是特征点法和直接法。而直接法中，光流则是其中的重点内容，比如LSD-SLAM中就使用到了光流的方法。本文将会就光流的理论原理、公式推导进行详细的剖析，以帮助读者深刻地理解。

光流算法

光流是关于视域中的物体运动检测中的概念，它用来描述相对于观察者的运动所造成的观测目标、表面或边缘。

简单来讲，光流描述了场景中物体运动在视觉中的变化。光流的概念，由Gibson在1950年提出，其通过相邻帧之间像素点的对应关系计算像素点的瞬时速度，从而描述物体信息。

为了应用光流计算物体运动法，相邻帧需要满足两个假设：

• 灰度不变：两个间隔帧之间描述相同目标的像素点的灰度需要保持恒定；
• 运动较小：两个间隔帧之间描述相同目标的像素点并没有发生剧烈变化；

假设图像上一个像素点$(x,y)$，它在时刻$t$的亮度为$I(x,y,t)$，用$u(x,y)$和$v(x,y)$表示该点光流在水平和垂直方向上的速度分量：

$$u=\frac{\mathrm{d}x}{\mathrm{d}t}$$

$$u=\frac{\mathrm{d}y}{\mathrm{d}t}$$

在经过时间间隔$dt$之后，该点的对应点的亮度变为$I(x+dx,y+dy,t+dt)$.

根据假设一：

$$I(x+dx,y+dy,t+dt)=I(x,y,t)$$

根据假设二，利用泰勒公式展开：

$$I(x+dx,y+dy,t+dt)\approx I(x,y,t)+\frac{\partial I}{\partial x}dx+\frac{\partial I}{\partial y}dy+\frac{\partial I}{\partial t}dt$$

根据上面两个式子，可以得到：

$$\frac{\partial I}{\partial x}dx+\frac{\partial I}{\partial y}dy+\frac{\partial I}{\partial t}dt=0$$

当$dt$足够小，趋向于0时：

$$\frac{\partial I}{\partial x}\frac{dx}{dt}+\frac{\partial I}{\partial y}\frac{dy}{dt}+\frac{\partial I}{\partial t}=0$$

$$-\frac{\partial I}{\partial t}=\frac{\partial I}{\partial x}\frac{dx}{dt}+\frac{\partial I}{\partial y}\frac{dy}{dt}=\frac{\partial I}{\partial x}u+\frac{\partial I}{\partial y}v$$

$$-I_t=I_x\frac{dx}{dt}+I_y\frac{dy}{dt}=I_{x}u+I_{y}v$$

$$-I_t=\left[\begin{array}{cc}{I}_x& I_y\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right]$$

其中，$I_x$和$I_y$为图像在$(x,y)$处的梯度，$I_t$为图像在$(x,y)$处关于时间$t$的导数，$I_t\approx I(x,y,t)-I(x,y,t-1)$。

这就是基本的光流约束方程。

不过对于某一个像素点而言，存在两个未知量$(u,v)$，但只有一个方程是无法求解的。

先小结一下光流的两个假设：

1. 灰度不变：图像场景中的目标的像素在帧间运动时灰度保持不变。用于得到光流法基本方程；
2. 运动较小：保证灰度可以对位置求偏导（换句话说，小运动情况下我们才能用前后顿之间单位位置变化引起的灰度变化去近似灰度对位置的偏导数)，这也是光流法不可或缺的假定。

LK光流

LK光流在基本的光流算法基础上增加了一个假设：

• 邻域相同：每个像素点的局部邻域内的所有像素点的光流保持一致。

LK光流证法一

由于基本的光流约束方程中，一个方程两个未知量没法求解。LK光流算法考虑到了像素点的邻域，将问题转变成了计算某些点集的光流，联立多个方程，从而解决了这个问题。

假设在邻域$R$内$n$个像素点，光流保持一致，可以得到：

$$-I_{tk}=\left[\begin{array}{cc}{I}_{xk}& I_{yk}\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right],k\in R$$

即：

$$\left[\begin{array}{cc}{I}_{x1}& I_{y1}\\ I_{x2}& I_{y2}\\ ...& ...\\ I_{xn}& I_{yn}\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right]=\left[\begin{array}{c}-I_{t1}\\ -I_{t2}\\ ...\\ -I_{tn}\end{array}\right]$$

为了简单起见，记作：

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

这是一个超正定方程，一种解法就是计算其最小二乘解。即，当$A^{T}A$可逆时：

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

带入：

$$\left[\begin{array}{c}u\\ v\end{array}\right]={\left[\begin{array}{cc}{\sum }_{i=1}^n{I}_{xi}^2& {\sum }_{i=1}^n{I}_{xi}{I}_{yi}\\ {\sum }_{i=1}^n{I}_{xi}{I}_{yi}& {\sum }_{i=1}^n{I}_{yi}^2\end{array}\right]}^{-1}\left[\begin{array}{c}-{\sum }_{i=1}^n{I}_{xi}{I}_t\\ -{\sum }_{i=1}^n{I}_{yi}{I}_t\end{array}\right]$$

一般来说，在图像中沿两个方向都有像素变化的区域（例如角点），其对应的$A^{T}A$是可逆的；相反的，在灰度变化很小的区域，一般对应的$A^{T}A$是不可逆的。

这限制了LK光流法的使用范围，这也是被称为稀疏光流法的主要原因。

LK光流证法二

假设相邻的两幅图像$I$和$J$，对于$I$中像素点$u={\left[\begin{array}{cc}{u}_x& u_y\end{array}\right]}^T$，需要在图像$J$中找到像素点$v=u+d={\left[\begin{array}{cc}{u}_x+d_x& u_y+d_y\end{array}\right]}^T$，使其最为相似（灰度值差別最小）。我们把$d={\left[\begin{array}{cc}{d}_x& d_y\end{array}\right]}^T$称为在$u$处的光流。

当然为了求解，需要引入一个具有相同光流的邻域，其大小使用$w_x$、$w_y$两个参数表示。则求解$d$转换为使下述目标函数取最小值的优化问题：

$$\epsilon (d)=\epsilon (d_x,d_y)=\sum _{x=u_x-w_x}^{u_x+w_x}\sum _{y=u_y-w_y}^{u_y+w_y}(I(x,y)-J(x+d_x,y+d_y){)}^2$$

在最优解处导数为0，则得到：

$$\frac{\partial \epsilon (d)}{\partial d}=-2\sum _{x=u_x-w_x}^{u_x+w_x}\sum _{y=u_y-w_y}^{u_y+w_y}(I(x,y)-J(x+d_x,y+d_y))\left[\begin{array}{cc}\frac{\partial J}{\partial x}& \frac{\partial J}{\partial y}\end{array}\right]$$

使用一阶泰勒展开，得到：

$$\frac{\partial \epsilon (d)}{\partial d}\approx -2\sum _{x=u_x-w_x}^{u_x+w_x}\sum _{y=u_y-w_y}^{u_y+w_y}(I(x,y)-J(x,y)-\left[\begin{array}{cc}\frac{\partial J}{\partial x}& \frac{\partial J}{\partial y}\end{array}\right]\left[\begin{array}{c}{d}_x\\ d_y\end{array}\right])\left[\begin{array}{cc}\frac{\partial J}{\partial x}& \frac{\partial J}{\partial y}\end{array}\right]$$

其中，$I(x,y)-J(x,y)$是像素点$(x,y)$关于时间的图像导数，记作$\delta I$；$\left[\begin{array}{cc}\frac{\partial J}{\partial x}& \frac{\partial J}{\partial y}\end{array}\right]$是图像$J$在像素点$(x,y)$处的梯度，将其近似为$\left[\begin{array}{cc}\frac{\partial I}{\partial x}& \frac{\partial I}{\partial y}\end{array}\right]$，并记$\nabla I={\left[\begin{array}{cc}\frac{\partial I}{\partial x}& \frac{\partial I}{\partial y}\end{array}\right]}^T$。

$$\frac{1}{2}\frac{\partial \epsilon (d)}{\partial d}\approx \sum _{x=u_x-w_x}^{u_x+w_x}\sum _{y=u_y-w_y}^{u_y+w_y}(\nabla I^{T}d-\delta I)\nabla I^T$$

需要注意，$\nabla I^{T}d-\delta I$为标量，所以：

$$\frac{1}{2}{\left[\begin{array}{c}\frac{\partial \epsilon (d)}{\partial d}\end{array}\right]}^T\approx \sum _{x=u_x-w_x}^{u_x+w_x}\sum _{y=u_y-w_y}^{u_y+w_y}(\nabla I^{T}d-\delta I)\nabla I$$

$$\frac{1}{2}{\left[\begin{array}{c}\frac{\partial \epsilon (d)}{\partial d}\end{array}\right]}^T\approx \sum _{x=u_x-w_x}^{u_x+w_x}\sum _{y=u_y-w_y}^{u_y+w_y}\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]d-\sum _{x=u_x-w_x}^{u_x+w_x}\sum _{y=u_y-w_y}^{u_y+w_y}\left[\begin{array}{c}\delta I\cdot I_x\\ \delta I\cdot I_y\end{array}\right]$$

记：

$$G=\sum _{x=u_x-w_x}^{u_x+w_x}\sum _{y=u_y-w_y}^{u_y+w_y}\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]$$

$$b=\sum _{x=u_x-w_x}^{u_x+w_x}\sum _{y=u_y-w_y}^{u_y+w_y}\left[\begin{array}{c}\delta I\cdot I_x\\ \delta I\cdot I_y\end{array}\right]$$

则：

$$d=G^{-1}b$$

两种证法

仔细观察两种证法的结果，可以看出最终得到的结果是一样的。现在梳理一下光流计算中，需要用到的值：

• $I_x$和$I_y$：图像在$(x,y)$处的梯度；

• $\delta I$或$I_t$：图像在$(x,y)$处对时间的变化。

迭代求解LK

有时候我们为了求得更加精确的$d$，可以使用迭代求解的方式更新$d$。

使用$K(K\ge 1)$代表迭代次数，对于第$k$次迭代，前面的第$k-1$次迭代用于提供初始值$d^{k-1}={\left[\begin{array}{cc}{d}_x^{k-1}& d_y^{k-1}\end{array}\right]}^T$，这里把它作为图像$J$移动的初始值，移动之后得到图像$J_k$，且$J_k(x,y)=J(x+d_x^{k-1},y+d_y^{k-1})$，之后的问题就是计算优化变量${\eta }^k={\left[\begin{array}{cc}{\eta }_x^k& {\eta }_y^k\end{array}\right]}^T$，使下面的目标函数最小：

$${\epsilon }^k({\eta }^k)=\epsilon ({\eta }_x^k,{\eta }_y^k)=\sum _{x=u_x-w_x}^{u_x+w_x}\sum _{y=u_y-w_y}^{u_y+w_y}(I(x,y)-J_k(x+{\eta }_x^k,y+{\eta }_y^k){)}^2$$

根据上文的第二种证法的推导形式，最终可以整理为：

$${\eta }^k=G^{-1}{b}_k$$

注意：这里的$G$在迭代计算中保持不变，因此只需要计算一次。$b_k$每次需要重新计算，之后第$k$次迭代的结果$d_k=d_{k-1}+{\eta }^k$。

当迭代次数达到设定次数，或者计算得到的${\eta }^k$的大小小于设定阀值时，迭代计算结束。

金字塔LK

这里重新把LK中的三条假设拿出来：

1. 同一个空间点的像素灰度值，在各个图像中固定不变；
2. 相邻帧间物体的运动足够小；
3. 邻域具有相同的光流。

这三条假设都是强假设，但实际场景中尤其是第二条很难满足。但是对于我们一直在用的一阶泰勒展开，只有在变量变化很小的情况下才会成立。但如果帧间像素点的运动比较大时，泰勒展开就不能适用了。

原理

金字塔LK主要就是用来解决帧间变化较大的问题。

举例来说，对于分辨率为$800\times 600$的图像，若两帧间某个像素点的运动为$40\times 40$，当把图像分辨率改为$400\times 300$时，运动变为$20\times 20$。以此类推，总能在图像像素降低到一定程度时，原本较大的像素间运动变得足够小，能够使用泰勒展开式进行计算。

把原始图像作为第$0$层，宽高缩小$2^L$倍的图像作为第$L$层，多层图像组织起来，就像是一个金字塔。

其中计算时需要注意一些小问题：

1. 在图像降采样之前通常使用低通滤波器滤波，防止降采样后发生锯齿现象，其中就会碰到对于图像边界像素如何滤波的问题；
2. 原始图像尺寸无法保证宽高严格满足整除$2^L$，因此降采样后新的图像的宽高如何确定，有不同的解决方法。可以对原始图像调整宽高满足整除$2^L$，或者用下面的公式确定新图像的宽高：

$$n_x^L\le \frac{n_x^{L-1}+1}{2}$$

$$n_y^L\le \frac{n_y^{L-1}+1}{2}$$

其中，$n_x^L$和$n_y^L$分别代表第$L$层图像的宽和高。

接下来，就是利用金字塔图像从最高的第$L$层图像，开始向下迭代计算每层中对应像素点的位移。其实这部分操作和迭代求解LK部分的原理相同。

使用$I^L$、$J^L$代表图像$I$和$J$的第$L$层金字塔图像，假设$u={\left[\begin{array}{cc}{u}_x& u_y\end{array}\right]}^T$为图像$I$中需要跟踪的像素点坐标，则在$I^L$中对应的坐标$u^L={\left[\begin{array}{cc}{u}_x^L& u_y^L\end{array}\right]}^T=u/2^L$。

首先在最高层$L_m$层计算光流大小，把它作为第$L_m-1$层光流的初始值，再计算该层的光流大小；之后把得到的光流大小作为$L_m-2$层光流大小的初始值。以此类推，直至第$0$层图像，得到光流的最终结果。

以第$L+1$层到第$L$层之间的计算为例，详细说明下计算流程。

假设通过第$L+1$层计算得到第$L$层光流大小的初始值为$g^L={\left[\begin{array}{cc}{g}_x^L& g_y^L\end{array}\right]}^T$，那么我们要做的就是求得$d^L=\left[\begin{array}{cc}{d}_x^L& d_y^L\end{array}\right]$，使
下面的优化函数最小：

$${\epsilon }^k(d^k)=\epsilon (d_x^k,d_y^k)=\sum _{x=u_x-w_x}^{u_x+w_x}\sum _{y=u_y-w_y}^{u_y+w_y}(I^L(x,y)-J^L(x+g_x^k+d_x^k,y+g_y^k+d_y^k){)}^2$$

由于$g^L={\left[\begin{array}{cc}{g}_x^L& g_y^L\end{array}\right]}^T$已经比较接近真实值了，所以这里我们很有底气地适用泰勒展开处理$J^L(x+g_x^k+d_x^k,y+g_y^k+d_y^k)$，后面的处理和常规LK类似，不再赘述。

计算得到$d^L=\left[\begin{array}{cc}{d}_x^L& d_y^L\end{array}\right]$后，定义$g^{L-1}=2(g^L+d^L)$作为第$L-1$层光流大小的初始值，重复上述计算过程逐层迭代计算即可。特别的，在最高层$L_m$处，定义$g^{L_m}={\left[\begin{array}{cc}0& 0\end{array}\right]}^T$。

通过观察我们可以发现，最终计算得到的光流$d={\sum }_{L=0}^{L_m}{2}^L{d}^L$。这个公式表明虽然原始图像中像素位移比较大，但对于每层图像来说对应的$d^L$都比较小，因此可以使用泰勒展开简化求解。

这里的金字塔LK的流程就有点像：准确值=估计值+残差。估计值由上一层金字塔提供，残差在本层中适用迭代法求出最佳的残差。估计值的初始值为顶层金字塔光流${\left[\begin{array}{cc}0& 0\end{array}\right]}^T$，残差在每层金字塔都会初始化为${\left[\begin{array}{cc}0& 0\end{array}\right]}^T$。

整体流程

现在把金字塔LK光流算法的整体流程使用伪代码梳理一遍。

目标：对于图像$I$中的像素点$u$，寻找在图像$J$中灰度最相似的像素点。

Step1

为图像$I$和$J$分别建立金字塔表达 { I L } L = 0 , 1 , . . . , L m \{I^{L}\}_{L=0,1,...,L_{m}} {IL}L=0,1,...,Lm​​和 { J L } L = 0 , 1 , . . . , L m \{J^{L}\}_{L=0,1,...,L_{m}} {JL}L=0,1,...,Lm​​

Step2

初始化最高层光流初始值$g^{L_m}={\left[\begin{array}{cc}0& 0\end{array}\right]}^T$

$forL=L_m;L\ge 0;L--$

        计算图像$I^L$中像素点$u$对应的位置：$u^L=\left[\begin{array}{cc}{u}_x^L& u_y^L\end{array}\right]=u/2^L$

        计算图像$I^L$在$x$方向地梯度：$I_x^L(x,y)=\frac{I^L(x+1,y)-I^L(x-1,y)}{2}$

        计算图像$I^L$在$y$方向地梯度：$I_y^L(x,y)=\frac{I^L(x,y+1)-I^L(x,y-1)}{2}$

        计算空间矩阵$G={\sum }_{x=u_x^L-w_x}^{u_x^L+w_x}{\sum }_{y=u_y^L-w_y}^{u_y^L+w_y}\left[\begin{array}{cc}{I}_x^{L^2}& I_x^L{I}_y^L\\ I_x^L{I}_y^L& I_y^{L^2}\end{array}\right]$

        初始化迭代LK初始值：$d^0={\left[\begin{array}{cc}0& 0\end{array}\right]}^T$

        $fork=1;k\le K||fabsf(\eta )<threshold;k++$

                计算图像差异：$\delta I^k(x,y)=I^L(x,y)-J^L(x+g_x^L+d_x^{k-1},y+g_y^L+d_y^{k-1})$

                计算匹配差异向量：$b_k={\sum }_{x=u_x^L-w_x}^{u_x^L+w_x}{\sum }_{y=u_y^L-w_y}^{u_y^L+w_y}\left[\begin{array}{c}\delta I^k(x,y)\cdot I_x(x,y)\\ \delta I^k(x,y)\cdot I_y(x,y)\end{array}\right]$

                计算光流：${\eta }^k=G^{-1}{b}_k$

                为下一次迭代提供初始值：$d^k=d^{k-1}+{\eta }^k$

         结束$k$上的迭代

         得到第$L$层图像上的光流优化值：$d^L=d^k$

         为第$L-1$层光流提供光流初始值：$g^{L-1}=2(g^L+d^L)$

结束$L$上的迭代

最终计算得到光流$d=g^0+d^0$

图像$J$中对应坐标$v=u+d$

补充

1. 实际计算中得到的类似于$(x+g_x^L+d_x^{k-1},y+g_y^L+d_y^{k-1})$往往不是整数点像素，所以其像素灰度值需要插值计算得到，一般可以用双线性插值计算；
2. 在对图像进行类似于卷积操作时（求梯度、高斯滤波等），需要注意位于图像边界处像素点的处理，这里会发生像素越界，需要采用一些方法处理；
3. 使用光流法计算匹配点时，当$v=u+d$超出图像$J$时，我们判定为匹配失败：不过即使$v=u+d$位于图像$J$中，也未必是正确的匹配点，毕竟光流法的前提假设过于理想。且该方法中未提供依据判断匹配是否成功，需要其他的方法来滤除该方法计算得到的误匹配，一种常用的方法是基于RANSAC计算两帧图像的变换矩阵。

应用：基于中值流的目标跟踪

TLD算法的跟踪模块借鉴了中值流算法，用于上一帧目标$bb$到当前帧目标$bb$的预测。

中值流算法的流程如下：

1. 初始化第$t$帧中的目标$bb$，在$bb$内均匀取点用于跟踪。TLD的作者在目标$bb$内均匀选取了100个像素点作为跟踪点；
2. 利用LK光流法确定第$t$帧与第$t+1$帧之间的光流情况。用预测到的第$t+1$帧中这些点的位置，反向利用LK光流法得到这些点的后向跟踪轨迹。这样，每个跟踪点都会在第一帧中由前向和后向得到两个位置，这两个位置的欧式距离就作为每个点的$FBerror$；
3. 计算每个跟踪点在两帧之间的匹配度（即计算$ncc$的值）；
4. 利用所有的跟踪点预测$bb$的变化有失妥当，因为有的跟踪点可能并不对$bb$的变化起作用，甚至会影响$bb$的预测。所以需要筛选得到一部分可靠的点。这也就是中值流算法的核心部分。将这些跟踪点在相邻两顿之间的$FBerror$和$ncc$值进行排序，得到中值，利用中值，得到同时小于$FBerror$中值与$ncc$中值的特征点（分别过滤$50\%$跟踪点），这样保留不到$50\%$的点作为跟踪成功的点；
5. 利用剩下的这些点来预测上一顿的目标$bb$在当前帧的位置和大小。

位置（位移估计）：用保留的这些点的$(x,y)$位移的中值作为目标$bb$的位移。

大小（尺度估计）：计算上一帧中保留的这些点之间的欧氏距离$d1$，当前帧中保留的这些点之间的欧氏距离$d2$。并以此得到伸缩比的中值$s$。$s=median(d2/d1)$，这个伸缩比的中值作为目标$bb$尺度伸缩值。

OpenCV接口

CV_EXPORTS_W void calcOpticalFlowPyrLK( InputArray prevImg, InputArray nextImg,
                           InputArray prevPts, CV_OUT InputOutputArray nextPts,
                           OutputArray status, OutputArray err,
                           Size winSize=Size(21,21), int maxLevel=3,
                           TermCriteria criteria=TermCriteria(
                            TermCriteria::COUNT+TermCriteria::EPS,
                            30, 0.01),
                           double derivLambda=0.5,
                           int flags=0,
                           double minEigThreshold=1e-4);

各个参数代表的含义如下：

• prevImg
  你的标定图像的灰度图
• nextImg
  你想搜寻的图像的灰度图
• prevPts
  输入的标定图像的特征点（可以是其他特征点检测方法找到的点）
• nextPts
  输出场景的特征点
• status
  输出状态向量（无符号char），如果在当前图像中能够光流得到标定的特征点位置改变，则设置status的对应位置为1，否则设置为0
• err
  输出错误向量；向量的每个元素被设为相应特征的一个错误，误差测量的类型可以在flags参数中设置；如果流不被发现然后错误未被定义（使用status（状态）参数找到此情形）。
• winSize
  在每个金字塔水平搜寻窗口的尺寸
• maxLevel
  金字塔的高度，初始为3层

相关阅读

• Pyramidal Implementation of the Lucas Kanade Feature Tracker Description of the algorithm
• 金字塔LK光流法数学原理学习笔记
• 总结：光流–LK光流–基于金字塔分层的LK光流–中值流
• 图像金字塔(高斯金字塔 与 拉普拉斯金字塔)
• 多元函数的泰勒(Taylor)展开式
• 光流金字塔calcOpticalFlowPyrLK进行特征点跟踪