利用图片光流算法扭曲图片

‘壹’ 运动目标检测——光流法与opencv代码实现

运动目标的检测的其主要目的是 获取目标对象的运动参数（位置、速度、加速度等）及运动轨迹 ，通过进一步分析处理，实现对目标行为更高层级上的理解。
运动目标检测技术目的是 从序列图像中将变化区域从背景图像中提取出来 ，常用于视频监视、图像压缩、三维重构、异常检测等。

运动目标检测主流方法有帧差法、背景差法、光流法等。光流法源于 仿生学 思想，更贴近于直觉，大量昆虫的视觉机理便是基于光流法。
二十世纪五十年代心理学家Gibson在他的着作“The Perception of Visual World”中首次提出了以心理学实验为基础的光流法基本概念，而直到八十年代才由Horn、Kanade、Lucash和Schunck创造性地将灰度与二维速度场相联系，引入光流约束方程的算法，对光流计算做了奠基性的工作。

光流（optical flow）：由于目标对象或者摄像机的移动造成的图像对象在连续两帧图像中的移动。

通俗说，对于一个图片序列，把每张图像每个像素在连续帧之间的运动速度和方向（ 某像素点在连续两帧上的位移矢量 ）找出来就是光流场。

第t帧的时A点的位置是(x1, y1)，第t+1帧时A点位置是(x2,y2)，则像素点A的位移矢量：(ux, vy) = (x2, y2) - (x1,y1)

如何知道第t+1帧的时候A点的位置涉及到不同的光流计算方法，主要有四种：基于梯度的方法、基于匹配的方法、基于能量的方法、基于相位的方法。

光流法依赖于三个假设：

根据所形成的光流场中 二维矢量的疏密程度 ，光流法可分为稠密光流与稀疏光流。

稀疏光流只对有 明显特征的组点 （如角点）进行跟踪，计算开销小。

http://www.opencv.org.cn/opencvdoc/2.3.2/html/moles/video/doc/motion_analysis_and_object_tracking.html#calcopticalflowfarneback
（1）calcOpticalFlowPyrLK
基于金字塔LK光流算法，计算某些点集的稀疏光流。
参考论文《Pyramidal Implementation of the Lucas Kanade Feature TrackerDescription of the algorithm》
（2）calcOpticalFlowFarneback
基于Gunnar Farneback 的算法计算稠密光流。
参考论文《Two-Frame Motion Estimation Based on PolynomialExpansion》
（3）CalcOpticalFlowBM
通过块匹配的方法来计算光流
（4）CalcOpticalFlowHS
基于Horn-Schunck 的算法计算稠密光流。
参考论文《Determining Optical Flow》
（5）calcOpticalFlowSF
论文《SimpleFlow: A Non-iterative, Sublinear Optical FlowAlgo》的实现

‘贰’ 视频检索的智能视频

