opencv 稀疏光流 稠密光流

opencv 稀疏光流 稠密光流

demo: http://download.csdn.net/detail/keen_zuxwang/9860696

参看、学习文档：
OpenCV学习笔记（七）Lucas-Kanade光流跟踪点的选择
http://blog.sina.com.cn/s/blog_674f0d390100i7bx.html
OpenCV之光流法跟踪运动目标
http://blog.csdn.net/u010684134/article/details/49185535

光流
是空间运动物体在观察成像平面上的像素运动的瞬时速度，是利用图像序列中像素在时间域上的变化以及相邻帧之间的相关性来找到上一帧跟当前帧之间存在的对应关系，
从而计算出相邻帧之间物体的运动信息的一种方法。一般而言，光流是由于场景中前景目标本身的移动、相机的运动，或者两者的共同运动所产生的。
（运动检测和图像分割，时间碰撞，运动补偿编码，三维立体视差，都是利用了这种边缘或表面运动的技术）

运动场
其实就是物体在三维真实世界中的运动；
光流场
是运动场在二维图像平面上（人的眼睛或者摄像头）的投影。通过一个图片序列，把每张图像中每个像素的运动速度和运动方向找出来就是光流场
第t帧的时候A点的位置是(x1, y1)，那么我们在第t+1帧的时候再找到A点，假如它的位置是(x2,y2)，那么我们就可以确定A点的运动了：
(ux, vy) = (x2, y2) - (x1,y1)

光流计算方法：基于匹配的方法、频域的方法、梯度的方法。
基于匹配的光流计算方法包括基于特征和基于区域两种
基于频域的方法，也称为基于能量的方法，利用速度可调的滤波组输出频率或相位信息。
基于梯度的方法利用图像序列亮度的时空微分计算2D速度场（光流）。

光流除了提供远近，还可以提供角度信息
光流通过不同目标的运动速度判断它们与我们的距离
光流就是感觉到的明显的视觉运动

光流算法假设
相邻帧之间的亮度恒定
相邻帧之间物体的运动比较微小
同一子图像的像素点具有相同的移动

金字塔Lucas-Kannade算法
在图像金字塔的最高层计算光流，用得到的运动估计结果作为下一层金字塔的起始点，重复这个过程直到到达金字塔的最底层。
这样就将不满足运动的假设可能性降到最小从而实现对更快和更长的运动的跟踪
L-K算法假设：
亮度恒定
即图像场景中目标的像素在帧间运动时外观上保持不变
时间连续或者运动是”小运动“
即图像的运动随时间的变化比较缓慢；
空间一致
即一个场景中同一表面上邻近的点具有相似的运动。

对大而连贯的运动情况，LK光流实际跟踪效果并不是很好：需要大的窗口来捕获大的运动，而大窗口往往违背运动连贯的假设！
而图像金字塔Lucas-Kannade算法可以解决这个问题。

Opencv实现
1、calcOpticalFlowPyrLK、cvCalcOpticalFlowPyrLK
基于图像金字塔的Lucas-Kanande算法，能够跟踪稀疏的光流，图像金字塔的作用是通过不断地下采样图片来减小两帧之间光流点的运动，
从而保证识别的准确性（L-K算法只在位移很小的时候才有作用）
而对于稀疏光流来说，计算时需要在被跟踪之前指定一组点（容易跟踪的点，例如角点），
因此在使用LK方法之前我们需要配合使用cvGoodFeatureToTrack()来寻找角点，然后利用金字塔LK光流算法，对运动进行跟踪

//Tracking.hpp
#define  CV_LKFLOW_PYR_A_READY       1
#define  CV_LKFLOW_PYR_B_READY       2
#define  CV_LKFLOW_INITIAL_GUESSES   4
#define  CV_LKFLOW_GET_MIN_EIGENVALS 8

//! computes sparse optical flow using multi-scale Lucas-Kanade algorithm
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),
                           int flags=0, double minEigThreshold=1e-4);

prevImg：
第一个8位输入图像或者通过 buildOpticalFlowPyramid()建立的金字塔
nextImg：
第二个输入图像或者和prevImg相同尺寸和类型的金字塔
prevPts：
二维点向量存储找到的光流，单精度浮点数
nextPts：
输出二维点向量（用单精度浮点坐标）包括第二幅图像中计算的输入特征的新点位置；当OPTFLOW_USE_INITIAL_FLOW 标志通过，向量必须有和输入一样的尺寸
status：
输出状态向量，若要相应的流特征被发现，向量的每个元素被设置为1，否则，被置为0
err：
输出错误向量；向量的每个元素被设为相应特征的一个错误，误差测量的类型可以在flags参数中设置；如果流不被发现然后错误未被定义（使用status（状态）参数找到此情形）。
winSize：
在每个金字塔水平搜寻窗口的尺寸，相当于一个均值滤波的作用，窗口大小决定了其噪声的抑制能力什么的
Count：
特征点的数目
level：
最大的金字塔层数。如果为 0 , 不使用金字塔 (即金字塔为单层)，如果为 1 , 使用两层，下面依次类推
criteria：
指定在每个金字塔层，为某点寻找光流的迭代过程的终止条件

/* It is Lucas & Kanade method, modified to use pyramids.
   Also it does several iterations to get optical flow for
   every point at every pyramid level.
   Calculates optical flow between two images for certain set of points (i.e.
   it is a "sparse" optical flow, which is opposite to the previous 3 methods) */
CVAPI(void)  cvCalcOpticalFlowPyrLK( const CvArr*  prev, const CvArr*  curr,
                                     CvArr*  prev_pyr, CvArr*  curr_pyr,
                                     const CvPoint2D32f* prev_features,
                                     CvPoint2D32f* curr_features,
                                     int       count,
                                     CvSize    win_size,
                                     int       level,
                                     char*     status,
                                     float*    track_error,
                                     CvTermCriteria criteria,
                                     int       flags )；       

功能：实现了金字塔中 Lucas-Kanade 光流计算的稀疏迭代版本，根据给出的前一帧特征点坐标计算当前视频帧上的特征点坐标。 函数寻找具有子象素精度的坐标值
prev
t 第一帧（跟踪图像的前一帧，一般是定位特征点）
curr
t+dt 第二帧/当前帧
prev_pyr
第一帧的金字塔缓存. 缓存必须有足够的空间来存储金字塔从层 1 到层 #level 的内容。尺寸 (image_width+8)*image_height/3 比特足够了
curr_pyr
第二帧的金字塔缓存,与prev_pyr类似
prev_features
第一帧特征点集（需要发现光流的特征点集，输入参数，如特征点为角点，则用cvGoodFeaturesToTrack（））
curr_features
新计算出来的光流位置的点集（第二帧光流特征点集，计算输出）
count
特征点的数目
win_size
每个金字塔层的搜索窗口尺寸
level
最大的金字塔层数。如果为0, 不使用金字塔 (即金字塔为单层), 如果为 1 , 使用两层，依次类推。
status
发现状态数组，如果对应特征的光流被发现，数组中的每一个元素都被设置为 1， 否则设置为 0。
track_error
双精度数组，包含原始图像碎片与移动点之间的差。为可选参数，可以是 NULL .
criteria
指定在每个金字塔层，为某点寻找光流的迭代过程的终止条件。
flags
其它选项：
CV_LKFLOW_PYR_A_READY , 在调用之前，第一帧的金字塔已经准备好
CV_LKFLOW_PYR_B_READY , 在调用之前，第二帧的金字塔已经准备好
CV_LKFLOW_INITIAL_GUESSES , 在调用之前，数组 B 包含特征的初始坐标 （Hunnish: 在本节中没有出现数组 B，不知是指的哪一个）

两个参数 prev_pyr 和 curr_pyr 都遵循下列规则：
flag根据我们的金字搭是否建立了，设置不同的值，一般在第一次使用0，在一帧处理完了保留当前金字塔为前一帧金字塔，下次处理时候直接可以使用（标识为 CV_LKFLOW_PYR_A_READY）
如果图像指针为 0, 函数在内部为其分配缓存空间，计算金字塔，然后再处理过后释放缓存。 否则，函数计算金字塔且存储它到缓存中，除非设置标识 CV_LKFLOW_PYR_A[B]_READY 。
图像应该足够大以便能够容纳 Gaussian 金字塔数据。调用函数以后，金字塔被计算而且相应图像的标识可以被设置(除了第一个图像的所有图像序列，标识 CV_LKFLOW_PYR_A_READY 被设置)，为下一次调用准备就绪。

