OpenCV3学习（11.2）LK光流法原理及opencv实现

光流的概念:(Optical flow or optic flow)
它是一种运动模式，这种运动模式指的是一个物体、表面、边缘在一个视角下由一个观察者（比如眼睛、摄像头等）和背景之间形成的明显移动。光流技术，如运动检测和图像分割，时间碰撞，运动补偿编码，三维立体视差，都是利用了这种边缘或表面运动的技术。光流，简单说也就是画面移动过程中，图像上每个像素的x，y位移量，比如第t帧的时候A点的位置是(x1, y1)，那么我们在第t+1帧的时候再找到A点，假如它的位置是(x2,y2)，那么我们就可以确定A点的运动了：(u, v) = (x2, y2) - (x1,y1)

光流算法：
它评估了两幅图像的之间的变形，它的基本假设是体素和图像像素守恒。它假设一个物体的颜色在前后两帧没有巨大而明显的变化。基于这个思路，我们可以得到图像约束方程。不同的光流算法解决了假定了不同附加条件的光流问题。

 

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

但是基本的光流约束方程的约束只有一个，无法求出下，x,y方向的速度u,v两个未知量，LK光流算法考虑到像素的邻域，将问题转变为计算某些点集的光流，通过联立多个方程，从而解决了这个问题。

LK光流算法

     在计算机视觉中，Lucas–Kanade光流算法是一种两帧差分的光流估计算法。这个算法是最常见，最流行的。它计算两帧在时间t 到t + δt之间每个像素点位置的移动。 由于它是基于图像信号的泰勒级数，这种方法称为差分，这就是对于空间和时间坐标使用偏导数。

LK光流算法的假设：

          

1、假设原图是I(x,y,z)  (这里是扩展到三维空间的，所以还有个z值)，移动后的图像是I(x+δx,y+δy,z+δz,t+δt),两者满足:
 

根据假设1得到：图像移动可以认为I (x ,y ,z ,t ) = I (x + δx ,y + δy ,z + δz ,t + δt )，H.O.T. 指更高阶，在移动足够小的情况下可以忽略。从这个方程中我们可以得到：

               

        或者

3、从这个方程中我们可以得到：

                                                 

其中= u， =v，Vx,Vy,Vz 分别是I(x,y,z,t)的光流向量中x，y，z的组成（二维图像没有z），   则   是图像在(x ,y,z ,t )这一点的梯度  （就是两帧图像块之间差值） 。

4、这个方程有三个未知量，尚不能被解决，这也就是所谓光流算法的光圈问题。那么要找到光流向量则需要另一套解决的方案。而Lucas-Kanade算法是一个非迭代的算法：假设流(Vx,Vy,Vz)在一个大小为m*m*m(m>1)的小窗中是一个常数，那么从像素1...n , n = m*m*m 中可以得到下列一组方程： 

三个未知数但是有多于三个的方程，这个方程组自然是个超定方程，也就是说方程组内有冗余，方程组可以表示为：

也就是：

采用最小二乘法：

得到：

 

其中的求和是从1到n。LK光流法是一种稀疏光流法，它假定光流在像素的局部邻域内是常数，利用最小二乘法对邻域所有像素求解基本光流方程，求解两帧图像在时间t和t+1之间像素点位置的移动。

这个算法的不足在于它不能产生一个密度很高的流向量，例如在运动的边缘和黑大的同质区域中的微小移动方面流信息会很快的褪去。它的优点在于有噪声存在的鲁棒性还是可以的。

Lucas-Kanade改进算法

为什么要用金字塔？因为lk算法的约束条件即：小速度，亮度不变以及区域一致性都是较强的假设，并不很容易得到满足。如当物体运动速度较快时，假设不成立，那么后续的假设就会有较大的偏差，使得最终求出的光流值有较大的误差。考虑物体的运动速度较大时，算法会出现较大的误差。那么就希望能减少图像中物体的运动速度。一个直观的方法就是，缩小图像的尺寸。假设当图像为400×400时，物体速度为[16 16],那么图像缩小为200×200时，速度变为[8,8]。缩小为100*100时，速度减少到[4,4]。所以在源图像缩放了很多以后，原算法又变得适用了。所以光流可以通过生成原图像的金字塔图像，逐层求解，不断精确来求得。简单来说上层金字塔（低分辨率）中的一个像素可以代表下层的两个。每一层的求解结果乘以2后加到下一层。

 