智能视频处理成为视频监控的“救命稻草”
智能视频源自计算机视觉技术，计算机视觉技术是人工智能研究的分支之一，它能够在图像及图像内容描述之间建立映射关系，从而使计算机能够通过数字图像处理和分析来有限理解视频画面中的内容。运用智能视频分析技术，当系统发现符合某种规则的行为（如定向运动、越界、游荡、遗留等）发生时，自动向监控系统发出报警信号（如声光报警），提示相关工作人员及时处理可疑事件。
智能视频算法的实现
目前，智能视频技术实现对移动目标的实时检测、识别、分类以及多目标跟踪等功能的主要算法分为以下五类：目标检测、目标跟踪、目标识别、行为分析、基于内容的视频检索和数据融合等。 目标检测（Object Detection）是按一定时间间隔从视频图像中抽取像素，采用软件技术来分析数字化的像素，将运动物体从视频序列中分离出来。运动目标检测技术是智能化分析的基础。常用的目标检测技术可以分为背景减除法（Background Subtraction）、时间差分法（Temporal Difference）和光流法（Optic Flow）三类。
背景减除法利用当前图像与背景图像的差分检测运动区域。背景减除法假设视频场景中有一个背景，而背景和前景并未给出严格定义，背景在实际使用中是变化的，所以背景建模是背景减除法中非常关键的一步。常用的背景建模方法有时间平均法、自适应更新法、高斯模型等。背景减除法能够提供相对来说比较完全的运动目标特征数据，但对于动态场景的变化，如光线照射情况、摄像机抖动和外来无关事件的干扰特别敏感。
时间差分法充分利用了视频图像的时域特征，利用相邻帧图像的相减来提取出前景移动目标的信息。该方法对于动态环境具有较强的自适应性，不对场景做任何假设，但一般不能完全提取出所有相关的特征像素点，在运动实体内部容易产生空洞现象，只能够检测到目标的边缘。当运动目标停止时，一般时间差分法便失效。 光流法通过比较连续帧为每个图像中的像素赋予一个运动矢量从而分割出运动物体。
光流法能够在摄像机运动的情况下检测出独立的运动目标，然而光流法运算复杂度高并且对噪声很敏感，所以在没有专门硬件支持下很难用于实时视频流检测中。 目标跟踪（Object Tracking）算法根据不同的分类标准，有着以下两种分类方法：根据目标跟踪与目标检测的时间关系分类和根据目标跟踪的策略分类。 根据目标跟踪与目标检测的时间关系的分类有三种：
一是先检测后跟踪(Detect before Track)，先检测每帧图像上的目标，然后将前后两帧图像上目标进行匹配，从而达到跟踪的目的。这种方法可以借助很多图像处理和数据处理的现有技术，但是检测过程没有充分利用跟踪过程提供的信息。
二是先跟踪后检测(Track before Detect)，先对目标下一帧所在的位置及其状态进行预测或假设，然后根据检测结果来矫正预测值。这一思路面临的难点是事先要知道目标的运动特性和规律。三是边检测边跟踪(Track while Detect)，图像序列中目标的检测和跟踪相结合，检测要利用跟踪来提供处理的对象区域，跟踪要利用检测来提供目标状态的观察数据。
根据目标跟踪的策略来分类，通常可分为3D方法和2D方法。相对3D方法而言，2D方法速度较快，但对于遮挡问题难以处理。基于运动估计的跟踪是最常用的方法之一。 目标识别（Object Recognize）利用物体颜色、速度、形状、尺寸等信息进行判别，区分人、交通工具和其他对象。目标识别常用人脸识别和车辆识别。
视频人脸识别的通常分为四个步骤：人脸检测、人脸跟踪、特征提取和比对。人脸检测指在动态的场景与复杂的背景中判断是否存在面像，并分离出这种面像。人脸跟踪指对被检测到的面貌进行动态目标跟踪。常用方法有基于模型的方法、基于运动与模型相结合的方法、肤色模型法等。
人脸特征提取方法归纳起来分为三类：第一类是基于边缘、直线和曲线的基本方法；第二类是基于特征模板的方法；第三类是考虑各种特征之间几何关系的结构匹配法。单一基于局部特征的提取方法在处理闭眼、眼镜和张嘴等情景时遇到困难，相对而言，基于整体特征统计的方法对于图像亮度和特征形变的鲁棒性更强。人脸比对是将抽取出的人脸特征与面像库中的特征进行比对，并找出最佳的匹配对象。
车辆识别主要分为车牌照识别、车型识别和车辆颜色识别等，应用最广泛和技术较成熟的是车牌照识别。 车牌照识别的步骤分别为：车牌定位、车牌字符分割、车牌字符特征提取和车牌字符识别。
车牌定位是指从车牌图像中找到车牌区域并把其分离出来。字符分割是将汉字、英文字母和数字字符从牌照中提取出来。车牌特征提取的基本任务是从众多特征中找出最有效的特征，常用的方法有逐像素特征提取法、骨架特征提取法、垂直水平方向数据统计特征提取法、特征点提取法和基于统计特征的提取法。车牌字符识别可以使用贝叶斯分离器、支持向量机(SVM)和神经网络分类器(NNC)等算法。 基于内容的图像检索技术是由用户提交检索样本，系统根据样本对象的底层物理特征生成特征集，然后在视频库中进行相似性匹配，得到检索结果的过程。现有基于内容的检索方法主要分为：基于颜色的检索方法、基于形状的检索方法和基于纹理的检索方法等。数据融合是将来自不同视频源的数据进行整合，以获得更丰富的数据分析结果。

‘叁’ 如何使用opencv实现金字塔光流lk跟踪算法

#include <stdio.h>
#include <windows.h>
#include "cv.h"
#include "cxcore.h"
#include "highgui.h"
#include <opencv2\opencv.hpp>
using namespace cv;

