基於MeanShift的目標跟蹤演算法、實現

馬衛飛發表於2017-12-10
原文網址 : https://blog.csdn.net/maweifei/article/details/78766784
演算法
這次將介紹基於MeanShift的目標跟蹤演算法，首先談談簡介，然後給出演算法實現流程，最後實現了一個單目標跟蹤的MeanShift演算法【matlab/c兩個版本】

      csdn貼公式比較煩，原諒我直接截圖了…

 

一、簡介

     首先扯扯無參密度估計理論，無參密度估計也叫做非引數估計，屬於數理統計的一個分支，和引數密度估計共同構成了概率密度估計方法。引數密度估計方法要求特徵空間服從一個已知的概率密度函式，在實際的應用中這個條件很難達到。而無引數密度估計方法對先驗知識要求最少，完全依靠訓練資料進行估計，並且可以用於任意形狀的密度估計。所以依靠無參密度估計方法，即不事先規定概率密度函式的結構形式，在某一連續點處的密度函式值可由該點鄰域中的若干樣本點估計得出。常用的無參密度估計方法有：直方圖法、最近鄰域法和核密度估計法。

     MeanShift演算法正是屬於核密度估計法，它不需要任何先驗知識而完全依靠特徵空間中樣本點的計算其密度函式值。對於一組取樣資料，直方圖法通常把資料的值域分成若干相等的區間，資料按區間分成若干組，每組資料的個數與總引數個數的比率就是每個單元的概率值；核密度估計法的原理相似於直方圖法，只是多了一個用於平滑資料的核函式。採用核函式估計法，在取樣充分的情況下，能夠漸進地收斂於任意的密度函式，即可以對服從任何分佈的資料進行密度估計。

     然後談談MeanShift的基本思想及物理含義：

    此外，從公式1中可以看到，只要是落入Sh的取樣點，無論其離中心x的遠近，對最終的Mh(x)計算的貢獻是一樣的。然而在現實跟蹤過程中，當跟蹤目標出現遮擋等影響時，由於外層的畫素值容易受遮擋或背景的影響，所以目標模型中心附近的畫素比靠外的畫素更可靠。因此，對於所有采樣點，每個樣本點的重要性應該是不同的，離中心點越遠，其權值應該越小。故引入核函式和權重係數來提高跟蹤演算法的魯棒性並增加搜尋跟蹤能力。

      接下來，談談核函式：

    核函式也叫視窗函式，在核估計中起到平滑的作用。常用的核函式有：Uniform，Epannechnikov，Gaussian等。本文演算法只用到了Epannechnikov，它數序定義如下：

二、基於MeanShift的目標跟蹤演算法

     基於均值漂移的目標跟蹤演算法通過分別計算目標區域和候選區域內畫素的特徵值概率得到關於目標模型和候選模型的描述，然後利用相似函式度量初始幀目標模型和當前幀的候選模版的相似性，選擇使相似函式最大的候選模型並得到關於目標模型的Meanshift向量，這個向量正是目標由初始位置向正確位置移動的向量。由於均值漂移演算法的快速收斂性，通過不斷迭代計算Meanshift向量，演算法最終將收斂到目標的真實位置，達到跟蹤的目的。

     下面通過圖示直觀的說明MeanShift跟蹤演算法的基本原理。如下圖所示：目標跟蹤開始於資料點xi0（空心圓點xi0，xi1，…，xiN表示的是中心點，上標表示的是的迭代次數，周圍的黑色圓點表示不斷移動中的視窗樣本點，虛線圓圈代表的是密度估計視窗的大小）。箭頭表示樣本點相對於核函式中心點的漂移向量，平均的漂移向量會指向樣本點最密集的方向，也就是梯度方向。因為 Meanshift 演算法是收斂的，因此在當前幀中通過反覆迭代搜尋特徵空間中樣本點最密集的區域，搜尋點沿著樣本點密度增加的方向“漂移”到區域性密度極大點點xiN，也就是被認為的目標位置，從而達到跟蹤的目的，MeanShift 跟蹤過程結束。

 

 

