【计算机视觉】目标跟踪| Meanshift均值漂移算法详细介绍|附代码

0、前言

在上篇文章【计算机视觉】目标跟踪| 光流算法详细介绍|附代码中，主要介绍了生成式目标追踪算法之一：光流法（optical flow)。在这篇文章中我们将介绍生成式算法的第二种算法：Meanshift均值漂移算法.

一、介绍&基本原理

背景： Meanshift算法是Fukunaga于1975年提出的，其基本思想是利用概率密度的梯度爬升来寻找局部最优。到了1995年，YizongCheng针对离x越近的采样点对x周围的统计特性越有效，定义了一族核函数，并根据所有样本点的重要性不同，设定了一个权重系数，扩大了MeanShift的使用范围。

啥是概率密度？？？

Wiki：“在数学中，连续型随机变量的概率密度函数（Probability density function，简写作PDF [1]），在不致于混淆时可简称为密度函数，是一个描述这个随机变量的输出值，在某个确定的取值点附近的可能性的函数。”

对于目标追踪算法简单来说呢，可以分开说

• 概率： 就是指这个目标在下一帧中出现的概率大小呗
• 密度： 这个区域可能是目标位置的概率大小呗，如果这片区域很有可能为目标位置，那当然这片区域的概率密度就很大。

⭐由下图可以得到以下信息：

• Meanshift：A tool for : Finding models in a set of data samples, manifesting an underlying probability density function(PDF) in $R^n$:
  1. 首先对目标对象进行特征建模得到model
  2. 在PDF ，也就是表示目标模型出现概率的函数处理的图片数据上进行PDF处理（也就是上面所说的利用概率密度的梯度爬升，找到目标出现可能性最大的区域）
  3. 使用MeanShit均值漂移算法，使得向目标出现可能性最大的位置进行漂移靠近（也就是下图中哪个坡度最大的地方）。
• PDF in feature space
  还是上面所说的，这个概率密度函数的这个概率，表示的是目标模型出现的概率，而这个目标模型是通过目标特征进行建模的（颜色直方图、纹理特征等）。当整个图片数据同样使用特征提取，PDF就负责表示与建模后特征相似的概率, 也就是说这个PDF是基于图像特征空间的。
  

基本追踪思路： 假设在视频的第一帧框定一个目标对象，对目标对象进行特征建模。在接下来的视频序列中，利用算法（不断向目标出现可能性最大的地方进行靠近），对该框定目标进行连续的跟踪，即给出当前目标框所需的偏移量。

• Motivation – to track non-rigid objects, (like a walking person), it is hard to specify an explicit 2D parametric motion model.

• Appearances of non-rigid objects can sometimes be modeled with color distributions

• The mean-shift algorithm is an efficient approach to tracking objects whose appearance is defined by color. (not limited to only color, however. Could also use edge orientations, texture, motion)

二、原理

• 参考https://github.com/xuweilin2014/learning/blob/master/mean_shift_tracking/mean_shift.pdf

2.1：直觉原理

直觉原理我定义为，就是不需要有任何公式，而是用语言逻辑和图像的形象化表示，让你直觉上对这个算法的原理和可行性有很好的理解。并且我始终不喜欢那种一上来讲原理就抛给你一堆公式的讲解方法。

其实认真阅读上面第一节的“介绍&基本原理”部分，我相信你已经达到了这一点，这里我就不再重复赘述。

下面再来看看理论&公式part

2.2：理论&公式part

 首先扯扯无参密度估计理论，无参密度估计也叫做非参数估计，属于数理统计的一个分支，和参数密度估计共同构成了概率密度估计方法。参数密度估计方法要求特征空间服从一个已知的概率密度函数，在实际的应用中这个条件很难达到。而无参数密度估计方法对先验知识要求最少，完全依靠训练数据进行估计，并且可以用于任意形状的密度估计（目标根据特征建模后的特征模型，就是“任意形状”，它不是确定的）。所以依靠无参密度估计方法，即不事先规定概率密度函数的结构形式，在某一连续点处的密度函数值可由该点邻域中的若干样本点估计得出。常用的无参密度估计方法有：直方图法、最近邻域法和核密度估计法。

