opencv 光流法详解

一.基本概念

光流的概念是Gibson于1950年提出的。所谓光流是指图像中模式运动的速度，光流场是一种二维(2D)瞬时速度场，其中二维速度向量是可见的三维速度向量在成像平面上的投影。光流法是把检测区域的图像变为速度的矢量场，每一个向量表示了景物中一个点在图像中位置的瞬时变化。因此，光流场携带了有关物体运动和景物三维结构的丰富信息，通过对速度场(光流场)的分析可以判断在检测区域内车辆的有无。

思路：求得整个图像检测区域的速度场，根据每个像素点的速度向量特征，可以对图像进行分析。如果是静止背景即无车通过时，则光流向量在整个检测区域是连续变化的；当有车通过时，光流向量必然和其邻域背景的光流向量不同，从而检测出车辆及其出现的位置。

光流法的前提假设：
（1）相邻帧之间的亮度恒定值，
（2）相邻视频帧的取帧的时间连续，或者相邻帧之间物体的运动比较“微小”;
（3）保持空间一致性，即，同一子图像的像素点具有相同的运动。

原理：
（1）对一个连续的视频帧序列进行处理；
（2）针对每一个视频序列，利用一定的目标检测方法，检测可能出现的前景目标；
（3）如果某一帧出现了前景目标，找到其具有代表性的关键特征点（如shi-Tomasi算法）；
（4）对之后的任意两个相邻视频帧而言，寻找上一帧中出现的关键特征点在当前帧中的最佳位置，从而得到前景目标在当前帧中的位置坐标；
（5）如此迭代进行，便可实现目标的跟踪；

二.程序中重要函数解释

      （1）void cvGoodFeaturesToTrack( const CvArr* image, CvArr* eig_image, CvArr* temp_image,CvPoint2D32f* corners, int* corner_count,double quality_level, double                   min_distance,const CvArr* mask=NULL );

                程序中此函数代码：cvGoodFeaturesToTrack(frame1_1C, eig_image, temp_image, frame1_features, &number_of_features, .01, .01, NULL);
                这是shi-Tomasi算法，该算法主要用于feature selection,即一张图中哪些是我 们感兴趣需要跟踪的点(interest point) 
                input:                     
                     * "frame1_1C" 输入图像. 
                     * "eig_image" and "temp_image" 只是给该算法提供可操作的内存区域
                     * 第一个".01" 规定了特征值的最小质量，因为该算法要得到好的特征点，哪就需要一个选择的阈值
                     * 第二个".01" 规定了像素之间最小的距离，用于减少运算复杂度，当然也一定程度降低了跟踪精度
                     * "NULL" 意味着处理整张图片，当然你也可以指定一块区域
                output: 
                     * "frame1_features" 将会包含fram1的特征值 
                     * "number_of_features" 将在该函数中自动填充上所找到特征值的真实数目

       （2）迭代算法的终止准则 

                 typedef struct CvTermCriteria  
                {  
   int    type;  /* CV_TERMCRIT_ITER 和CV_TERMCRIT_EPS二值之一，或者二者的组合 */  
  int    max_iter; /* 最大迭代次数 */  
           double epsilon; /* 结果的精确性 */  
                 } 

       （3） void cvCalcOpticalFlowPyrLK( const CvArr* prev, const CvArr* curr, CvArr* prev_pyr,CvArr* curr_pyr, const CvPoint2D32f* prev_features, CvPoint2D32f*curr_features, int count, CvSize win_size, int level, char* status,  float*track_error, CvTermCriteria criteria, int flags );

   计算一个稀疏特征集的光流，使用金字塔中的迭代 Lucas-Kanade 方法。参数：*prev在时间 t 的第一帧；*curr在时间 t + dt 的第二帧；*prev_pyr第一帧的金字塔缓存. 如果指针非 NULL , 则缓存必须有足够的空间来存储金字塔从层 1 到层 #level 的内容。尺寸 (image_width+8)*image_height/3 比特足够了;*curr_pyr与 prev_pyr 类似， 用于第二帧；*prev_features需要发现光流的点集；*curr_features包含新计算出来的位置的 点集；*count特征点的数目；*win_size每个金字塔层的搜索窗口尺寸；*level最大的金字塔层数。如果为 0 , 不使用金字塔 (即金字塔为单层), 如果为 1 , 使用两层，下面依次类推；*status数组。如果对应特征的光流被发现，数组中的每一个元素都被设置为 1，否则设置为 0；*error双精度数组，包含原始图像碎片与移动点之间的差。为可选参数，可以是 NULL ；*criteria准则，指定在每个金字塔层，为某点寻找光流的迭代过程的终止条件；*flags其它选项:(1)CV_LKFLOW_PYR_A_READY , 在调用之前，第一帧的金字塔已经准备好(2)CV_LKFLOW_PYR_B_READY , 在调用之前，第二帧的金字塔已经准备好(3)CV_LKFLOW_INITIAL_GUESSES , 在调用之前，数组 B 包含特征的初始坐标 （Hunnish: 在本节中没有出现数组 B，不知是指的哪一个）