2、calcOpticalFlowFarneback、 cvCalcOpticalFlowFarneback
利用Gunnar Farneback算法计算全局性的稠密光流算法（即图像上所有像素点的光流都计算出来），由于要计算图像上所有点的光流，故计算耗时，速度慢
稠密光流需要使用某种插值方法在比较容易跟踪的像素之间进行插值以解决那些运动不明确的像素

Two-Frame Motion Estimation Based on PolynomialExpansion
//! computes dense optical flow using Farneback algorithm
CV_EXPORTS_W void calcOpticalFlowFarneback( InputArray prev, InputArray next,
                           CV_OUT InputOutputArray flow, double pyr_scale, int levels, int winsize,
                           int iterations, int poly_n, double poly_sigma, int flags );

prev：
前一帧图像
next：
后一帧图像
flow：
输出的光流矩阵。矩阵大小同输入的图像一样大，但是矩阵中的每一个元素可不是一个值，而是两个值，分别表示这个点在x方向与y方向的运动量（偏移量）。
pyr_scale：
金字塔上下两层之间的尺度关系
levels：
金字塔层数
winsize：
均值窗口大小，越大越能denoise并且能够检测快速移动目标，但会引起模糊运动区域
iterations：
迭代次数
poly_n：
像素领域大小，一般为5，7等
poly_sigma：
高斯标注差，一般为1-1.5
flags：
计算方法。

主要包括OPTFLOW_USE_INITIAL_FLOW和OPTFLOW_FARNEBACK_GAUSSIAN

enum
{
    OPTFLOW_USE_INITIAL_FLOW = CV_LKFLOW_INITIAL_GUESSES,
    OPTFLOW_LK_GET_MIN_EIGENVALS = CV_LKFLOW_GET_MIN_EIGENVALS,
    OPTFLOW_FARNEBACK_GAUSSIAN = 256
};

/* Estimate optical flow for each pixel using the two-frame G. Farneback algorithm */
CVAPI(void) cvCalcOpticalFlowFarneback( const CvArr* prev, const CvArr* next,
                                        CvArr* flow, double pyr_scale, int levels,
                                        int winsize, int iterations, int poly_n,
                                        double poly_sigma, int flags );

3、calcOpticalFlowSF 
SimpleFlow: A Non-iterative, Sublinear Optical FlowAlgorithm
//! computes dense optical flow using Simple Flow algorithm
CV_EXPORTS_W void calcOpticalFlowSF(Mat& from,
                                    Mat& to,
                                    Mat& flow,
                                    int layers,
                                    int averaging_block_size,
                                    int max_flow);

CV_EXPORTS_W void calcOpticalFlowSF(Mat& from,
                                    Mat& to,
                                    Mat& flow,
                                    int layers,
                                    int averaging_block_size,
                                    int max_flow,
                                    double sigma_dist,
                                    double sigma_color,
                                    int postprocess_window,
                                    double sigma_dist_fix,
                                    double sigma_color_fix,
                                    double occ_thr,
                                    int upscale_averaging_radius,
                                    double upscale_sigma_dist,
                                    double upscale_sigma_color,
                                    double speed_up_thr); 

//! finds the strong enough corners where the cornerMinEigenVal() or cornerHarris() report the local maxima
CV_EXPORTS_W void goodFeaturesToTrack( InputArray image, OutputArray corners,
                                     int maxCorners, double qualityLevel, double minDistance,
                                     InputArray mask=noArray(), int blockSize=3,
                                     bool useHarrisDetector=false, double k=0.04 );

cvGoodFeaturesToTrack 确定图像的强角点
/* Finds a sparse set of points within the selected region
   that seem to be easy to track */
CVAPI(void)  cvGoodFeaturesToTrack( const CvArr* image, CvArr* eig_image,
                                    CvArr* temp_image, CvPoint2D32f* corners,
                                    int* corner_count, double  quality_level,
                                    double  min_distance,
                                    const CvArr* mask CV_DEFAULT(NULL),
                                    int block_size CV_DEFAULT(3),
                                    int use_harris CV_DEFAULT(0),
                                    double k CV_DEFAULT(0.04) );