对于Lucas-Kanade改进算法来说，主要的步骤有三步：建立金字塔，基于金字塔跟踪，迭代过程。

1、首先对每一帧图像建立图像金字塔，原始图像在底部，最低分辨率在顶层；

2、金字塔跟踪
      总体来讲，金字塔特征跟踪算法描述如下：首先，光流和仿射变换矩阵在最高一层的图像上计算出；将上一层的计算结果作为初始值传递给下一层图像，这一层的图像在这个初始值的基础上，计算这一层的光流和仿射变化矩阵；再将这一层的光流和仿射矩阵作为初始值传递给下一层图像，直到传递给最后一层，即原始图像层，这一层计算出来的光流和仿射变换矩阵作为最后的光流和仿射变换矩阵的结果。

 

                   

 
 

from：https://blog.csdn.net/sgfmby1994/article/details/68489944

from：https://blog.csdn.net/longlovefilm/article/details/79824723

在OpenCV中主要有个基于图像金字塔的LK算法的函数CalcOpticalFlowPyrLK（）：

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，则过滤掉相应的特征，相应的流也不进行处理。因此可以移除不好的点并提升性能。

实例：

#include<opencv2\opencv.hpp>
#include<iostream>
#include <ctype.h>
#include <opencv2\video\tracking.hpp>
#include <opencv2\highgui\highgui.hpp>
#include <opencv2\imgproc\imgproc.hpp>


using namespace std;
using namespace cv;

/*void tracking(Mat &frame, Mat &output);
bool addNewPoints();
bool acceptTrackedPoint(int i);

string window_name = "optical flow tracking";
Mat gray;	// 当前图片
Mat gray_prev;	// 前一张图
vector<Point2f> points[2];	// point0为特征点原位置，point1为新位置,[2]表示包含两个vector<Point2f>元素的数组
vector<Point2f> initial;	// 初始化跟踪点的位置
vector<Point2f> features;	// 检测的特征
int maxCount = 500;		// 检测的最大特征数
double qLevel = 0.01;	// 特征检测的等级
double minDist = 10.0;	// 两特征点之间的最小距离
vector<uchar> status;	// 跟踪特征的状态，特征的流发现为1，否则为0
vector<float> err;
int main()
{
	Mat frame;
	Mat result;

	VideoCapture capture("jianpan.avi");

	if (capture.isOpened())	// 摄像头读取文件开关
	{
		while (true)
		{
			capture >> frame;

			if (!frame.empty())
			{
				tracking(frame, result);
			}
			else
			{
				printf(" --(!) No captured frame -- Break!");
				break;
			}

			int c = waitKey(100);
			if ((char)c == 27)
			{
				break;
			}
		}
	}
	return 0;
	//destroyAllWindows();
}

//
// function: tracking
// brief: 跟踪
// parameter: frame	输入的视频帧
//			  output 有跟踪结果的视频帧
// return: void
//
void tracking(Mat &frame, Mat &output)
{
	cvtColor(frame, gray, CV_BGR2GRAY);
	frame.copyTo(output);
	// 添加特征点
	if (addNewPoints())
	{
		goodFeaturesToTrack(gray, features, maxCount, qLevel, minDist);
		points[0].insert(points[0].end(), features.begin(), features.end());
		initial.insert(initial.end(), features.begin(), features.end());
	}

	if (gray_prev.empty())
	{
		gray.copyTo(gray_prev);
	}
	// l-k光流法运动估计
	calcOpticalFlowPyrLK(gray_prev, gray, points[0], points[1], status, err);//status标志位
	// 去掉一些不好的特征点
	int k = 0;
	for (size_t i = 0; i<points[1].size(); i++)
	{
		if (acceptTrackedPoint(i))
		{
			initial[k] = initial[i];
			points[1][k++] = points[1][i];
		}
	}
	points[1].resize(k);//取前k个元素
	initial.resize(k);
	// 显示特征点和运动轨迹
	for (size_t i = 0; i<points[1].size(); i++)
	{
		line(output, initial[i], points[1][i], Scalar(0, 0, 255));
		circle(output, points[1][i], 3, Scalar(255, 0, 0), -1);
	}

	// 把当前跟踪结果作为下一此参考
	swap(points[1], points[0]);//把跟踪到的点points[1]放到初始点points【0】
	swap(gray_prev, gray);//把当前放到gray_prev,作为前一张图，与新的当前图gray比较

	imshow(window_name, output);
}

//
// function: addNewPoints
// brief: 检测新点是否应该被添加
// parameter:
// return: 是否被添加标志
//
bool addNewPoints()
{
	return points[0].size() <= 10;
}

//
// function: acceptTrackedPoint
// brief: 决定哪些跟踪点被接受
// parameter:
// return:
//
bool acceptTrackedPoint(int i)
{
	return status[i] && ((abs(points[0][i].x - points[1][i].x) + abs(points[0][i].y - points[1][i].y)) > 2);
}*/