運動目標的實現過程【具體演算法】：

 

三、程式碼實現

說明：

1.       RGB顏色空間刨分，採用16*16*16的直方圖

2.       目標模型和候選模型的概率密度計算公式參照上文

3.       opencv版本執行：按P停止，擷取目標，再按P，進行單目標跟蹤

4.       Matlab版本，將視訊改為圖片序列，第一幀停止，手工標定目標，雙擊目標區域，進行單目標跟蹤。

 

matlab版本：

 

[plain] view plain copy
  1. function [] = select()  
  2. close all;  
  3. clear all;  
  4. %%%%%%%%%%%%%%%%%%根據一幅目標全可見的影象圈定跟蹤目標%%%%%%%%%%%%%%%%%%%%%%%  
  5. I=imread('result72.jpg');  
  6. figure(1);  
  7. imshow(I);  
  8.   
  9.   
  10. [temp,rect]=imcrop(I);  
  11. [a,b,c]=size(temp);         %a:row,b:col  
  12.   
  13.   
  14. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%計算目標影象的權值矩陣%%%%%%%%%%%%%%%%%%%%%%%  
  15. y(1)=a/2;  
  16. y(2)=b/2;  
  17. tic_x=rect(1)+rect(3)/2;  
  18. tic_y=rect(2)+rect(4)/2;  
  19. m_wei=zeros(a,b);%權值矩陣  
  20. h=y(1)^2+y(2)^2 ;%頻寬  
  21.   
  22.   
  23. for i=1:a  
  24.     for j=1:b  
  25.         dist=(i-y(1))^2+(j-y(2))^2;  
  26.         m_wei(i,j)=1-dist/h; %epanechnikov profile  
  27.     end  
  28. end  
  29. C=1/sum(sum(m_wei));%歸一化係數  
  30.   
  31.   
  32. %計算目標權值直方圖qu  
  33. %hist1=C*wei_hist(temp,m_wei,a,b);%target model  
  34. hist1=zeros(1,4096);  
  35. for i=1:a  
  36.     for j=1:b  
  37.         %rgb顏色空間量化為16*16*16 bins  
  38.         q_r=fix(double(temp(i,j,1))/16);  %fix為趨近0取整函式  
  39.         q_g=fix(double(temp(i,j,2))/16);  
  40.         q_b=fix(double(temp(i,j,3))/16);  
  41.         q_temp=q_r*256+q_g*16+q_b;            %設定每個畫素點紅色、綠色、藍色分量所佔比重  
  42.         hist1(q_temp+1)= hist1(q_temp+1)+m_wei(i,j);    %計算直方圖統計中每個畫素點佔的權重  
  43.     end  
  44. end  
  45. hist1=hist1*C;  
  46. rect(3)=ceil(rect(3));  
  47. rect(4)=ceil(rect(4));  
  48.   
  49.   
  50.   
  51.   
  52. %%%%%%%%%%%%%%%%%%%%%%%%%讀取序列影象  
  53. myfile=dir('D:\matlab7\work\mean shift\image\*.jpg');  
  54. lengthfile=length(myfile);  
  55.   
  56.   
  57. for l=1:lengthfile  
  58.     Im=imread(myfile(l).name);  
  59.     num=0;  
  60.     Y=[2,2];  
  61.       
  62.       
  63.     %%%%%%%mean shift迭代  
  64.     while((Y(1)^2+Y(2)^2>0.5)&num<20)   %迭代條件  
  65.         num=num+1;  
  66.         temp1=imcrop(Im,rect);  
  67.         %計算侯選區域直方圖  
  68.         %hist2=C*wei_hist(temp1,m_wei,a,b);%target candidates pu  
  69.         hist2=zeros(1,4096);  
  70.         for i=1:a  
  71.             for j=1:b  
  72.                 q_r=fix(double(temp1(i,j,1))/16);  
  73.                 q_g=fix(double(temp1(i,j,2))/16);  
  74.                 q_b=fix(double(temp1(i,j,3))/16);  
  75.                 q_temp1(i,j)=q_r*256+q_g*16+q_b;  
  76.                 hist2(q_temp1(i,j)+1)= hist2(q_temp1(i,j)+1)+m_wei(i,j);  
  77.             end  
  78.         end  
  79.         hist2=hist2*C;  
  80.         figure(2);  
  81.         subplot(1,2,1);  
  82.         plot(hist2);  
  83.         hold on;  
  84.           
  85.         w=zeros(1,4096);  
  86.         for i=1:4096  
  87.             if(hist2(i)~=0) %不等於  
  88.                 w(i)=sqrt(hist1(i)/hist2(i));  
  89.             else  
  90.                 w(i)=0;  
  91.             end  
  92.         end  
  93.           
  94.           
  95.           
  96.         %變數初始化  
  97.         sum_w=0;  
  98.         xw=[0,0];  
  99.         for i=1:a;  
  100.             for j=1:b  
  101.                 sum_w=sum_w+w(uint32(q_temp1(i,j))+1);  
  102.                 xw=xw+w(uint32(q_temp1(i,j))+1)*[i-y(1)-0.5,j-y(2)-0.5];  
  103.             end  
  104.         end  
  105.         Y=xw/sum_w;  
  106.         %中心點位置更新  
  107.         rect(1)=rect(1)+Y(2);  
  108.         rect(2)=rect(2)+Y(1);  
  109.     end  
  110.       
  111.       
  112.     %%%跟蹤軌跡矩陣%%%  
  113.     tic_x=[tic_x;rect(1)+rect(3)/2];  
  114.     tic_y=[tic_y;rect(2)+rect(4)/2];  
  115.       
  116.     v1=rect(1);  
  117.     v2=rect(2);  
  118.     v3=rect(3);  
  119.     v4=rect(4);  
  120.     %%%顯示跟蹤結果%%%  
  121.     subplot(1,2,2);  
  122.     imshow(uint8(Im));  
  123.     title('目標跟蹤結果及其運動軌跡');  
  124.     hold on;  
  125.     plot([v1,v1+v3],[v2,v2],[v1,v1],[v2,v2+v4],[v1,v1+v3],[v2+v4,v2+v4],[v1+v3,v1+v3],[v2,v2+v4],'LineWidth',2,'Color','r');  
  126.     plot(tic_x,tic_y,'LineWidth',2,'Color','b');  
  127.       
  128.       
  129. end  


 執行結果：

 

 

 

