|
光流法是比較經(jīng)典的運(yùn)動(dòng)估計(jì)方法，本文不僅敘述簡(jiǎn)單明了，而且附代碼，故收藏．
在空間中，運(yùn)動(dòng)可以用運(yùn)動(dòng)場(chǎng)描述。而在一個(gè)圖像平面上，物體的運(yùn)動(dòng)往往是通過(guò)圖像序列中不同圖象灰度分布的不同體現(xiàn)的。從而，空間中的運(yùn)動(dòng)場(chǎng)轉(zhuǎn)移到圖像上就表示為光流場(chǎng)，光流場(chǎng)反映了圖像上每一點(diǎn)灰度的變化趨勢(shì)。
光流可以看作帶有灰度的像素點(diǎn)在圖像平面運(yùn)動(dòng)產(chǎn)生的瞬時(shí)速度場(chǎng)。下面我們推導(dǎo)光流方程：
假設(shè)E(x,y,t)為(x,y)點(diǎn)在時(shí)刻t的灰度（照度）。設(shè)t+dt時(shí)刻該點(diǎn)運(yùn)動(dòng)到(x+dx,y+dy)點(diǎn)，他的照度為E(x+dx,y+dy,t+dt)。我們認(rèn)為，由于對(duì)應(yīng)同一個(gè)點(diǎn)，所以
E(x,y,t) =
E(x+dx,y+dy,t+dt) ——
光流約束方程
將上式右邊做泰勒展開(kāi)，并令dt->0，則得到：Exu+Eyv+Et
= 0，其中：
Ex =
dE/dx Ey =
dE/dy Et =
dE/dt u =
dx/dt v = dy/dt
上面的Ex,Ey,Et的計(jì)算都很簡(jiǎn)單，用離散的差分代替導(dǎo)數(shù)就可以了。光流法的主要任務(wù)就是通過(guò)求解光流約束方程求出u,v。但是由于只有一個(gè)方程，所以這是個(gè)病態(tài)問(wèn)題。所以人們提出了各種其他的約束方程以聯(lián)立求解。但是由于我們用于攝像機(jī)固定的這一特定情況，所以問(wèn)題可以大大簡(jiǎn)化。
攝像機(jī)固定的情形
在攝像機(jī)固定的情形下，運(yùn)動(dòng)物體的檢測(cè)其實(shí)就是分離前景和背景的問(wèn)題。我們知道對(duì)于背景，理想情況下，其光流應(yīng)當(dāng)為0，只有前景才有光流。所以我們并不要求通過(guò)求解光流約束方程求出u,v。我么只要求出亮度梯度方向的速率就可以了，即求出sqrt(u*u+v*v)。
而由光流約束方程可以很容易求到梯度方向的光流速率為 V =
abs(Et/sqrt(Ex*Ex+Ey*Ey))。這樣我們?cè)O(shè)定一個(gè)閾值T。
V(x,y)
> T 則(x,y)是前景 ，反之是背景
C++實(shí)現(xiàn)
在實(shí)現(xiàn)攝像機(jī)固定情況的光流法時(shí)，需要有兩幀連續(xù)的圖像，下面的算法針對(duì)RGB24格式的圖像計(jì)算光流：
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);
}
//////////////////////////////////////////////////////////////////////////
/另一個(gè)代碼
/
/////////////////////////////////////////////////////////////////////////////
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
;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; //計(jì)算光流的窗口大小
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 );
}
}
}
}
cvShowImage( "Contour",
abs_img);
cvShowImage( "Contour2", vely);
cvReleaseImage(&velx);
cvReleaseImage(&vely);
cvWaitKey(20);
return WW_OK;
}
|
