碼上快樂

相關內容    簡體    繁體

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

本文轉載自   查看原文   2015-04-17 12:52   10840    圖像處理

這次將介紹基於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 plaincopy在CODE上查看代碼片派生到我的代碼片
  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 plaincopy在CODE上查看代碼片派生到我的代碼片
  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. }  

運行結果：

 

初始目標直方圖：

候選目標直方圖：


免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。



猜您在找 目標跟蹤算法meanshift優缺點 目標跟蹤學習筆記_1(opencv中meanshift和camshift例子的應用) 目標跟蹤算法綜述 單目標跟蹤算法SiamRPN TLD目標跟蹤算法 目標跟蹤算法 目標跟蹤算法分類（一） 【目標跟蹤+光流算法】運動跟蹤合集 運動目標跟蹤算法綜述 多目標跟蹤算法——SORT
 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM