Camshift算法原理及其Opencv实现

Camshift算法原理及其Opencv实现
Camshift原理
camshift利用目标的颜色直方图模型将图像转换为颜色概率分布图，初始化一个搜索窗的大小和位置，并根据上一帧得到的结果自适应调整搜索窗口的位置和大小，从而定位出当前图像中目标的中心位置。

分为三个部分：
1--色彩投影图（反向投影）：
(1).RGB 颜色空间对光照亮度变化较为敏感，为了减少此变化对跟踪效果的影响，首先将图像从RGB空间转换到HSV空间。

(2).然后对其中的H分量作直方图，在直方图中代表了不同H分量值出现的概率或者像素个数，就是说可以查找出H分量大小为h的概率或者像素个数，即得到了颜色概率查找表。

(3).将图像中每个像素的值用其颜色出现的概率对替换，就得到了颜色概率分布图。

这个过程就叫反向投影，颜色概率分布图是一个灰度图像。

2--meanshift
meanshift算法是一种密度函数梯度估计的非参数方法，通过迭代寻优找到概率分布的极值来定位目标。

算法过程为：
(1).在颜色概率分布图中选取搜索窗W
(2).计算零阶距：
计算一阶距：
计算搜索窗的质心：
(3).调整搜索窗大小
宽度为；长度为1.2s；
(4).移动搜索窗的中心到质心，如果移动距离大于预设的固定阈值，则重复2)3)4)，直到搜索窗的中心与质心间的移动距离小于预设的固定阈值，或者循环运算的次数达到某一最大值，停止计算。

关于meanshift的收敛性证明可以google相关文献。

3--camshift
将meanshift算法扩展到连续图像序列，就是camshift算法。

它将视频的所有帧做meanshift运算，并将上一帧的结果，即搜索窗的大小和中心，作为下一帧meanshift算法搜索窗的初始值。

如此迭代下去，就可以实现对目标的跟踪。

算法过程为：
(1).初始化搜索窗
(2).计算搜索窗的颜色概率分布（反向投影）
(3).运行meanshift算法，获得搜索窗新的大小和位置。

(4).在下一帧视频图像中用(3)中的值重新初始化搜索窗的大小和位置，再跳转到(2)继续进行。

camshift 能有效解决目标变形和遮挡的问题，对系统资源要求不高，时间复杂度低，在简单背景下能够取得良好的跟踪效果。

但当背景较为复杂，或者有许多与目标颜色相似像素干扰的情况下，会导致跟踪失败。

因为它单纯的考虑颜色直方图，忽略了目标的空间分布特性，所以这种情况下需加入对跟踪目标的预测算法。

Camshift的opencv实现
原文/houdy/archive/2004/11/10/175739.aspx
1--Back Projection
计算Back Projection的OpenCV代码。

(1).准备一张只包含被跟踪目标的图片，将色彩空间转化到HSI空间，获得其中的H分量：
IplImage* target=cvLoadImage("target.bmp",-1); //装载图片
IplImage* target_hsv=cvCreateImage( cvGetSize(target), IPL_DEPTH_8U, 3 );
IplImage* target_hue=cvCreateImage( cvGetSize(target), IPL_DEPTH_8U, 3 );
cvCvtColor(target,target_hsv,CV_BGR2HSV); //转化到HSV空间cvSplit( target_hsv, target_hue, NULL, NULL, NULL ); //获得H分量
(2).计算H分量的直方图，即1D直方图：
IplImage*
h_plane=cvCreateImage( cvGetSize(target_hsv),IPL_DEPTH_8U,1 );
int hist_size[]={255}; //将H分量的值量化到[0,255]
float* ranges[]={ {0,360} }; //H分量的取值范围是[0,360) CvHistogram* hist=cvCreateHist(1, hist_size, ranges, 1);
cvCalcHist(&target_hue, hist, 0, NULL);
在这里需要考虑H分量的取值范围的问题，H分量的取值范围是[0,360)，这个取值范围的值不能用一个byte来表示，为了能用一个byte表示，需要将H值做适当的量化处理，在这里我们将H分量的范围量化到[0,255]。

