光流法

光流法的介绍（含Ｃ＋＋代码）(转载)

默认分类 2009-07-11 17:41:48 阅读30 评论0   字号：大中小 订阅

光流法是比较经典的运动估计方法，本文不仅叙述简单明了，而且附代码，故收藏．

在空间中，运动可以用运动场描述。而在一个图像平面上，物体的运动往往是通过图像序列中不同图象灰度分布的不同体现的。从而，空间中的运动场转移到图像上就表示为光流场，光流场反映了图像上每一点灰度的变化趋势。

光流可以看作带有灰度的像素点在图像平面运动产生的瞬时速度场。下面我们推导光流方程：

假设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);

    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)

{

 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 );

      

     }

    }

   }

  }

 cvShowImage( "Contour", abs_img);

 cvShowImage( "Contour2", vely);

 cvReleaseImage(&velx);

 cvReleaseImage(&vely);

 cvWaitKey(20);

 

 return WW_OK;

}