光流法c++完整代码

1、光流法是比较经典的运动估计方法，本文不仅叙述简单明了，而且附代码，故收藏.在空间中，运动可以用运动场描述。而在一个图像平面上，物体的运动往往是通过图像序列中不同图象灰 度分布的不同体现的。从而，空间中的运动场转移到图像上就表示为光流场，光流场反映了图像上每一点 灰度的变化趋势。光流可以看作带有灰度的像素点在图像平面运动产生的瞬时速度场。下面我们推导光流方程：假设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+

2、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。但是由于只有一个方程，所以这是个病态问题。所以人们提出了各种其他的约束方程以 联立求解。但是由于我们用于摄像机固定的这一特定情况，所以问题可以大大简化。摄像机固定的情形在摄像机固定的情形下，运动物体的检测其实就是分离前景和背景的问题。我们知道对于背景，理想情况 下，其光流应当

3、为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;mems

4、et(opticalflow,0,width*height*sizeof(int);memset(output,255,size);for(i=2;iheight-2;i+)(for(j=2;j10)(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:Descrip

5、tion:光流法计算运动速度与方向Date:2006-6-14Author: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,

6、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.val0;

7、float y_shift= (float)yc.val0;const int winsize=5; 计算光流的窗口大小if(x%(winsize*2)=0) & (y%(winsize*2)=0)(if(x_shift!=0 | y_shift!=0)(if(xwinsize & ywinsize & x width-winsize) & yheight-winsize) (cvSetImageROI( velx, cvRect( x-winsize, y-winsize, 2*winsize, 2*winsize);CvScalar total_x = cvSum(velx);float

8、xx = (float)total_x.val0;cvResetlmageROI(velx);cvSetImageROI( vely, cvRect( x-winsize, y-winsize, 2*winsize, 2*winsize);CvScalar total_y = cvSum(vely);float yy = (float)total_y.val0;cvResetImageROI(vely);cvSetImageROI( abs_img, cvRect( x-winsize, y-winsize, 2*winsize, 2*winsize);CvScalar total_speed

9、 = cvSum(abs_img);float ss = (float)total_speed.val0/(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;elsealpha_angle = abs(atan(yy/xx);if(xx0) alpha_angle = pi - alpha_angle ;if(xx0 & yy0 & yy0

10、) 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

11、_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;cvR