光流法是比较经典的运动估计方法，本文不仅叙述简单明了，而且附代码，故收藏．
光流是空间运动物体在观测成像面上的像素运动的瞬时速度。光流的研究是利用图像序列中的像素强度数据的时域变化和相关性来确定各自像素位置的“运动”，即研究图像灰度在时间上的变化与景象中物体结构及其运动的关系。一般情况下，光流由相机运动、场景中目标运动或两者的共同运动产生。光流计算方法大致可分为三类：基于匹配的、频域的和梯度的方法。
(1) 基于匹配的光流计算方法包括基于特征和基于区域两种。基于特征的方法不断地对目标主要特征进行定位和跟踪，对大目标的运动和亮度变化具有鲁棒性。存在的问题是光流通常很稀疏，而且特征提取和精确匹配也十分困难。基于区域的方法先对类似的区域进行定位，然后通过相似区域的位移计算光流。这种方法在视频编码中得到了广泛的应用。然而，它计算的光流仍不稠密。
(2) 基于频域的方法利用速度可调的滤波组输出频率或相位信息。虽然能获得高精度的初始光流估计，但往往涉及复杂的计算。另外，进行可靠性评价也十分困难。
(3) 基于梯度的方法利用图像序列的时空微分计算2D速度场（光流）。由于计算简单和较好的效果，基于梯度的方法得到了广泛的研究。虽然很多基于梯度的光流估计方法取得了较好的光流估计，但由于在计算光流时涉及到可调参数的人工选取、可靠性评价因子的选择困难，以及预处理对光流计算结果的影响，在应用光流对目标进行实时监测与自动跟踪时仍存在很多问题。
光流法检测运动物体的基本原理是：给图像中的每一个像素点赋予一个速度矢量，这就形成了一个图像运动场，在运动的一个特定时刻，图像上的点与三维物体上的点一一对应，这种对应关系可由投影关系得到，根据各个像素点的速度矢量特征，可以对图像进行动态分析。如果图像中没有运动物体，则光流矢量在整个图像区域是连续变化的。当图像中有运动物体时，目标和图像背景存在相对运动，运动物体所形成的速度矢量必然和邻域背景速度矢量不同，从而检测出运动物体及位置。采用光流法进行运动物体检测的问题主要在于大多数光流法计算耗时，实时性和实用性都较差。但是光流法的优点在于光流不仅携带了运动物体的运动信息，而且还携带了有关景物三维结构的丰富信息，它能够在不知道场景的任何信息的情况下，检测出运动对象。
对于视频监控系统来说，所用的图像基本都是摄像机静止状态下摄取得，所以对有实时性和准确性要求的系统来说，纯粹使用光流法来检测目标不太实际。更多的是利用光流计算方法与其它方法相结合来实现对目标检测和运动估计。
然而，在实际应用中，由于遮挡性、多光源、透明性和噪声等原因，使得光流场基本方程的灰度守恒假设条件不能满足，不能求解出正确的光流场，同时大多数的光流计算方法相当复杂，计算量巨大，不能满足实时的要求，因此，一般不被对精度和实时性要求比较高的监控系统所采用。
在空间中，运动可以用运动场描述。而在一个图像平面上，物体的运动往往是通过图像序列中不同图象灰度分布的不同体现的。从而，空间中的运动场转移到图像上就表示为光流场，光流场反映了图像上每一点灰度的变化趋势。
光流可以看作带有灰度的像素点在图像平面运动产生的瞬时速度场。下面我们推导光流方程：
假设E(x,y,t)为(x,y)点在时刻t的灰度（照度）。设t+dt时刻该点运动到(x+dx,y+dy)点，他的照度为E(x+dx,y+dy,t+dt)。我们认为，由于对应同一个点，所以
E(x,y,t) = E(x+dx,y+dy,t+dt) —— 光流约束方程
将上式右边做泰勒展开，并令dt->0，则得到：Exu+Eyv+Et = 0，其中：
Ex = dE/dx Ey = dE/dy Et = dE/dt u = dx/dt v = dy/dt
上面的Ex,Ey,Et的计算都很简单，用离散的差分代替导数就可以了。光流法的主要任务就是通过求解光流约束方程求出u,v。但是由于只有一个方程，所以这是个病态问题。所以人们提出了各种其他的约束方程以联立求解。但是由于我们用于摄像机固定的这一特定情况，所以问题可以大大简化。
摄像机固定的情形
在摄像机固定的情形下，运动物体的检测其实就是分离前景和背景的问题。我们知道对于背景，理想情况下，其光流应当为0，只有前景才有光流。所以我们并不要求通过求解光流约束方程求出u,v。我么只要求出亮度梯度方向的速率就可以了，即求出sqrt(u*u+v*v)。
而由光流约束方程可以很容易求到梯度方向的光流速率为 V = abs(Et/sqrt(Ex*Ex+Ey*Ey))。这样我们设定一个阈值T。
V(x,y) > T 则(x,y)是前景 ，反之是背景
C++实现
在实现摄像机固定情况的光流法时，需要有两帧连续的图像，下面的算法针对RGB24格式的图像计算光流：
void calculate(unsigned char* buf)
{
int Ex,Ey,Et;
int gray1,gray2;
int u;
int i,j;
memset(opticalflow,0,width*height*sizeof(int));
memset(output,255,size);
for(i=2;i<height-2;i++){
for(j=2;j<width-2;j++){
gray1 = int(((int)(buf[(i*width+j)*3])
+(int)(buf[(i*width+j)*3+1])
+(int)(buf[(i*width+j)*3+2]))*1.0/3);
gray2 = int(((int)(prevframe[(i*width+j)*3])
+(int)(prevframe[(i*width+j)*3+1])
+(int)(prevframe[(i*width+j)*3+2]))*1.0/3);
Et = gray1 - gray2;
gray2 = int(((int)(buf[(i*width+j+1)*3])
+(int)(buf[(i*width+j+1)*3+1])
+(int)(buf[(i*width+j+1)*3+2]))*1.0/3);
Ex = gray2 - gray1;
gray2 = int(((int)(buf[((i+1)*width+j)*3])
+(int)(buf[((i+1)*width+j)*3+1])
+(int)(buf[((i+1)*width+j)*3+2]))*1.0/3);
Ey = gray2 - gray1;
Ex = ((int)(Ex/10))*10;
Ey = ((int)(Ey/10))*10;
Et = ((int)(Et/10))*10;
u = (int)((Et*10.0)/(sqrt((Ex*Ex+Ey*Ey)*1.0))+0.1);
opticalflow[i*width+j] = u;
if(abs(u)>10){
output[(i*width+j)*3] = 0;
output[(i*width+j)*3+1] = 0;
output[(i*width+j)*3+2] = 0;
}
}
}
memcpy(prevframe,buf,size);
}
//////////////////////////////////////////////////////////////////////////
/另一个代码
/
/////////////////////////////////////////////////////////////////////////////
WW_RETURN HumanMotion::ImgOpticalFlow(IplImage *pre_grey,IplImage *grey)
/*************************************************
Function:
Description: 光流法计算运动速度与方向
Date: 2006－6－14
Author:
Input:
Output:
Return:
Others:
*************************************************/
{
IplImage *velx = cvCreateImage( cvSize(grey->width ,grey->height),IPL_DEPTH_32F, 1 );
IplImage *vely = cvCreateImage( cvSize(grey->width ,grey->height),IPL_DEPTH_32F, 1 );
velx->origin = vely->origin = grey->origin;
CvSize winSize = cvSize(5,5);
cvCalcOpticalFlowLK( prev_grey, grey, winSize, velx, vely );
cvAbsDiff( grey,prev_grey, abs_img );
cvThreshold( abs_img, abs_img, 29, 255, CV_THRESH_BINARY);
CvScalar xc,yc;
for(int y =0 ;yheight; y++)
for(int x =0;xwidth;x++ )
{
xc = cvGetAt(velx,y,x);
yc = cvGetAt(vely,y,x);
float x_shift= (float)xc.val[0];
float y_shift= (float)yc.val[0];
const int winsize=5; //计算光流的窗口大小
if((x%(winsize*2)==0) && (y%(winsize*2)==0) )
{
if(x_shift!=0 || y_shift!=0)
{
if(x>winsize && y>winsize && x <(velx->width-winsize) && y<(velx->height-winsize) )
{
cvSetImageROI( velx, cvRect( x-winsize, y-winsize, 2*winsize, 2*winsize));
CvScalar total_x = cvSum(velx);
float xx = (float)total_x.val[0];
cvResetImageROI(velx);
cvSetImageROI( vely, cvRect( x-winsize, y-winsize, 2*winsize, 2*winsize));
CvScalar total_y = cvSum(vely);
float yy = (float)total_y.val[0];
cvResetImageROI(vely);
cvSetImageROI( abs_img, cvRect( x-winsize, y-winsize, 2*winsize, 2*winsize));
CvScalar total_speed = cvSum(abs_img);
float ss = (float)total_speed.val[0]/(4*winsize*winsize)/255;
cvResetImageROI(abs_img);
const double ZERO = 0.000001;
const double pi = 3.1415926;
double alpha_angle;
if(xx-ZERO)
alpha_angle = pi/2;
else
alpha_angle = abs(atan(yy/xx));
if(xx<0 && yy>0) alpha_angle = pi - alpha_angle ;
if(xx<0 && yy<0) alpha_angle = pi + alpha_angle ;
if(xx>0 && yy<0) alpha_angle = 2*pi - alpha_angle ;
CvScalar line_color;
float scale_factor = ss*100;
line_color = CV_RGB(255,0,0);
CvPoint pt1,pt2;
pt1.x = x;
pt1.y = y;
pt2.x = static_cast(x + scale_factor*cos(alpha_angle));
pt2.y = static_cast(y + scale_factor*sin(alpha_angle));
cvLine( image, pt1, pt2 , line_color, 1, CV_AA, 0 );
CvPoint p;
p.x = (int) (pt2.x + 6 * cos(alpha_angle - pi / 4*3));
p.y = (int) (pt2.y + 6 * sin(alpha_angle - pi / 4*3));
cvLine( image, p, pt2, line_color, 1, CV_AA, 0 );
p.x = (int) (pt2.x + 6 * cos(alpha_angle + pi / 4*3));
p.y = (int) (pt2.y + 6 * sin(alpha_angle + pi / 4*3));
cvLine( image, p, pt2, line_color, 1, CV_AA, 0 );
/*
line_color = CV_RGB(255,255,0);
pt1.x = x-winsize;
pt1.y = y-winsize;
pt2.x = x+winsize;
pt2.y = y+winsize;
cvRectangle(image, pt1,pt2,line_color,1,CV_AA,0);
*/
}
}
}
}
cvShowImage( "Contour", abs_img);
cvShowImage( "Contour2", vely);
cvReleaseImage(&velx);
cvReleaseImage(&vely);
cvWaitKey(20);
return WW_OK;