光流场法的基本思想：在空间中，运动可以用运动场描述，而在一个图像平面上，物体的运动往往是通过图像序列中不同图像灰度分布的不同体现的，从而，空间中的运动场转移到图像上就表示为光流场(Optical Flow Field)。光流场反映了图像上每一点灰度的变化趋势，可看成是带有灰度的像素点在图像平面上运动而产生的瞬时速度场，也是一种对真实运动场的近似估计。
在比较理想的情况下，它能够检测独立运动的对象，不需要预先知道场景的任何信息，可以很精确地计算出运动物体的速度，并且可用于摄像机运动的情况。但光流法存在下面的缺点：有时即使没有发生运动，在外部照明发生变化时，也可以观测到光流；另外，在缺乏足够的灰度等级变化的区域，实际运动也往往观测不到。三维物体的运动投影到二维图像的亮度变化，本身由于部分信息的丢失而使光流法存在孔径问题和遮挡问题，用光流法估算二维运动场是不确定的，需要附加的假设模型来模拟二维运动场的结构；在准确分割时，光流法还需要利用颜色、灰度、边缘等空域特征来提高分割精度；同时由于光流法采用迭代的方法，计算复杂耗时，如果没有特殊的硬件支持，很难应用于视频序列的实时检测。
推导光流方程过程：
假设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 for(j=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);
//this postion the continual pix add together and get the mean 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);
//last pic do the same processing
Et = gray1 - gray2; //Et means that the difference between the two pic at the same point
//next point in the horizontal directiongray2 = 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; //x direction (next point ) the difference in the current pic
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; //y direction (next point ) the difference in the current pic
Ex = ((int)(Ex/10))*10; //int 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 mean that 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;
}
