粒子濾波初探（二）利用粒子濾波實現視訊目標跟蹤-程式碼部分（C++&&opencv2.49）
阿新 • • 發佈：2018-12-24
利用粒子濾波實現視訊目標跟蹤工程實戰
放在最前：致謝taotao1233、yangyangcv、yang_xian521 以及先驅 Rob Hess 所開源的程式碼和思路。
本篇：基本為工程翻譯，以及對上面版本的一些修正，使用的是opencv2.49，以Mat型別代替了容易導致記憶體洩漏的Iplimage*型別，總的來說，就是從舊版的opencv資料結構翻譯到opencv2.49下。yang_xian521 也做過類似工作，但是不知道是否我操作有誤，執行其程式，處理速度很慢，而且跟蹤效果不好，既然如此，何不自己動手，豐衣足食？
程式的流程可以大致參考： 粒子濾波初探（一）利用粒子濾波實現視訊目標跟蹤的大致流程 ，會有一點點不同和簡化
- 我嘗試使用opencv2.49自帶的隨機數生成器RNG，效果不如gsl1.8的隨機數生成器，因此本工程依舊採用gsl1.8【因此需要配置gsl1.8，這裡不贅述】，如果哪位爺研究出來了，請下方留言。
- 本文不再贅述通過滑鼠選取目標的回撥部分，詳細請參見：ROI選取
直接獲取完整程式碼見這裡
既然是粒子濾波，首先我們構造粒子的結構體：
/*********************結構體************************/ // 粒子結構體 typedef struct particle { double x; // 當前x座標 double y; // 當前y座標 double scale; // 視窗比例係數 double xPre; // x座標預測位置 double yPre; // y座標預測位置 double scalePre; // 視窗預測比例係數 double xOri; // 原始x座標 double yOri; // 原始y座標 int width; // 原始區域寬度 int height; // 原始區域高度 //Rect rect; // 原始區域大小 MatND hist; // 粒子區域的特徵直方圖 double weight; // 該粒子的權重 } PARTICLE;
零、建立粒子陣列：
int num_particles = NUM_PARTICLE; //粒子數
PARTICLE particles[NUM_PARTICLE];
PARTICLE new_particles[NUM_PARTICLE];一、t-1時刻：視訊讀取【不展開】
二、等待滑鼠選取目標【不展開】
三、選取目標完成，則進行如下操作
【這裡對應部落格一的步驟一到步驟五，並作出了一些簡化（直接初始化粒子到目標區域，省去了一些步驟）】
（3.1）將選取的目標區域轉換到hsv空間下
（3.2）計算目標區域的hsv直方圖，只計算h和s兩個通道即可
（3.3）歸一化上面計算的直方圖到（0，1）
（3.4）初始化粒子
涉及函式如下：
/*************************粒子初始化******************************************/
void particle_init(particle* particles,int _num_particle,MatND hist)
{
for (int i = 0; i<_num_particle; i++)
{
//所有粒子初始化到框中的目標中心
particles[i].x = roiRect.x + 0.5 * roiRect.width;
particles[i].y = roiRect.y + 0.5 * roiRect.height;
particles[i].xPre = particles[i].x;
particles[i].yPre = particles[i].y;
particles[i].xOri = particles[i].x;
particles[i].yOri = particles[i].y;
//pParticles->rect = roiRect;
particles[i].width = roiRect.width;
particles[i].height = roiRect.height;
particles[i].scale = 1.0;
particles[i].scalePre = 1.0;
particles[i].hist = hist;
//權重全部為0？
particles[i].weight = 0;
}
}
/******************************************************************************/第三步程式碼如下：
if (getTargetFlag == true){
//目標區域轉換到hsv空間
cvtColor(roiImage, hsv_roiImage, COLOR_BGR2HSV);
//計算目標區域的直方圖
calcHist(&hsv_roiImage, 1, channels, Mat(), hist, 2, histSize, ranges);
normalize(hist, hist,0,1,NORM_MINMAX,-1,Mat()); // 歸一化L2
//粒子初始化
particle_init(particles, num_particles, hist);
}四、選取目標完成&&初始化完成，進入迴圈：【進入t時刻】
- 讀取新的一幀frame
4.1、粒子狀態生成（部落格一的第六、七步）&&重要性取樣（部落格一的第八步）
- 由上一時刻粒子群的狀態和權重，利用gsl庫生成高斯分佈隨機數，進而生成出新的粒子群的位置和狀態（函式transition()）
- 對每個粒子進行操作：
- （1）以粒子為中心，以粒子結構體中的.width和.height以及縮放因子scale，作為粒子所在區域的範圍，在frame中擷取該範圍，轉換到hsv空間下
- （2）計算上面擷取的hsv區域，計算其直方圖，並歸一化到（0，1）
- （3）這個新計算的直方圖與第三步初始化時計算的目標直方圖進行比較（cv::compareHist()）【計算BHATTACHARYYA距離】，將結果作為該粒子的權重
涉及函式如下：
/************************粒子狀態轉移（新位置生成預測）***********************/
//相關定義
/* standard deviations for gaussian sampling in transition model */
#define TRANS_X_STD 1.0
#define TRANS_Y_STD 0.5
#define TRANS_S_STD 0.001
/* autoregressive dynamics parameters for transition model */
#define A1 2.0//2.0
#define A2 -1.0//-1.0
#define B0 1.0000
particle transition(particle p, int w, int h, gsl_rng* rng)
{
//double rng_nu_x = rng.uniform(0., 1.);
//double rng_nu_y = rng.uniform(0., 0.5);
float x, y, s;
particle pn;
/* sample new state using second-order autoregressive dynamics */
x = A1 * (p.x - p.xOri) + A2 * (p.xPre - p.xOri) +
B0 * gsl_ran_gaussian(rng, TRANS_X_STD)/*rng.gaussian(TRANS_X_STD)*/ + p.xOri; //計算該粒子下一時刻的x
pn.x = MAX(0.0, MIN((float)w - 1.0, x));
y = A1 * (p.y - p.yOri) + A2 * (p.yPre - p.yOri) +
B0 * gsl_ran_gaussian(rng, TRANS_Y_STD)/*rng.gaussian(TRANS_Y_STD)*/ + p.yOri;
pn.y = MAX(0.0, MIN((float)h - 1.0, y));
s = A1 * (p.scale - 1.0) + A2 * (p.scalePre - 1.0) +
B0 * gsl_ran_gaussian(rng, TRANS_S_STD)/*rng.gaussian(TRANS_S_STD)*/ + 1.0;
pn.scale = MAX(0.1, s);
pn.xPre = p.x;
pn.yPre = p.y;
pn.scalePre = p.scale;
pn.xOri = p.xOri;
pn.yOri = p.yOri;
pn.width = p.width;
pn.height = p.height;
//pn.hist = p.hist;
pn.weight = 0;
return pn;
}對每個粒子的操作如下：
//對每個粒子的操作：
for (j = 0; j < num_particles; j++){
//rng = rng.next();
//這裡利用高斯分佈的隨機數來生成每個粒子下一次的位置以及範圍
particles[j] = transition(particles[j], frame.cols, frame.rows, rng);
s = particles[j].scale;
//根據新生成的粒子資訊擷取對應frame上的區域
Rect imgParticleRect = Rect(std::max(0, std::min(cvRound(particles[j].x - 0.5*particles[j].width), cvRound(frame.cols - particles[j].width*s))),
std::max(0, std::min(cvRound(particles[j].y - 0.5*particles[j].height), cvRound(frame.rows - particles[j].height*s))),
cvRound(particles[j].width*s),
cvRound(particles[j].height*s));
Mat imgParticle = current_frame(imgParticleRect).clone();
//上述區域轉換到hsv空間
cvtColor(imgParticle, imgParticle, CV_BGR2HSV);
//計算區域的直方圖
calcHist(&imgParticle, 1, channels, Mat(), particles[j].hist, 2, histSize, ranges);
//直方圖歸一化到（0，1）
normalize(particles[j].hist, particles[j].hist, 0, 1, NORM_MINMAX, -1, Mat()); // 歸一化L2
//畫出藍色的粒子框
rectangle(frame, imgParticleRect, Scalar(255, 0, 0), 1, 8);
imshow("particle", imgParticle);
//比較目標的直方圖和上述計算的區域直方圖,更新particle權重
particles[j].weight = exp(-100 * (compareHist(hist, particles[j].hist, CV_COMP_BHATTACHARYYA))); //CV_COMP_CORREL
int jj = 0;
}4.2 歸一化權重
/*************************粒子權重歸一化****************************/
void normalize_weights(particle* particles, int n)
{
float sum = 0;
int i;
for (i = 0; i < n; i++)
sum += particles[i].weight;
for (i = 0; i < n; i++)
particles[i].weight /= sum;
}
//歸一化權重
normalize_weights(particles, num_particles);4.3 重取樣（部落格一的第九步） 【原理參見：白巧克力亦唯心的總結：Particle Filter Tutorial 粒子濾波：從推導到應用（三）】
- （1）將粒子按權重從高到低排序
- （2）重取樣（函式resample()）
涉及函式如下：
/*************************粒子重取樣********************************/
void resample(particle* particles,particle* new_particles,int num_particles)
{
//計算每個粒子的概率累計和
double sum[NUM_PARTICLE], temp_sum = 0;
int k = 0;
for (int j = num_particles - 1; j >= 0; j--){
temp_sum += particles[j].weight;
sum[j] = temp_sum;
}
//為每個粒子生成一個均勻分佈【0，1】的隨機數
RNG sum_rng(time(NULL));
double Ran[NUM_PARTICLE];
for (int j = 0; j < num_particles; j++){
sum_rng = sum_rng.next();
Ran[j] = sum_rng.uniform(0., 1.);
}
//在粒子概率累積和陣列中找到最小的大於給定隨機數的索引，複製該索引的粒子一次到新的粒子陣列中 【從權重高的粒子開始】
for (int j = 0; j <num_particles; j++){
int copy_index = get_min_index(sum, num_particles, Ran[j]);
new_particles[k++] = particles[copy_index];
if (k == num_particles)
break;
}
//如果上面的操作完成，新粒子陣列的數量仍少於原給定粒子數量，則複製權重最高的粒子，直到粒子數相等
while (k < num_particles)
{
new_particles[k++] = particles[0]; //複製權值最高的粒子
}
//以新粒子陣列覆蓋久的粒子陣列
for (int i = 0; i<num_particles; i++)
{
particles[i] = new_particles[i]; //複製新粒子到particles
}
}操作步驟如下：
int np, k = 0;
//將粒子按權重從高到低排序
qsort(particles, num_particles, sizeof(particle), &particle_cmp);
//重取樣
resample(particles, new_particles, num_particles);4.4 求期望，獲取跟蹤位置 【這裡筆者偷了個懶，直接以權重最高的粒子作為目標了，正統做法應當是按粒子權重求期望】
- （1）將粒子按權重從高到低排序
- （2）提取目標位置
//排序
qsort(particles, num_particles, sizeof(particle), &particle_cmp);
//這裡直接取權重最高的作為目標了，標準做法應該是按加權平均來計算目標位置
s = particles[0].scale;
Rect rectTrackingTemp = Rect(std::max(0, std::min(cvRound(particles[0].x - 0.5*particles[0].width), cvRound(frame.cols - particles[0].width*s))),
std::max(0, std::min(cvRound(particles[0].y - 0.5*particles[0].height), cvRound(frame.rows - particles[0].height*s))),
cvRound(particles[0].width*s),
cvRound(particles[0].height*s));
rectangle(frame, rectTrackingTemp, Scalar(0, 0, 255), 1, 8, 0);完整程式碼：
依個人來講，初學者大可不必看上面的程式碼，直接執行一次程式碼再說。在看程式碼的過程中有不明白的，再慢慢翻看。
這裡就直接給出可執行的程式碼了，基於【vs2013+opencv2.49+gsl1.8】：
如果你覺得值：可以從csdn下載，給小弟一個賺積分的機會
（1）csdn渠道：https://download.csdn.net/download/dieju8330/10858046
如果你沒有csdn積分，沒關係，
（2）百度wanpan：連結：https://pan.baidu.com/s/1VR_UYAqKk0fBTwZB_8ehPw 提取碼：6xvc
效果圖如下：
