MeanShift 均值漂移跟踪算法原理

概述

整个暑假都在做Tracking，其中最为重要的核心就是均值漂移了，均值漂移算法指的是一个迭代的步骤，即先算出当前点的漂移均值，移动该点到其漂移均值，然后以此为新的起始点，继续移动，知道满足一定条件结束。所以，均值漂移实际上是一种在一种数据的密度分布中寻找局部极值的稳定的方法。如果分布连续，那么处理则变得非常容易，在这种情况下本质上只需要对数据的密度直方图应用爬山算法。更加准确的说，均值漂移是与核密度估计的规则有关的算法。而所谓“核”，实际上是一个如同高斯分布的局部函数。如果在充足的点处拥有足够合适的带权重和尺度的核，数据的分布则可以完全根据这些核来表示。然而又与核密度的估计不同，均值漂移仅仅估计数据分布的梯度。如果变化为0的地方则表示是这个分布的峰值（当然，也有可能是局部的）。当然在附近或其他尺度上还是有可能有峰值的。

定义

给定$d$维空间$R^{d}$中$n$个样本点$x_{i},i=1,2,…,n$，$n$在$x$点的均值漂移向量的基本形式定义为：

$$ \begin{equation} M_{h}(x)=\frac{1}{k}\sum_{x_{i}\in S_{k}}\left(x_{i}-x\right) \end{equation} $$

目标模型

均值漂移采用的是特征值的加权概率分布来描述目标模型，属于模式识别中主要描述目标的模型，不同于自动控制理论中采用的状态方程。目标模型一共具有$m$个特征值（可以理解为像素的灰度值），于是对于序列$q={q_{n}}$，而$u\in {1,…,m}$有

$$ \begin{equation} q\left(u\right)=C\sum_{i=1}^{n}k\left( \left\Vert \frac{X_{i}-X_{0}}{H} \right\Vert ^{2}\right) \end{equation} $$

其中，$X_{0}$为窗口中心点向量值（可能为RGB向量或者灰度值），$X_{i}$是窗口内第$i$点向量值。$C$为归一化常数，保障$\sum_{i=1}^{m}q_{i}=1$ ，$H$为核函数的带宽向量。$M$为特征值的个数，对于图像处理可以理解为灰度等级划分的个数，从而特征值$u$为对应的灰度等级。 $\delta$函数为脉冲函数，保证只具有$u$特征值的像素才对概率分布作出贡献。从而$k$函数可以理解为$u$灰度值的一个加权频数。

匹配对象

同样采用的是特征值加权概率分布：

$$ \begin{equation} P_{u}(Y)=C_{h} \sum_{i=1}^{n_{k}}k\left( \left\Vert \frac{X_{i}-Y}{H_{h}} \right\Vert ^{2}\right) \delta\left(b(X_{i})-u\right) \end{equation} $$

其中$i\in [1,…,n_{h}]$，$Y$为匹配对象的中心， $X_{i}$是匹配窗口内第$i$点向量值。 $H_{h}$为匹配窗口的核函数带宽向量，$C_{h}$为匹配窗口特征向量的归一化常数。

匹配相似

匹配对象与目标模型的相似程度，相似函数采用的是Bhattacharyya函数

$$ \begin{equation} \rho\left(p(Y),q\right)=\sum_{u=1}^{m}\sqrt{P_{u}(Y)q_{u}} \end{equation} $$

匹配过程

均值漂移采用梯度下降法，首先$$\rho(Y)$$在\rho(Y_{0})附近进行泰勒展开，去前两项，得到

$$ \begin{equation} \rho(Y) \approx \rho(Y_{0})+\frac{d\rho}{dp}\left(p(Y)-p(Y_{0})\right) \end{equation} $$

定义

$$ \begin{equation} \rho_{u}(Y)=\sqrt{p_{u}(Y)q_{u}} \end{equation} $$

从而

$ \begin{equation} \rho_{u}(Y)=\rho_{u}(Y_{0})+\frac{q_{u}}{ 2 \sqrt{p_{u}(Y_{0})q_{u}} }\left( p_{u}(Y) - p_{u}(Y_{0}) \right) \end{equation} $

即

