MeanShift知识整理

第一个问题：无参数密度估计

给定任意一组观测数据或数据采样值，估计出样本的分布。

无参数密度估计，它对数据分布规律没有附加任何假设，而是直接从数据样本本身出发研究数据分布特征，对先验知识要求最少，完全依靠训练数据进行估计，而且能处理任意的概率分布。

例如：直方图法，最近邻域法，核密度估计方法。

有参数密度估计：高斯统计模型。

更形象的说明：

已经有N个点，它们的坐标分布如下面所示，如何求出这个区域中，哪个位置的样本分布密度最大，换句话说，如果有第N+1个点，它最大的可能出现在哪个位置。

下面我们就用Mean Shift的原理来解决这类的问题。

第二个问题：核函数的密度估计

给定维空间中的样本集合，则点关于核函数和带宽矩阵的核函数密度估计表示为：

其中

其中表示带宽矩阵。我们用比例单位矩阵来表示，将上面的核密度方程写成下面的典型的表达式：

由于，可将密度估计写成关于核函数的轮廓函数形式：

由上面的式子可以看出，核函数是一种权值函数，其作用是将每个样本点按到中心点距离的远近进行加权，距离中心点近的样本点概率密度估计影响大，赋予大的权值，反之，赋予小的权值。

第三个问题：Mean Shift 向量

我们在第二个问题中，已经由核函数来估计出来了一个样本集合的概率密度，现在我们要根据这个概率密度方程，来分析数据集合中密度最大数据的分布位置，我们可以对密度函数求导，分析梯度变化。

令上面的导数，则可以得到

也就是说上面的点就是概率密度最大的点。

第二个中括号内的内容就是Mean Shift向量，表示为：

为了更直观的说明上面公式的意义，我们把核函数，这样上面的式子就可以写为：

中间的实心黑点表示x点，也是核函数的中心。周围的空心白点是样本点。箭头表示样本点相对于核函数中心点的偏移向量，平均的偏移量会指向样本点最密的方向，也就是梯度方向。因此，Mean Shift向量应该转移到样本点相对点变化最多的地方，其方向也就是密度梯度的方向。

但是上面中所有样本点对的贡献是一样的。一般而言，离越近的采样点对估计周围的统计特性越重要，因此引入核函数的概念，式中就是对每个采样点的权值。

第四个问题：基于颜色的Mean Shift目标跟踪算法

4.1 直方图

灰度直方图

直方图反映了图像中每一种颜色出现的次数。

那么彩色直方图的话就是每一种（R,G,B）这样的组合出现的次数。

灰度图像直方图横坐标0-255

而彩色直方图，共有256*256*256种色彩，就是直方图的横坐标0-16777215

我们通过降维来简化彩色颜色直方图

我们把R，G，B三个通道0-255区间内的值都分为16个区间，即将0-255转换为0-15

比如（35，122，200）按上面将转换为（2，7，12），这样整个颜色空间就可以用16*16*16种组合来表示。

将三维向量映射为一维的

（R，G，B）== 256*R+16*G+B

4.2 目标模型

设是目标模板区域的像素位置，且目标区域的中心为，定义函数为像素处的颜色索引函数，为处像素点有相应直方图中颜色索引值。则目标模型中第个特征值估计的概率密度为：

其中，为目标特征值。为核函数的轮廓函数，是一种加权函数，其作用是给目标区域中的像素点设置权值，对靠近目标模板中心像素赋予较大的权值，而对于远离目标模板中心的像素赋予较小的权值。目标区域中心的权值最大，区域边缘的权值小，对于远离目标模板中心的像素容易受到边缘噪声，干扰物的影响，加权后进行密度估计时增加了它们的鲁棒性，提高了跟踪的抗干扰能力。

由于，即所有特征值的概率和为1。因为可以推导出：

颜色索引函数的值为：

4.3 候选目标模型的表示