函数 cvCalcOpticalFlowPyrLK 实现了金字塔中 Lucas-Kanade 光流计算的稀疏迭代版本([Bouguet00])。它根据给出的前一帧特征点坐标计算当前视频帧上的特征点坐标。函数寻找具有子象素精度的坐标值。

两个参数 prev_pyr 和 curr_pyr 都遵循下列规则：如果图像指针为 0, 函数在内部为其分配缓存空间，计算金字塔，然后再处理过后释放缓存。否则，函数计算金字塔且存储它到缓存中，除非设置标识 CV_LKFLOW_PYR_A[B]_READY 。 图像应该足够大以便能够容纳Gaussian 金字塔数据。调用函数以后，金字塔被计算而且相应图像的标识可以被设置，为下一次调用准备就绪 (比如：对除了第一个图像的所有图像序列，标识 CV_LKFLOW_PYR_A_READY 被设置).

三.程序源代码

[cpp] view plaincopy

1. #include <stdio.h>  
2. #include <cv.h>  
3. #include <highgui.h>  
4. #include <math.h>  
5. static const double pi = 3.14159265358979323846;  
6.   
7.   
8. inline static double square(int a)  
9. {  
10.     return a * a;  
11. }  
12.   
13.   
14. /*此函数主要目的：给img分配内存空间（除非此图像已经非NULL），并设定图像大小、位深以及通道数。 
15.   即使该图像的大小、深度或信道数与要求的不同，也会创建一个非NULL图像。*/  
16. inline static void allocateOnDemand( IplImage **img, CvSize size, int depth, int channels)  
17. {  
18.     if ( *img != NULL ) return;  
19.         *img = cvCreateImage( size, depth, channels );  
20.     if ( *img == NULL )  
21.     {  
22.         fprintf(stderr, "Error: Couldn't allocate image. Out of memory?\n");  
23.         exit(-1);  
24.     }  
25. }  
26.   
27.   
28. int main(void)  
29. {  
30.     CvCapture *input_video = cvCaptureFromFile("video1.avi");  //创建一个对象，读取avi视频。  
31.     if (input_video == NULL)  
32.     {  
33.         fprintf(stderr, "Error: Can't open video.\n");  //视频不存在或格式不支持时显示  
34.         return -1;  
35.     }  
36.   
37.     cvQueryFrame( input_video );  // 读取一帧是为了获得帧的长宽。  
38.                                       
39.     CvSize frame_size;  
40.     frame_size.height = (int) cvGetCaptureProperty( input_video, CV_CAP_PROP_FRAME_HEIGHT );  
41.     frame_size.width =  (int) cvGetCaptureProperty( input_video, CV_CAP_PROP_FRAME_WIDTH );  
42.   
43.     long number_of_frames;  //视频帧的长度  
44.   
45.     cvSetCaptureProperty( input_video, CV_CAP_PROP_POS_AVI_RATIO, 1. );  //跳到视频结束，以便获取视频长度（帧数）  
46.     number_of_frames = (int) cvGetCaptureProperty( input_video, CV_CAP_PROP_POS_FRAMES );  //得到帧数  
47.   
48.     cvSetCaptureProperty( input_video, CV_CAP_PROP_POS_FRAMES, 0. );  //重新回到视频开始，以便下续步骤进行  
49.   
50.       
51. /*开始进行光流法*/  
52.     cvNamedWindow("Optical Flow", CV_WINDOW_AUTOSIZE);  //创建窗口，显示输出，大小自动调节  
53.     long current_frame = 0;  
54.     while(true)  
55.     {  
56.         static IplImage *frame = NULL, *frame1 = NULL, *frame1_1C = NULL, *frame2_1C =  
57.         NULL, *eig_image = NULL, *temp_image = NULL, *pyramid1 = NULL, *pyramid2 = NULL;  
58.   
59.           
60.         cvSetCaptureProperty( input_video, CV_CAP_PROP_POS_FRAMES, current_frame );  
61.       
62.         frame = cvQueryFrame( input_video );  //读取第一帧  
63.         if (frame == NULL)  
64.         {  
65.             fprintf(stderr, "Error: Hmm. The end came sooner than we thought.\n");   
66.             return -1;  
67.         }  
68.   
69.         allocateOnDemand( &frame1_1C, frame_size, IPL_DEPTH_8U, 1 );  //分配给frame1_1C内存空间  
70.         cvConvertImage(frame, frame1_1C, CV_CVTIMG_FLIP);  //将帧数据赋给frame1_1C  
71.   
72.         allocateOnDemand( &frame1, frame_size, IPL_DEPTH_8U, 3 );  //把具有全部颜色信息的原帧保存，以备最后在屏幕上显示用  
73.         cvConvertImage(frame, frame1, CV_CVTIMG_FLIP);  
74.   
75.         frame = cvQueryFrame( input_video );  //读取第二帧  
76.         if (frame == NULL)  
77.         {  
78.             fprintf(stderr, "Error: Hmm. The end came sooner than we thought.\n");  
79.             return -1;  
80.         }  
81.         allocateOnDemand( &frame2_1C, frame_size, IPL_DEPTH_8U, 1 );  
82.         cvConvertImage(frame, frame2_1C, CV_CVTIMG_FLIP);  
83.   
84.     /* 施和托马斯特征跟踪! */  
85.   
86.         allocateOnDemand( &eig_image, frame_size, IPL_DEPTH_32F, 1 );  //分配需要的内存  
87.         allocateOnDemand( &temp_image, frame_size, IPL_DEPTH_32F, 1 );  //分配需要的内存  
88.   
89.         CvPoint2D32f frame1_features[400];  //创建数组，存放一帧的特征  
90.   
91.         int number_of_features;  //函数运行前先设定特征数的最大值，运行后将是找到的特征数真正数量  
92.         number_of_features = 400;  //可以改变此值，折中精确性  
93.   
94.         cvGoodFeaturesToTrack(frame1_1C, eig_image, temp_image, frame1_features, &number_of_features, .01, .01, NULL);  //施和托马斯算法  
95.   
96.     /* 金字塔的 Lucas Kanade 光流法! */  
97.       
98.         CvPoint2D32f frame2_features[400];  //数组用于存放帧二中来自帧一的特征。  
99.         char optical_flow_found_feature[400];  //当且仅当frame1中的特征值在frame2中找到时，对应值为非零。  
100.         float optical_flow_feature_error[400];  //数组第i个元素表对应点光流误差  
101.           
102.         CvSize optical_flow_window = cvSize(3,3);  //lucas-kanade光流法运算窗口，可以去其他大小，不过运算量加大  
103.   
104.         CvTermCriteria optical_flow_termination_criteria = cvTermCriteria( CV_TERMCRIT_ITER | CV_TERMCRIT_EPS, 20, .3 ); //终止条件  
105.   
106.         allocateOnDemand( &pyramid1, frame_size, IPL_DEPTH_8U, 1 );  
107.         allocateOnDemand( &pyramid2, frame_size, IPL_DEPTH_8U, 1 );  
108.   
109.         //跑Lucas Kanade算法  
110.         cvCalcOpticalFlowPyrLK(frame1_1C, frame2_1C, pyramid1, pyramid2, frame1_features,  
111.                                frame2_features, number_of_features, optical_flow_window, 5,  
112.                                optical_flow_found_feature, optical_flow_feature_error,  
113.                                optical_flow_termination_criteria, 0 );    
114.   
115.     /* 画光流场，画图是依据两帧对应的特征值，这个特征值就是图像上我们感兴趣的点，如边缘上的点P(x,y)*/  
116.         for(int i = 0; i < number_of_features; i++)  
117.         {  
118.             if ( optical_flow_found_feature[i] == 0 ) continue;  //如果Lucas Kanade算法没找到特征点，跳过  
119.               
120.             int line_thickness;   
121.             line_thickness = 1;  
122.               
123.             CvScalar line_color;   
124.             line_color = CV_RGB(255,0,0);  
125.   
126.             /* 画箭头,因为帧间的运动很小，所以需要缩放，不然看不见箭头，缩放因子为3 */  
127.             CvPoint p,q;  
128.             p.x = (int) frame1_features[i].x;  
129.             p.y = (int) frame1_features[i].y;  
130.             q.x = (int) frame2_features[i].x;  
131.             q.y = (int) frame2_features[i].y;  
132.   
133.             double angle;  
134.             angle = atan2( (double) p.y - q.y, (double) p.x - q.x ); //方向  
135.             double hypotenuse;   
136.             hypotenuse = sqrt( square(p.y - q.y) + square(p.x - q.x) );  //长度  
137.   
138.             q.x = (int) (p.x - 3 * hypotenuse * cos(angle)); //放大三倍  
139.             q.y = (int) (p.y - 3 * hypotenuse * sin(angle));  
140.   
141.             cvLine( frame1, p, q, line_color, line_thickness, CV_AA, 0 );  //画箭头主体  
142.             //画箭的头部  
143.             p.x = (int) (q.x + 9 * cos(angle + pi / 4));  
144.             p.y = (int) (q.y + 9 * sin(angle + pi / 4));  
145.             cvLine( frame1, p, q, line_color, line_thickness, CV_AA, 0 );  
146.             p.x = (int) (q.x + 9 * cos(angle - pi / 4));  
147.             p.y = (int) (q.y + 9 * sin(angle - pi / 4));  
148.             cvLine( frame1, p, q, line_color, line_thickness, CV_AA, 0 );  
149.         }  
150.   
151.         cvShowImage("Optical Flow", frame1);  
152.   
153.         int key_pressed;  
154.         key_pressed = cvWaitKey(33);  //隔33ms显示一帧  
155.   
156.         /* 触发B键，视频回一帧 */  
157.         if (key_pressed == 'b' || key_pressed == 'B') current_frame--;  
158.         else current_frame++;  
159.   
160.         /* 不要超出视频的头和尾 */  
161.         if (current_frame < 0) current_frame = 0;  
162.         if (current_frame >= number_of_frames - 1) current_frame = number_of_frames - 2;  
163.     }  
164. }  

【整合】
【链接】http://download.csdn.net/detail/u012756029/6518397（《Optical Flow》 Written by David Stavens有ppt和注释，详细的英文解释）