源代碼代碼解析

1、Mean Shift，我們翻譯為“均值飄移”。其在聚類，圖像平滑。圖像分割和跟蹤方面得到了比較廣泛的應用。由于本人目前研究跟蹤方面的東西，故此主要介紹利用Mean Shift方法進行目標跟蹤，從而對MeanShift 有一個比較全面的介紹。(以下某些部分轉載常峰學長的“Mean Shift概述”)Mean Shift這個概念最早是由Fukunaga等人于1975 年在一篇關于概率密度梯度函數(shù)的估計(The Estimation of the Gradient of a Density Function, with Applications in Pattern Recognition)中提出來

2、的,其最初含義正如其名，就是偏移的均值向量,在這里Mean Shift是一個名詞，它指代的是一個向量,但隨著Mean Shift理論的發(fā)展,Mean Shift的含義也發(fā)生了變化，如果 我們說Mean Shift算法,一般是指一個迭代的步驟,即先算出當前點的偏移均值，移動該點到其偏移均值，然后 以此為新的起始點，繼續(xù)移動,直到滿足一定的條件結束.然而在以后的很長一段時間內Mean Shift并沒有引起人們的注意，直到20年 以后,也就是1995年,另外一篇關于Mean Shift的重要文獻(Mean shift, mode seeking, and clustering)才發(fā)表.在這篇重要的文

3、獻中,Yizong Cheng對基本 的Mean Shift算法在以下兩個方面做了推廣，首先Yizong Cheng定義了一族核函 數(shù),使得隨著樣本與被偏移點的距離不同，其偏移量對均值偏移向量的貢獻也不同, 其次Yizong Cheng還設定了一個權重系數(shù)，使得不同的樣本點重要性不一樣，這大 大擴大了 Mean Shift的適用范圍.另外Yizong Cheng指出了 Mean Shift可能應用 的領域，并給出了具體的例子。Comaniciu 等人在還(Mean-shift Blob Tracking through Scale Space) 中把非剛體的跟蹤問題近似為一個Mean Shif