同理，设（n为候选区域中像素点的个数）是候选目标区域的像素位置，在当前帧以y为中心，候选模型的特征值为的颜色概率估计密度为：

其中归一化常数为：

4.4 相似性度量

选用Bhattacharrya系数相似性函数来表示目标模型和目标候选模型的相似程度。

用Bhattacharrya系数表示和的相似程度如下式：

和越相似，的值越大，在理想情况下和的概率分布是完全相同的，值为1，否则的值介于0~1之间

使用上式，两个模型的距离就可以表示为：

4.5 目标定位

当有了目标和候选目标的模型以及度量它们似度的准则后，目标跟踪的任务就是在当前帧中找到目标的新位置。

我们可以通过最小化的距离值来获得目标当前帧的最可能的位置。

在当前帧对寻找目标的新位置是根据前一帧目标位置的估计开始的，并在周围邻域内寻找。于是，必须先计算出在当前帧中位于处的候选目标的颜色概率，对在处进行泰勒展开，可以得到下面的式子：

其中，

我们只分析后面的第二项，

可以看出第二项就是在当前帧中利用核函数和图像像素加权值计算得到的概率密度估计。

根据Mean Shift向量，我们就可以得到候选区域中心移向真实目标区域的向量

4.6 跟踪目标窗口自适应大小

采用的带宽来计算三次候选目标的颜色直方图，然后按下面的公式分别计算相似性：

最其中最小的带宽作为目标窗口的大小

4.7 代码实现

说明：

1.       RGB颜色空间刨分，采用16*16*16的直方图

2.       目标模型和候选模型的概率密度计算公式参照上文

3.       opencv版本运行：按P停止，截取目标，再按P，进行单目标跟踪

4.       Matlab版本，将视频改为图片序列，第一帧停止，手工标定目标，双击目标区域，进行单目标跟踪。

opencv版本:

  
#include "cv.h"  
#include "highgui.h"  
#define  u_char unsigned char  
#define  DIST 0.5  
#define  NUM 20  

//全局变量  
bool pause = false;
bool is_tracking = false;
CvRect drawing_box;
IplImage *current;
double *hist1, *hist2;
double *m_wei;                                                                  //权值矩阵  
double C = 0.0;                                                                //归一化系数  

void init_target(double *hist1, double *m_wei, IplImage *current)
{
	IplImage *pic_hist = 0;
	int t_h, t_w, t_x, t_y;
	double h, dist;
	int i, j;
	int q_r, q_g, q_b, q_temp;

	t_h = drawing_box.height;
	t_w = drawing_box.width;
	t_x = drawing_box.x;
	t_y = drawing_box.y;

	h = pow(((double)t_w) / 2, 2) + pow(((double)t_h) / 2, 2);            //带宽  
	pic_hist = cvCreateImage(cvSize(300, 200), IPL_DEPTH_8U, 3);     //生成直方图图像  

	//初始化权值矩阵和目标直方图  
	for (i = 0; i < t_w*t_h; i++)
	{
		m_wei[i] = 0.0;
	}

	for (i = 0; i<4096; i++)
	{
		hist1[i] = 0.0;
	}

	for (i = 0; i < t_h; i++)
	{
		for (j = 0; j < t_w; j++)
		{
			dist = pow(i - (double)t_h / 2, 2) + pow(j - (double)t_w / 2, 2);
			m_wei[i * t_w + j] = 1 - dist / h;
			//printf("%f\n",m_wei[i * t_w + j]);  
			C += m_wei[i * t_w + j];
		}
	}

	//计算目标权值直方  
	for (i = t_y; i < t_y + t_h; i++)
	{
		for (j = t_x; j < t_x + t_w; j++)
		{
			//rgb颜色空间量化为16*16*16 bins  
			q_r = ((u_char)current->imageData[i * current->widthStep + j * 3 + 2]) / 16;
			q_g = ((u_char)current->imageData[i * current->widthStep + j * 3 + 1]) / 16;
			q_b = ((u_char)current->imageData[i * current->widthStep + j * 3 + 0]) / 16;
			q_temp = q_r * 256 + q_g * 16 + q_b;
			hist1[q_temp] = hist1[q_temp] + m_wei[(i - t_y) * t_w + (j - t_x)];
		}
	}

	//归一化直方图  
	for (i = 0; i<4096; i++)
	{
		hist1[i] = hist1[i] / C;
		//printf("%f\n",hist1[i]);  
	}

	//生成目标直方图  
	double temp_max = 0.0;

	for (i = 0; i < 4096; i++)         //求直方图最大值，为了归一化  
	{
		//printf("%f\n",val_hist[i]);  
		if (temp_max < hist1[i])
		{
			temp_max = hist1[i];
		}
	}
	//画直方图  
	CvPoint p1, p2;
	double bin_width = (double)pic_hist->width / 4096;
	double bin_unith = (double)pic_hist->height / temp_max;

	for (i = 0; i < 4096; i++)
	{
		p1.x = i * bin_width;
		p1.y = pic_hist->height;
		p2.x = (i + 1)*bin_width;
		p2.y = pic_hist->height - hist1[i] * bin_unith;
		//printf("%d,%d,%d,%d\n",p1.x,p1.y,p2.x,p2.y);  
		cvRectangle(pic_hist, p1, p2, cvScalar(0, 255, 0), -1, 8, 0);
	}
	cvSaveImage("hist1.jpg", pic_hist);
	cvReleaseImage(&pic_hist);
}