static const double pi = 3.14159265358979323846;
inline static double square(int a)
{
return a * a;
}
/*该函数目的：给img分配内存空间，并设定format，如位深以及channel数*/
inline static void allocateOnDemand(IplImage **img, CvSize size, int depth, int channels)
{
if (*img != NULL) return;
*img = cvCreateImage(size, depth, channels);
if (*img == NULL)
{
fprintf(stderr, "Error: Couldn't allocate image. Out of memory?\n");
exit(-1);
}
}
/*主函数，原程序是读取avi视频文件，然后处理，我简单改成从摄像头直接读取数据*/
int main(int argc, char *argv[])
{

//读取摄像头
VideoCapture cap(0);
//读取视频文件

//VideoCapture cap; cap.open("optical_flow_input.avi");
if (!cap.isOpened())
{
return -1;
}
Mat frame;

/*
bool stop = false;
while (!stop)
{
cap >> frame;
// cvtColor(frame, edges, CV_RGB2GRAY);
// GaussianBlur(edges, edges, Size(7, 7), 1.5, 1.5);
// Canny(edges, edges, 0, 30, 3);
// imshow("当前视频", edges);
imshow("当前视频", frame);
if (waitKey(30) >= 0)
stop = true;
}
*/

//CvCapture *input_video = cvCaptureFromFile( "optical_flow_input.avi" );
//cv::VideoCapture cap = *(cv::VideoCapture *) userdata;

//if (input_video == NULL)
// {
// fprintf(stderr, "Error: Can't open video device.\n");
// return -1;
// }

/*先读取一帧，以便得到帧的属性，如长、宽等*/
//cvQueryFrame(input_video);

/*读取帧的属性*/
CvSize frame_size;
frame_size.height = cap.get(CV_CAP_PROP_FRAME_HEIGHT);
frame_size.width = cap.get(CV_CAP_PROP_FRAME_WIDTH);

/*********************************************************/

/*用于把结果写到文件中去,非必要
int frameW = frame_size.height; // 744 for firewire cameras
int frameH = frame_size.width; // 480 for firewire cameras
VideoWriter writer("VideoTest.avi", -1, 25.0, cvSize(frameW, frameH), true);

/*开始光流法*/
//VideoWriter writer("VideoTest.avi", CV_FOURCC('D', 'I', 'V', 'X'), 25.0, Size(640, 480), true);

while (true)
{
static IplImage *frame = NULL, *frame1 = NULL, *frame1_1C = NULL,
*frame2_1C = NULL, *eig_image = NULL, *temp_image = NULL,
*pyramid1 = NULL, *pyramid2 = NULL;

Mat framet;
/*获取第一帧*/
// cap >> framet;
cap.read(framet);
Mat edges;
//黑白抽象滤镜模式
// cvtColor(framet, edges, CV_RGB2GRAY);
// GaussianBlur(edges, edges, Size(7, 7), 1.5, 1.5);
// Canny(edges, edges, 0, 30, 3);

//转换mat格式到lpiimage格式
frame = &IplImage(framet);
if (frame == NULL)
{
fprintf(stderr, "Error: Hmm. The end came sooner than we thought.\n");
return -1;
}

/*由于opencv的光流函数处理的是8位的灰度图，所以需要创建一个同样格式的
IplImage的对象*/
allocateOnDemand(&frame1_1C, frame_size, IPL_DEPTH_8U, 1);

/* 把摄像头图像格式转换成OpenCV惯常处理的图像格式*/
cvConvertImage(frame, frame1_1C, 0);

/* 我们需要把具有全部颜色信息的原帧保存，以备最后在屏幕上显示用*/
allocateOnDemand(&frame1, frame_size, IPL_DEPTH_8U, 3);
cvConvertImage(frame, frame1, 0);

/* 获取第二帧 */
//cap >> framet;
cap.read(framet);
// cvtColor(framet, edges, CV_RGB2GRAY);
// GaussianBlur(edges, edges, Size(7, 7), 1.5, 1.5);
// Canny(edges, edges, 0, 30, 3);
frame = &IplImage(framet);
if (frame == NULL)
{
fprintf(stderr, "Error: Hmm. The end came sooner than we thought.\n");
return -1;
}

/*原理同上*/
allocateOnDemand(&frame2_1C, frame_size, IPL_DEPTH_8U, 1);
cvConvertImage(frame, frame2_1C, 0);

/*********************************************************
开始shi-Tomasi算法，该算法主要用于feature selection,即一张图中哪些是我
们感兴趣需要跟踪的点(interest point)
input:
* "frame1_1C" 输入图像.
* "eig_image" and "temp_image" 只是给该算法提供可操作的内存区域.
* 第一个".01" 规定了特征值的最小质量，因为该算法要得到好的特征点，哪就
需要一个选择的阈值
* 第二个".01" 规定了像素之间最小的距离，用于减少运算复杂度，当然也一定
程度降低了跟踪精度
* "NULL" 意味着处理整张图片，当然你也可以指定一块区域
output:
* "frame1_features" 将会包含fram1的特征值
* "number_of_features" 将在该函数中自动填充上所找到特征值的真实数目,
该值<= 400
**********************************************************/

/*开始准备该算法需要的输入*/

/* 给eig_image,temp_image分配空间*/
allocateOnDemand(&eig_image, frame_size, IPL_DEPTH_32F, 1);
allocateOnDemand(&temp_image, frame_size, IPL_DEPTH_32F, 1);

/* 定义存放frame1特征值的数组，400只是定义一个上限 */
CvPoint2D32f frame1_features[400];
int number_of_features = 400;

/*开始跑shi-tomasi函数*/
cvGoodFeaturesToTrack(frame1_1C, eig_image, temp_image,
frame1_features, &number_of_features, .01, .01, NULL);

/**********************************************************
开始金字塔Lucas Kanade光流法，该算法主要用于feature tracking,即是算出
光流，并跟踪目标。
input:
* "frame1_1C" 输入图像，即8位灰色的第一帧
* "frame2_1C" 第二帧，我们要在其上找出第一帧我们发现的特征点在第二帧
的什么位置
* "pyramid1" and "pyramid2" 是提供给该算法可操作的内存区域，计算中间
数据
* "frame1_features" 由shi-tomasi算法得到的第一帧的特征点.
* "number_of_features" 第一帧特征点的数目
* "optical_flow_termination_criteria" 该算法中迭代终止的判别，这里是
epsilon<0.3，epsilon是两帧中对应特征窗口的光度之差的平方，这个以后的文
章会讲
* "0" 这个我不知道啥意思，反正改成1就出不来光流了，就用作者原话解释把
means disable enhancements. (For example, the second array isn't
pre-initialized with guesses.)
output:
* "frame2_features" 根据第一帧的特征点，在第二帧上所找到的对应点
* "optical_flow_window" lucas-kanade光流算法的运算窗口,具体lucas-kanade
会在下一篇详述
* "5" 指示最大的金字塔层数，0表示只有一层，那就是没用金字塔算法
* "optical_flow_found_feature" 用于指示在第二帧中是否找到对应特征值，
若找到，其值为非零
* "optical_flow_feature_error" 用于存放光流误差
**********************************************************/

/*开始为pyramid lucas kanade光流算法输入做准备*/
CvPoint2D32f frame2_features[400];

/* 该数组相应位置的值为非零，如果frame1中的特征值在frame2中找到 */
char optical_flow_found_feature[400];

/* 数组第i个元素表对应点光流误差*/
float optical_flow_feature_error[400];

/*lucas-kanade光流法运算窗口,这里取3*3的窗口,可以尝试下5*5,区别就是5*5
出现aperture problem的几率较小,3*3运算量小，对于feature selection即shi-tomasi算法来说足够了*/
CvSize optical_flow_window = cvSize(5, 5);
// CvSize optical_flow_window = cvSize(5, 5);
/* 终止规则，当完成20次迭代或者当epsilon<=0.3，迭代终止，可以尝试下别的值*/
CvTermCriteria optical_flow_termination_criteria= cvTermCriteria(CV_TERMCRIT_ITER | CV_TERMCRIT_EPS, 20, .3);

/*分配工作区域*/
allocateOnDemand(&pyramid1, frame_size, IPL_DEPTH_8U, 1);
allocateOnDemand(&pyramid2, frame_size, IPL_DEPTH_8U, 1);

/*开始跑该算法*/
cvCalcOpticalFlowPyrLK(frame1_1C, frame2_1C, pyramid1, pyramid2,frame1_features, frame2_features, number_of_features,
optical_flow_window, 5, optical_flow_found_feature,optical_flow_feature_error, optical_flow_termination_criteria, 0);

/*画光流场，画图是依据两帧对应的特征值，
这个特征值就是图像上我们感兴趣的点，如边缘上的点P(x,y)*/
for (int i = 0; i< number_of_features; i++)
{
/* 如果没找到对应特征点 */
if (optical_flow_found_feature[i] == 0)
continue;
int line_thickness;
line_thickness = 1;

/* CV_RGB(red, green, blue) is the red, green, and blue components
* of the color you want, each out of 255.
*/
CvScalar line_color;
line_color = CV_RGB(255, 0, 0);

/*画箭头,因为帧间的运动很小，所以需要缩放，不然看不见箭头，缩放因子为3*/
CvPoint p, q;
p.x = (int)frame1_features[i].x;
p.y = (int)frame1_features[i].y;
q.x = (int)frame2_features[i].x;
q.y = (int)frame2_features[i].y;

double angle;
angle = atan2((double)p.y - q.y, (double)p.x - q.x);
double hypotenuse;
hypotenuse = sqrt(square(p.y - q.y) + square(p.x - q.x));

/*执行缩放*/
q.x = (int)(p.x - 5 * hypotenuse * cos(angle));
q.y = (int)(p.y - 5 * hypotenuse * sin(angle));

/*画箭头主线*/
/* "frame1"要在frame1上作画.
* "p" 线的开始点.
* "q" 线的终止点.
* "CV_AA" 反锯齿.
* "0" 没有小数位.
*/
cvLine(frame1, p, q, line_color, line_thickness, CV_AA, 0);

/* 画箭的头部*/
p.x = (int)(q.x + 9 * cos(angle + pi / 4));
p.y = (int)(q.y + 9 * sin(angle + pi / 4));
cvLine(frame1, p, q, line_color, line_thickness, CV_AA, 0);
p.x = (int)(q.x + 9 * cos(angle - pi / 4));
p.y = (int)(q.y + 9 * sin(angle - pi / 4));
cvLine(frame1, p, q, line_color, line_thickness, CV_AA, 0);
}
/*显示图像*/

/*创建一个名为optical flow的窗口，大小自动改变*/
cvNamedWindow("Optical Flow", CV_WINDOW_NORMAL);
cvFlip(frame1, NULL, 2);
cvShowImage("Optical Flow", frame1);

/*延时，要不放不了*/
cvWaitKey(33);

/*写入到文件中去*/

// cv::Mat m = cv::cvarrToMat(frame1);//转换lpimgae到mat格式
// writer << m;//opencv3.0 version writer

}
cap.release();
cvWaitKey(33);
system("pause");
}

