0 简介

在学习vins-mono过程中，算法前端采用基于KLT光流的跟踪方法，对光流和KLT光流进行简单的总结。包括LK光流跟踪算法原理，基于金字塔改进的LK光流跟踪，KLT光流跟踪算法，以及在光流跟踪时使用的角点检测算法Harris角点检测和Shi-Tomasi角点检测。

光流跟踪分为稠密光流和稀疏光流，稠密光流即对图像中所有的像素点进行跟踪，稀疏光流则是对图像中的特征点进行跟踪。KLT光流跟踪是一种稀疏光流跟踪算法。

1 光流的三大假设条件

光流法跟踪的的主要原理就是通过计算像素点在相邻帧图像之间的速度，可以得到像素在下一帧的位置，具体计算方法则是根据光流法的三大假设联立方程组，基于最小二乘对像素在两帧之间的运动进行求解。

在vins系列中，算法前端采用KLT光流跟踪算法，函数入口为calcOpticalFlowPyrLK

光流跟踪法的三大假设条件如下所示：

• 亮度恒定即灰度不变性：同一个空间点的像素灰度值，在各个图像中是固定不变的
• 小运动即时间连续性：相邻帧间保持微小的运动，时间的变化不会引起像素的剧烈变化
• 空间一致性：特征点及其周围邻域的像素均具有相同的运动

2 LK光流跟踪算法原理

(此处参考14讲第8章)
对于t时刻的某个像素(x,y)，当相机运动到t+1时刻时，基于灰度不变性假设可以得到如下关系

对公式左边进行一阶泰勒展开，并根据灰度不变性可以得到以下公式：

首先dx/dt和dy/dt分别为像素点在x轴和y轴方向上的导数，分别记为u和v；
然后∂I/∂x和∂I/∂y分别为像素值在(x,y)处的偏导数，记为在(x,y)处的梯度；
则可以将上述公式写成矩阵的形式：

为了求解像素在相邻两帧的速度(u,v)，根据空间一致性假设同一像素块中的像素具有相同的运动，则可以在特征点像素周围选择一个邻域范围，该范围内的像素和特征点具有相同的运动，则可以联立方程组然后根据最小二乘法对变量(u,v)进行求解。

3 LK-金字塔光流

这部分的参考链接：https://blog.csdn.net/App_12062011/article/details/51880258
当相机运动速度过快时，相邻帧图像之间就会出现较大的位移，这时小运动这一较强的假设就不再成立，由于同一特征点在相邻图像中的位移Δ_pixel较大，在进行迭代求解时很容易陷入局部最优的情况。

金字塔光流则是根据这一问题提出的改进措施。
通过对图像进行resize，在进行跟踪时首先在分辨率较小的金字塔图像之间进行跟踪。因为对图像进行了下采样，如果说在原图像中像素点在800x800的图像上的运动速度为[8,8]，那么在400x400的图像上的运动速度可能会接近[4x4]，在200x200的图像上速度接近[2x2]，此时在分辨率为200x200又可以使用LK光流进行求解。在求解得到像素在金字塔顶层图像的速度后，会将该速度和在像素点处的梯度信息作为初值继续计算特征点在下一层的速度和梯度信息，直到求解出特征点在原图像的速度信息。

现在假设某一特征点，分别在图像I和I时刻的坐标为x_I(x,y)、x_J(x+dx,y+dy)，dx=(dx,dy)即为图像在x处的速度。实际上，I时刻的图像到J时刻的图像之间经历的一个仿射变换，因此对于光流跟踪问题可以转化为一个求解一个速度dx=[dx,dy]和一个仿射变换矩阵A，使得特征点邻域内的像素灰度值之差最小，仿射变换矩阵的定义如下：

3.1 金字塔创建

一般金字塔为4层，这里涉及到一个精度和实时性之间的权衡，金字塔层数越高则需要进行光流计算的次数越多，实时性越差。

3.2 金字塔跟踪

这部分的参考链接：https://blog.csdn.net/App_12062011/article/details/51880258

• 建立金字塔后，首先对于最高层的图像计算出其上的光流；
• 通过上一层的计算结果对下一层的的图像进预平移，并在L层的基础上计算出该层的残余光流向量dL，简单来说下一层像素点的位置=初始图像下一层的位置+上一层光流+上一层与下一层之间的位置差值
• 对于每一层迭代光流计算，最终得到的光流也即是所有层光流的叠加。