(3).计算Back Projection：
IplImage* rawImage;
//get from video frame,unsigned byte,one channel
IplImage* result=cvCreateImage(cvGetSize(rawImage),IPL_DEPTH_8U,1); cvCalcBackProject(&rawImage,result,hist);
(4). result即为我们需要的.
2--Mean Shift算法
质心可以通过以下公式来计算：
(1).计算区域内0阶矩
for(int i=0;i< height;i++)
for(int j=0;j< width;j++)
M00+=I(i,j)
(2).区域内1阶矩：
for(int i=0;i< height;i++)
for(int j=0;j< width;j++)
{
M10+=i*I(i,j);
M01+=j*I(i,j);
}
(3).则Mass Center为：
Xc=M10/M00; Yc=M01/M00
在OpenCV中，提供Mean Shift算法的函数，函数的原型是：
int cvMeanShift(IplImage* imgprob,CvRect windowIn,
CvTermCriteria criteria,CvConnectedComp* out);
需要的参数为：
(1).IplImage* imgprob：2D概率分布图像，传入；
(2).CvRect windowIn：初始的窗口，传入；
(3).CvTermCriteria criteria：停止迭代的标准，传入；
(4).CvConnectedComp* out:查询结果，传出。

注：构造CvTermCriteria变量需要三个参数，一个是类型，另一个是迭代的最大次数，最后一个表示特定的阈值。

例如可以这样构造 criteria：
criteria=cvTermCriteria(CV_TERMCRIT_ITER|CV_TERMCRIT_EPS,10,0.1)。

3--CamShift算法
整个算法的具体步骤分5步：
Step 1：将整个图像设为搜寻区域。

Step 2：初始话Search Window的大小和位置。

Step 3：计算Search Window内的彩色概率分布，此区域的大小比Search Window 要稍微大一点。

Step 4：运行MeanShift。

获得Search Window新的位置和大小。

Step 5：在下一帧视频图像中，用Step 3获得的值初始化Search Window的位置和大小。

跳转到Step 3继续运行。

OpenCV代码：
在OpenCV中，有实现CamShift算法的函数，此函数的原型是：
cvCamShift(IplImage* imgprob, CvRect windowIn,
CvTermCriteria criteria,
CvConnectedComp* out, CvBox2D* box=0);
其中：
imgprob：色彩概率分布图像。

windowIn：Search Window的初始值。

Criteria：用来判断搜寻是否停止的一个标准。

out：保存运算结果,包括新的Search Window的位置和面积。