‘肆’ 用C#可以写光流法吗

自己实现的一个光流算法，通过模式搜索匹配的方式计算相邻两张图片的平移量。

模式匹配：选择方块模式或者X形模式，在两个图片中，将该模式像素灰度值做差并求和，相差最小的认为上最匹配的。

多模式匹配：在图片中选择多个位置，检索最符合的位置，最后将多个位置的匹配结果作平均值。

经过测试，在草地、柏油路面、地毯等非规则图形的粗糙表面上表现良好。

// optical flow use Multi-Pattern-Match algrithm, use the big inner diff pattern to do multi pattern match, then averge the result// already implemented: square-pattern, X-patternclass COpticalFlow_MPM
{public:
COpticalFlow_MPM(){} virtual ~COpticalFlow_MPM(){} static bool AddImplementation(COpticalFlow_MPM* imp)
{ if(m_impNum < c_maxImpNum){
m_impTbl[m_impNum++] = imp; return true;
} return false;
} static void SetImageDimesion(int width, int height, int lineBytes)
{ for(int i = 0; i < m_impNum; ++i){
m_impTbl[i]->m_width = width;
m_impTbl[i]->m_height = height;
m_impTbl[i]->m_lineBytes = lineBytes;
m_impTbl[i]->GenerateSearchTable();
m_impTbl[i]->GeneratePatternTable();
}
} // auto choose the pattern to do optical flow
static void AutoOpticalFlow(uint8_t* image1, uint8_t* image2)
{
m_impTbl[m_impCurr]->calcOpticalFlow(image1, image2); // check if need switch pattern
static int s_goodCount = 0; static int s_badCount = 0; if(m_quality > 0){
s_goodCount++;
}else{
s_badCount++;
} if(s_goodCount + s_badCount > 30){ if(s_badCount * 2 > s_goodCount){
m_impCurr = m_impCurr < (m_impNum - 1) ? m_impCurr + 1 : 0;
}
s_goodCount = s_badCount = 0;
}
} // the result
static uint8_t m_quality; // 0 ~ 255, 0 means the optical flow is invalid.
static float m_offset_x; // unit is pixel
static float m_offset_y;protected: virtual const char* Name() = 0; virtual void GeneratePatternTable() = 0; // prepare the address offset tables, that can make the calculation simple and fast.
void GenerateSearchTable()
{ // generate the search offset from corresponding location to the max distance
int index = 0; int yNum, ay[2]; for (int dist = 1; dist <= c_searchD; ++dist){ for (int x = -dist; x <= dist; ++x){ // for each x, only have 1 or 2 dy choices.
ay[0] = dist - abs(x); if (ay[0] == 0){
yNum = 1;
} else{
yNum = 2;
ay[1] = -ay[0];
} for (int iy = 0; iy < yNum; ++iy){
m_searchOffsets[index++] = ay[iy] * m_lineBytes + x;
}
}
} // generate the watch points.
index = 0; int center = m_width * m_height / 2 + m_width / 2; for (int y = -c_watchN; y <= c_watchN; ++y){ for (int x = -c_watchN; x <= c_watchN; ++x){
m_watchPoints[index++] = center + y * c_watchG * m_lineBytes + x * c_watchG * m_width / m_height;
}
}
} void ResetResult()
{
m_quality = 0;
m_offset_x = 0;
m_offset_y = 0;
} void calcOpticalFlow(uint8_t* image1, uint8_t* image2)
{
ResetResult(); int betterStart; int matchedOffset; int x1, y1, x2, y2; int matchedCount = 0; int offset_x[c_watchS]; int offset_y[c_watchS]; for (int i = 0; i < c_watchS; ++i){ if (SearchMaxInnerDiff(image1, m_watchPoints[i], betterStart)){
int32_t minDiff = SearchBestMatch(image1 + betterStart, m_patternOffsets, c_patternS, image2, betterStart, matchedOffset); if (minDiff < c_patternS * c_rejectDiff){
x1 = betterStart % m_lineBytes; y1 = betterStart / m_lineBytes;
x2 = matchedOffset % m_lineBytes; y2 = matchedOffset / m_lineBytes;
m_offset_x += (x2 - x1);
m_offset_y += (y2 - y1);
offset_x[matchedCount] = (x2 - x1);
offset_y[matchedCount] = (y2 - y1);
matchedCount++;
}
}
} if (matchedCount >= 4){
m_offset_x /= matchedCount;
m_offset_y /= matchedCount; // calculate the variance, and use the variance to get the quality.
float varX = 0, varY = 0; for (int i = 0; i < matchedCount; ++i){
varX += (offset_x[i] - m_offset_x) * (offset_x[i] - m_offset_x);
varY += (offset_y[i] - m_offset_y) * (offset_y[i] - m_offset_y);
}
varX /= (matchedCount - 1);
varY /= (matchedCount - 1); float varMax = varX > varY ? varX : varY;
m_quality = (uint8_t)(varMax > 2 ? 0 : (2-varMax) * 255 / 2); if(m_quality == 0){
ResetResult();
}
}
} // get the pattern inner diff, the pattern is center of the area.
inline int32_t InnerDiff(const uint8_t* center, const int* patternPoints, const int patternSize)
{
int32_t sum = 0;
int32_t mean = 0; for (int i = 0; i < patternSize; ++i){
sum += center[patternPoints[i]];
}
mean = sum / patternSize;

int32_t sumDiff = 0; for (int i = 0; i < patternSize; ++i){
sumDiff += abs(center[patternPoints[i]] - mean);
} return sumDiff;
} // get the sum diff between two pattern, the pattern is the center of the area.
inline int32_t PatternDiff(const uint8_t* center1, const uint8_t* center2, const int* patternPoints, const int patternSize)
{
int32_t sumDiff = 0; for (int i = 0; i < patternSize; ++i){
sumDiff += abs(center1[patternPoints[i]] - center2[patternPoints[i]]);
} return sumDiff;
} // search the max inner diff location, image is the full image begining, the return value searchOffset is base on the image begining.
inline bool SearchMaxInnerDiff(const uint8_t* image, int searchStart, int& betterStart)
{ // if the inner diff is less than this number, cannot use this pattern to do search.
const int c_minInnerDiff = c_patternS * 4; const int c_acceptInnerDiff = c_patternS * 12; const uint8_t* searchCenter = image + searchStart;
int32_t currDiff = InnerDiff(searchCenter, m_patternOffsets, c_patternS);
int32_t maxDiff = currDiff;
betterStart = 0; for (int i = 0; i < c_searchS; ++i){
currDiff = InnerDiff(searchCenter + m_searchOffsets[i], m_patternOffsets, c_patternS); if (currDiff > maxDiff){
maxDiff = currDiff;
betterStart = m_searchOffsets[i];
} if (maxDiff > c_acceptInnerDiff){ break;
}
} if (maxDiff < c_minInnerDiff){ return false;
}

betterStart += searchStart; return true;
} // get the minnmum pattern diff with the 8 neighbors.
inline int32_t MinNeighborDiff(const uint8_t* pattern)
{ const int32_t threshDiff = c_patternS * c_acceptDiff; // eight neighbors of a pattern
const int neighborOffsets[8] = { -1, 1, -m_lineBytes, m_lineBytes, -m_lineBytes - 1, -m_lineBytes + 1, m_lineBytes - 1, m_lineBytes + 1 }; int minDiff = PatternDiff(pattern, pattern + neighborOffsets[0], m_patternOffsets, c_patternS); if (minDiff < threshDiff){ return minDiff;
} int diff; for (int i = 1; i < 8; ++i){
diff = PatternDiff(pattern, pattern + neighborOffsets[i], m_patternOffsets, c_patternS); if (diff < minDiff){
minDiff = diff; if (minDiff < threshDiff){ return minDiff;
}
}
} return minDiff;
} // search the pattern that have max min_diff with neighbors, image is the full image begining, the return value betterStart is base on the image begining.
inline bool SearchMaxNeighborDiff(const uint8_t* image, int searchStart, int& betterStart)
{ const uint8_t* searchCenter = image + searchStart;
int32_t currDiff = MinNeighborDiff(searchCenter);
int32_t maxDiff = currDiff;
betterStart = 0; for (int i = 0; i < c_searchS; ++i){
currDiff = MinNeighborDiff(searchCenter + m_searchOffsets[i]); if (currDiff > maxDiff){
maxDiff = currDiff;
betterStart = m_searchOffsets[i];
}
} if (maxDiff <= c_patternS * c_acceptDiff){ return false;
}

betterStart += searchStart; return true;
} // match the target pattern in the image, return the best match quality and matched offset; the pattern is the center, image is the full image begining.
inline int32_t SearchBestMatch(const uint8_t* target, const int* patternPoints, const int patternSize, const uint8_t* image, int searchStart, int& matchedOffset)
{ const int thinkMatchedDiff = patternSize * c_acceptDiff; const uint8_t* searchCenter = image + searchStart; const uint8_t* matched = searchCenter;
int32_t currDiff = PatternDiff(target, matched, patternPoints, patternSize);
int32_t minDiff = currDiff; for (int i = 0; i < c_searchS; ++i){
currDiff = PatternDiff(target, searchCenter + m_searchOffsets[i], patternPoints, patternSize); if (currDiff < minDiff){
minDiff = currDiff;
matched = searchCenter + m_searchOffsets[i];
} if (minDiff < thinkMatchedDiff){ break;
}
}

matchedOffset = matched - image; return minDiff;
} int m_width, m_height, m_lineBytes; static const int c_acceptDiff = 2; // if the average pixel error is less than this number, think already matched
static const int c_rejectDiff = 8; // if the average pixel error is larger than this number, think it's not matched // all address offset to the pattern key location, the size is according to the square pattern.
static const int c_patternN = 3; static const int c_patternS = (2 * c_patternN + 1) * (2 * c_patternN + 1); int m_patternOffsets[c_patternS]; // the offsets to the image start for each seed point, the match is around these seed points.
static const int c_watchN = 2; static const int c_watchS = (2 * c_watchN + 1) * (2 * c_watchN + 1); static const int c_watchG = 30; // The gap of the watch grid in height direction
int m_watchPoints[c_watchS]; // the search offset to the search center, match the pattern from the corresponding location to the max distance. (not include distance 0.)
static const int c_searchD = 10; // search street-distance from the key location
static const int c_searchS = 2 * c_searchD * c_searchD + 2 * c_searchD; int m_searchOffsets[c_searchS]; // The implements table that use various pattern
static int m_impCurr; static int m_impNum; static const int c_maxImpNum = 16; static COpticalFlow_MPM* m_impTbl[c_maxImpNum];
};// save the optical flow resultuint8_t COpticalFlow_MPM::m_quality; // 0 ~ 255, 0 means the optical flow is invalid.float COpticalFlow_MPM::m_offset_x; // unit is pixelfloat COpticalFlow_MPM::m_offset_y;// the implements that use different patternint COpticalFlow_MPM::m_impCurr = 0;int COpticalFlow_MPM::m_impNum = 0;
COpticalFlow_MPM* COpticalFlow_MPM::m_impTbl[COpticalFlow_MPM::c_maxImpNum];// Multi-Pattern-Match-Squareclass COpticalFlow_MPMS : public COpticalFlow_MPM
{public:
COpticalFlow_MPMS(){} virtual ~COpticalFlow_MPMS(){} virtual const char* Name() { return "Square"; }protected: // prepare the address offset tables, that can make the calculation simple and fast.
virtual void GeneratePatternTable()
{ // generate the address offset of the match area to the center of the area.
int index = 0; for (int y = -c_patternN; y <= c_patternN; ++y){ for (int x = -c_patternN; x <= c_patternN; ++x){
m_patternOffsets[index++] = y * m_lineBytes + x;
}
}
}
};// Multi-Pattern-Match-Xclass COpticalFlow_MPMX : public COpticalFlow_MPM
{public:
COpticalFlow_MPMX(){} virtual ~COpticalFlow_MPMX(){} virtual const char* Name() { return "X"; }protected: // prepare the address offset tables, that can make the calculation simple and fast.
virtual void GeneratePatternTable()
{ // generate the address offset of the match area to the center of the area.
int index = 0; int armLen = (c_patternS - 1) / 4; for (int y = -armLen; y <= armLen; ++y){ if(y == 0){
m_patternOffsets[index++] = 0;
}else{
m_patternOffsets[index++] = y * m_lineBytes - y;
m_patternOffsets[index++] = y * m_lineBytes + y;
}
}
}
};static COpticalFlow_MPMS of_mpms;static COpticalFlow_MPMX of_mpmx;void OpticalFlow::init()
{ // set the optical flow implementation table
COpticalFlow_MPM::AddImplementation(&of_mpms);
COpticalFlow_MPM::AddImplementation(&of_mpmx);
COpticalFlow_MPM::SetImageDimesion(m_width, m_height, m_lineBytes);
}

