欢迎来到网际学院,让您的头脑满载而归!

【目标跟踪】基于粒子滤波算法的目标跟踪:Rob Hess源码分析

发布日期:2018-03-12 16:53:58 作者:管理员 阅读:471

前言:粒子滤波广泛的应用于目标跟踪,粒子滤波器是一种序列蒙特卡罗滤波方法,其实质是利用一系列随机抽取的样本(即粒子)来替代状态的后验概率分布。在此不打算介绍和推理繁杂的概率公式,我们来分析Rob Hess源码从而深入理解粒子滤波算法。试验平

【目标跟踪】基于粒子滤波算法的目标跟踪:Rob Hess源码分析

前言:

粒子滤波广泛的应用于目标跟踪,粒子滤波器是一种序列蒙特卡罗滤波方法,其实质是利用一系列随机抽取的样本(即粒子)来替代状态的后验概率分布。在此不打算介绍和推理繁杂的概率公式,我们来分析Rob Hess源码从而深入理解粒子滤波算法。


试验平台:

VS2010 + opencv2.4.10 + gsl1.8库 + RobHess粒子滤波源码

相关资料:

       Rob Hess粒子滤波的相关代码:http://blogs.oregonstate.edu/hess/code/particles/

       Rob Hess有关粒子滤波的文章:http://web.engr.oregonstate.edu/~afern/papers/cvpr09.pdf

       Rob Hess多目标跟踪效果视频:http://v.youku.com/v_show/id_XOTQ5NDA5ODgw.html

       yangyangcv博主的分析及代码(建议先看):http://www.cnblogs.com/yangyangcv/archive/2010/05/23/1742263.html

下载 yangyangcv提供的代码,配置好后即可使用;本博文也提供下载,见文章末尾


源码分析:

下面用粒子滤波算法来实现单目标的跟踪,在初始时需要手动选取待跟踪目标区域。

步骤一:粒子初始化

1、目标的选取:

在起始帧画面中标记待跟踪的区域(目标物体),在以后的连续帧视频中跟踪目标物体。

2、粒子的定义:

粒子即样本,Rob Hess源码中默认粒子数目PARTICLES=100,下面是粒子的属性(可以看出一个粒子代表的就是一个矩形区域):


[cpp] view plain copy

  1. typedef struct particle {  

  2.   float x;          /**< 当前x坐标 */  

  3.   float y;          /**< 当前y坐标 */  

  4.   float s;          /**< scale */  

  5.   float xp;         /**< 之前x坐标 */  

  6.   float yp;         /**< 之前y坐标 */  

  7.   float sp;         /**< previous scale */  

  8.   float x0;         /**< x0 */  

  9.   float y0;         /**< y0 */  

  10.   int width;        /**< 粒子宽度 */  

  11.   int height;       /**< 粒子高度 */  

  12.   histogram* histo; /**< 跟踪区域的参考直方图 */  

  13.   float w;          /**< 权重 */  

  14. } particle;  



3、特征提取:

目标跟踪最重要的就是特征,应该选取好的特征(拥有各种不变性的特征当然是最好的);另外考虑算法的效率,目标跟踪一般是实时跟踪,所以对算法实时性有一定的要求。Rob Hess源码提取的是目标的颜色特征(颜色特征对图像本身的尺寸、方向、视角的依赖性较小,从而具有较高的鲁棒性),粒子与目标的直方图越相似,则说明越有可能是目标。

捕捉到一帧图像,标记待跟踪的区域,将其从BGR转化到HSV空间,提取感兴趣部分(目标)的HSV,进行直方图统计并归一化直方图


[cpp] view plain copy

  1. /* increment appropriate histogram bin for each pixel */  

  2. //计算直方图  

  3. for( r = 0; r < img->height; r++ )  

  4.     for( c = 0; c < img->width; c++ )  

  5.     {  

  6.         bin = histo_bin( /*pixval32f( h, r, c )*/((float*)(h->imageData + h->widthStep*r) )[c],  

  7.                     ((float*)(s->imageData + s->widthStep*r) )[c],  

  8.                     ((float*)(v->imageData + v->widthStep*r) )[c] );  

  9.         hist[bin] += 1;  

  10.     }  






