光流法C++源代码

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