个人博客：http://www.chenjianqu.com/
原文链接：http://www.chenjianqu.com/show-89.html
光流(Optical Flow)是一种描述像素随时间在图像之间运动的方法。随着时间的流逝，同一个像素会在图像中运动，我们希望追踪它的运动过程。其中，计算部分像素运动的称为稀疏光流，计算所有像素的称为稠密光流。稀疏光流以Lucas Kanade光流为代表，并可以在SLAM中用于跟踪特征点位置。
Lucas-Kanade光流
Lucas–Kanade光流算法是一种两帧差分的光流估计算法。在LK光流中，认为来自相机的图像是随时间变化的。图像可以看作时间的函数: I(t)，在t时刻，位于(x,y)处的像素的灰度可以写成 I(x,y,t) 。这种方式把图像看成了关于位置与时间的函数，它的值域就是图像中像素的灰度。现在考虑某个固定的空间点，它在t时刻的像索坐标为x,y。由于相机的运动，它的图像坐标将发生变化。我们希望估计这个空间点在其他时刻图像中位置。LK光流法有三个基本假设：
1. 灰度不变：一个像素点随着时间的变化，其亮度值（像素灰度值）是恒定不变的。这是光流法的基本设定。所有光流法都必须满足。灰度不变假设是很强的假设， 实际当中很可能不成立。事实上，由于物体的材质不同，像素会出现高光和阴影部分；有时，相机会自动调整曝光参数，使得图像整体变亮或变暗。这此时候灰度不变假设都是不成立的，因此光流的结果也不一定可靠。
2. 小运动： 时间的变化不会引起位置的剧烈变化。这样才能利用相邻帧之间的位置变化引起的灰度值变化，去求取灰度对位置的偏导数。所有光流法必须满足。
3. 空间一致：假设某一个窗口内的像素具有相同的运动。这是LK光流法独有的假定。因为为了求取x,y方向的速度，需要建立多个方程联立求解。而空间一致假设就可以利用邻域n个像素点来建立n个方程。
对于t时刻位于(x,y)处的像素，我们设t+dt时刻它运动到(x+dx,y+dy)处。由于灰度不变，有:
I (x+dx, y+dy, t+dt) = I (x,y,t).
对上式的左边进行泰勒展开，保留一阶项，得：
I (x+dx, y+dy, t+dt) ≈ I (x,y,t) + (∂I/∂x)dx + (∂I/∂y)dy + (∂I/∂t)dt
由灰度不变，即下一个时刻的灰度等于之前的灰度，从而：
(∂I/∂x)dx + (∂I/∂y)dy + (∂I/∂t)dt = 0
两边除以dt，得：
(∂I/∂x)(dx/dt) + (∂I/∂y)(dy/dt) = - ∂I/∂t
dx/dt是像素在x轴上的运动速度，而dy/dt为y轴上的速度，记为u,v。这也是我们想知道的变量。
∂I/∂x 是图像在该点处的x方向的梯度，∂I/∂y是y方向的梯度，记为Ix,Iy。图像梯度的计算。
∂I/∂t 是图像灰度对时间的变化量，记为It ，即两帧之间的该点的灰度变化量。
原式写成矩阵形式为：
我们想要计算的是像素的运动u,v，上式有两个变量但只有一个方程，因此引入空间一致假设。考虑一个大小为w*w的窗口，该窗口内的像素具有同样的运动，因此我们共有w*w个方程：
于是整个方程为：
这是一个关于u,v的超定线性方程，可用最小二乘法解得：
这样就得到了像素在图像间的运动速度u,v。t取离散时刻，可以估计某块像素在若干图像中出现的位置。
改进的LK光流
原始的LK光流假设:灰度不变、小运动、空间一致都是较强的假设，并不容易得到满足。考虑物体的运动速度较大时，算法会出现较大的误差，那么我们希望能减少图像中物体的运动速度。假设当图像为400×400时，物体速度为[16 16],那么图像缩小为200×200时，速度变为[8,8]。缩小为100*100时，速度减少到[4,4]。在源图像缩放后，原算法又变得适用了。所以光流可以通过生成原图像的金字塔图像，逐层求解，不断精确来求得。简单来说上层金字塔（低分辨率）中的一个像素可以代表下层的两个。每一层的求解结果乘以2后加到下一层。主要的步骤有三步：建立金字塔，基于金字塔跟踪，迭代过程。
OpenCV代码
OpenCV中使用CalcOpticalFlowPyrLK()这个函数实现LK光流，参数如下：
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), int flags=0, double minEigThreshold=1e-4); prevImg – 第一个8位输入图像或者通过 buildOpticalFlowPyramid()建立的金字塔 nextImg – 第二个输入图像或者和prevImg相同尺寸和类型的金字塔 prevPts – 二维点向量存储找到的光流；点坐标必须是单精度浮点数 nextPts – 输出二维点向量（用单精度浮点坐标）包括第二幅图像中计算的输入特征的新点位置；当OPTFLOW_USE_INITIAL_FLOW 标志通过，向量必须有和输入一样的尺寸。 status – 输出状态向量（无符号char）；如果相应的流特征被发现，向量的每个元素被设置为1，否则，被置为0. err – 输出错误向量；向量的每个元素被设为相应特征的一个错误，误差测量的类型可以在flags参数中设置；如果流不被发现然后错误未被定义（使用status（状态）参数找到此情形）。 winSize – 在每个金字塔水平搜寻窗口的尺寸。 maxLevel – 基于最大金字塔层次数。如果设置为0，则不使用金字塔（单级）；如果设置为1，则使用两个级别，等等。如果金字塔被传递到input，那么算法使用的级别与金字塔同级别但不大于MaxLevel。 criteria – 指定迭代搜索算法的终止准则（在指定的最大迭代次数标准值（criteria.maxCount）之后，或者当搜索窗口移动小于criteria.epsilon。） flags – 操作标志，可选参数： OPTFLOW_USE_INITIAL_FLOW – 使用初始估计，存储在nextPts中；如果未设置标志，则将prevPts复制到nextPts并被视为初始估计。 OPTFLOW_LK_GET_MIN_EIGENVALS – 使用最小本征值作为误差度量（见minEigThreshold描述）；如果未设置标志，则将原始周围的一小部分和移动的点之间的 L1 距离除以窗口中的像素数，作为误差度量。 minEigThreshold – 算法所计算的光流方程的2x2标准矩阵的最小本征值（该矩阵称为[Bouguet00]中的空间梯度矩阵）÷ 窗口中的像素数。如果该值小于MinEigThreshold，则过滤掉相应的特征，相应的流也不进行处理。因此可以移除不好的点并提升性能。
C++代码如下：
CMakeLists.txt
cmake_minimum_required(VERSION 2.6)
project(lk_flow)
set( CMAKE_BUILD_TYPE Release )
set( CMAKE_CXX_FLAGS "-std=c++11 -O3" )
find_package( OpenCV )
include_directories( ${OpenCV_INCLUDE_DIRS} )
add_executable(lk_flow main.cpp)
target_link_libraries( lk_flow ${OpenCV_LIBS} )
install(TARGETS lk_flow RUNTIME DESTINATION bin)
main.cpp
#include <iostream>
#include <fstream>
#include <list>
#include <vector>
using namespace std;
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/features2d/features2d.hpp>
#include <opencv2/video/tracking.hpp>
using namespace cv;
int main( int argc, char** argv )
{
//特征点 因为要删除跟踪失败的点，使用list
list<Point2f> keypoints;
Mat color, last_color;
color = imread("1.png");
Mat showFlowImg = imread("9.png");
// 对第一帧提取FAST特征点
vector<KeyPoint> kps;
Ptr<FastFeatureDetector> detector = FastFeatureDetector::create();
detector->detect(color, kps);
for(auto kp:kps)
keypoints.push_back( kp.pt );
last_color = color;
//进行光流跟踪
for ( int index=2; index<10; index++ ){
//读取图片文件
color = imread(to_string(index)+".png");
vector<Point2f> next_keypoints;
vector<Point2f> prev_keypoints;
for(auto kp:keypoints)
prev_keypoints.push_back(kp);
//用LK跟踪特征点
vector<unsigned char> status;
vector<float> error;
calcOpticalFlowPyrLK(last_color, color, prev_keypoints, next_keypoints, status, error );
// 把跟丢的点删掉，把新的位置赋给keypoints
int i=0;
for(auto iter=keypoints.begin(); iter!=keypoints.end(); i++){
if(status[i] == 0 ){
iter = keypoints.erase(iter);
continue;
}
//画跟踪的线条
Point2f beforePoint=*iter;
Point2f afterPoint=next_keypoints[i];
line(showFlowImg, beforePoint,afterPoint, Scalar(0,240,0), 1);
*iter = next_keypoints[i];
iter++;
}
cout<<"tracked keypoints: "<<keypoints.size()<<endl;
//如果全部跟丢了
if (keypoints.size() == 0){
cout<<"all keypoints are lost."<<endl;
break;
}
// 画出 keypoints
Mat img_show = color.clone();
for (auto kp:keypoints)
circle(img_show, kp, 5, Scalar(0, 240, 0), 1);
imshow("corners", img_show);
waitKey(0);
last_color = color;
}
imshow("showFlowImg", showFlowImg);
waitKey(0);
return 0;
}
特征点：
特征点的轨迹：
程序运行结果显示，图像中大部分特征点能够顺利跟踪到，但某些特征点会丢失。丢失的特征点或是被移出视野外，或是被其它物体挡住，如果不提取新的特征点，那么光流的跟踪会越来越少。
位于物体角点处的特征更加稳定，边缘处的特征会沿着边缘滑动，因为沿着边缘移动时特征块的内容基本不变，因此被程序认为是同一个地方。而其它地方的特征点则会频繁跳动。因此最好提取的是角点，其次是边缘点。
LK光流点不需要计算和匹配描述子，但是本身也需要一定的计算量。另外LK光流跟踪能直接得到特征点的对应关系，不太会误匹配，但是光流必需要求相机的运动是微小的。
参考文献
[0]高翔.视觉SLAM14讲
[1]一度逍遥.光流法详解之一(LK光流）.https://www.cnblogs.com/riddick/p/10586662.html .2019-03-25
[2]菜鸟知识搬运工.OpenCV3学习（11.2）LK光流法原理及opencv实现. https://blog.****.net/qq_30815237/article/details/87208319 .2019-02-13