源代码代码解析

1、Mean Shift，我们翻译为“均值飘移”。其在聚类，图像平滑。图像分割和跟踪方面得到了比较广泛的应用。由于本人目前研究跟踪方面的东西，故此主要介绍利用Mean Shift方法进行目标跟踪，从而对MeanShift 有一个比较全面的介绍。(以下某些部分转载常峰学长的“Mean Shift概述”)Mean Shift这个概念最早是由Fukunaga等人于1975 年在一篇关于概率密度梯度函数的估计(The Estimation of the Gradient of a Density Function, with Applications in Pattern Recognition)中提出来

2、的,其最初含义正如其名，就是偏移的均值向量,在这里Mean Shift是一个名词，它指代的是一个向量,但随着Mean Shift理论的发展,Mean Shift的含义也发生了变化，如果 我们说Mean Shift算法,一般是指一个迭代的步骤,即先算出当前点的偏移均值，移动该点到其偏移均值，然后 以此为新的起始点，继续移动,直到满足一定的条件结束.然而在以后的很长一段时间内Mean Shift并没有引起人们的注意，直到20年 以后,也就是1995年,另外一篇关于Mean Shift的重要文献(Mean shift, mode seeking, and clustering)才发表.在这篇重要的文

3、献中,Yizong Cheng对基本 的Mean Shift算法在以下两个方面做了推广，首先Yizong Cheng定义了一族核函 数,使得随着样本与被偏移点的距离不同，其偏移量对均值偏移向量的贡献也不同, 其次Yizong Cheng还设定了一个权重系数，使得不同的样本点重要性不一样，这大 大扩大了 Mean Shift的适用范围.另外Yizong Cheng指出了 Mean Shift可能应用 的领域，并给出了具体的例子。Comaniciu 等人在还(Mean-shift Blob Tracking through Scale Space) 中把非刚体的跟踪问题近似为一个Mean Shif

4、t最优化问题,使得跟踪可以实时的 进行。目前，利用Mean Shift进行跟踪已经相当成熟。目标跟踪不是一个新的问题，目前在计算机视觉领域内有不少人在研究。所 谓跟踪，就是通过已知的图像帧中的目标位置找到目标在下一帧中的位置。下面主要以代码形式展现Mean Shift在跟踪中的应用。void CObjectTracker:ObjeckTrackerHandlerByUser(IplImage *frame)/ 跟踪函数 (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的原理用数学公式说一下吧。 1、目标模型，算法采用的是特征值的加权概率分布来描述目标模型。这应该是模式识别中主要描述目标的 模型，不同于自动控制理论中采用的状态方程。目标模型共m个特征值（可以理解为像素灰度值）g =/曰1,.辨其中X0是窗口中心点向量值（可能为RBG向量或者灰度值），Xi是窗口内第i点向量值。C为归一 化常数，保障q1+q2+q3+qm=1 ，H为核函数的带宽向量。M为特征值的个数，对应于图像处理可以 理解为灰度等级划分的个数，从而特征值u为对应的灰度等级。d函数为脉冲函数，保证只有具

10、有u特 征值的像素才对概率分布作出贡献。从而k函数可以理解为u灰度值的一个加权频数。2、匹配对象，也采用特征值加权概率分布5(Z?(A_) -：/) j e爵 I M)二 G其中，Y为匹配对象的中心，Xi5(Z?(A_) -：/) j e爵 I M)二 G其中，Y为匹配对象的中心，Xi是匹配窗口内第i点向量值，Hh为匹配窗口的核函数带宽向量。Ch为匹配窗口特征向量的归一化常数。3、匹配对象与目标模型的相似程度，相似函数可采用Bhattacharyya函数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、匹配过程就是寻找相似函数最大值的寻优过程，Mean-Shift采用的是梯度下降法。首先将(Y)在扩(Y0)附近进行泰勒级数展开，取前两项。即:Vp()=独I旗叫g( 1=1m I n涉次”)为权值。因此如果如下确定Y1,那么Y1-Y0将与梯度方向一致。以上为mean shift的数学原理。有关文字的叙述已经在上一篇中提到了。用mean shift来跟踪属于确定性 算法，粒子滤波器属于统计学方法。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输入新图像，迭代直到收敛：计算图像上对应区域的新直方图；新直方图与目标直方图比较，计算权重；根据权重,计算图像上对应区域的形心/质心；根据形心，修正目标位置；直方图分为两部分，每部分大小4096,RGB的256*256*256种组合，缩减为16*16*16=4096种组合.如果目标区域

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;初始化一个目标的参数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种组合如果目标区域的点是边缘点，则计入直方图的后一部分，否则计入直方图的前一部分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; 根据边缘信息，累计直方图/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