功能：查找图像中的角点（满足quality_level，min_distance距离等参数）
image：
输入图像，8-位或浮点32-比特，单通道
eig_image：
浮点32位图像，尺寸与输入图像一致
temp_image：
另外一个临时图像，格式与尺寸与eig_image一致
corners：
检测到的角点
corner_count：
检测到的角点数目
quality_level：
最大最小特征值的乘法因子。定义可接受图像角点的最小质量因子，控制特征点的选取。
min_distance：
限制因子。得到的角点间最小距离。使用 Euclidian 距离
mask：
ROI:感兴趣区域。函数在ROI中计算角点，如果 mask 为 NULL，则选择整个图像

cvFindCornerSubPix 用于发现亚像素精度的角点位置

/* Adjust corner position using some sort of gradient search */
CVAPI(void)  cvFindCornerSubPix( const CvArr* image, CvPoint2D32f* corners,
                                 int count, CvSize win, CvSize zero_zone,
                                 CvTermCriteria  criteria );

image
灰度图像
corners
为整数值的像素位置，corners设定了角点的初始位置。
count
计算角点数目
win
指定产生的窗口尺寸
zero_zone
定义一个禁区(与win相似，但通常比win小),该区域在方程组以及自相关矩阵中不被考虑。如果不需要这样一个禁区，则zero_zone应设置为cvSize(-1,-1)
criteria
迭代终止条件，最大迭代次数CV_TERMCRIT_ITER类型、精度CV_TERMCRIT_EPS类型(或者是两者的组合)
终止条件的设置在极大程度上影响最终得到的亚像素值的精度

#include <string>
#include <stdio.h>
#include <stdlib.h>
#include <jni.h>
#include <android/log.h>
#include <iostream>
#include <fstream>
#include <opencv2/core/core.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/objdetect/objdetect.hpp>
#include "opencv2/features2d/features2d.hpp"
#include "opencv2/calib3d/calib3d.hpp"
#include "opencv2/imgproc/types_c.h"
#include "opencv2/imgproc/imgproc_c.h"
#include "opencv2/video/tracking.hpp"
#include "opencv2/video/video.hpp"
#include <opencv2/ml/ml.hpp>

using namespace std;
using namespace cv;

#define LOG_TAG "FeatureTest"
#define LOGD(...) ((void)__android_log_print(ANDROID_LOG_DEBUG, LOG_TAG, __VA_ARGS__))