4、t最優(yōu)化問題,使得跟蹤可以實時的 進行。目前，利用Mean Shift進行跟蹤已經(jīng)相當成熟。目標跟蹤不是一個新的問題，目前在計算機視覺領域內有不少人在研究。所 謂跟蹤，就是通過已知的圖像幀中的目標位置找到目標在下一幀中的位置。下面主要以代碼形式展現(xiàn)Mean Shift在跟蹤中的應用。void CObjectTracker:ObjeckTrackerHandlerByUser(IplImage *frame)/ 跟蹤函數(shù) (m_cActiveObject = 0;if (m_sTrackingObjectTablem_cActiveObject.Status)(if (!m_sTrackingO

5、bjectTablem_cActiveObject.assignedAnObject)(FindHistogram(frame,m_sTrackingObjectTablem_cActiveObject.initHisto gram);m_sTrackingObjectTablem_cActiveObject.assignedAnObject = true;else(FindNextLocation(frame);/利用mean shift迭代找出目標下一個位置點DrawObjectBox(frame);void CObjectTracker:FindNextLocation(IplImage

6、 *frame)(int i, j, opti, optj;SINT16 scale3 = -3, 3, 0;FLOAT32 dist, optdist;SINT16 h, w, optX, optY;/try no-scalingFindNextFixScale(frame);/找出目標的下一個大致范圍optdist二LastDist;optX=m_sTrackingObjectTablem_cActiveObject.X;optY=m_sTrackingObjectTablem_cActiveObject.Y;/try one of the 9 possible scalingi=rand

7、()*2/RAND_MAX;j=rand()*2/RAND_MAX;h=m_sTrackingObjectTablem_cActiveObject.H;w=m_sTrackingObjectTablem_cActiveObject.W;if(h+scalei10 & w+scalej10 & h+scaleim_nImageHeight/2 &w+scalejm_nImageWidth/2)(m_sTrackingObjectTablem_cActiveObject.H=h+scalei;m_sTrackingObjectTablem_cActiveObject.W=w+scalej;Find

8、NextFixScale(frame);if( (dist二LastDist) n， dist);else /no scaling is better(m_sTrackingObjectTablem_cActiveObject.X=optX;m_sTrackingObjectTablem_cActiveObject.Y=optY;m_sTrackingObjectTablem_cActiveObject.H=h;m_sTrackingObjectTablem_cActiveObject.W=w;TotalDist+=optdist; /the latest distance/ printf(n

9、)；接著上面，這里仍然在跟蹤的基礎上講解mean shift。首先還是把mean shift的原理用數(shù)學公式說一下吧。 1、目標模型，算法采用的是特征值的加權概率分布來描述目標模型。這應該是模式識別中主要描述目標的 模型，不同于自動控制理論中采用的狀態(tài)方程。目標模型共m個特征值（可以理解為像素灰度值）g =/曰1,.辨其中X0是窗口中心點向量值（可能為RBG向量或者灰度值），Xi是窗口內第i點向量值。C為歸一 化常數(shù)，保障q1+q2+q3+qm=1 ，H為核函數(shù)的帶寬向量。M為特征值的個數(shù)，對應于圖像處理可以 理解為灰度等級劃分的個數(shù)，從而特征值u為對應的灰度等級。d函數(shù)為脈沖函數(shù)，保證只有具

10、有u特 征值的像素才對概率分布作出貢獻。從而k函數(shù)可以理解為u灰度值的一個加權頻數(shù)。2、匹配對象，也采用特征值加權概率分布5(Z?(A_) -：/) j e爵 I M)二 G其中，Y為匹配對象的中心，Xi5(Z?(A_) -：/) j e爵 I M)二 G其中，Y為匹配對象的中心，Xi是匹配窗口內第i點向量值，Hh為匹配窗口的核函數(shù)帶寬向量。Ch為匹配窗口特征向量的歸一化常數(shù)。3、匹配對象與目標模型的相似程度，相似函數(shù)可采用Bhattacharyya函數(shù)P(YP(Y.) + (P(Y)-p(Y定義，5WD = WL)+ 氣耳20(約以2 2十1 所 1 所因此， p(n=-zVK(W + ?Z

11、()Al.仃)要使得人(Y)向最大值迭代，只要Y的搜索方向與梯度方向一致即可，通過求導可得到Y0的梯度方向為:ffl pSY).q) = 2 如(4、匹配過程就是尋找相似函數(shù)最大值的尋優(yōu)過程，Mean-Shift采用的是梯度下降法。首先將(Y)在擴(Y0)附近進行泰勒級數(shù)展開，取前兩項。即:Vp()=獨I旗叫g( 1=1m I n涉次”)為權值。因此如果如下確定Y1,那么Y1-Y0將與梯度方向一致。以上為mean shift的數(shù)學原理。有關文字的敘述已經(jīng)在上一篇中提到了。用mean shift來跟蹤屬于確定性 算法，粒子濾波器屬于統(tǒng)計學方法。meanshift跟蹤算法相對于粒子濾波器來說可能實

12、時性更好一些，但是 跟蹤的準確性在理論上還是略遜于粒子濾波器的。mean shift跟蹤的的實質就是通過對應的模板來確定目 標的下一個位置。通過迭代找到新的中心點(即是目標的新的位置點)。有關跟蹤的code如下所示：/*Bilkent University:Mean-shift Tracker based Moving Object Tracker in VideoVersion: 1.0Compiler: Microsoft Visual C+ 6.0 (tested in both debug and releasemode)Modified by Mr Zhou*/#include Ob

13、jectTracker.h#include utils.h#include #include #include /* #define GetRValue(rgb) (UBYTE8) (rgb)#define GetGValue(rgb) (UBYTE8) (ULONG_32) (rgb) 8)#define GetBValue(rgb) (UBYTE8) (rgb) 16)*/#define RGB(r, g ,b) (ULONG_32) (UBYTE8) (r) | (UBYTE8) (g) 8) | (ULONG_32) (UBYTE8)(b) 16)#define min(a, b) (

14、a) (b) ? (a) : (b)#define MEANSHIFT_ITARATION_NO 5#define DISTANCE_ITARATION_NO 1#define ALPHA 1#define EDGE_DETECT_TRESHOLD 32/*1給定目標的初始位置和尺寸，計算目標在圖像中的直方圖；2輸入新圖像，迭代直到收斂：計算圖像上對應區(qū)域的新直方圖；新直方圖與目標直方圖比較，計算權重；根據(jù)權重,計算圖像上對應區(qū)域的形心/質心；根據(jù)形心，修正目標位置；直方圖分為兩部分，每部分大小4096,RGB的256*256*256種組合，縮減為16*16*16=4096種組合.如果目標區(qū)域

15、的點是邊緣點，則計入直方圖的后一部分，否則計入直方圖的前一部分.*/CObjectTracker:CObjectTracker(INT32 imW，INT32 imH，IMAGE_TYPE eImageType)(m_nImageWidth = imW;m_nImageHeight = imH;m_eIMAGE_TYPE = elmageType;m_cSkipValue = 0;for (UBYTE8 i=0;iMAX_OBJECT_TRACK_NUMBER;i+)/ 初始化各個目標(m_sTrackingObjectTablei.Status = false;for(SINT16 j=0;

16、jHISTOGRAM_LENGTH;j+)m_sTrackingObjectTablei.initHistogramj = 0;m_nFrameCtr = 0;m_uTotalTime = 0;m_nMaxEstimationTime = 0;m_cActiveObject = 0;TotalDist=0.0;LastDist=0.0;switch (eImageType)(case MD_RGBA:m_cSkipValue = 4 ;break ;case MD_RGB:m_cSkipValue = 3 ;break ;CObjectTracker:CObjectTracker()/retu

17、rns pixel values in format |0|B|G|R| wrt to (x.y)/*ULONG_32 CObjectTracker:GetPixelValues(UBYTE8 *frame,SINT16 x,SINT16 y)(ULONG_32 pixelValues = 0;pixelValues = *(frame+(y*m_nImageWidth+x)*m_cSkipValue+2)|/0BGR*(frame+(y*m_nImageWidth+x)*m_cSkipValue+1) 8| *(frame+(y*m_nImageWidth+x)*m_cSkipValue)

18、8) & 0 xFF);/ *(frame+(y*m_nImageWidth+x)*m_cSkipValue) = UBYTE8(pixelValues 16) & 0 xFF);/setpix32fsetpix8c(r, y, x, UBYTE8(pixelValues & 0 xFF);setpix8c(g, y, x, UBYTE8(pixelValues 8) & 0 xFF);setpix8c(b, y, x, UBYTE8(pixelValues 16) & 0 xFF);/ returns box colorULONG_32 CObjectTracker:GetBoxColor(

19、)(ULONG_32 pixelValues = 0;switch(m_cActiveObject)(case 0:pixelValues = RGB(255,0,0);break;case 1:pixelValues = RGB(0,255,0);break;case 2:pixelValues = RGB(0,0,255);break;case 3:pixelValues = RGB(255,255,0);break;case 4:pixelValues = RGB(255,0,255);break;case 5:pixelValues = RGB(0,255,255);break;cas

20、e 6:pixelValues = RGB(255,255,255);break;case 7:pixelValues = RGB(128,0,128);break;case 8:pixelValues = RGB(128,128,0);break;case 9:pixelValues = RGB(128,128,128);break;case 10:pixelValues = RGB(255,128,0);break;case 11:pixelValues = RGB(0,128,128);break;case 12:pixelValues = RGB(123,50,10);break;ca

21、se 13:pixelValues = RGB(10,240,126);break;case 14:pixelValues = RGB(0,128,255);break;case 15:pixelValues = RGB(128,200,20);break;default:break;初始化一個目標的參數(shù)void CObjectTracker:ObjectTrackerInitObjectParameters(SINT16 x,SINT16 y,SINT16 Width,SINT16Height)m_cActiveObject = 0;m_sTrackingObjectTablem_cActi

22、veObject.X = x;m_sTrackingObjectTablem_cActiveObject.Y = y;m_sTrackingObjectTablem_cActiveObject.W = Width;m_sTrackingObjectTablem_cActiveObject.H = Height;m_sTrackingObjectTablem_cActiveObject.vectorX = 0;m_sTrackingObjectTablem_cActiveObject.vectorY = 0;m_sTrackingObjectTablem_cActiveObject.Status

23、 = true;m_sTrackingObjectTablem_cActiveObject.assignedAnObject = false;進行一次跟蹤void CObjectTracker:ObjeckTrackerHandlerByUser(IplImage *frame)(m_cActiveObject = 0;if (m_sTrackingObjectTablem_cActiveObject.Status)(if (!m_sTrackingObjectTablem_cActiveObject.assignedAnObject)(計算目標的初始直方圖FindHistogram(fram

24、e,m_sTrackingObjectTablem_cActiveObject.initHistogram);m_sTrackingObjectTablem_cActiveObject.assignedAnObject = true;else(在圖像上搜索目標FindNextLocation(frame);DrawObjectBox(frame); /Extracts the histogram of box/frame.圖像/histogram:直方圖在圖像frame中計算當前目標的直方圖histogram直方圖分為兩部分，每部分大小4096,/RGB的256*256*256種組合，縮減為1

25、6*16*16=4096種組合如果目標區(qū)域的點是邊緣點，則計入直方圖的后一部分，否則計入直方圖的前一部分void CObjectTracker:FindHistogram(IplImage *frame, FLOAT32 (*histogram)(SINT16 i = 0;SINT16 x = 0;SINT16 y = 0;UBYTE8 E = 0;UBYTE8 qR = 0,qG = 0,qB = 0;/ ULONG_32 pixelValues = 0;UINT32 numberOfPixel = 0;IplImage* r, * g, * b;r = cvCreateImage( cvG

26、etSize(frame), frame-depth, 1 );g = cvCreateImage( cvGetSize(frame), frame-depth, 1 );b = cvCreateImage( cvGetSize(frame), frame-depth, 1 );cvCvtPixToPlane( frame, b, g, r, NULL ); /divide color image into separate planes r, g, b. The exact sequence doesnt matter.for (i=0;iHISTOGRAM_LENGTH;i+) /rese

27、t all histogram histogrami = 0.0;/for all the pixels in the regionfor(y=max(m_sTrackingObjectTablem_cActiveObject.Y-m_sTrackingObjectTablem_cActiveObject.H/2,0);y=min(m_sTrackingObjectTablem_cActiveObject.Y+m_sTrackingObjectTablem_cActiveObject.H/2,m_nImageHeight-1);y+)for (x=max(m_sTrackingObjectTa

28、blem_cActiveObject.X-m_sTrackingObjectTablem_cActiveObject.W/2,0);x=min(m_sTrackingObjectTablem_cActiveObject.X+m_sTrackingObjectTablem_cActiveObject.W /2,m_nImageWidth-1);x+)(邊緣信息：當前點與上下左右4點灰度差異是否超過閾值E = CheckEdgeExistance(r, g, b,x,y);qR = (UBYTE8)pixval8c( r, y, x )/16;/quantize R componentqG = (

29、UBYTE8)pixval8c( g, y, x )/16;/quantize G componentqB = (UBYTE8)pixval8c( b, y, x )/16;/quantize B componenthistogram4096*E+256*qR+16*qG+qB += 1; 根據(jù)邊緣信息，累計直方圖/HISTOGRAM_LENGTH=8192numberOfPixel+;for (i=0;iHISTOGRAM_LENGTH;i+) /normalizehistogrami = histogrami/numberOfPixel;/for (i=0;idepth, 1 );g =

30、cvCreateImage( cvGetSize(frame), frame-depth, 1 );b = cvCreateImage( cvGetSize(frame), frame-depth, 1 );cvCvtPixToPlane( frame, b, g, r, NULL );pixelValues = GetBoxColor();/the x left and right boundsx_sum =min(m_sTrackingObjectTablem_cActiveObject.X+m_sTrackingObjectTablem_cActiveObject.W/2+1,m_nIm

31、ageWidth-1);/ 右邊界x_diff =max(m_sTrackingObjectTablem_cActiveObject.X-m_sTrackingObjectTablem_cActiveObject.W/2,0);/左邊界/the y upper and lower boundsy_sum =min(m_sTrackingObjectTablem_cActiveObject.Y+m_sTrackingObjectTablem_cActiveObject.H/2+1,m_nImageHeight-1);/ 下邊界y_diff =max(m_sTrackingObjectTablem

32、_cActiveObject.Y-m_sTrackingObjectTablem_cActiveObject.H/2,0);/上邊界for (y=y_diff;y=y_sum;y+)(SetPixelValues(r, g, b,pixelValues,x_diff,y);SetPixelValues(r, g, b,pixelValues,x_diff+1,y);SetPixelValues(r, g, b,pixelValues,x_sum-1,y);SetPixelValues(r, g, b,pixelValues,x_sum,y);for (x=x_diff;x=x_sum;x+)S

33、etPixelValues(r, g, b,pixelValues,x,y_diff);SetPixelValues(r, g, b,pixelValues,x,y_diff+1);SetPixelValues(r, g, b,pixelValues,x,y_sum-1);SetPixelValues(r, g, b,pixelValues,x,y_sum);cvCvtPlaneToPix(b, g, r, NULL, frame);cvReleaseImage(&r);cvReleaseImage(&g);cvReleaseImage(&b);/ Computes weights and d

34、rives the new location of object in the next frame/frame.圖像/histogram:直方圖計算權重，更新目標的坐標void CObjectTracker:FindWightsAndCOM(IplImage *frame, FLOAT32 (*histogram)(SINT16 i = 0;SINT16 x = 0;SINT16 y = 0;UBYTE8 E = 0;FLOAT32 sumOfWeights = 0;SINT16 ptr = 0;UBYTE8 qR = 0,qG = 0,qB = 0;FLOAT32 newX = 0.0;F

35、LOAT32 newY = 0.0;/ ULONG_32 pixelValues = 0;IplImage* r, * g, * b;FLOAT32 *weights = new FLOAT32HISTOGRAM_LENGTH;for (i=0;i0.0 )weightsi = m_sTrackingObjectTablem_cActiveObject.initHistogrami/histogrami; /qu/pu(y0)elseweightsi = 0.0;r = cvCreateImage( cvGetSize(frame), frame-depth, 1 );g = cvCreate

36、Image( cvGetSize(frame), frame-depth, 1 );b = cvCreateImage( cvGetSize(frame), frame-depth, 1 );cvCvtPixToPlane( frame, b, g, r, NULL ); /divide color image into separate planes r, g, b. The exact sequence doesnt matter.for(y=max(m_sTrackingObjectTablem_cActiveObject.Y-m_sTrackingObjectTablem_cActiv

37、eObject.H/2,0);y=min(m_sTrackingObjectTablem_cActiveObject.Y+m_sTrackingObjectTablem_cActiveObject.H/2,m_nImageHeight-1);y+)for(x=max(m_sTrackingObjectTablem_cActiveObject.X-m_sTrackingObjectTablem_cActiveObject.W/2,0);x0)(m_sTrackingObjectTablem_cActiveObject.X = SINT16(newX/sumOfWeights) + 0.5); /

38、update locationm_sTrackingObjectTablem_cActiveObject.Y = SINT16(newY/sumOfWeights) + 0.5);cvReleaseImage(&r);cvReleaseImage(&g);cvReleaseImage(&b);delete weights, weights = 0;/ Returns the distance between two histograms.FLOAT32 CObjectTracker:FindDistance(FLOAT32 (*histogram)(SINT16 i = 0;FLOAT32 d

39、istance = 0;for(i=0;iHISTOGRAM_LENGTH;i+)distance += FLOAT32(sqrt(DOUBLE64(m_sTrackingObjectTablem_cActiveObject.initHistogrami*histogrami);return(sqrt(1-distance);/An alternative distance measurementFLOAT32 CObjectTracker:CompareHistogram(UBYTE8 (*histogram)(SINT16 i = 0;FLOAT32 distance = 0.0;FLOA

40、T32 difference = 0.0;for (i=0;i0)distance += difference;elsedistance -= difference;return(distance);/ Returns the edge insformation of a pixel at (x,y), assume a large jump of value around edge pixelsUBYTE8 CObjectTracker:CheckEdgeExistance(IplImage *r, IplImage *g, IplImage *b, SINT16_x,SINT16 _y)(

41、UBYTE8 E = 0;SINT16 GrayCenter = 0;SINT16 GrayLeft = 0;SINT16 GrayRight = 0;SINT16 GrayUp = 0;SINT16 GrayDown = 0;/ ULONG_32 pixelValues = 0;/ pixelValues = GetPixelValues(frame,_x,_y);GrayCenter = SINT16(3*pixval8c( r, _y, _x )+6*pixval8c( g, _y, _x )+pixval8c( b, _y, _x );if (_x0)(/ pixelValues =

42、GetPixelValues(frame,_x-1,_y);GrayLeft = SINT16(3*pixval8c( r, _y, _x-1 )+6*pixval8c( g, _y, _x-1 )+pixval8c( b, _y, _x-1 );if (_x 0)(/ pixelValues = GetPixelValues(frame,_x,_y-1);GrayUp = SINT16(3*pixval8c( r, _y-1, _x )+6*pixval8c( g, _y-1, _x )+pixval8c( b, _y-1, _x );if (_yEDGE_DETECT_TRESHOLD)E

43、 = 1;if (abs(GrayCenter-GrayRight)/10)EDGE_DETECT_TRESHOLD)E = 1;if (abs(GrayCenter-GrayUp)/10)EDGE_DETECT_TRESHOLD)E = 1;if (abs(GrayCenter-GrayDown)/10)EDGE_DETECT_TRESHOLD)E = 1;return(E);/ Alpha blending: used to update initial histogram by the current histogramvoid CObjectTracker:UpdateInitialHistogram(UBYTE8 (*histogram)(SINT16 i = 0;for (i=0; i10 & w+scalej10 & h+scaleim_nImageHeight/2 & w+scalejm_nImageWidth/2) (m