MeanShift算法正是属于 核密度估计法 ，它不需要任何先验知识而完全依靠 特征空间中样本点的计算其密度函数值。对于一组采样数据，直方图法通常把数据的值域分成若干相等的区间，数据按区间分成若干组，每组数据的个数与总参数个数的比率就是每个单元的概率值；核密度估计法的原理相似于直方图法，只是多了一个用于 平滑数据的核函数 。采用核函数估计法，在采样充分的情况下，能够渐进地收敛于任意的密度函数，即可以对服从任何分布的数据进行密度估计。

然后谈谈MeanShift的基本思想及物理含义：

此外，从公式1中可以看到，只要是落入Sh的采样点，无论其离中心x的远近，对最终的Mh(x)计算的贡献是一样的。然而在现实跟踪过程中，当跟踪目标出现遮挡等影响时，由于外层的像素值容易受遮挡或背景的影响，所以目标模型中心附近的像素比靠外的像素更可靠。因此，对于所有采样点，每个样本点的重要性应该是不同的，离中心点越远，其权值应该越小。故引入 核函数和权重系数来提高跟踪算法的鲁棒性并增加搜索跟踪能力。

接下来，谈谈核函数：⭐


核函数也叫窗口函数，在核估计中起到平滑的作用。常用的核函数有：Uniform，Epannechnikov，Gaussian等。本文算法只用到了Epannechnikov，它数序定义如下：

2.3：Meanshift目标跟踪算法流程

基于均值漂移的目标跟踪算法通过分别计算目标区域和候选区域内像素的特征值概率得到关于目标模型和候选模型的描述，然后利用 相似函数度量初始帧目标模型和当前帧的候选模版的相似性，选择使相似函数最大的候选模型 并得到关于 目标模型的Meanshift向量，这个向量正是目标 由初始位置向正确位置移动的向量。由于均值漂移算法的 快速收敛性，通过不断迭代计算Meanshift向量，算法最终将收敛到目标的真实位置，达到跟踪的目的。

下面通过图示直观的说明MeanShift跟踪算法的基本原理。

如下图所示：

• 目标跟踪开始于数据点$x{i}^0$（空心圆点$x{i}^0$，$x{i}^1$，…，$x{i}^N$表示的是中心点，上标表示的是的迭代次数;
• 周围的黑色圆点表示不断移动中的窗口样本点
• 虚线圆圈代表的是密度估计窗口的大小）
• 箭头表示样本点相对于核函数中心点的漂移向量，平均的漂移向量会指向样本点最密集的方向，也就是梯度方向。
• 因为 Meanshift 算法是收敛的，因此在当前帧中通过反复迭代搜索特征空间中样本点最密集的区域，搜索点沿着样本点密度增加的方向“漂移”到局部密度极大点点$x{i}^N$，也就是被认为的目标位置，从而达到跟踪的目的，MeanShift 跟踪过程结束。

运动目标的实现过程【具体算法】：

2.4：MeanShift 算法的优缺点

优点：

1. 算法计算量不大，在目标区域已知的情况下完全可以做到实时跟踪；

2. 采用核函数直方图模型，对边缘遮挡、目标旋转、变形和背景运动不敏感。

缺点：

1. 跟踪过程中由于窗口宽度大小保持不变，框出的区域不会随着目标的扩大（或缩小）而扩大（或缩小）；

2. 当目标速度较快时，跟踪效果不好；

3. 直方图特征在目标颜色特征描述方面略显匮乏，缺少空间信息；

三、python示例代码

说明：

首先第一帧中选择目标区域后，按下enter键后即可对视频中的目标进行追踪

import math
import numpy as np
import cv2


def get_tr(img):
    # 定义需要返回的参数
    mouse_params = {'x': None, 'width': None, 'height': None,
                    'y': None, 'temp': None}
    cv2.namedWindow('image')
    # 鼠标框选操作函数
    cv2.setMouseCallback('image', on_mouse, mouse_params)
    cv2.imshow('image', img)
    cv2.waitKey(0)
    return [mouse_params['x'], mouse_params['y'], mouse_params['width'],
            mouse_params['height']], mouse_params['temp']