$$ \begin{equation} \rho_{u}(Y)= \frac{1}{2}\rho_{u}(Y_{0}) + \frac{1}{2}\sqrt{\frac{p_{u}(Y)q_{u}}{p_{u}(Y_{0})}} \end{equation} $$

因此

$$ \begin{equation} \rho(Y)=\frac{1}{2}\left(\sum_{u=1}^{m}\sqrt{p_{u}(Y_{0})q_{u}}+\sum_{u=1}^{m}\left(p_{u}(Y_0)\sqrt{\frac{q_{u}}{p_{u}(Y_{0})}}\right) \right) \end{equation} $$

要使得$\rho(Y)$向最大值迭代，只要$Y$的搜索方向与梯度方向一致即可，通过求导得到$Y_{0}$的梯度方向为：

$$ \begin{equation} \nabla\rho(Y_{0}) = \frac{C_{h}}{H_{h}^{2}} \left[ \sum_{i=1}^{n_{k}}w_{i}g\left( \left\Vert \frac{Y_{0}-X_{i}}{H_{h}} \right\Vert ^{2} \right) \right] \left[ \frac{\sum_{i=1}^{n_{k}}X_{i}w_{i}g\left( \left\Vert \frac{Y_{0}-X_{i}}{H_{h}} \right\Vert ^{2} \right)}{\sum_{i=1}^{n_{k}}w_{i}g\left( \left\Vert \frac{Y_{0}-X_{i}}{H_{h}} \right\Vert ^{2} \right)} - Y_{0} \right] \end{equation} $$

其中，

$$ \begin{equation} w_{i} = \sum_{u=1}{m}\sqrt{\frac{q_{u}}{p_{u}(Y_{0})}}\delta(b(X_{i})-u) \end{equation} $$

为权值。设

$$ \begin{equation} Y_{1} = \frac{\sum_{i=1}^{n_{k}}X_{i}w_{i}g\left( \left\Vert \frac{Y_{0}-X_{i}}{H_{h}} \right\Vert ^{2} \right)}{\sum_{i=1}^{n_{k}}w_{i}g\left( \left\Vert \frac{Y_{0}-X_{i}}{H_{h}} \right\Vert ^{2} \right)} \end{equation} $$

因此如果如下确定$Y_{1}$，那么$Y_{1}-Y_{0}$（此即为均值漂移向量）将于梯度方向一致。

总结

均值漂移跟踪算法的思路很清晰，我们从一个核

$$ \begin{equation} K(X-X_{i}) = ck\left( \left\Vert \frac{X-X_{i}}{H} \right\Vert ^2 \right) \end{equation} $$

和一个近似概率分布

$$ \begin{equation} P(x) = \frac{1}{n} \sum_{i=1}^{n} K(x-x_{i}) \end{equation} $$

出发，重点关心$P(x)$的梯度

$$ \begin{equation} \nabla P(x) = \frac{1}{n}\nabla K(x-x_{i}) \end{equation} $$

令$g(x)=-k'(x)$，并得出了这样的式子：

$$ \begin{equation} \nabla P(x) = \frac{c}{n} \left[ \sum_{i=1}^{n}g_{i}\left( \left\Vert \frac{x-x_{i}}{H} \right\Vert ^{2} \right) \right] \left[ \frac{\sum_{i=1}^{n}x_{i}g_{i}\left( \left\Vert \frac{x-x_{i}}{H} \right\Vert ^{2} \right)}{\sum_{i=1}^{n}g_{i}\left( \left\Vert \frac{x-x_{i}}{H} \right\Vert ^{2} \right)} - x \right] \end{equation} $$

对于$x$向量通过上面的公式得到新的向量为

$$ \begin{equation} \left[ \frac{\sum_{i=1}^{n}x_{i}g_{i}\left( \left\Vert \frac{x-x_{i}}{H} \right\Vert ^{2} \right)}{\sum_{i=1}^{n}g_{i}\left( \left\Vert \frac{x-x_{i}}{H} \right\Vert ^{2} \right)} - x \right] \end{equation} $$

而$H$则是对于以$x$为心的半径。矩形核并不随着到中心的距离下降，而是一个突然变成零的突然转换。这个与高斯核的指数衰减不同，与Epanechnikov核的随着到中心的距离的开放衰减也不同。我们用自然语言来描述整个过程：