opencv版本：

[cpp] view plain copy
  1. #include "stdafx.h"  
  2. #include "cv.h"  
  3. #include "highgui.h"  
  4. #define  u_char unsigned char  
  5. #define  DIST 0.5  
  6. #define  NUM 20  
  7.   
  8. //全域性變數  
  9. bool pause = false;  
  10. bool is_tracking = false;  
  11. CvRect drawing_box;  
  12. IplImage *current;  
  13. double *hist1, *hist2;  
  14. double *m_wei;                                                                  //權值矩陣  
  15. double C = 0.0;                                                                //歸一化係數  
  16.   
  17. void init_target(double *hist1, double *m_wei, IplImage *current)  
  18. {  
  19.     IplImage *pic_hist = 0;  
  20.     int t_h, t_w, t_x, t_y;  
  21.     double h, dist;  
  22.     int i, j;  
  23.     int q_r, q_g, q_b, q_temp;  
  24.       
  25.     t_h = drawing_box.height;  
  26.     t_w = drawing_box.width;  
  27.     t_x = drawing_box.x;  
  28.     t_y = drawing_box.y;  
  29.   
  30.     h = pow(((double)t_w)/2,2) + pow(((double)t_h)/2,2);            //頻寬  
  31.     pic_hist = cvCreateImage(cvSize(300,200),IPL_DEPTH_8U,3);     //生成直方圖影象  
  32.   
  33.     //初始化權值矩陣和目標直方圖  
  34.     for (i = 0;i < t_w*t_h;i++)  
  35.     {  
  36.         m_wei[i] = 0.0;  
  37.     }  
  38.   
  39.     for (i=0;i<4096;i++)  
  40.     {  
  41.         hist1[i] = 0.0;  
  42.     }  
  43.   
  44.     for (i = 0;i < t_h; i++)  
  45.     {  
  46.         for (j = 0;j < t_w; j++)  
  47.         {  
  48.             dist = pow(i - (double)t_h/2,2) + pow(j - (double)t_w/2,2);  
  49.             m_wei[i * t_w + j] = 1 - dist / h;   
  50.             //printf("%f\n",m_wei[i * t_w + j]);  
  51.             C += m_wei[i * t_w + j] ;  
  52.         }  
  53.     }  
  54.   
  55.     //計算目標權值直方  
  56.     for (i = t_y;i < t_y + t_h; i++)  
  57.     {  
  58.         for (j = t_x;j < t_x + t_w; j++)  
  59.         {  
  60.             //rgb顏色空間量化為16*16*16 bins  
  61.             q_r = ((u_char)current->imageData[i * current->widthStep + j * 3 + 2]) / 16;  
  62.             q_g = ((u_char)current->imageData[i * current->widthStep + j * 3 + 1]) / 16;  
  63.             q_b = ((u_char)current->imageData[i * current->widthStep + j * 3 + 0]) / 16;  
  64.             q_temp = q_r * 256 + q_g * 16 + q_b;  
  65.             hist1[q_temp] =  hist1[q_temp] +  m_wei[(i - t_y) * t_w + (j - t_x)] ;  
  66.         }  
  67.     }  
  68.   
  69.     //歸一化直方圖  
  70.     for (i=0;i<4096;i++)  
  71.     {  
  72.         hist1[i] = hist1[i] / C;  
  73.         //printf("%f\n",hist1[i]);  
  74.     }  
  75.   
  76.     //生成目標直方圖  
  77.     double temp_max=0.0;  
  78.   
  79.     for (i = 0;i < 4096;i++)         //求直方圖最大值，為了歸一化  
  80.     {  
  81.         //printf("%f\n",val_hist[i]);  
  82.         if (temp_max < hist1[i])  
  83.         {  
  84.             temp_max = hist1[i];  
  85.         }  
  86.     }  
  87.     //畫直方圖  
  88.     CvPoint p1,p2;  
  89.     double bin_width=(double)pic_hist->width/4096;  
  90.     double bin_unith=(double)pic_hist->height/temp_max;  
  91.   
  92.     for (i = 0;i < 4096; i++)  
  93.     {  
  94.         p1.x = i * bin_width;  
  95.         p1.y = pic_hist->height;  
  96.         p2.x = (i + 1)*bin_width;  
  97.         p2.y = pic_hist->height - hist1[i] * bin_unith;  
  98.         //printf("%d,%d,%d,%d\n",p1.x,p1.y,p2.x,p2.y);  
  99.         cvRectangle(pic_hist,p1,p2,cvScalar(0,255,0),-1,8,0);  
  100.     }  
  101.     cvSaveImage("hist1.jpg",pic_hist);  
  102.     cvReleaseImage(&pic_hist);  
  103. }  
  104.   
  105. void MeanShift_Tracking(IplImage *current)  
  106. {  
  107.     int num = 0, i = 0, j = 0;  
  108.     int t_w = 0, t_h = 0, t_x = 0, t_y = 0;  
  109.     double *w = 0, *hist2 = 0;  
  110.     double sum_w = 0, x1 = 0, x2 = 0,y1 = 2.0, y2 = 2.0;  
  111.     int q_r, q_g, q_b;  
  112.     int *q_temp;  
  113.     IplImage *pic_hist = 0;  
  114.   
  115.     t_w = drawing_box.width;  
  116.     t_h = drawing_box.height;  
  117.       
  118.     pic_hist = cvCreateImage(cvSize(300,200),IPL_DEPTH_8U,3);     //生成直方圖影象  
  119.     hist2 = (double *)malloc(sizeof(double)*4096);  
  120.     w = (double *)malloc(sizeof(double)*4096);  
  121.     q_temp = (int *)malloc(sizeof(int)*t_w*t_h);  
  122.   
  123.     while ((pow(y2,2) + pow(y1,2) > 0.5)&& (num < NUM))  
  124.     {  
  125.         num++;  
  126.         t_x = drawing_box.x;  
  127.         t_y = drawing_box.y;  
  128.         memset(q_temp,0,sizeof(int)*t_w*t_h);  
  129.         for (i = 0;i<4096;i++)  
  130.         {  
  131.             w[i] = 0.0;  
  132.             hist2[i] = 0.0;  
  133.         }  
  134.   
  135.         for (i = t_y;i < t_h + t_y;i++)  
  136.         {  
  137.             for (j = t_x;j < t_w + t_x;j++)  
  138.             {  
  139.                 //rgb顏色空間量化為16*16*16 bins  
  140.                 q_r = ((u_char)current->imageData[i * current->widthStep + j * 3 + 2]) / 16;  
  141.                 q_g = ((u_char)current->imageData[i * current->widthStep + j * 3 + 1]) / 16;  
  142.                 q_b = ((u_char)current->imageData[i * current->widthStep + j * 3 + 0]) / 16;  
  143.                 q_temp[(i - t_y) *t_w + j - t_x] = q_r * 256 + q_g * 16 + q_b;  
  144.                 hist2[q_temp[(i - t_y) *t_w + j - t_x]] =  hist2[q_temp[(i - t_y) *t_w + j - t_x]] +  m_wei[(i - t_y) * t_w + j - t_x] ;  
  145.             }  
  146.         }  
  147.   
  148.         //歸一化直方圖  
  149.         for (i=0;i<4096;i++)  
  150.         {  
  151.             hist2[i] = hist2[i] / C;  
  152.             //printf("%f\n",hist2[i]);  
  153.         }  
  154.         //生成目標直方圖  
  155.         double temp_max=0.0;  
  156.   
  157.         for (i=0;i<4096;i++)         //求直方圖最大值，為了歸一化  
  158.         {  
  159.             if (temp_max < hist2[i])  
  160.             {  
  161.                 temp_max = hist2[i];  
  162.             }  
  163.         }  
  164.         //畫直方圖  
  165.         CvPoint p1,p2;  
  166.         double bin_width=(double)pic_hist->width/(4368);  
  167.         double bin_unith=(double)pic_hist->height/temp_max;  
  168.   
  169.         for (i = 0;i < 4096; i++)  
  170.         {  
  171.             p1.x = i * bin_width;  
  172.             p1.y = pic_hist->height;  
  173.             p2.x = (i + 1)*bin_width;  
  174.             p2.y = pic_hist->height - hist2[i] * bin_unith;  
  175.             cvRectangle(pic_hist,p1,p2,cvScalar(0,255,0),-1,8,0);  
  176.         }  
  177.         cvSaveImage("hist2.jpg",pic_hist);  
  178.       
  179.         for (i = 0;i < 4096;i++)  
  180.         {  
  181.             if (hist2[i] != 0)  
  182.             {  
  183.                 w[i] = sqrt(hist1[i]/hist2[i]);  
  184.             }else  
  185.             {  
  186.                 w[i] = 0;  
  187.             }  
  188.         }  
  189.               
  190.         sum_w = 0.0;  
  191.         x1 = 0.0;  
  192.         x2 = 0.0;  
  193.   
  194.         for (i = 0;i < t_h; i++)  
  195.         {  
  196.             for (j = 0;j < t_w; j++)  
  197.             {  
  198.                 //printf("%d\n",q_temp[i * t_w + j]);  
  199.                 sum_w = sum_w + w[q_temp[i * t_w + j]];  
  200.                 x1 = x1 + w[q_temp[i * t_w + j]] * (i - t_h/2);  
  201.                 x2 = x2 + w[q_temp[i * t_w + j]] * (j - t_w/2);  
  202.             }  
  203.         }  
  204.         y1 = x1 / sum_w;  
  205.         y2 = x2 / sum_w;  
  206.           
  207.         //中心點位置更新  
  208.         drawing_box.x += y2;  
  209.         drawing_box.y += y1;  
  210.   
  211.         //printf("%d,%d\n",drawing_box.x,drawing_box.y);  
  212.     }  
  213.     free(hist2);  
  214.     free(w);  
  215.     free(q_temp);  
  216.     //顯示跟蹤結果  
  217.     cvRectangle(current,cvPoint(drawing_box.x,drawing_box.y),cvPoint(drawing_box.x+drawing_box.width,drawing_box.y+drawing_box.height),CV_RGB(255,0,0),2);  
  218.     cvShowImage("Meanshift",current);  
  219.     //cvSaveImage("result.jpg",current);  
  220.     cvReleaseImage(&pic_hist);  
  221. }  
  222.   
  223. void onMouse( int event, int x, int y, int flags, void *param )  
  224. {  
  225.     if (pause)  
  226.     {  
  227.         switch(event)  
  228.         {  
  229.         case CV_EVENT_LBUTTONDOWN:   
  230.             //the left up point of the rect  
  231.             drawing_box.x=x;  
  232.             drawing_box.y=y;  
  233.             break;  
  234.         case CV_EVENT_LBUTTONUP:  
  235.             //finish drawing the rect (use color green for finish)  
  236.             drawing_box.width=x-drawing_box.x;  
  237.             drawing_box.height=y-drawing_box.y;  
  238.             cvRectangle(current,cvPoint(drawing_box.x,drawing_box.y),cvPoint(drawing_box.x+drawing_box.width,drawing_box.y+drawing_box.height),CV_RGB(255,0,0),2);  
  239.             cvShowImage("Meanshift",current);  
  240.               
  241.             //目標初始化  
  242.             hist1 = (double *)malloc(sizeof(double)*16*16*16);  
  243.             m_wei =  (double *)malloc(sizeof(double)*drawing_box.height*drawing_box.width);  
  244.             init_target(hist1, m_wei, current);  
  245.             is_tracking = true;  
  246.             break;  
  247.         }  
  248.         return;  
  249.     }  
  250. }  
  251.   
  252.   
  253.   
  254. int _tmain(int argc, _TCHAR* argv[])  
  255. {  
  256.     CvCapture *capture=cvCreateFileCapture("test.avi");  
  257.     current = cvQueryFrame(capture);  
  258.     char res[20];  
  259.     int nframe = 0;  
  260.   
  261.     while (1)  
  262.     {     
  263.     /*  sprintf(res,"result%d.jpg",nframe); 
  264.         cvSaveImage(res,current); 
  265.         nframe++;*/  
  266.         if(is_tracking)  
  267.         {  
  268.             MeanShift_Tracking(current);  
  269.         }  
  270.   
  271.         int c=cvWaitKey(1);  
  272.         //暫停  
  273.         if(c == 'p')   
  274.         {  
  275.             pause = true;  
  276.             cvSetMouseCallback( "Meanshift", onMouse, 0 );  
  277.         }  
  278.         while(pause){  
  279.             if(cvWaitKey(0) == 'p')  
  280.                 pause = false;  
  281.         }  
  282.         cvShowImage("Meanshift",current);  
  283.         current = cvQueryFrame(capture); //抓取一幀  
  284.     }  
  285.   
  286.     cvNamedWindow("Meanshift",1);  
  287.     cvReleaseCapture(&capture);  
  288.     cvDestroyWindow("Meanshift");  
  289.     return 0;  
  290. }  

執行結果：

 

初始目標直方圖：


候選目標直方圖：


原始碼及素材的下載地址我過會在評論裡給出。

相關文章