#ifdef __cplusplus
extern "C" {
#endif

#define BYTE unsigned char

void EqualizeGray(IplImage* pImg)  //直方图均衡化
{
    IplImage* pGray=cvCreateImage(cvSize(pImg->width,pImg->height),IPL_DEPTH_8U,1);
    cvCvtColor(pImg,pGray,CV_BGR2GRAY);
    pGray->origin=pImg->origin;
    //cvShowImage("Gray Image",pGray);

    IplImage* pEqualizedImg=cvCreateImage(cvSize(pImg->width,pImg->height),IPL_DEPTH_8U,1);
    cvEqualizeHist(pGray, pEqualizedImg);
    //cvShowImage("Equalize Image",pEqualizedImg);
    cvReleaseImage(&pGray);
}

void cvFindContours2(Mat pImg)   //颜色识别 轮廓检测  IplImage* pImg
{
    Mat pGray;
    cvtColor(pImg, pGray, CV_BGR2GRAY); // 1、灰度图

    blur(pGray, pGray, Size(5,5)); //2、 均值滤波

    //Finde vertical edges. Car plates have high density of vertical lines
    Mat img_sobel;
    Sobel(pGray, img_sobel, CV_8U, 1, 1, 3, 1, 0, BORDER_DEFAULT);// 3、sobel 边缘提取 xorder=1,yorder=0,kernelsize=3

    Mat img_threshold;
    threshold(img_sobel, img_threshold, 0, 255, CV_THRESH_OTSU); // 4、二值化  img_sobel CV_THRESH_BINARY

    //Morphplogic operation close:remove blank spaces and connect all regions that have a high number of edges
    Mat element = getStructuringElement(MORPH_RECT, Size(17, 3) ); // 5、高级形态学变化，闭操作
    morphologyEx(img_threshold, img_threshold, CV_MOP_CLOSE, element);

    //Find 轮廓 of possibles plates
    vector< vector< Point> > contours; // a vector of contours
    findContours(img_threshold, contours,
                CV_RETR_EXTERNAL, // 提取外部轮廓
                CV_CHAIN_APPROX_NONE); // all pixels of each contours

    LOGD("         contours.size() = %d      !", contours.size());
    drawContours(pImg, contours, -1, Scalar(0,0,255), 2);
}

static IplImage *image=0;
IplImage *grey=0,*prev_grey=0,*pyramid=0,*prev_pyramid=0,*swap_temp;   //创建各类指针，以便存储
CvPoint2D32f *points[2]={0,0};
CvPoint2D32f *swap_points;
char* status=0;   //对变量进行初始化
//int count=0;
int need_to_init=1;
int night_mode=0;
int flags=0;

int runLKOpenCV(IplImage* imgA)
{
    int i=0,k=0,c=0;
    int win_size=15;//10;
    int r=3;
    static int count=0;
    const int MAX_COUNT=500; // 500

    if (!image){
        image=cvCreateImage(cvGetSize(imgA),8,3);
        image->origin=imgA->origin;

        grey=cvCreateImage(cvGetSize(imgA),8,1);
        prev_grey=cvCreateImage(cvGetSize(imgA),8,1);

        pyramid=cvCreateImage(cvGetSize(imgA),8,1);
        prev_pyramid=cvCreateImage(cvGetSize(imgA),8,1);

        points[0]=(CvPoint2D32f*)cvAlloc(MAX_COUNT*sizeof(points[0][0]));
        points[1]=(CvPoint2D32f*)cvAlloc(MAX_COUNT*sizeof(points[0][0]));
        status=(char*)cvAlloc(MAX_COUNT);
        flags=0;
    }

    cvCvtColor(imgA, image, CV_BGRA2BGR);
    cvCvtColor(image,grey,CV_BGR2GRAY);

    if (need_to_init)
    {
       //automatic initialization
       IplImage* eig=cvCreateImage(cvGetSize(grey),32,1);
       IplImage* temp=cvCreateImage(cvGetSize(grey),32,1);
       double quality=0.01;
       double min_distance= 15;//10;

       count=MAX_COUNT;
       cvGoodFeaturesToTrack(grey, eig, temp,
                points[1],
                &count,
                quality,
                min_distance,
                0,
                3,
                0,
                0.04);   //读取第一帧影像

       //能够将角点位置精确到亚像素级精度, 提取易于跟踪的特征点，特征点精确描述
       cvFindCornerSubPix(grey, points[1], count,
                  cvSize(win_size,win_size),
                  cvSize(-1,-1),  // cvSize(1,1)就表示成忽略掉相邻1个像素
                  cvTermCriteria(CV_TERMCRIT_ITER|CV_TERMCRIT_EPS,20, 0.03)); //迭代次数iteration最小精度epsilon

       for (i=0;i<count;i++){
           cvRectangle(image, cvPoint(points[1][i].x-r,points[1][i].y-r), cvPoint(points[1][i].x+r, points[1][i].y+r), cvScalar(0,0,255), 2);
       }

       cvReleaseImage(&eig);
       cvReleaseImage(&temp);
    }

    if (count>0)
    {
        //金字塔LK算法（cvCalcOpticalFlowPyrLK）的两个特点，一、金字塔原理；二、偏导方法求解位移。
        cvCalcOpticalFlowPyrLK(prev_grey,grey,prev_pyramid,pyramid,
            points[0],
            points[1],
            count,
            cvSize(win_size,win_size),
            3,
            status,
            0,
            cvTermCriteria(CV_TERMCRIT_ITER|CV_TERMCRIT_EPS,20,0.03),
            flags);  //读取下一帧影像，进行KLT跟踪，特征点筛选

        flags|=CV_LKFLOW_PYR_A_READY;

        for (i=k=0;i<count;i++)
        {
           if (!status[i])
             continue;
           points[1][k++]=points[1][i];
           cvCircle(image,cvPointFrom32f(points[1][i]),3,CV_RGB(255,0,0),-1,8,0);
         }
        count=k;
    }

    CV_SWAP(prev_grey,grey,swap_temp);
    CV_SWAP(prev_pyramid,pyramid,swap_temp);
    CV_SWAP(points[0],points[1],swap_points);

    need_to_init=0;
}

JNIEXPORT jlong JNICALL Java_com_example_optical_MainActivity_doPyrLKConfig(JNIEnv *env, jclass clz)
{
    LOGD("          doPyrLKConfig         !");
    need_to_init = 1;
}

JNIEXPORT jlong JNICALL Java_com_example_optical_MainActivity_doCofPyrLK0(JNIEnv *env, jclass clz, jlong imageGray)
{
    LOGD("          doCofPyrLK0         !");
    Mat imageMat = Mat(*(Mat*)imageGray);

    IplImage  temp_src = imageMat;
    IplImage* imageg = &temp_src;

    runLKOpenCV(imageg);

    Mat *hist = new Mat(image);
    return (jlong) hist;
}

Mat gray;
Mat gray0;                  // 当前图片
Mat gray_prev;              // 预测图片
vector<Point2f> points0[2]; // point0为特征点的原来位置，point1为特征点的新位置
vector<Point2f> initial;    // 初始化跟踪点的位置
vector<Point2f> features;   // 检测的特征
int maxCount   = 500;       // 检测的最大特征数
double qLevel  = 0.01;      // 特征检测的等级
double minDist = 10.0;      // 两特征点之间的最小距离
vector<uchar> status0;      // 跟踪特征的状态，特征的流发现为1，否则为0
vector<float> err;

//检测新点是否应该被添加
bool addNewPoints(){
    return points0[0].size() <= 10;
}

//决定哪些跟踪点被接受
bool acceptTrackedPoint(int i){
    return status0[i] && ((abs(points0[0][i].x - points0[1][i].x)
            + abs(points0[0][i].y - points0[1][i].y)) > 2);
}

void tracking(Mat &frame, Mat &output)
{
    cvtColor(frame, gray, CV_BGR2GRAY);
    frame.copyTo(output);
    // 添加特征点
    if (addNewPoints()) {
        /*
        //! finds the strong enough corners where the cornerMinEigenVal() or cornerHarris() report the local maxima
        CV_EXPORTS_W void goodFeaturesToTrack( InputArray image, OutputArray corners,
                                             int maxCorners, double qualityLevel, double minDistance,
                                             InputArray mask=noArray(), int blockSize=3,
                                             bool useHarrisDetector=false, double k=0.04 );*/
        goodFeaturesToTrack(gray, features, maxCount, qLevel, minDist);
        points0[0].insert(points0[0].end(), features.begin(), features.end());
        initial.insert(initial.end(), features.begin(), features.end());
    }

    if (gray_prev.empty()) {
        gray.copyTo(gray_prev);
    }

    // l-k光流法运动估计
    /*
    //! computes sparse optical flow using multi-scale Lucas-Kanade algorithm
    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),
                               int flags=0, double minEigThreshold=1e-4);*/
    calcOpticalFlowPyrLK(gray_prev, gray, points0[0], points0[1], status0, err);

    // 去掉一些不好的特征点
    int k = 0;
    for (size_t i=0; i<points0[1].size(); i++) {
        if (acceptTrackedPoint(i)) {
            initial[k] = initial[i];
            points0[1][k++] = points0[1][i];
        }
    }
    points0[1].resize(k);
    initial.resize(k);

    // 显示特征点和运动轨迹
    for (size_t i=0; i<points0[1].size(); i++) {
        line(output, initial[i], points0[1][i], Scalar(0, 0, 255));
        circle(output, points0[1][i], 3, Scalar(255, 255, 0), -1);
    }

    // 把当前跟踪结果作为下一次参考
    swap(points0[1], points0[0]);
    swap(gray_prev, gray);
}

JNIEXPORT jlong JNICALL Java_com_example_optical_MainActivity_doCofPyrLK(JNIEnv *env, jclass clz, jlong imageGray, jlong imageGray0)
{
    LOGD("          doCofPyrLK         !");

    Mat frame = Mat(*(Mat*)imageGray);
    Mat frame0 = Mat(*(Mat*)imageGray0);
    //Mat result;
    //tracking(frame, result);

    cvtColor(frame, gray_prev, CV_BGR2GRAY);
    cvtColor(frame0, gray, CV_BGR2GRAY);

    // 添加特征点
    goodFeaturesToTrack(gray, features, maxCount, qLevel, minDist);
    points0[0].insert(points0[0].end(), features.begin(), features.end());
    initial.insert(initial.end(), features.begin(), features.end());

    // l-k光流法运动估计
    calcOpticalFlowPyrLK(gray_prev, gray, points0[0], points0[1], status0, err);

    // 去掉一些不好的特征点
    int k = 0;
    for (size_t i=0; i<points0[1].size(); i++) {
        if (acceptTrackedPoint(i)) {
            initial[k] = initial[i];
            points0[1][k++] = points0[1][i];
        }
    }
    points0[1].resize(k);
    initial.resize(k);

    // 显示特征点和运动轨迹
    for (size_t i=0; i<points0[1].size(); i++) {
        line(frame0, initial[i], points0[1][i], Scalar(0, 0, 255));
        circle(frame0, points0[1][i], 3, Scalar(255, 255, 0), -1);
    }

    //Mat mtx(image,0);
    Mat *hist = new Mat(frame0); // result
    return (jlong) hist;
}

#define UNKNOWN_FLOW_THRESH 1e9
void makecolorwheel(vector<Scalar> &colorwheel)
{
    int RY = 15;
    int YG = 6;
    int GC = 4;
    int CB = 11;
    int BM = 13;
    int MR = 6;
    int i;
    for (i = 0; i < RY; i++) colorwheel.push_back(Scalar(255,       255*i/RY,     0));
    for (i = 0; i < YG; i++) colorwheel.push_back(Scalar(255-255*i/YG, 255,       0));
    for (i = 0; i < GC; i++) colorwheel.push_back(Scalar(0,         255,      255*i/GC));
    for (i = 0; i < CB; i++) colorwheel.push_back(Scalar(0,         255-255*i/CB, 255));
    for (i = 0; i < BM; i++) colorwheel.push_back(Scalar(255*i/BM,      0,        255));
    for (i = 0; i < MR; i++) colorwheel.push_back(Scalar(255,       0,        255-255*i/MR));
}

void motionToColor(Mat flow, Mat &color)
{
    if (color.empty())
        color.create(flow.rows, flow.cols, CV_8UC3);

    static vector<Scalar> colorwheel; //Scalar r,g,b
    if (colorwheel.empty())
        makecolorwheel(colorwheel);

    // determine motion range:
    float maxrad = -1;

    // Find max flow to normalize fx and fy
    for (int i= 0; i < flow.rows; ++i) {
        for (int j = 0; j < flow.cols; ++j) {
            Vec2f flow_at_point = flow.at<Vec2f>(i, j);
            float fx = flow_at_point[0];
            float fy = flow_at_point[1];
            if ((fabs(fx) >  UNKNOWN_FLOW_THRESH) || (fabs(fy) >  UNKNOWN_FLOW_THRESH))
                continue;
            float rad = sqrt(fx * fx + fy * fy);
            maxrad = maxrad > rad ? maxrad : rad;
        }
    }

    for (int i= 0; i < flow.rows; ++i)
    {
        for (int j = 0; j < flow.cols; ++j){
            uchar *data = color.data + color.step[0] * i + color.step[1] * j;
            Vec2f flow_at_point = flow.at<Vec2f>(i, j);

            float fx = flow_at_point[0] / maxrad;
            float fy = flow_at_point[1] / maxrad;
            if ((fabs(fx) >  UNKNOWN_FLOW_THRESH) || (fabs(fy) >  UNKNOWN_FLOW_THRESH))
            {
                data[0] = data[1] = data[2] = 0;
                continue;
            }
            float rad = sqrt(fx * fx + fy * fy);

            float angle = atan2(-fy, -fx) / CV_PI;
            float fk = (angle + 1.0) / 2.0 * (colorwheel.size()-1);
            int k0 = (int)fk;
            int k1 = (k0 + 1) % colorwheel.size();
            float f = fk - k0;
            //f = 0; // uncomment to see original color wheel

            for (int b = 0; b < 3; b++)
            {
                float col0 = colorwheel[k0][b] / 255.0;
                float col1 = colorwheel[k1][b] / 255.0;
                float col = (1 - f) * col0 + f * col1;
                if (rad <= 1)
                    col = 1 - rad * (1 - col); // increase saturation with radius
                else
                    col *= .75; // out of range
                data[2 - b] = (int)(255.0 * col);
            }
        }
    }
}

JNIEXPORT jlong JNICALL Java_com_example_optical_MainActivity_doCofFarneback(JNIEnv *env, jclass clz, jlong imageGray, jlong imageGray0)
{
    LOGD("            doCofFarneback        !");
    //Mat& img_display_0  = *(Mat*)imageGray;
    Mat input1;
    Mat input2;
    Mat flow;
    Mat I1 = Mat(*(Mat*)imageGray);
    Mat I2 = Mat(*(Mat*)imageGray0);

    input1.create(I1.size(),CV_8UC1);//CV_32FC1
    cvtColor(I1,input1,CV_BGR2GRAY);

    input2.create(I2.size(),CV_8UC1);//CV_32FC1
    cvtColor(I2,input2,CV_BGR2GRAY);

    calcOpticalFlowFarneback(input1,input2,flow,0.5,3,20,3,5,1.2,0);
    LOGD("      calcOpticalFlowFarneback      !");

    //flow = abs(flow);
    Mat motion2color;
    motionToColor(flow, motion2color);

    LOGD("      motionToColor      !");

    Mat *hist = new Mat(motion2color);
    return (jlong) hist;
}

void detect_and_draw_objects( IplImage* image, CvHaarClassifierCascade* cascade,int do_pyramids )
{
    IplImage* small_image = image;
    CvMemStorage* storage = cvCreateMemStorage(0);
    CvSeq* faces;
    int i, scale = 1;

    //金字塔分解，对输入图像向下采样，提升计算性能
    if( do_pyramids ) {
        small_image = cvCreateImage(cvSize(image->width/2,image->height/2), IPL_DEPTH_8U, 3); //sceneMat.type()
        cvPyrDown( image, small_image, CV_GAUSSIAN_5x5 );
        scale = 2;
    }

    faces = cvHaarDetectObjects(small_image, cascade, storage, 1.1, 2, CV_HAAR_DO_CANNY_PRUNING);

    /* draw all the rectangles */
    for(i = 0; i < faces->total; i++) {
        CvRect face_rect = *(CvRect*)cvGetSeqElem(faces, i);
        cvRectangle(image, cvPoint(face_rect.x*scale,face_rect.y*scale),
                     cvPoint((face_rect.x+face_rect.width)*scale,(face_rect.y+face_rect.height)*scale),
                     CV_RGB(255,0,0), 3);
    }

    if(small_image != image)
        cvReleaseImage(&small_image);
    cvReleaseMemStorage(&storage);
}

JNIEXPORT jlong JNICALL Java_com_example_optical_MainActivity_doHaarClassifier(JNIEnv *env, jclass clz, jlong imageGray)
{
    LOGD("          doHaarClassifier         !");
    Mat imageMat = Mat(*(Mat*)imageGray);
    IplImage  temp_src = imageMat;
    IplImage* image = &temp_src;

    //const char* classPath = "/storage/emulated/0/OCR.xml";
    CvHaarClassifierCascade* cascade=(CvHaarClassifierCascade*)cvLoad("/storage/emulated/0/haarcascade_frontalface_default.xml");
    //haarcascade_fullbody.xml

    detect_and_draw_objects(image, cascade, 0 );
    cvReleaseHaarClassifierCascade( &cascade );

    Mat *hist = new Mat(image);
    return (jlong) hist;
}

#ifdef __cplusplus
}
#endif

cvCalcOpticalFlowPyrLK（）：

calcOpticalFlowFarneback（）：

CvHaarClassifierCascade：