[cpp] view plain copy

  1. int histo_bin( float h, float s, float v )  

  2. {  

  3.   int hd, sd, vd;  

  4.   

  5.   /* if S or V is less than its threshold, return a "colorless" bin */  

  6.   vd = MIN( (int)(v * NV / V_MAX), NV-1 );  

  7.   if( s < S_THRESH  ||  v < V_THRESH )  

  8.     return NH * NS + vd;  

  9.     

  10.   /* otherwise determine "colorful" bin */  

  11.   hd = MIN( (int)(h * NH / H_MAX), NH-1 );  

  12.   sd = MIN( (int)(s * NS / S_MAX), NS-1 );  

  13.   return sd * NH + hd;  

  14. }  


img是目标矩形区域,h、s、v是个分量,hist就是直方图统计;NV、V_MAX....等是宏定义的固定值。

[cpp] view plain copy

  1. void normalize_histogram( histogram* histo )//直方图归一化  

  2. {  

  3.   float* hist;  

  4.   float sum = 0, inv_sum;  

  5.   int i, n;  

  6.   

  7.   hist = histo->histo;  

  8.   n = histo->n;  

  9.   

  10.   /* compute sum of all bins and multiply each bin by the sum's inverse */  

  11.   for( i = 0; i < n; i++ )  

  12.     sum += hist[i];  

  13.   inv_sum = 1.0 / sum;  

  14.   for( i = 0; i < n; i++ )  

  15.     hist[i] *= inv_sum;  

  16. }  

4、粒子初始化

根据选定的目标区域来初始化粒子,初始时所有粒子都为等权重,具有同样的属性。

[cpp] view plain copy

  1. /* create particles at the centers of each of n regions */  

  2. for( i = 0; i < n; i++ )  

  3. {  

  4.     width = regions[i].width;  

  5.     height = regions[i].height;  

  6.     x = regions[i].x + width / 2;<span style="white-space:pre">   </span>//粒子中心  

  7.     y = regions[i].y + height / 2;  

  8.     for( j = 0; j < np; j++ )  

  9.     {  

  10.         particles[k].x0 = particles[k].xp = particles[k].x = x;  

  11.         particles[k].y0 = particles[k].yp = particles[k].y = y;  

  12.         particles[k].sp = particles[k].s = 1.0;  

  13.         particles[k].width = width;  

  14.         particles[k].height = height;  

  15.         particles[k].histo = histos[i];  

  16.         particles[k++].w = 0;  

  17.     }  

  18. }  

步骤二、粒子相似度搜索、计算

5、粒子搜索

在步骤一中,初始化了100个粒子,由于初始帧中指定了目标区域,而该目标会在下一帧中发生偏移,但是相邻帧目标移动得不是太远,所以在目标区域附近随机撒出100个粒子。(此处使用的是二阶动态回归来估计偏移后的粒子位置)

[cpp] view plain copy

  1. particle transition( particle p, int w, int h, gsl_rng* rng )//随机撒出一个粒子  

  2. {  

  3.   float x, y, s;  

  4.   particle pn;  

  5.   //回归模型的参数即A1、A2、B0等Rob Hess在代码中已设定(我不知道是怎么来的?)  

  6.   /* sample new state using second-order autoregressive dynamics (使用二阶动态回归来自动更新粒子状态)*/  

  7.   x = A1 * ( p.x - p.x0 ) + A2 * ( p.xp - p.x0 ) +  

  8.     B0 * gsl_ran_gaussian( rng, TRANS_X_STD ) + p.x0;  

  9.   pn.x = MAX( 0.0, MIN( (float)w - 1.0, x ) );  

  10.   y = A1 * ( p.y - p.y0 ) + A2 * ( p.yp - p.y0 ) +  

  11.     B0 * gsl_ran_gaussian( rng, TRANS_Y_STD ) + p.y0;  

  12.   pn.y = MAX( 0.0, MIN( (float)h - 1.0, y ) );  

  13.   s = A1 * ( p.s - 1.0 ) + A2 * ( p.sp - 1.0 ) +  

  14.     B0 * gsl_ran_gaussian( rng, TRANS_S_STD ) + 1.0;  

  15.   pn.s = MAX( 0.1, s );  

  16.   pn.xp = p.x;  

  17.   pn.yp = p.y;  

  18.   pn.sp = p.s;  

  19.   pn.x0 = p.x0;  

  20.   pn.y0 = p.y0;  

  21.   pn.width = p.width;  

  22.   pn.height = p.height;  

  23.   pn.histo = p.histo;  

  24.   pn.w = 0;  

  25.   

  26.   return pn;  

  27. }  

6、相似度计算