void MeanShift_Tracking(IplImage *current)
{
	int num = 0, i = 0, j = 0;
	int t_w = 0, t_h = 0, t_x = 0, t_y = 0;
	double *w = 0, *hist2 = 0;
	double sum_w = 0, x1 = 0, x2 = 0, y1 = 2.0, y2 = 2.0;
	int q_r, q_g, q_b;
	int *q_temp;
	IplImage *pic_hist = 0;

	t_w = drawing_box.width;
	t_h = drawing_box.height;

	pic_hist = cvCreateImage(cvSize(300, 200), IPL_DEPTH_8U, 3);     //生成直方图图像  
	hist2 = (double *)malloc(sizeof(double)* 4096);
	w = (double *)malloc(sizeof(double)* 4096);
	q_temp = (int *)malloc(sizeof(int)*t_w*t_h);

	while ((pow(y2, 2) + pow(y1, 2) > 0.5) && (num < NUM))
	{
		num++;
		t_x = drawing_box.x;
		t_y = drawing_box.y;
		memset(q_temp, 0, sizeof(int)*t_w*t_h);
		for (i = 0; i<4096; i++)
		{
			w[i] = 0.0;
			hist2[i] = 0.0;
		}

		for (i = t_y; i < t_h + t_y; i++)
		{
			for (j = t_x; j < t_w + t_x; j++)
			{
				//rgb颜色空间量化为16*16*16 bins  
 				q_r = ((u_char)current->imageData[i * current->widthStep + j * 3 + 2]) / 16;
				q_g = ((u_char)current->imageData[i * current->widthStep + j * 3 + 1]) / 16;
				q_b = ((u_char)current->imageData[i * current->widthStep + j * 3 + 0]) / 16;
				q_temp[(i - t_y) *t_w + j - t_x] = q_r * 256 + q_g * 16 + q_b;
				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];
			}
		}

		//归一化直方图  
		for (i = 0; i<4096; i++)
		{
			hist2[i] = hist2[i] / C;
			//printf("%f\n",hist2[i]);  
		}
		//生成目标直方图  
		double temp_max = 0.0;

		for (i = 0; i<4096; i++)         //求直方图最大值，为了归一化  
		{
			if (temp_max < hist2[i])
			{
				temp_max = hist2[i];
			}
		}
		//画直方图  
		CvPoint p1, p2;
		double bin_width = (double)pic_hist->width / (4368);
		double bin_unith = (double)pic_hist->height / temp_max;

		for (i = 0; i < 4096; i++)
		{
			p1.x = i * bin_width;
			p1.y = pic_hist->height;
			p2.x = (i + 1)*bin_width;
			p2.y = pic_hist->height - hist2[i] * bin_unith;
			cvRectangle(pic_hist, p1, p2, cvScalar(0, 255, 0), -1, 8, 0);
		}
		cvSaveImage("hist2.jpg", pic_hist);

		for (i = 0; i < 4096; i++)
		{
			if (hist2[i] != 0)
			{
				w[i] = sqrt(hist1[i] / hist2[i]);
			}
			else
			{
				w[i] = 0;
			}
		}

		sum_w = 0.0;
		x1 = 0.0;
		x2 = 0.0;

		for (i = 0; i < t_h; i++)
		{
			for (j = 0; j < t_w; j++)
			{
				//printf("%d\n",q_temp[i * t_w + j]);  
				sum_w = sum_w + w[q_temp[i * t_w + j]];
				x1 = x1 + w[q_temp[i * t_w + j]] * (i - t_h / 2);
				x2 = x2 + w[q_temp[i * t_w + j]] * (j - t_w / 2);
			}
		}
		y1 = x1 / sum_w;
		y2 = x2 / sum_w;

		//中心点位置更新  
		drawing_box.x += y2;
		drawing_box.y += y1;

		//printf("%d,%d\n",drawing_box.x,drawing_box.y);  
	}
	free(hist2);
	free(w);
	free(q_temp);
	//显示跟踪结果  
	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);
	cvShowImage("Meanshift", current);
	//cvSaveImage("result.jpg",current);  
	cvReleaseImage(&pic_hist);
}