box：包含被跟踪物体的最小矩形
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ opencv中自带camshift代码
2011-12-02 14:48:23| 分类：opencv |举报 |字号订阅
cvCamShift( const void* imgProb, CvRect windowIn,
CvTermCriteria criteria,
CvConnectedComp* _comp,
CvBox2D* box )
{
const int TOLERANCE = 10; // 容差，计算色度直方图的窗口比搜索窗口稍大CvMoments moments; // 矩
double m00 = 0, m10, m01, mu20, mu11, mu02, inv_m00; // 几个特征矩double a, b, c, xc, yc;
double rotate_a, rotate_c; // 旋转
double theta = 0, square; // 角度，
double cs, sn;
double length = 0, width = 0;
int itersUsed = 0; //
CvConnectedComp comp;
CvMat cur_win, stub, *mat = (CvMat*)imgProb;
CV_FUNCNAME( "cvCamShift" );
comp.rect = windowIn; // 连接部件矩形区域初始化-输入窗口
__BEGIN__;
CV_CALL( mat = cvGetMat( mat, &stub ));
CV_CALL( itersUsed = cvMeanShift( mat, windowIn, criteria, &comp )); // 执行cvMeanShift，返回迭代次数
windowIn = comp.rect; //收敛的搜索窗口坐标 (comp->rect 字段)
//设置计算色度直方图窗口，它的大小是在搜索窗口的基础上四边加上TOLERANCE，同时还要保证该窗口在图像大小范围内
windowIn.x -= TOLERANCE;
if( windowIn.x < 0 )
windowIn.x = 0;
windowIn.y -= TOLERANCE;
if( windowIn.y < 0 )
windowIn.y = 0;
windowIn.width += 2 * TOLERANCE;
if( windowIn.x + windowIn.width > mat->width )
windowIn.width = mat->width - windowIn.x;
windowIn.height += 2 * TOLERANCE;
if( windowIn.y + windowIn.height > mat->height )
windowIn.height = mat->height - windowIn.y;
CV_CALL( cvGetSubRect( mat, &cur_win, windowIn ));//取 mat子集
/* Calculating moments in new center mass */
CV_CALL( cvMoments( &cur_win, &moments ));
//空间矩
m00 = moments.m00;
m10 = moments.m10;
m01 = moments.m01;
//中心距
mu11 = moments.mu11; //mu11=m11-m10*(m10/m00),ok
mu20 = moments.mu20; //mu20=m20-m10*(m10/m00),ok
mu02 = moments.mu02; //mu02=m02-m01*(m01/m00),ok
if( fabs(m00) < DBL_EPSILON )
EXIT;
inv_m00 = 1. / m00;
//(m10 * inv_m00,m01 * inv_m00)是窗口质心坐标
xc = cvRound( m10 * inv_m00 + windowIn.x ); // 图像质心坐标，x分量
yc = cvRound( m01 * inv_m00 + windowIn.y ); // 图像质心坐标，y分量
a = mu20 * inv_m00;
b = mu11 * inv_m00;
c = mu02 * inv_m00;
/* Calculating width & height *///=====================================================从这里
square = sqrt( 4 * b * b + (a - c) * (a - c) );
/* Calculating orientation */
theta = atan2( 2 * b, a - c + square );
/* Calculating width & length of figure */
cs = cos( theta );
sn = sin( theta );
rotate_a = cs * cs * mu20 + 2 * cs * sn * mu11 + sn * sn * mu02; rotate_c = sn * sn * mu20 - 2 * cs * sn * mu11 + cs * cs * mu02; length = sqrt( rotate_a * inv_m00 ) * 4;
width = sqrt( rotate_c * inv_m00 ) * 4;//================================================到这里
/* In case, when tetta is 0 or 1.57... the Length & Width may be exchanged */
if( length < width )
{
double t;
CV_SWAP( length, width, t );
CV_SWAP( cs, sn, t );
theta = CV_PI*0.5 - theta;
}
/* Saving results */
if( _comp || box )
{
int t0, t1;
int _xc = cvRound( xc ); //坐标取整
int _yc = cvRound( yc );
t0 = cvRound( fabs( length * cs ));
t1 = cvRound( fabs( width * sn ));
t0 = MAX( t0, t1 ) + 2;
comp.rect.width = MIN( t0, (mat->width - _xc) * 2 );
t0 = cvRound( fabs( length * sn ));
t1 = cvRound( fabs( width * cs ));
t0 = MAX( t0, t1 ) + 2;
comp.rect.height = MIN( t0, (mat->height - _yc) * 2 );
// comp.rect顶点和高宽的设置，保证其不超出图像范围
comp.rect.x = MAX( 0, _xc - comp.rect.width / 2 );
comp.rect.y = MAX( 0, _yc - comp.rect.height / 2 );
comp.rect.width = MIN( mat->width - comp.rect.x, comp.rect.width ); comp.rect.height = MIN( mat->height - comp.rect.y, comp.rect.height ); // end 设置
comp.area = (float) m00; //
}
__END__;
if( _comp )
*_comp = comp;
if( box )
{
box->size.height = (float)length; //目标边带的高
box->size.width = (float)width; //目标边带的宽
box->angle = (float)(theta*180./CV_PI); //方向
box->center = cvPoint2D32f( comp.rect.x + comp.rect.width*0.5f, comp.rect.y + comp.rect.height*0.5f);// box的中心就是收敛窗口的中心
}
return itersUsed; // 返回迭代次数
}。