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

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.14159
struct 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*16
q[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]);//计算每个特征值的权值
else
w[j] = 0;
}

//求新的中心c
float 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;//带宽h
for(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);

}
}

//求目标模板的核概率密度p
unsigned 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
= p
/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;
}
}