视频跟踪之均值漂移算法的实现

meanshift(均值漂移算法)实际上是一种基于梯度的搜索算法，先看看meanshift简单的推导:给定d维空间Rd中的几个样本点，任选一点x0,,(n维)

meanshift向量的基本形式定义为:

                                                                                   

其中Sk是以点x0为球心，半径为h的一个高维球，若是二维，那自然就是一个圆了。k为落入高维球区域内样本的个数.而且这个区域的所有点y满足以下关系:

                                                                                       

                                                                                       

如下图，假设样本点处在一个二维空间内

                                                                                               

如图蓝色点是我们选定的迭代的初始点即x0，将蓝色圆(圆的半径为h)内所有向量相加，相加的结果为如图为黑色向量，其终点指向的是如图所示的红色点，然后以红色点为圆心，h为半径再画一个圆，再求这个圆内以圆心为起点所有向量的和。如此迭代下去，圆的中心点为收敛于一个固定的点，也就是概率密度最大的地方。

以上只是meanshift大概的解释，还没有涉及到核函数.

我们知道密度估计理论分为两种：一种是参数密度估计，如高斯建模，即对模型的参数进行估计，其估计方法有两个：最大似然估计和贝叶斯估计。另一种是无参密度估计：如直方图，最近邻域法，核密度估计等。

核密度估计有一点类似于直方图，它也把数据值域分成若干相等的区间，每个区间称为一个bin,将样本数据按照这些区间划分到不同的bin，每个bin中样本的个数与总样本数量的比值就是每个bin的值，如:假设有1,2,3,4,5,6,7,8,9这9个样本，1到3之间的数分到一个bin（称为bin1）,4到6的分到一个bin(bin2)，7到9的分到一个bin(bin3)，则bin1的值为1/3,bin2=1/3,bin3=1/3;

相对于直方图法，核密度估计多了一个平滑数据的核函数，如高斯核函数。之所以要用到核函数，主要是对样本点进行加权，我们知道对于高维球内所有的样本点，他们对求解的贡献是不一样的，那么怎样去衡量这个贡献呢，为此便引入了核函数，说白了，核函数就是用来求每个样本点的权重的。

例如我们选择高斯函数作为核函数  ，引入核函数后的meanshift向量为

                                                    

  其中n为搜索区域内样本的个数.                                                                                      

假设我们搜素的区域是二维，且上一次迭代得到的搜索区域的中心坐标为x0,那么我们把x0带入上面的函数，即可得到漂移向量m，得到漂移向量之后就可以得到本次迭代搜索区域的中心点坐标了。直到整个迭代的误差小于我们设定的值(一般小于一个像素就可以了)，或者迭代次数大于我们设定的次数。

下面是结合opencv实现的meanshift算法

kernel_mean.h主要实现带核函数的meanshift算法

#include "cv/cv.h"#include "cv/cxerror.h"#include <iostream>using namespace std;#define e 2.71828#define PI 3.14159struct center{float x;float y;};CvRect kernel_meanshift(IplImage* image, CvRect windowIn,//计算概率密度图的重心windowIn为当前帧的初始搜索窗口int criteria,  float p[], float C, float h)//{int eps = 0.5;int height = windowIn.height;int width = windowIn.width;float *w = new float[4096];center c_center;center n_center;//迭代之后目标中心for(int times = 0; times <  criteria; times++){for(int wn = 0; wn < 4096; wn++)w[wn] = 0;float q[4096] = {0};//存放当前搜索窗口计算出来的核概率密度c_center.x = windowIn.x + windowIn.width/2 ;//获得搜索窗口的中心c_center.y = windowIn.y + windowIn.height/2 ;unsigned char *pdata = (unsigned char *)image->imageData;for(int j = windowIn.y; j < windowIn.height + windowIn.y; j++){for(int k =windowIn.x; k < windowIn.width + windowIn.x; k++){int delta_x = k - windowIn.width/2 - windowIn.x;int delta_y = j - windowIn.height/2 - windowIn.y;float d = sqrt(float(delta_x*delta_x) + float(delta_y*delta_y))/h;int r = pdata[j * image->widthStep + k * 3 + 2];int g = pdata[j * image->widthStep + k * 3 + 1];int b = pdata[j * image->widthStep + k * 3 + 0];int index = (r/16)*256 +(g/16)*16 +b/16;//三维直方图实际上是一个三维数组16*16*16q[index] = q[index] + 1.57*(1 - d*d);}}////////////////for(int m = 0; m < 4096; m++){q[m] = q[m]/C;}//for(int j =0; j < 4096; j++){if(q[j] != 0) w[j] = sqrt(p[j]/q[j]);//计算每个特征值的权值elsew[j] = 0;}//求新的中心cfloat x = 0;float y = 0;float div = 0;for(int  j = windowIn.y; j < windowIn.height + windowIn.y; j++){for(int k =windowIn.x; k < windowIn.width + windowIn.x; k++){int xi = k;int yi = j;int r = pdata[j * image->widthStep + k * 3 + 2];int g = pdata[j * image->widthStep + k * 3 + 1];int b = pdata[j * image->widthStep + k * 3 + 0];int index = (r/16)*256 +(g/16)*16 +b/16;x += (xi - c_center.x) * w[index];y += (yi - c_center.y) * w[index];div += w[index];}}x = x / div;y = y / div;n_center.x = c_center.x + x ;n_center.y = c_center.y + y;windowIn.x = n_center.x - width/2;windowIn.y = n_center.y - height/2;float error = sqrt(x*x+y*y);if(  error < 0.5){windowIn.x = n_center.x - width/2;windowIn.y = n_center.y - height/2;delete []w;return windowIn;}}delete []w;return windowIn;}/* End of file. */

