作者：gnuhpc
出處：http://www.cnblogs.com/gnuhpc/
1.引言
這個項目是由俄亥俄州立大學（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.
大致翻譯如下：
物體追蹤是一個棘手的問題。一個普適的，通用的物體追蹤算法必須應對諸如攝像頭運動、不穩定物體的追蹤、復雜的背景、以及存在其他移動物體等困難的狀況。粒子濾波算法是一個采用蒙特卡羅采樣進行貝葉斯濾波的方法。這種方法有許多的用途，但它已經成為進行物體追蹤最好的方法。從概念上講，一個粒子濾波算法包含一個被監視系統的狀態的概率分布。在本項目中，狀態就是指被追蹤物體的位置、大小等等。在許多情況下，非線性和非高斯型在物體的運動和相似性建模上會得到一個難以處理的濾波分布。粒子濾波采用將這個分布重新表示為一組加權值，或稱為粒子的方法克服了這個困難。每個粒子表示一個可能的系統狀態實例。換句話說，每個粒子描述了被追蹤物體可能處於的一個方位。一個粒子集包含了被追蹤物體最有可能處於的方位。因此，我們可以通過尋找在粒子濾波分布中最大的權重來確定物體最有可能處於的狀態。
2.程序流程總述
1.命令行參數處理 －>
2.設置隨機數生成器環境，創建隨機數生成器、並且對其初始化。->
3.初始化視頻句柄 －>
4.取視頻中的一幀進行處理 －>
1）GRB－>HSV
2）保存當前幀在frames
3) 判斷是否為第一幀，
若是則，
（1）忙等用戶選定欲跟蹤的區域
（2）計算相關區域直方圖
（3）得到跟蹤粒子
若不是則，
（1）對每個粒子作變換，並計算每個粒子的權重
（2）對粒子集合進行歸一化
（3）重新采樣粒子
4）畫出粒子所代表的區域
5.釋放圖像
3.命令行參數處理
void arg_parse( int argc, char** argv )
{
int i = 0;
/*extract program name from command line (remove path, if present) */
pname = remove_path( argv[0] );
/*parse commandline options */
while( TRUE )
{
char* arg_check;
int arg = getopt( argc, argv, OPTIONS );
if( arg == -1 )
break;
switch( arg )
{
/* user asked for help */
case 'h':
usage( pname );
exit(0);
break;
case 'a':
show_all = TRUE;
break;
/* user wants to output tracking sequence */
case 'o':
export = TRUE;
break;
/* user wants to set number of particles */
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;
/* catch invalid arguments */
default:
fatal_error( "-%c: invalid option/nTry '%s -h' for help.",
optopt, pname );
}
}
/* make sure input and output files are specified */
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 );
/* record video file name */
vid_file = argv[optind];
}
作者使用Getopt這個系統函數對命令行進行解析，－h表示顯示幫助，－a表示將所有粒子所代表的位置都顯示出來，－o表示輸出tracking的幀，－p number進行粒子數的設定，然后再最后指定要處理的視頻文件。
4.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空間。
5.設定隨機數
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部分所述進行隨機數生成有三個步驟：
隨機數生成器環境建立，隨機數生成器的創建，隨機數生成器的初始化。
6.計算選定區域直方圖
/*
Computes a reference histogram for each of the object regions defined by
the user
@param frame video frame in which to compute histograms; should have been
converted to hsv using bgr2hsv in observation.h
@param regions regions of /a frame over which histograms should be computed
@param n number of regions in /a regions
@param export if TRUE, object region images are exported
@return Returns an /a n element array of normalized histograms corresponding
to regions of /a frame specified in /a regions.
*/
histogram** compute_ref_histos( IplImage* frame, CvRect* regions, int n )
{
histogram** histos = malloc( n * sizeof( histogram* ) );
IplImage* tmp;
int i;
/* extract each region from frame and compute its histogram */
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對直方圖進行歸一化。
計算直方圖的函數詳解如下：
/*
Calculates a cumulative histogram as defined above for a given array
of images
@param img an array of images over which to compute a cumulative histogram;
each must have been converted to HSV colorspace using bgr2hsv()
@param n the number of images in imgs
@return Returns an un-normalized HSV histogram for /a imgs
*/
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++ )
{
/* extract individual HSV planes from image */
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 );
/* increment appropriate histogram bin for each pixel */
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
7.初始化粒子集
/*
Creates an initial distribution of particles at specified locations
@param regions an array of regions describing player locations around
which particles are to be sampled
@param histos array of histograms describing regions in /a regions
@param n the number of regions in /a regions
@param p the total number of particles to be assigned
@return Returns an array of /a p particles sampled from around regions in
/a regions
*/
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;
/* create particles at the centers of each of n regions */
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;
}
}
/* make sure to create exactly p particles */
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個粒子被初始化。
8.位置可能性確定
/*
Computes the likelihood of there being a player at a given location in
an image
@param img image that has been converted to HSV colorspace using bgr2hsv()
@param r row location of center of window around which to compute likelihood
@param c col location of center of window around which to compute likelihood
@param w width of region over which to compute likelihood
@param h height of region over which to compute likelihood
@param ref_histo reference histogram for a player; must have been
normalized with normalize_histogram()
@return Returns the likelihood of there being a player at location
(/a r, /a c) in /a img
*/
float likelihood( IplImage* img, int r, int c,
int w, int h, histogram* ref_histo )
{
IplImage* tmp;
histogram* histo;
float d_sq;
/* extract region around (r,c) and compute and normalize its histogram */
cvSetImageROI( img, cvRect( c - w / 2, r - h / 2, w, h ) );
tmp = cvCreateImage( cvGetSize(img), IPL_DEPTH_32F, 3 );
cvCopy( img, tmp, NULL );
cvResetImageROI( img );
histo = calc_histogram( &tmp, 1 );
cvReleaseImage( &tmp );
normalize_histogram( histo );
/* compute likelihood as e^{/lambda D^2(h, h^*)} */
d_sq = histo_dist_sq( histo, ref_histo );
free( histo );
return exp( -LAMBDA * d_sq );
}
程序首先取出對相關粒子表示的區域，然后計算其直方圖，並且歸一化。將這個直方圖和原來用戶選定區域的直方圖傳入函數histo_dist_sq進行比較，最后返回e^(-Lambda*d_sq)返回，成為這個粒子的權重。
函數histo_dist_sq的實現如下：
/*
Computes squared distance metric based on the Battacharyya similarity
coefficient between histograms.
@param h1 first histogram; should be normalized
@param h2 second histogram; should be normalized
@return Returns a squared distance based on the Battacharyya similarity
coefficient between /a h1 and /a h2
*/
float histo_dist_sq( histogram* h1, histogram* h2 )
{
float* hist1, * hist2;
float sum = 0;
int i, n;
n = h1->n;
hist1 = h1->histo;
hist2 = h2->histo;
/*
According the the Battacharyya similarity coefficient,
D = /sqrt{ 1 - /sum_1^n{ /sqrt{ h_1(i) * h_2(i) } } }
*/
for( i = 0; i < n; i++ )
sum += sqrt( hist1[i]*hist2[i] );
return 1.0 - sum;
}
采用統計學上的巴氏距離Bhattacharyya distance，根據wiki的描述， Bhattacharyya distance 描述的是兩個離散概率分布的相似性，它通常在分類操作中被用來度量不同類型的可分離性，也就是說這個距離算式就是評定相似度的。嚴格定義為：
For discrete probability distributions p and q over the same domain X, it is defined as:
where:
is the Bhattacharyya coefficient .
該程序中的算式和這個式子略有差別。
9.粒子集合變換
/*
Samples a transition model for a given particle
@param p a particle to be transitioned
@param w video frame width
@param h video frame height
@param rng a random number generator from which to sample
@return Returns a new particle sampled based on <EM>p</EM>'s transition
model
*/
particle transition( particle p, int w, int h, gsl_rng* rng )
{
float x, y, s;
particle pn;
/* sample new state using second-order autoregressive dynamics */
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之內的條件。
10.粒子集重新采樣
/*
Re-samples a set of particles according to their weights to produce a
new set of unweighted particles
@param particles an old set of weighted particles whose weights have been
normalized with normalize_weights()
@param n the number of particles in /a particles
@return Returns a new set of unweighted particles sampled from /a particles
*/
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排序函數，按照權重，由大到小將原粒子集排序。然后將權重大的在新的粒子集中分配的多一點。
11.權重歸一化
/*
Normalizes particle weights so they sum to 1
@param particles an array of particles whose weights are to be normalized
@param n the number of particles in /a particles
*/
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;
}