【opencv】光流法

简介：在计算机视觉中，Lucas–Kanade光流算法是一种两帧差分的光流估计算法。它由Bruce D. Lucas 和 Takeo Kanade提出。

光流的概念:(Optical flow or optic flow)
它是一种运动模式，这种运动模式指的是一个物体、表面、边缘在一个视角下由一个观察者（比如眼睛、摄像头等）和背景之间形成的明显移动。光流技术，如运动检测和图像分割，时间碰撞，运动补偿编码，三维立体视差，都是利用了这种边缘或表面运动的技术。

二维图像的移动相对于观察者而言是三维物体移动的在图像平面的投影。
有序的图像可以估计出二维图像的瞬时图像速率或离散图像转移。

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

Lucas–Kanade算法：
这个算法是最常见，最流行的。它计算两帧在时间t 到t + δt之间每个每个像素点位置的移动。 由于它是基于图像信号的泰勒级数，这种方法称为差分，这就是对于空间和时间坐标使用偏导数。
图像约束方程可以写为I (x ,y ,z ,t ) = I (x + δx ,y + δy ,z + δz ,t + δt )
I(x, y,z, t) 为在（x,y,z）位置的体素。
我们假设移动足够的小，那么对图像约束方程使用泰勒公式，我们可以得到：

H.O.T. 指更高阶，在移动足够小的情况下可以忽略。从这个方程中我们可以得到：

或者

我们得到：

V x ,V y ,V z 分别是I(x,y,z,t)的光流向量中x，y，z的组成。 , , 和 则是图像在(x ,y ,z ,t )这一点向相应方向的差分 。
所以

I x V x + I y V y + I z V z = − I t。

写做：

这个方程有三个未知量，尚不能被解决，这也就是所谓光流算法的光圈问题。那么要找到光流向量则需要另一套解决的方案。而Lucas-Kanade算法是一个非迭代的算法：

假设流(Vx,Vy,Vz)在一个大小为m*m*m(m>1)的小窗中是一个常数，那么从像素1…n , n = m 3 中可以得到下列一组方程：

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

记作：

为了解决这个超定问题，我们采用最小二乘法：


得到：

其中的求和是从1到n。

这也就是说寻找光流可以通过在四维上图像导数的分别累加得出。我们还需要一个权重函数W(i, j,k) ， 来突出窗口中心点的坐标。高斯函数做这项工作是非常合适的，

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

补充：opencv里实现的看上去蛮复杂，现在还不是太明白。其中LK经典算法也是迭代法，是由高斯迭代法解线性方程组进行迭代的。

参考文献：Lucas B and Kanade T. An Iterative Image Registration Technique with an Application to Stereo Vision. Proc. Of 7th International Joint Conference on Artificial Intelligence (IJCAI), pp.674-679.是81年发表的。

opencv中实现的是84年发表的算法，

参考文献：Bruce D. Lucas , “Generalized Image Matching by the Method of Differences,” doctoral dissertation, tech. report , Robotics Institute, Carnegie Mellon University, July, 1984

这篇论文是收费的，一直都没下到还。

#include <opencv2/video/video.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/core/core.hpp>

#include <iostream>
#include <cstdio>


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为特征点的新位置

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;


//         CvCapture* capture = cvCaptureFromCAM( -1 );        // 摄像头读取文件开关

        VideoCapture capture(0);


        if(capture.isOpened()/*capture*/)        // 摄像头读取文件开关

        {

                while(true)

                {

//                         frame = cvQueryFrame( capture );        // 摄像头读取文件开关

                        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;

}


//////////////////////////////////////////////////////////////////////////

// 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);

        // 去掉一些不好的特征点

        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);

        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]);

        swap(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);

}