static void help()
{
	// print a welcome message, and the OpenCV version
	cout << "\nThis is a demo of Lukas-Kanade optical flow lkdemo(),\n"
		"Using OpenCV version " << CV_VERSION << endl;
	cout << "\nIt uses camera by default, but you can provide a path to video as an argument.\n";
	cout << "\nHot keys: \n"
		"\tESC - quit the program\n"
		"\tr - auto-initialize tracking\n"
		"\tc - delete all the points\n"
		"\tn - switch the \"night\" mode on/off\n"
		"To add/remove a feature point click it\n" << endl;
}

Point2f point;
bool addRemovePt = false;

static void onMouse(int event, int x, int y, int /*flags*/, void* /*param*/)
{
	if (event == CV_EVENT_LBUTTONDOWN)
	{
		point = Point2f((float)x, (float)y);
		addRemovePt = true;
	}
}

int main(int argc, char** argv)
{
	help();

	VideoCapture cap;
	TermCriteria termcrit(CV_TERMCRIT_ITER | CV_TERMCRIT_EPS, 20, 0.03);
	Size subPixWinSize(10, 10), winSize(31, 31);

	const int MAX_COUNT = 500;
	bool needToInit = false;
	bool nightMode = false;

	if (argc == 1 || (argc == 2 && strlen(argv[1]) == 1 && isdigit(argv[1][0])))
		cap.open(argc == 2 ? argv[1][0] - '0' : 0);
	else if (argc == 2)
		cap.open(argv[1]);

	if (!cap.isOpened())
	{
		cout << "Could not initialize capturing...\n";
		return 0;
	}

	namedWindow("LK Demo", 1);
	setMouseCallback("LK Demo", onMouse, 0);

	Mat gray, prevGray, image, frame;
	vector<Point2f> points[2];

	for (;;)
	{
		cap >> frame;
		if (frame.empty())
			break;

		frame.copyTo(image);
		cvtColor(image, gray, COLOR_BGR2GRAY);

		if (nightMode)
			image = Scalar::all(0);

		if (needToInit)
		{
			// automatic initialization
			goodFeaturesToTrack(gray, points[1], MAX_COUNT, 0.01, 10, Mat(), 3, 0, 0.04);
			cornerSubPix(gray, points[1], subPixWinSize, Size(-1, -1), termcrit);
			addRemovePt = false;
		}
		else if (!points[0].empty())
		{
			vector<uchar> status;
			vector<float> err;
			if (prevGray.empty())
				gray.copyTo(prevGray);
			calcOpticalFlowPyrLK(prevGray, gray, points[0], points[1], status, err, winSize,
				3, termcrit, 0, 0.001);
			size_t i, k;
			for (i = k = 0; i < points[1].size(); i++)
			{
				if (addRemovePt)
				{
					if (norm(point - points[1][i]) <= 5)
					{
						addRemovePt = false;
						continue;
					}
				}

				if (!status[i])
					continue;

				points[1][k++] = points[1][i];
				circle(image, points[1][i], 3, Scalar(0, 255, 0), -1, 8);
			}
			points[1].resize(k);
		}

		if (addRemovePt && points[1].size() < (size_t)MAX_COUNT)//当鼠标选中跟踪点
		{
			vector<Point2f> tmp;
			tmp.push_back(point);
			cornerSubPix(gray, tmp, winSize, cvSize(-1, -1), termcrit);
			points[1].push_back(tmp[0]);//将跟踪点存入
			addRemovePt = false;
		}

		needToInit = false;
		imshow("LK Demo", image);

		char c = (char)waitKey(10);
		if (c == 27)
			break;
		switch (c)
		{
		case 'r':
			needToInit = true;
			break;
		case 'c':
			points[0].clear();
			points[1].clear();
			break;
		case 'n':
			nightMode = !nightMode;
			break;
		}

		std::swap(points[1], points[0]);
		cv::swap(prevGray, gray);
	}

	return 0;
}

 运行前需要配置argv参数，很简单可参见：https://blog.csdn.net/qq_30815237/article/details/87006586