main.cpp参照了opencv自带的camshift视频跟踪的例子

/*****************************************在RGB颜色空间计算的核概率密度*跟踪算法只是实现了传统的meanshift算法，比较粗糙*搜索区域的大小是固定的有待改进*带宽h也是固定的有待改进*bin的个数是16*16*16********************************************/#include "cv/cv.h"#include "cv/highgui.h"#include "cv/cxcore.h"#include "kernel_mean.h"#include <iostream>#include <windows.h>using namespace std;#pragma comment (lib, "cv/cv.lib")#pragma comment (lib, "cv/highgui.lib")#pragma comment (lib, "cv/cxcore.lib")CvRect selection;CvPoint origin;int select_object = 0;IplImage *image1 = 0;IplImage *image = 0;int track_object = 0;int initial_pCalc = 0;float p[4096] = {0};//初始的目标概率直方图float C =0;//概率密度直方图的归一化系数void on_mouse( int event, int x, int y, int flags, void* param );float h = 0;CvRect track_window ;void main(){int iname = 0;CvPoint rightup;IplImage *frame = 0;CvCapture* capture = 0;capture = cvCaptureFromCAM( 0 );//capture = cvCaptureFromAVI( "1.avi" );cvNamedWindow( "CamShiftDemo", 1 );cvSetMouseCallback( "CamShiftDemo", on_mouse, 0 );for(;;){frame = cvQueryFrame( capture );if( !image1){image1 = cvCreateImage(cvGetSize(frame), 8, 3);image = cvCreateImage(cvGetSize(frame), 8, 3);}cvCopy(frame, image1, 0);image1->origin = frame->origin;image->origin = frame->origin;cvSmooth(image1, image, CV_GAUSSIAN, 3, 0, 0, 0);//计算人为选取的目标模板的概率直方图if( track_object )//track_object为1表示目标选取完毕可以对目标模板进行计算{if( track_object < 0){h = sqrt(pow((double)selection.height, 2.0) + pow((double)selection.width, 2.0));h = h /2;//带宽hfor(int i =0; i < selection.height; i++)//求归一化系数C{for(int j = 0; j < selection.width; j++){int deltax=  j - selection.width/2;int deltay = i - selection.height/2;float d = sqrt(float(deltax*deltax) + float(deltay*deltay))/h;C += 1.57*(1 - d*d);}}//求目标模板的核概率密度punsigned char *pdata = (unsigned char *)image->imageData;for(int i =selection.y; i < selection.height+selection.y; i++){for(int j = selection.x; j < selection.width+selection.x; j++){int deltax=  j - selection.width/2-selection.x;int deltay = i - selection.height/2-selection.y;float d = sqrt(float(deltax*deltax) + float(deltay*deltay))/h;int r = pdata[i * image->widthStep + j * 3 + 2];int g = pdata[i * image->widthStep + j * 3 + 1];int b = pdata[i * image->widthStep + j * 3 + 0];int index = (r/16)*256 +(g/16)*16 +b/16;p[index] = p[index] + 1.57*(1 - d*d);}}for(int n = 0; n < 4096; n++){p[n] = p[n]/C;//对核概率密度进行归一化}track_window = selection;track_object = 1;}//long time1 = GetTickCount();CvRect object_window = kernel_meanshift(image, track_window,40, p, C, h);//long time2 = GetTickCount();//cout<<time2-time1<<endl;track_window = object_window;rightup.x = object_window.x + object_window.width -1;rightup.y = object_window.y + object_window.height -1;cvRectangle(image, cvPoint(object_window.x, object_window.y), rightup, cvScalar(0,0,255), 2, 8, 0);}if( select_object && selection.width > 0 && selection.height > 0 ){cvSetImageROI( image, selection );cvXorS( image, cvScalarAll(255), image, 0 );//矩阵与给定值进行异或操作cvResetImageROI( image );}cvShowImage("CamShiftDemo", image);int c = cvWaitKey(40);if( (char) c == 27 )break;}cvReleaseCapture( &capture );cvDestroyWindow("CamShiftDemo");}void on_mouse( int event, int x, int y, int flags, void* param ){if( !image )return;if( image->origin )y = image->height - y;if( select_object ){selection.x = MIN(x,origin.x);selection.y = MIN(y,origin.y);selection.width = selection.x + CV_IABS(x - origin.x);selection.height = selection.y + CV_IABS(y - origin.y);selection.x = MAX( selection.x, 0 );selection.y = MAX( selection.y, 0 );selection.width = MIN( selection.width, image->width );selection.height = MIN( selection.height, image->height );selection.width -= selection.x;selection.height -= selection.y;}switch( event ){case CV_EVENT_LBUTTONDOWN:origin = cvPoint(x,y);selection = cvRect(x,y,0,0);select_object = 1;break;case CV_EVENT_LBUTTONUP:select_object = 0;if( selection.width > 0 && selection.height > 0 )track_object = -1;break;}}