选择搜索域（包括域的初始位置，域的类型[均匀、多项式、指数或者高斯类型]，域的形状，域的大小）。
计算域（可能是带权重的）的重心。
将域的重心设置在计算出的重心处。
返回（2），直到域的位置不再变化（通常会，迭代过程由最大迭代次数或者两次迭代中心变化的程度进行限制。虽然如此，迭代过程最后还是会收敛。）

实现

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92


int cvMeanShift( const void* imgProb, CvRect windowIn, CvTermCriteria criteria, CvConnectedComp* comp )
{
	// CvMoments计算矩形的形状特征（重心、面积等）
	CvMoments moments;
	int    i = 0, eps;
	CvMat  stub, *mat = (CvMat*)imgProb;
	CvMat  cur_win;
	CvRect cur_rect = windowIn;

	// 初始化跟踪窗口
	if( comp )
		comp-&gt;rect = windowIn;

	// 把0阶矩和1阶矩初始化置0
	moments.m00 = moments.m10 = moments.m01 = 0;

	mat = cvGetMat( mat, &amp;stub );

	if( CV_MAT_CN( mat-&gt;type ) &gt; 1 )
		CV_Error( CV_BadNumChannels, cvUnsupportedFormat );

	if( windowIn.height &lt;= 0 || windowIn.width &lt;= 0 )
 		CV_Error( CV_StsBadArg, &quot;输入窗口大小错误&quot; );
 	windowIn = cv::Rect(windowIn) &amp; cv::Rect(0, 0, mat-&gt;cols, mat-&gt;rows);

	// 迭代的标准，精度为1.0，迭代次数为100。
	criteria = cvCheckTermCriteria( criteria, 1., 100 );
	// 精度eps = 1
	eps = cvRound( criteria.epsilon * criteria.epsilon );

	// 最大循环次数 = 最大迭代次数criteria.max_iter = 100
	for( i = 0; i &lt; criteria.max_iter; i++ )
 	{
 		int dx, dy, nx, ny;
 		double inv_m00;
 		cur_rect = cv::Rect(cur_rect) &amp; cv::Rect(0, 0, mat-&gt;cols, mat-&gt;rows);
		if( cv::Rect(cur_rect) == cv::Rect() )
		{
			cur_rect.x = mat-&gt;cols/2;
			cur_rect.y = mat-&gt;rows/2;
		}
		cur_rect.width = MAX(cur_rect.width, 1);
		cur_rect.height = MAX(cur_rect.height, 1);

		// 选取搜索区域，对改矩形区域计算它的0,1阶矩
		cvGetSubRect( mat, &amp;cur_win, cur_rect );
		cvMoments( &amp;cur_win, &amp;moments );

		// 计算质心
		if( fabs(moments.m00) &lt; DBL_EPSILON )
			break;

		// 搜索区域的质量m00
		inv_m00 = moments.inv_sqrt_m00*moments.inv_sqrt_m00;
		// 搜索区域的水平重心偏移量为dx
		dx = cvRound( moments.m10 * inv_m00 - windowIn.width*0.5 );
		// 搜索区域的垂直重心偏移量为dy
		dy = cvRound( moments.m01 * inv_m00 - windowIn.height*0.5 );

		// 搜索区域的重心坐标(nx,ny)
		nx = cur_rect.x + dx;
		ny = cur_rect.y + dy;

		// 跟踪目标处于图像边缘情况下的处理
		if( nx &lt; 0 )
 			nx = 0;
 		else if( nx + cur_rect.width &gt; mat-&gt;cols )
			nx = mat-&gt;cols - cur_rect.width;

		if( ny &lt; 0 )
 			ny = 0;
 		else if( ny + cur_rect.height &gt; mat-&gt;rows )
			ny = mat-&gt;rows - cur_rect.height;

		dx = nx - cur_rect.x;
		dy = ny - cur_rect.y;
		cur_rect.x = nx;
		cur_rect.y = ny;

		// 覆盖中心质量和窗口的检查，当精度达到要求则跳出循环
		if( dx*dx + dy*dy &lt; eps )
 			break;
 	}
 	// 返回值
 	if( comp )
 	{
 		comp-&gt;rect = cur_rect;
		comp-&gt;area = (float)moments.m00;
	}

	return i;
}