轉自：http://blog.163.com/prevBlogPerma.do?host=quanhuatang&srl=222682420080137347226&mode=prev
光流法是比較經典的運動估計方法，本文不僅敘述簡單明了，而且附代碼，故收藏．
在空間中，運動可以用運動場描述。而在一個圖像平面上，物體的運動往往是通過圖像序列中不同圖象灰度分布的不同體現的。從而，空間中的運動場轉移到圖像上就表示為光流場，光流場反映了圖像上每一點灰度的變化趨勢。
光流可以看作帶有灰度的像素點在圖像平面運動產生的瞬時速度場。下面我們推導光流方程：
假設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 ;y<velx->height; y++)
for(int x =0;x<velx->width;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 && 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<int>(x + scale_factor*cos(alpha_angle));
pt2.y = static_cast<int>(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;}