然后计算这100个粒子hsv空间直方图与目标hsv空间直方图相似成度,马氏距离(Battacharyya)来度量两个粒子的相似度系数。


[cpp] view plain copy

  1. //计算粒子与跟踪区域直方图相似程度,越相似则权值越大  

  2. particles[j].w = likelihood( hsv_frame, cvRound(particles[j].y),  

  3.             cvRound( particles[j].x ),  

  4.             cvRound( particles[j].width * s ),  

  5.             cvRound( particles[j].height * s ),  

  6.             particles[j].histo );  



[cpp] view plain copy

  1. float histo_dist_sq( histogram* h1, histogram* h2 )//两个粒子的  

  2. {  

  3.   float* hist1, * hist2;  

  4.   float sum = 0;  

  5.   int i, n;  

  6.   

  7.   n = h1->n;  

  8.   hist1 = h1->histo;  

  9.   hist2 = h2->histo;  

  10.   

  11.   /* 

  12.     According the the Battacharyya similarity coefficient, 

  13.      

  14.     D = \sqrt{ 1 - \sum_1^n{ \sqrt{ h_1(i) * h_2(i) } } }//马氏距离公式 

  15.   */  

  16.   for( i = 0; i < n; i++ )  

  17.     sum += sqrt( hist1[i]*hist2[i] );  

  18.   return 1.0 - sum;  

  19. }  

步骤三、重采样


7、重采样

由步骤二知,经过一次搜索后,粒子的权重会发生改变,离目标距离远的粒子权重小,距离近的权重大。那么经过多帧的跟踪后,有的粒子权重会变得相当的小,也就是与目标不相似了,即粒子退化现象,这种粒子我们要抛弃,那么什么时候该抛弃它呢?就得设一个权重阈值,凡是权重低于阈值的粒子就抛弃。OK,原来有100个粒子,然后总会有被抛弃的粒子,似的粒子总数不满100个,此时就要找一些新的粒子来补充,那么就用最大权值来填充。(比如现在抛弃了20个粒子,我们就复制20个权值最大的粒子,放到里面填满100个。)——这就是重采样的过程。

[cpp] view plain copy

  1. particle* resample( particle* particles, int n )  

  2. {  

  3.   particle* new_particles;  

  4.   int i, j, np, k = 0;  

  5.   

  6.   qsort( particles, n, sizeof( particle ), &particle_cmp );//根据权重进行排序  

  7.   new_particles = malloc( n * sizeof( particle ) );  

  8.     for( i = 0; i < n; i++ )  

  9.     {  

  10.         np = cvRound( particles[i].w * n );//淘汰弱权值样本,保留阈值以上样本  

  11.         for( j = 0; j < np; j++ )  

  12.         {  

  13.             new_particles[k++] = particles[i];  

  14.             if( k == n )  

  15.             goto exit;  

  16.         }  

  17.     }  

  18.   while( k < n )  

  19.     new_particles[k++] = particles[0];//复制大权值样本以填充满样本空间  

  20.   

  21.  exit:  

  22.   return new_particles;  

  23. }  

8、更新粒子

将重采样后的100个粒子更新到步骤一4中。至此完成了一次粒子滤波。(需要注意的是,经过一次迭代后,步骤一4中的粒子权值已经不是等权值了)

9、小结

先初始化(完成1、2、3、4),再不停的迭代5、6、7、8过程(即:5—>6—>7—>8—>5......)。


以上是粒子滤波的全部过程,以及一些核心源码,可以认为权重最大的粒子是目标物体。值得关注的是,重采样会降低粒子的多样性(因为是许多粒子是直接复制过来的),这样也会对目标跟踪产生影响,有兴趣的可以继续研究改进算法以保证粒子的多样性。

over!

以上是对粒子滤波的个人见解,至于理论上看着比较复杂的理论公式推导没有体现出来(其实有些概率论知识我还是没看懂),再次膜拜一下Rob Hess大婶,牛的一逼(他也实现了Lowe的SIFT算法开源代码出来)。

本文所需材料都已打包上传,点击此处下载(配置gsl,请百度):http://download.csdn.net/detail/hujingshuang/8669945


欢迎转载,转载请注明出处


Copyright oneie ©2014-2017 All Rights Reserved. 所有资料来源于互联网对相关版权责任概不负责。如发现侵犯了您的版权请与我们联系。 网际学院 版权所有
免责声明  商务合作及投稿请联系 QQ:86662817