void onMouse(int event, int x, int y, int flags, void *param)
{
	if (pause)
	{
		switch (event)
		{
		case CV_EVENT_LBUTTONDOWN:
			//the left up point of the rect  
			drawing_box.x = x;
			drawing_box.y = y;
			break;
		case CV_EVENT_LBUTTONUP:
			//finish drawing the rect (use color green for finish)  
			drawing_box.width = x - drawing_box.x;
			drawing_box.height = y - drawing_box.y;
			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);
			cvShowImage("Meanshift", current);

			//目标初始化  
			hist1 = (double *)malloc(sizeof(double)* 16 * 16 * 16);
			m_wei = (double *)malloc(sizeof(double)*drawing_box.height*drawing_box.width);
			init_target(hist1, m_wei, current);
			is_tracking = true;
			break;
		}
		return;
	}
}



int main(int argc, char* argv[])
{
	//******************usb camera**********************//
	//CvCapture* capture = cvCreateCameraCapture(0);
	//current = cvQueryFrame(capture);



	//******************video**********************//
	CvCapture *capture = cvCreateFileCapture("car.mpg");
	current = cvQueryFrame(capture);


	char res[20];
	int nframe = 0;

	while (1)
	{
		/*  sprintf(res,"result%d.jpg",nframe);
		cvSaveImage(res,current);
		nframe++;*/
		if (is_tracking)
		{
			MeanShift_Tracking(current);
		}

		int c = cvWaitKey(1);
		//暂停  
		if (c == 'p')
		{
			pause = true;
			cvSetMouseCallback("Meanshift", onMouse, 0);
		}
		while (pause){
			if (cvWaitKey(0) == 'p')
				pause = false;
		}
		cvShowImage("Meanshift", current);
		current = cvQueryFrame(capture); //抓取一帧  
	}

	cvNamedWindow("Meanshift", 1);
	cvReleaseCapture(&capture);
	cvDestroyWindow("Meanshift");
	return 0;
}

运行结果：

 

初始目标直方图：

候选目标直方图：

源码及素材的下载地址：http://download.csdn.net/detail/jinshengtao/7489501