4 KLT光流跟踪

LK光流法旨在说明对图像像素的一种光流跟踪算法，根据三大假设联立方程组最终将问题转换为求解最小二乘的问题，即AX=b的形式。由于Hessian矩阵的特征值都比较大时才能满足空间一致性的假设条件，即在选择角点作为关键点时符合这一情况。

实际上，一方面像素移动并不会那么简单，另一方面窗口内像素并不都是同样的移动方式，因为这样的近似必然会带来误差。而现在的问题就是如何去选择合适的窗口，或者特征点，从而获得最为精确的跟踪。KLT角点检测方法就是为了选择一个适合跟踪的特征点，它认为一个好的特征点的定义应当就是能被好的跟踪。

在opencv中集成的goodFeaturesToTrack函数即使用的Shi-Tomasi角点检测。

5 Harris角点检测

参考链接：https://zhuanlan.zhihu.com/p/150411153

算法的核心是利用局部窗口在图像上进行移动，判断灰度是否发生较大的变化。如果窗口内的灰度值（在梯度图上）都有较大的变化，那么这个窗口所在区域就存在角点。

1.Harris角点检测步骤:

• 当窗口（局部区域）同时向 x （水平）和 y（垂直） 两个方向移动时，计算窗口内部的像素值变化量 E(x,y) ；
• 对于每个窗口，都计算其对应的一个角点响应函数 R；
• 然后对该函数进行阈值处理，如果 R>threshold，表示该窗口对应一个角点特征。

2.Harris角点数学模型
建立数学模型的目的，是为了确定哪些窗口会引起较大的灰度值变化。对于窗口中某个点(x,y)其灰度值为I(x,y)，当其移动一个(u,v)的距离后，其对应灰度值为I(x+u,y+v)，I(x+u,y+v) - I(x,y)即为窗口在移动是该像素灰度值的变化。

为了更为清晰的确定窗口内的像素灰度变化，选择计算窗口内所有像素灰度值的变化。如果窗口中心点像素是角点，那么窗口移动前后，中心点的灰度值变化非常强烈，所以该点权重系数应该设大一点，表示该点对灰度变化的贡献较大；而离窗口中心（角点）较远的点，这些点的灰度变化比较小，于是将权重系数设小一点，表示该点对灰度变化的贡献较小。

因此这里引入窗口函数w(x,y)(其实就是一个高斯核)，以计算整个窗口中像素灰度值得变化。设窗口在移动前后的灰度值变化为E(u,v)，则当窗口中心像素为角点时，E(u,v)会变得比较大。

为简化计算，将I(x+u,y+v)进行一阶泰勒展开。

其中对应的矩阵M为：

3.Harris角点求解
现在定义了Harris角点的数学模型，下一步就是如何去寻找角点的坐标以及存在角点的对应的窗口。

（1）通过M矩阵的特征值判断
如果特征值在两个方向都比较大就是角点；如果特征值在一个方向比较大是边；如果特征值都很小就是平面。

可以把R看成旋转因子，其不影响两个正交方向的变化分量。经对角化处理后，将两个正交方向的变化分量提取出来，就是 λ1 和 λ2（特征值）。

（2）通过R响应判断
本质上也是根据M矩阵的特征值计算一个分数R，根据分数进行判断:R大于某个阈值，认为是角点；R大于0但小于阈值是图像块；R小于0是边缘。
R的计算方式如下：

Harris 角点检测的结果是带有这些分数R的灰度图像，设定一个阈值，分数大于这个阈值的像素就对应角点。

6 Shi-Tomasi角点检测

参考链接：https://blog.csdn.net/qq_40374812/article/details/117016214
Shi-Tomasi角点是Harris角点的改进版。

Harris角点检测算法引入角点响应值R，根据R值的大小来判断是否存在角点和边界。

Shi-Tomasi角点检测算法就是针对R值的计算公式进行了优化和修改。由于 Harris角点检测算法的稳定性和α值有关，而α是个经验值，需要我们根据自身经验去设置，不好设定最佳值。

Shi-Tomasi发现，角点的稳定性其实和矩阵M的较小特征值有关(即λ1，λ2)，于是直接用较小的那个特征值作为分数。这样就不用调整α值了。

所以 Shi-Tomasi 将R值的计算公式改为如下形式：R = min(λ1，λ2)

和 Harris 一样，如果该分数大于设定的阈值，我们就认为它是一个角点。