uint32_t OpticalFlow::flow_image_in(const uint8_t *buf, int len, uint8_t *quality, int32_t *centi_pixel_x, int32_t *centi_pixel_y)
{ static uint8_t s_imageBuff1[m_pixelNum]; static uint8_t s_imageBuff2[m_pixelNum]; static uint8_t* s_imagePre = NULL; static uint8_t* s_imageCurr = s_imageBuff1; *quality = 0; *centi_pixel_x = 0; *centi_pixel_y = 0;

memcpy(s_imageCurr, buf, len); // first image
if(s_imagePre == NULL){
s_imagePre = s_imageCurr;
s_imageCurr = s_imageCurr == s_imageBuff1 ? s_imageBuff2 : s_imageBuff1; // switch image buffer
return 0;
}

COpticalFlow_MPM::AutoOpticalFlow(s_imagePre, s_imageCurr); if(COpticalFlow_MPM::m_quality > 0){ *quality = COpticalFlow_MPM::m_quality; *centi_pixel_x = (int32_t)(COpticalFlow_MPM::m_offset_x * 100); *centi_pixel_y = (int32_t)(COpticalFlow_MPM::m_offset_y * 100);
}

s_imagePre = s_imageCurr;
s_imageCurr = s_imageCurr == s_imageBuff1 ? s_imageBuff2 : s_imageBuff1; // switch image buffer
return 0;
}