( 注：本文為這學期一個作業，關於粒子濾波的介紹在前面2篇博客中已提到過，即：
目標跟蹤學習筆記_2(particle filter初探1)
目標跟蹤學習筆記_3(particle filter初探2)
前面2篇博客已經提到當粒子數增加時會內存報錯，后面又仔細查了下程序，是代碼方面的問題。所以本次的代碼與前幾次改變比較小。當然這些code基本也是參考網上的。代碼寫得很不規范，時間不夠，等以后有機會將其優化並整理成類的形式。)
Opencv實現粒子濾波算法
摘要
本文通過opencv實現了一種目標跟蹤算法——粒子濾波算法，算法的思想來源於文獻[1][2],且在其思想上稍微做了些修改。其大概過程是：首先手動用鼠標框出一個目標區域，計算其直方圖特征值作為模板，然后在該目標中心周圍撒粒子，根據所撒粒子為中心的矩形框內計算其直方圖特征，並與目標相比較，最后根據比較出的結果重復上面過程，即重采樣的方法撒粒子，粒子擴散，狀態觀察，目標預測。最后通過實驗證明，取得了較好的效果。
關鍵字：目標跟蹤，粒子濾波，opencv
前言
目標跟蹤過程分為2部分，即目標特征提取和目標跟蹤算法。
其中目標特征提取又包括以下幾種：1. 各種色彩空間直方圖，利用色彩空間的直方圖分布作為目標跟蹤的特征，可以減少物體遠近距離的影響，因為其顏色分布大致相同。2.輪廓特征，提取目標的輪廓特征，可以加快算法的速度，且可以在目標有小部分影響的情況下同樣有效果。3. 紋理特征，如果被跟蹤目標是有紋理的，則根據其紋理特征來跟蹤效果會有所改善。
目標跟蹤算法目前大概分為以下4種：1. 基於meanshift算法，即利用meanshift算法可以快速找到領域目標最相似的地方，效果還不錯，但是其只能找到局部最大值，且不能解決遮擋問題以及不能自適應跟蹤目標的形狀，方向等。其后面有學者對其做了改進，比如說camshift，就可以自適應物體的大小，方向，具有較好的跟蹤效果。2. Kalman濾波的思想，該思想是利用物體的運動模型來，即服從高斯模型，來對目標狀態進行預測，然后與觀察模型進行比較，根據2者之間的誤差來尋找運動目標的狀態，但是該算法的精度不高，因為其高斯運動模型在現實生活中很多條件下並得不到滿足，並且該算法對雜亂的背景也很敏感。3. 基於粒子濾波的思想，每次通過實驗可以重采樣粒子的分布，根據該分布對粒子進行擴散，然后通過擴散的結果來觀察目標的狀態，最后更新目標的狀態。該算法最大的特點是跟蹤速度快，且能解決一部分遮擋問題，在實際應用過程中越來越多。4.基於目標建模的方法。該方法具有一定的針對性，需要提前知道所需跟蹤的目標是什么，比如說車輛，人臉，行人等。由於已經知道了跟蹤目標，所以必須對目標進行建模，然后利用該模型來進行跟蹤。該方法的局限性是必須提前知道所跟蹤的目標是什么，因而其推廣性比較差。
本文通過學習文獻[1],初步從理論上了解了粒子濾波的幾個步驟，然后參考文獻[2]中基於顏色直方圖的跟蹤算法，自己實現了一個粒子濾波目標跟蹤器。本文實現的算法只能進行單目標跟蹤，后面有學者將粒子濾波的思想擴展到多目標跟蹤，比如說文獻[3].
實現過程
1. 在攝像頭采集到的視頻序列中手動標注一個目標區域，用矩形框表示，並計算該目標區域內的直方圖，作為匹配模板。
2. 在選中的目標區域中心處撒預定值為100的粒子數目，並對這100個粒子的結構體都初始化為同樣的值，比如說粒子位置為目標區域中心，粒子矩形長寬為目標矩形長寬，粒子的初始尺寸為1.等等。
其粒子的結構體和注釋如下：
/****定義粒子結構體****/
typedef struct particle
{
int orix,oriy;//原始粒子坐標
int x,y;//當前粒子的坐標
double scale;//當前粒子窗口的尺寸
int prex,prey;//上一幀粒子的坐標
double prescale;//上一幀粒子窗口的尺寸
Rect rect;//當前粒子矩形窗口
Mat hist;//當前粒子窗口直方圖特征
double weight;//當前粒子權值
}PARTICLE;
3. 利用二階動態模型對這100個粒子進行隨機擴散，每個粒子根據二階隨機動態模型擴散到一個新的位置。
4. 計算以新位置處每個粒子為中心的矩形框內圖像的直方圖特征，然后每個特征都與目標模板直方圖進行比較，計算其相似度。
5. 根據計算得到的相似度算出每個粒子的權值，即相似度偉大的權值越大，並且對權值進行歸一化。
6. 取最大權值處的粒子中心為跟蹤目標中心，並根據粒子結構體中的尺寸成員算出目標的矩形框。
7. 根據上一個狀態的粒子分布情況，按照其權值乘以粒子總數100進行重采用，即權值大的粒子其采樣得到的粒子數也大。這樣可以重采用重新得到100個粒子。
8. 根據這100個粒子的情況，重新進入步驟3進行循環，經過對粒子權值分布的重采樣，根據二階隨機運動模型進行粒子擴散，根據直方圖的巴氏距離推斷粒子權值，然后取權值最大的那個作為目標，如此循環。
程序流程圖
實驗代碼
實驗代碼如下：
// particle_tracking.cpp : 定義控制台應用程序的入口點。 // #include "stdafx.h" #include <opencv2/core/core.hpp> #include "opencv2/imgproc/imgproc.hpp" #include <opencv2/highgui/highgui.hpp> #include <stdio.h> #include <iostream> using namespace cv; using namespace std; Rect select; bool select_flag=false; bool tracking=false;//跟蹤標志位 bool select_show=false; Point origin; Mat frame,hsv; int after_select_frames=0;//選擇矩形區域完后的幀計數 /****rgb空間用到的變量****/ //int hist_size[]={16,16,16};//rgb空間各維度的bin個數 //float rrange[]={0,255.0}; //float grange[]={0,255.0}; //float brange[]={0,255.0}; //const float *ranges[] ={rrange,grange,brange};//range相當於一個二維數組指針 /****hsv空間用到的變量****/ int hist_size[]={16,16,16}; float hrange[]={0,180.0}; float srange[]={0,256.0}; float vrange[]={0,256.0}; //int hist_size[]={32,32,32}; //float hrange[]={0,359.0.0}; //float srange[]={0,1.0}; //float vrange[]={0,1.0}; const float *ranges[]={hrange,srange,vrange}; int channels[]={0,1,2}; /****有關粒子窗口變化用到的相關變量****/ int A1=2; int A2=-1; int B0=1; double sigmax=1.0; double sigmay=0.5; double sigmas=0.001; /****定義使用粒子數目宏****/ #define PARTICLE_NUMBER 100 //如果這個數設定太大，經測試這個數字超過25就會報錯，則在運行時將會出現錯誤 /****定義粒子結構體****/ typedef struct particle { int orix,oriy;//原始粒子坐標 int x,y;//當前粒子的坐標 double scale;//當前粒子窗口的尺寸 int prex,prey;//上一幀粒子的坐標 double prescale;//上一幀粒子窗口的尺寸 Rect rect;//當前粒子矩形窗口 Mat hist;//當前粒子窗口直方圖特征 double weight;//當前粒子權值 }PARTICLE; PARTICLE particles[PARTICLE_NUMBER]; /************************************************************************************************************************/ /**** 如果采用這個onMouse()函數的話，則可以畫出鼠標拖動矩形框的4種情形 ****/ /************************************************************************************************************************/ void onMouse(int event,int x,int y,int,void*) { //Point origin;//不能在這個地方進行定義，因為這是基於消息響應的函數，執行完后origin就釋放了，所以達不到效果。 if(select_flag) { select.x=MIN(origin.x,x);//不一定要等鼠標彈起才計算矩形框，而應該在鼠標按下開始到彈起這段時間實時計算所選矩形框 select.y=MIN(origin.y,y); select.width=abs(x-origin.x);//算矩形寬度和高度 select.height=abs(y-origin.y); select&=Rect(0,0,frame.cols,frame.rows);//保證所選矩形框在視頻顯示區域之內 // rectangle(frame,select,Scalar(0,0,255),3,8,0);//顯示手動選擇的矩形框 } if(event==CV_EVENT_LBUTTONDOWN) { select_flag=true;//鼠標按下的標志賦真值 tracking=false; select_show=true; after_select_frames=0;//還沒開始選擇，或者重新開始選擇，計數為0 origin=Point(x,y);//保存下來單擊是捕捉到的點 select=Rect(x,y,0,0);//這里一定要初始化，因為在opencv中Rect矩形框類內的點是包含左上角那個點的，但是不含右下角那個點。 } else if(event==CV_EVENT_LBUTTONUP) { select_flag=false; tracking=true; select_show=false; after_select_frames=1;//選擇完后的那一幀當做第1幀 } } /****粒子權值降序排列函數****/ int particle_decrease(const void *p1,const void *p2) { PARTICLE* _p1=(PARTICLE*)p1; PARTICLE* _p2=(PARTICLE*)p2; if(_p1->weight<_p2->weight) return 1; else if(_p1->weight>_p2->weight) return -1; return 0;//相等的情況下返回0 } int main(int argc, unsigned char* argv[]) { char c; Mat target_img,track_img; Mat target_hist,track_hist; PARTICLE *pParticle; /***打開攝像頭****/ VideoCapture cam(0); if (!cam.isOpened()) return -1; /****讀取一幀圖像****/ cam>>frame; if(frame.empty()) return -1; VideoWriter output_dst( "demo.avi", CV_FOURCC('M', 'J', 'P', 'G'), 10, frame.size(), 1 ); /****建立窗口****/ namedWindow("camera",1);//顯示視頻原圖像的窗口 /****捕捉鼠標****/ setMouseCallback("camera",onMouse,0); while(1) { /****讀取一幀圖像****/ cam>>frame; if(frame.empty()) return -1; /****將rgb空間轉換為hsv空間****/ cvtColor(frame,hsv,CV_BGR2HSV); if(tracking) { if(1==after_select_frames)//選擇完目標區域后 { /****計算目標模板的直方圖特征****/ target_img=Mat(hsv,select);//在此之前先定義好target_img,然后這樣賦值也行，要學會Mat的這個操作 calcHist(&target_img,1,channels,Mat(),target_hist,3,hist_size,ranges); normalize(target_hist,target_hist); /****初始化目標粒子****/ pParticle=particles;//指針初始化指向particles數組 for(int x=0;x<PARTICLE_NUMBER;x++) { pParticle->x=cvRound(select.x+0.5*select.width);//選定目標矩形框中心為初始粒子窗口中心 pParticle->y=cvRound(select.y+0.5*select.height); pParticle->orix=pParticle->x;//粒子的原始坐標為選定矩形框(即目標)的中心 pParticle->oriy=pParticle->y; pParticle->prex=pParticle->x;//更新上一次的粒子位置 pParticle->prey=pParticle->y; pParticle->rect=select; pParticle->prescale=1; pParticle->scale=1; pParticle->hist=target_hist; pParticle->weight=0; pParticle++; } } else if(2==after_select_frames)//從第二幀開始就可以開始跟蹤了 { double sum=0.0; pParticle=particles; RNG rng;//隨機數產生器 /****更新粒子結構體的大部分參數****/ for(int i=0;i<PARTICLE_NUMBER;i++) { int x,y; int xpre,ypre; double s,pres; xpre=pParticle->x; ypre=pParticle->y; pres=pParticle->scale; /****更新粒子的矩形區域即粒子中心****/ x=cvRound(A1*(pParticle->x-pParticle->orix)+A2*(pParticle->prex-pParticle->orix)+ B0*rng.gaussian(sigmax)+pParticle->orix); pParticle->x=max(0,min(x,frame.cols-1)); y=cvRound(A1*(pParticle->y-pParticle->oriy)+A2*(pParticle->prey-pParticle->oriy)+ B0*rng.gaussian(sigmay)+pParticle->oriy); pParticle->y=max(0,min(y,frame.rows-1)); s=A1*(pParticle->scale-1)+A2*(pParticle->prescale-1)+B0*(rng.gaussian(sigmas))+1.0; pParticle->scale=max(1.0,min(s,3.0)); pParticle->prex=xpre; pParticle->prey=ypre; pParticle->prescale=pres; // pParticle->orix=pParticle->orix; // pParticle->oriy=pParticle->oriy; //注意在c語言中，x-1.0，如果x是int型，則這句語法有錯誤,但如果前面加了cvRound(x-0.5)則是正確的 pParticle->rect.x=max(0,min(cvRound(pParticle->x-0.5*pParticle->scale*pParticle->rect.width),frame.cols)); pParticle->rect.y=max(0,min(cvRound(pParticle->y-0.5*pParticle->scale*pParticle->rect.height),frame.rows)); pParticle->rect.width=min(cvRound(pParticle->rect.width),frame.cols-pParticle->rect.x); pParticle->rect.height=min(cvRound(pParticle->rect.height),frame.rows-pParticle->rect.y); // pParticle->rect.width=min(cvRound(pParticle->scale*pParticle->rect.width),frame.cols-pParticle->rect.x); // pParticle->rect.height=min(cvRound(pParticle->scale*pParticle->rect.height),frame.rows-pParticle->rect.y); /****計算粒子區域的新的直方圖特征****/ track_img=Mat(hsv,pParticle->rect); calcHist(&track_img,1,channels,Mat(),track_hist,3,hist_size,ranges); normalize(track_hist,track_hist); /****更新粒子的權值****/ // pParticle->weight=compareHist(target_hist,track_hist,CV_COMP_INTERSECT); //采用巴氏系數計算相似度,永遠與最開始的那一目標幀相比較 pParticle->weight=1.0-compareHist(target_hist,track_hist,CV_COMP_BHATTACHARYYA); /****累加粒子權值****/ sum+=pParticle->weight; pParticle++; } /****歸一化粒子權重****/ pParticle=particles; for(int i=0;i<PARTICLE_NUMBER;i++) { pParticle->weight/=sum; pParticle++; } /****根據粒子的權值降序排列****/ pParticle=particles; qsort(pParticle,PARTICLE_NUMBER,sizeof(PARTICLE),&particle_decrease); /****根據粒子權重重采樣粒子****/ PARTICLE newParticle[PARTICLE_NUMBER]; int np=0,k=0; for(int i=0;i<PARTICLE_NUMBER;i++) { np=cvRound(pParticle->weight*PARTICLE_NUMBER); for(int j=0;j<np;j++) { newParticle[k++]=particles[i]; if(k==PARTICLE_NUMBER) goto EXITOUT; } } while(k<PARTICLE_NUMBER) newParticle[k++]=particles[0]; EXITOUT: for(int i=0;i<PARTICLE_NUMBER;i++) particles[i]=newParticle[i]; }//end else //????????這個排序很慢，粒子數一多就卡 // qsort(pParticle,PARTICLE_NUMBER,sizeof(PARTICLE),&particle_decrease); /****計算粒子期望，采用所有粒子位置的期望值做為跟蹤結果****/ /*Rect_<double> rectTrackingTemp(0.0,0.0,0.0,0.0); pParticle=particles; for(int i=0;i<PARTICLE_NUMBER;i++) { rectTrackingTemp.x+=pParticle->rect.x*pParticle->weight; rectTrackingTemp.y+=pParticle->rect.y*pParticle->weight; rectTrackingTemp.width+=pParticle->rect.width*pParticle->weight; rectTrackingTemp.height+=pParticle->rect.height*pParticle->weight; pParticle++; }*/ /****計算最大權重目標的期望位置，作為跟蹤結果****/ Rect rectTrackingTemp(0,0,0,0); pParticle=particles; rectTrackingTemp.x=pParticle->x-0.5*pParticle->rect.width; rectTrackingTemp.y=pParticle->y-0.5*pParticle->rect.height; rectTrackingTemp.width=pParticle->rect.width; rectTrackingTemp.height=pParticle->rect.height; /****計算最大權重目標的期望位置，采用權值最大的1/4個粒子數作為跟蹤結果****/ /*Rect rectTrackingTemp(0,0,0,0); double weight_temp=0.0; pParticle=particles; for(int i=0;i<PARTICLE_NUMBER/4;i++) { weight_temp+=pParticle->weight; pParticle++; } pParticle=particles; for(int i=0;i<PARTICLE_NUMBER/4;i++) { pParticle->weight/=weight_temp; pParticle++; } pParticle=particles; for(int i=0;i<PARTICLE_NUMBER/4;i++) { rectTrackingTemp.x+=pParticle->rect.x*pParticle->weight; rectTrackingTemp.y+=pParticle->rect.y*pParticle->weight; rectTrackingTemp.width+=pParticle->rect.width*pParticle->weight; rectTrackingTemp.height+=pParticle->rect.height*pParticle->weight; pParticle++; }*/ /****計算最大權重目標的期望位置，采用所有粒子數作為跟蹤結果****/ /*Rect rectTrackingTemp(0,0,0,0); pParticle=particles; for(int i=0;i<PARTICLE_NUMBER;i++) { rectTrackingTemp.x+=cvRound(pParticle->rect.x*pParticle->weight); rectTrackingTemp.y+=cvRound(pParticle->rect.y*pParticle->weight); pParticle++; } pParticle=particles; rectTrackingTemp.width = pParticle->rect.width; rectTrackingTemp.height = pParticle->rect.height;*/ //創建目標矩形區域 Rect tracking_rect(rectTrackingTemp); pParticle=particles; /****顯示各粒子運動結果****/ for(int m=0;m<PARTICLE_NUMBER;m++) { rectangle(frame,pParticle->rect,Scalar(255,0,0),1,8,0); pParticle++; } /****顯示跟蹤結果****/ rectangle(frame,tracking_rect,Scalar(0,0,255),3,8,0); after_select_frames++;//總循環每循環一次，計數加1 if(after_select_frames>2)//防止跟蹤太長，after_select_frames計數溢出 after_select_frames=2; } if(select_show) rectangle(frame,select,Scalar(0,0,255),3,8,0);//顯示手動選擇的矩形框 output_dst<<frame; //顯示視頻圖片到窗口 imshow("camera",frame); // select.zeros(); //鍵盤響應 c=(char)waitKey(20); if(27==c)//ESC鍵 return -1; } return 0; }
實驗結果
實驗結果如下圖所示，其中紅色框為跟蹤到的目標，綠色框為跟蹤過程中粒子所在位置的矩形框。分為4組實驗對比，分別手動標注4個不同的部位進行跟蹤，達到了較好的效果。（此處圖略，圖片是用攝像頭跟蹤本人多個部位的結果，不好貼圖）
目標1的跟蹤效果：(略)
目標2的跟蹤效果：(略)
目標3的跟蹤效果：(略)
目標4的跟蹤效果：(略)
實驗總結
通過本次實驗可以更加深刻的理解了粒子濾波的算法思想，熟悉了粒子濾波算法的幾個步驟，並且對opencv這個工具有了更深一步的認識，編程功底也得到了相應的提高，體會 到了科研的樂趣。
參考文獻：
1. Isard, M. and A. Blake (1998). "Condensation—conditional density propagation for visual tracking." International journal of computer vision 29(1): 5-28.
2. Nummiaro, K., E. Koller-Meier, et al. (2003). "An adaptive color-based particle filter." Image and Vision Computing 21(1): 99-110.
3. Okuma, K., A. Taleghani, et al. (2004). "A boosted particle filter: Multitarget detection and tracking." Computer Vision-ECCV 2004: 28-39.