这个项目是由俄亥俄州立大学（OSU）一位博士生所写，http://web.engr.oregonstate.edu/~hess/，这位博士在其个人主页上对该项目进行了如下描述：
Object tracking is a tricky problem. A general, all-purpose object tracking algorithm must deal with difficulties like camera motion, erratic object motion, cluttered backgrounds, and other moving objects. Such hurdles render general image processing techniques an inadequate solution to the object tracking problem.
Particle filtering is a Monte Carlo sampling approach to Bayesian filtering. It has many uses but has become the state of the art in object tracking. Conceptually, a particle filtering algorithm maintains a probability distribution over the state of the system it is monitoring, in this case, the state -- location, scale, etc. -- of the object being tracked. In most cases, non-linearity and non-Gaussianity in the object's motion and likelihood models yields an intractable filtering distribution. Particle filtering overcomes this intractability by representing the distribution as a set of weighted samples, or particles. Each particle represents a possible instantiation of the state of the system. In other words, each particle describes one possible location of the object being tracked. The set of particles contains more weight at locations where the object being tracked is more likely to be. We can thus determine the most probable state of the object by finding the location in the particle filtering distribution with the highest weight.
大致翻译如下：
物体追踪是一个棘手的问题。一个普适的，通用的物体追踪算法必须应对诸如摄像头运动、不稳定物体的追踪、复杂的背景、以及存在其他移动物体等困难的状况。粒子滤波算法是一个采用蒙特卡罗采样进行贝叶斯滤波的方法。这种方法有许多的用途，但它已经成为进行物体追踪最好的方法。从概念上讲，一个粒子滤波算法包含一个被监视系统的状态的概率分布。在本项目中，状态就是指被追踪物体的位置、大小等等。在许多情况下，非线性和非高斯型在物体的运动和相似性建模上会得到一个难以处理的滤波分布。粒子滤波采用将这个分布重新表示为一组加权值，或称为粒子的方法克服了这个困难。每个粒子表示一个可能的系统状态实例。换句话说，每个粒子描述了被追踪物体可能处于的一个方位。一个粒子集包含了被追踪物体最有可能处于的方位。因此，我们可以通过寻找在粒子滤波分布中最大的权重来确定物体最有可能处于的状态。
程序流程：
1.命令行参数处理 －>
2.设置随机数生成器环境，创建随机数生成器、并且对其初始化。->
3.初始化视频句柄 －>
4.取视频中的一帧进行处理 －>
1）GRB－>HSV
2）保存当前帧在frames
3) 判断是否为第一帧，
若是则，
（1）忙等用户选定欲跟踪的区域
（2）计算相关区域直方图
（3）得到跟踪粒子
若不是则，
（1）对每个粒子作变换，并计算每个粒子的权重
（2）对粒子集合进行归一化
（3）重新采样粒子
4）画出粒子所代表的区域
5.释放图像
OpenCV学习——物体跟踪的粒子滤波算法实现之命令行参数处理
void arg_parse( int argc, char** argv )
{
int i = 0;
pname = remove_path( argv[0] );
while( TRUE )
{
char* arg_check;
int arg = getopt( argc, argv, OPTIONS );
if( arg == -1 )
break;
switch( arg )
{
case 'h':
usage( pname );
exit(0);
break;
case 'a':
show_all = TRUE;
break;
case 'o':
export = TRUE;
break;
case 'p':
if( ! optarg )
fatal_error( "error parsing arguments at -%c\n" \
"Try '%s -h' for help.", arg, pname );
num_particles = strtol( optarg, &arg_check, 10 );
if( arg_check == optarg || *arg_check != '\0' )
fatal_error( "-%c option requires an integer argument\n" \
"Try '%s -h' for help.", arg, pname );
break;
default:
fatal_error( "-%c: invalid option\nTry '%s -h' for help.",
optopt, pname );
}
}
if( argc - optind < 1 )
fatal_error( "no input image specified.\nTry '%s -h' for help.", pname );
if( argc - optind > 2 )
fatal_error( "too many arguments.\nTry '%s -h' for help.", pname );
vid_file = argv[optind];
}
作者使用Getopt这个系统函数对命令行进行解析，－h表示显示帮助，－a表示将所有粒子所代表的位置都显示出来，－o表示输出tracking的帧，－p number进行粒子数的设定，然后再最后指定要处理的视频文件。
OpenCV学习——物体跟踪的粒子滤波算法实现之RGB->HSV
IplImage* bgr2hsv( IplImage* bgr )
{
IplImage* bgr32f, * hsv;
bgr32f = cvCreateImage( cvGetSize(bgr), IPL_DEPTH_32F, 3 );
hsv = cvCreateImage( cvGetSize(bgr), IPL_DEPTH_32F, 3 );
cvConvertScale( bgr, bgr32f, 1.0 / 255.0, 0 );
cvCvtColor( bgr32f, hsv, CV_BGR2HSV );
cvReleaseImage( &bgr32f );
return hsv;
}
程序现将图像的像素值归一化，然后使用OpenCV中的cvCvtcolor函数将图像从RGB空间转换到HSV空间。
OpenCV学习——物体跟踪的粒子滤波算法实现之设定随机数
gsl_rng_env_setup();//setup the enviorment of random number generator
rng = gsl_rng_alloc( gsl_rng_mt19937 );//create a random number generator
gsl_rng_set( rng, time(NULL) );//initializes the random number generator.
作者使用GSL库进行随机数的产生，GSL是GNU的科学计算库，其中手册中random部分所述进行随机数生成有三个步骤：
随机数生成器环境建立，随机数生成器的创建，随机数生成器的初始化。
OpenCV学习——物体跟踪的粒子滤波算法实现之计算选定区域直方图
histogram** compute_ref_histos( IplImage* frame, CvRect* regions, int n )
{
histogram** histos = malloc( n * sizeof( histogram* ) );
IplImage* tmp;
int i;
for( i = 0; i < n; i++ )
{
cvSetImageROI( frame, regions[i] );//set the region of interest
tmp = cvCreateImage( cvGetSize( frame ), IPL_DEPTH_32F, 3 );
cvCopy( frame, tmp, NULL );
cvResetImageROI( frame );//free the ROI
histos[i] = calc_histogram( &tmp, 1 );//calculate the hisrogram
normalize_histogram( histos[i] );//Normalizes a histogram so all bins sum to 1.0
cvReleaseImage( &tmp );
}
return histos;
}
程序中先设置了一个类型为histogram的指向指针的指针，是histogram指针数组的指针，这个数组是多个选定区域的直方图数据存放的位置。然后对于每一个用户指定的区域，在第一帧中都进行了ROI区域设置，通过对ROI区域的设置取出选定区域，交给函数calc_histogram计算出直方图，并使用normalize_histogram对直方图进行归一化。
计算直方图的函数详解如下：
histogram* calc_histogram( IplImage** imgs, int n )
{
IplImage* img;
histogram* histo;
IplImage* h, * s, * v;
float* hist;
int i, r, c, bin;
histo = malloc( sizeof(histogram) );
histo->n = NH*NS + NV;
hist = histo->histo;
memset( hist, 0, histo->n * sizeof(float) );
for( i = 0; i < n; i++ )
{
img = imgs[i];
h = cvCreateImage( cvGetSize(img), IPL_DEPTH_32F, 1 );
s = cvCreateImage( cvGetSize(img), IPL_DEPTH_32F, 1 );
v = cvCreateImage( cvGetSize(img), IPL_DEPTH_32F, 1 );
cvCvtPixToPlane( img, h, s, v, NULL );
for( r = 0; r < img->height; r++ )
for( c = 0; c < img->width; c++ )
{
bin = histo_bin( pixval32f( h, r, c ),
pixval32f( s, r, c ),
pixval32f( v, r, c ) );
hist[bin] += 1;
}
cvReleaseImage( &h );
cvReleaseImage( &s );
cvReleaseImage( &v );
}
return histo;
}
这个函数将h、s、 v分别取出，然后以从上到下，从左到右的方式遍历以函数histo_bin的评判规则放入相应的bin中（很形象的）。函数histo_bin的评判规则详见下图：
｜－－－－｜－－－－｜－－－－｜。。。。｜－－－－｜－－－－－－｜－－－－－－｜。。。。｜－－－－－－－｜
1NH 2NH 3NH NS＊NH NS＊NH＋1 NS*NH+2 NS＊NH＋NV
OpenCV学习——物体跟踪的粒子滤波算法实现之初始化粒子集
particle* init_distribution( CvRect* regions, histogram** histos, int n, int p)
{
particle* particles;
int np;
float x, y;
int i, j, width, height, k = 0;
particles = malloc( p * sizeof( particle ) );
np = p / n;
for( i = 0; i < n; i++ )
{
width = regions[i].width;
height = regions[i].height;
x = regions[i].x + width / 2;
y = regions[i].y + height / 2;
for( j = 0; j < np; j++ )
{
particles[k].x0 = particles[k].xp = particles[k].x = x;
particles[k].y0 = particles[k].yp = particles[k].y = y;
particles[k].sp = particles[k].s = 1.0;
particles[k].width = width;
particles[k].height = height;
particles[k].histo = histos[i];
particles[k++].w = 0;
}
}
i = 0;
while( k < p )
{
width = regions[i].width;
height = regions[i].height;
x = regions[i].x + width / 2;
y = regions[i].y + height / 2;
particles[k].x0 = particles[k].xp = particles[k].x = x;
particles[k].y0 = particles[k].yp = particles[k].y = y;
particles[k].sp = particles[k].s = 1.0;
particles[k].width = width;
particles[k].height = height;
particles[k].histo = histos[i];
particles[k++].w = 0;
i = ( i + 1 ) % n;
}
return particles;
}
程序中的变量np是指若有多个区域n，则一个区域内的粒子数为p/n，这样粒子的总数为p。然后程序对每个区域（n个）中p/n个粒子进行初始化，三个位置坐标都为选定区域的中点，比例都为1，宽度和高度为选定区域的高度。然后又跑了个循环确定p个粒子被初始化。
OpenCV学习——物体跟踪的粒子滤波算法实现之粒子集合变换
particle transition( particle p, int w, int h, gsl_rng* rng )
{
float x, y, s;
particle pn;
x = A1 * ( p.x - p.x0 ) + A2 * ( p.xp - p.x0 ) +
B0 * gsl_ran_gaussian( rng, TRANS_X_STD ) + p.x0;
pn.x = MAX( 0.0, MIN( (float)w - 1.0, x ) );
y = A1 * ( p.y - p.y0 ) + A2 * ( p.yp - p.y0 ) +
B0 * gsl_ran_gaussian( rng, TRANS_Y_STD ) + p.y0;
pn.y = MAX( 0.0, MIN( (float)h - 1.0, y ) );
s = A1 * ( p.s - 1.0 ) + A2 * ( p.sp - 1.0 ) +
B0 * gsl_ran_gaussian( rng, TRANS_S_STD ) + 1.0;
pn.s = MAX( 0.1, s );
pn.xp = p.x;
pn.yp = p.y;
pn.sp = p.s;
pn.x0 = p.x0;
pn.y0 = p.y0;
pn.width = p.width;
pn.height = p.height;
pn.histo = p.histo;
pn.w = 0;
return pn;
}
程序使用动态二阶自回归模型作为基本变换思路，变换的对象有坐标x，坐标y，和比例s。变换的x和y要符合在width和height之内的条件。
OpenCV学习——物体跟踪的粒子滤波算法实现之粒子集重新采样
particle* resample( particle* particles, int n )
{
particle* new_particles;
int i, j, np, k = 0;
qsort( particles, n, sizeof( particle ), &particle_cmp );
new_particles = malloc( n * sizeof( particle ) );
for( i = 0; i < n; i++ )
{
np = cvRound( particles[i].w * n );
for( j = 0; j < np; j++ )
{
new_particles[k++] = particles[i];
if( k == n )
goto exit;
}
}
while( k < n )
new_particles[k++] = particles[0];
exit:
return new_particles;
}
程序先使用C标准库中的qsort排序函数，按照权重，由大到小将原粒子集排序。然后将权重大的在新的粒子集中分配的多一点。
OpenCV学习——物体跟踪的粒子滤波算法实现之权重归一化
void normalize_weights( particle* particles, int n )
{
float sum = 0;
int i;
for( i = 0; i < n; i++ )
sum += particles[i].w;
for( i = 0; i < n; i++ )
particles[i].w /= sum;
}
转自：http://blog.csdn.net/gnuhpc/category/549384.aspx?PageNumber=3
OpenCV中实现了粒子滤波的代码，位置在c:\program files\opencv\cv\src\cvcondens.cpp文件，通过分析这个文件，可以知道库函数中如何实现粒子滤波过程的。
首先是从手册上拷贝的粒子滤波跟踪器的数据结构：
typedef struct CvConDensation
{
int MP; // 测量向量的维数： Dimension of measurement vector
int DP; // 状态向量的维数： Dimension of state vector
float* DynamMatr; // 线性动态系统矩阵：Matrix of the linear Dynamics system
float* State; // 状态向量： Vector of State
int SamplesNum; // 粒子数： Number of the Samples
float** flSamples; // 粒子向量数组： array of the Sample Vectors
float** flNewSamples; // 粒子向量临时数组： temporary array of the Sample Vectors
float* flConfidence; // 每个粒子的置信度(译者注：也就是粒子的权值)：Confidence for each Sample
float* flCumulative; // 权值的累计： Cumulative confidence
float* Temp; // 临时向量：Temporary vector
float* RandomSample; // 用来更新粒子集的随机向量： RandomVector to update sample set
CvRandState* RandS; // 产生随机向量的结构数组： Array of structures to generate random vectors
} CvConDensation;
与粒子滤波相关的几个函数：
cvCreateConDensation：用于构造上述滤波器数据结构
cvReleaseConDensation：释放滤波器
cvConDensInitSampleSet：初始化粒子集
cvConDensUpdateByTime：更新粒子集
下面着重对这几个函数进行分析。
CV_IMPL CvConDensation* cvCreateConDensation( int DP, int MP, int SamplesNum )
{
int i;
CvConDensation *CD = 0;
CV_FUNCNAME( "cvCreateConDensation" );
__BEGIN__;
if( DP < 0 || MP < 0 || SamplesNum < 0 )
CV_ERROR( CV_StsOutOfRange, "" );
/* allocating memory for the structure */
CV_CALL( CD = (CvConDensation *) cvAlloc( sizeof( CvConDensation )));
/* setting structure params */
CD->SamplesNum = SamplesNum;
CD->DP = DP;
CD->MP = MP;
/* allocating memory for structure fields */
CV_CALL( CD->flSamples = (float **) cvAlloc( sizeof( float * ) * SamplesNum ));
CV_CALL( CD->flNewSamples = (float **) cvAlloc( sizeof( float * ) * SamplesNum ));
CV_CALL( CD->flSamples[0] = (float *) cvAlloc( sizeof( float ) * SamplesNum * DP ));
CV_CALL( CD->flNewSamples[0] = (float *) cvAlloc( sizeof( float ) * SamplesNum * DP ));
/* setting pointers in pointer's arrays */
for( i = 1; i < SamplesNum; i++ )
{
CD->flSamples[i] = CD->flSamples[i - 1] + DP;
CD->flNewSamples[i] = CD->flNewSamples[i - 1] + DP;
}
CV_CALL( CD->State = (float *) cvAlloc( sizeof( float ) * DP ));
CV_CALL( CD->DynamMatr = (float *) cvAlloc( sizeof( float ) * DP * DP ));
CV_CALL( CD->flConfidence = (float *) cvAlloc( sizeof( float ) * SamplesNum ));
CV_CALL( CD->flCumulative = (float *) cvAlloc( sizeof( float ) * SamplesNum ));
CV_CALL( CD->RandS = (CvRandState *) cvAlloc( sizeof( CvRandState ) * DP ));
CV_CALL( CD->Temp = (float *) cvAlloc( sizeof( float ) * DP ));
CV_CALL( CD->RandomSample = (float *) cvAlloc( sizeof( float ) * DP ));
/* Returning created structure */
__END__;
return CD;
}
输入参数分别是系统状态向量维数、测量向量维数以及粒子个数，然后根据这些参数为滤波器相应结构分配空间。
CV_IMPL void
cvReleaseConDensation( CvConDensation ** ConDensation )
{
CV_FUNCNAME( "cvReleaseConDensation" );
__BEGIN__;
CvConDensation *CD = *ConDensation;
if( !ConDensation )
CV_ERROR( CV_StsNullPtr, "" );
if( !CD )
EXIT;
/* freeing the memory */
cvFree( &CD->State );
cvFree( &CD->DynamMatr);
cvFree( &CD->flConfidence );
cvFree( &CD->flCumulative );
cvFree( &CD->flSamples[0] );
cvFree( &CD->flNewSamples[0] );
cvFree( &CD->flSamples );
cvFree( &CD->flNewSamples );
cvFree( &CD->Temp );
cvFree( &CD->RandS );
cvFree( &CD->RandomSample );
/* release structure */
cvFree( ConDensation );
__END__;
}
释放滤波器占用的空间，没什么好说的
CV_IMPL void
cvConDensInitSampleSet( CvConDensation * conDens, CvMat * lowerBound, CvMat * upperBound )
{
int i, j;
float *LBound;
float *UBound;
float Prob = 1.f / conDens->SamplesNum;
CV_FUNCNAME( "cvConDensInitSampleSet" );
__BEGIN__;
if( !conDens || !lowerBound || !upperBound )
CV_ERROR( CV_StsNullPtr, "" );
if( CV_MAT_TYPE(lowerBound->type) != CV_32FC1 ||
!CV_ARE_TYPES_EQ(lowerBound,upperBound) )
CV_ERROR( CV_StsBadArg, "source has not appropriate format" );
if( (lowerBound->cols != 1) || (upperBound->cols != 1) )
CV_ERROR( CV_StsBadArg, "source has not appropriate size" );
if( (lowerBound->rows != conDens->DP) || (upperBound->rows != conDens->DP) )
CV_ERROR( CV_StsBadArg, "source has not appropriate size" );
LBound = lowerBound->data.fl;
UBound = upperBound->data.fl;
/* Initializing the structures to create initial Sample set */
这里根据输入的动态范围给每个系统状态分配一个产生随即数的结构
for( i = 0; i < conDens->DP; i++ )
{
cvRandInit( &(conDens->RandS[i]),
LBound[i],
UBound[i],
i );
}
/* Generating the samples */
根据产生的随即数，为每个粒子的每个系统状态分配初始值，并将每个粒子的置信度设置为相同的1/n
for( j = 0; j < conDens->SamplesNum; j++ )
{
for( i = 0; i < conDens->DP; i++ )
{
cvbRand( conDens->RandS + i, conDens->flSamples[j] + i, 1 );
}
conDens->flConfidence[j] = Prob;
}
/* Reinitializes the structures to update samples randomly */
产生以后更新粒子系统状态的随即结构，采样范围为原来初始范围的-1/5到1/5
for( i = 0; i < conDens->DP; i++ )
{
cvRandInit( &(conDens->RandS[i]),
(LBound[i] - UBound[i]) / 5,
(UBound[i] - LBound[i]) / 5,
i);
}
__END__;
}
CV_IMPL void
cvConDensUpdateByTime( CvConDensation * ConDens )
{
int i, j;
float Sum = 0;
CV_FUNCNAME( "cvConDensUpdateByTime" );
__BEGIN__;
if( !ConDens )
CV_ERROR( CV_StsNullPtr, "" );
/* Sets Temp to Zero */
icvSetZero_32f( ConDens->Temp, ConDens->DP, 1 );
/* Calculating the Mean */
icvScaleVector_32f是内联函数，表示flConfidence乘flSamples，放到State里，即求系统状态的过程，系统状态保存在temp中
for( i = 0; i < ConDens->SamplesNum; i++ )
{
icvScaleVector_32f( ConDens->flSamples[i], ConDens->State, ConDens->DP,
ConDens->flConfidence[i] );
icvAddVector_32f( ConDens->Temp, ConDens->State, ConDens->Temp, ConDens->DP );
Sum += ConDens->flConfidence[i];
ConDens->flCumulative[i] = Sum;
}
/* Taking the new vector from transformation of mean by dynamics matrix */
icvScaleVector_32f( ConDens->Temp, ConDens->Temp, ConDens->DP, 1.f / Sum );
icvTransformVector_32f( ConDens->DynamMatr, ConDens->Temp, ConDens->State, ConDens->DP,
ConDens->DP );
Sum = Sum / ConDens->SamplesNum;
/* Updating the set of random samples */
重采样，将新采样出来的粒子保存在flNewSamples中
for( i = 0; i < ConDens->SamplesNum; i++ )
{
j = 0;
while( (ConDens->flCumulative[j] <= (float) i * Sum)&&(j<ConDens->SamplesNum-1))
{
j++;
}
icvCopyVector_32f( ConDens->flSamples[j], ConDens->DP, ConDens->flNewSamples[i] );
}
/* Adding the random-generated vector to every vector in sample set */
for( i = 0; i < ConDens->SamplesNum; i++ )
{
为每个新粒子产生一个随即量
for( j = 0; j < ConDens->DP; j++ )
{
cvbRand( ConDens->RandS + j, ConDens->RandomSample + j, 1 );
}
将flNewSamples加一个随即量，保存入flSamples中
icvTransformVector_32f( ConDens->DynamMatr, ConDens->flNewSamples[i],
ConDens->flSamples[i], ConDens->DP, ConDens->DP );
icvAddVector_32f( ConDens->flSamples[i], ConDens->RandomSample, ConDens->flSamples[i],
ConDens->DP );
}
__END__;
}
运行完上述函数之后，经重采样后的粒子保存在flSamples中，同时上次迭代的加权系统状态保存在State中