def on_mouse(event, x, y, flags, param):
    global img, point1
    img2 = img.copy()
    if event == cv2.EVENT_LBUTTONDOWN:  # 左键点击
        point1 = (x, y)
        cv2.circle(img2, point1, 10, (0, 255, 0), 5)
        cv2.imshow('image', img2)
    elif event == cv2.EVENT_MOUSEMOVE and (flags & cv2.EVENT_FLAG_LBUTTON):  # 按住左键拖曳
        cv2.rectangle(img2, point1, (x, y), (255, 0, 0), 5)
        cv2.imshow('image', img2)
    elif event == cv2.EVENT_LBUTTONUP:  # 左键释放
        point2 = (x, y)
        cv2.rectangle(img2, point1, point2, (0, 0, 255), 5)
        cv2.imshow('image', img2)
        # 返回框选矩形左上角点的坐标、矩形宽度、高度以及矩形包含的图像
        param['x'] = min(point1[0], point2[0])
        param['y'] = min(point1[1], point2[1])
        param['width'] = abs(point1[0] - point2[0])
        param['height'] = abs(point1[1] - point2[1])
        param['temp'] = img[param['y']:param['y'] + param['height'],
                        param['x']:param['x'] + param['width']]


def main():
    global img
    cap = cv2.VideoCapture("D:/CollegeStudy/AI/windElectProject/VS2019_cpp/data/pump/pump2.mp4")
    # 获取视频第一帧
    ret, frame = cap.read()
    img = frame
    # 框选目标并返回相应信息：rect为四个信息，temp为框选出来的图像
    rect, temp = get_tr(img)
    (a, b, c) = temp.shape
    y = [a / 2, b / 2]

    # 计算目标图像的权值矩阵
    m_wei = np.zeros((a, b))
    for i in range(a):
        for j in range(b):
            z = (i - y[0]) ** 2 + (j - y[1]) ** 2
            m_wei[i, j] = 1 - z / (y[0] ** 2 + y[1] ** 2)

    # 计算目标权值直方图
    C = 1 / sum(sum(m_wei))
    hist1 = np.zeros(16 ** 3)
    for i in range(a):
        for j in range(b):
            q_b = math.floor(float(temp[i, j, 0]) / 16)
            q_g = math.floor(float(temp[i, j, 1]) / 16)
            q_r = math.floor(float(temp[i, j, 2]) / 16)
            q_temp1 = q_r * 256 + q_g * 16 + q_b
            hist1[int(q_temp1)] = hist1[int(q_temp1)] + m_wei[i, j]
    hist1 = hist1 * C

    # 接着读取视频并进行目标跟踪
    while (1):
        ret, frame = cap.read()

        if ret == True:
            Img = frame
            num = 0
            Y = [1, 1]

            # mean shift迭代
            while (np.sqrt(Y[0] ** 2 + Y[1] ** 2) > 0.5) & (num < 20):
                num = num + 1

                # 计算候选区域直方图
                temp2 = Img[int(rect[1]):int(rect[1] + rect[3]), int(rect[0]):int(rect[0] + rect[2])]
                hist2 = np.zeros(16 ** 3)
                q_temp2 = np.zeros((a, b))
                for i in range(a):
                    for j in range(b):
                        q_b = math.floor(float(temp2[i, j, 0]) / 16)
                        q_g = math.floor(float(temp2[i, j, 1]) / 16)
                        q_r = math.floor(float(temp2[i, j, 2]) / 16)
                        q_temp2[i, j] = q_r * 256 + q_g * 16 + q_b
                        hist2[int(q_temp2[i, j])] = hist2[int(q_temp2[i, j])] + m_wei[i, j]
                hist2 = hist2 * C

                w = np.zeros(16 ** 3)
                for i in range(16 ** 3):
                    if hist2[i] != 0:
                        w[i] = math.sqrt(hist1[i] / hist2[i])
                    else:
                        w[i] = 0

                sum_w = 0
                sum_xw = [0, 0]
                for i in range(a):
                    for j in range(b):
                        sum_w = sum_w + w[int(q_temp2[i, j])]
                        sum_xw = sum_xw + w[int(q_temp2[i, j])] * np.array([i - y[0], j - y[1]])
                Y = sum_xw / sum_w

                # 位置更新
                rect[0] = rect[0] + Y[1]
                rect[1] = rect[1] + Y[0]

            v0 = int(rect[0])
            v1 = int(rect[1])
            v2 = int(rect[2])
            v3 = int(rect[3])
            pt1 = (v0, v1)
            pt2 = (v0 + v2, v1 + v3)

            # 画矩形
            IMG = cv2.rectangle(Img, pt1, pt2, (0, 0, 255), 2)
            cv2.imshow('IMG', IMG)
            k = cv2.waitKey(60) & 0xff
            if k == 27:
                break
        else:
            break


if __name__ == '__main__':
    main()