Tracking-KCF Algorithm

注：本文涉及的算法的代码实践已上传至GitHub，恳求大佬们指点！^ _ ^

1. 目标检测跟踪与算法背景概述

目标跟踪任务在许多的计算机视觉系统中都是极为关键的一个组成部分。对于任意给定的一个初始图像的Patch（Filter滑过的区域），目标跟踪任务的目的在于训练一个分类器来将待跟踪的目标与它所处的环境区分开，为了能够在后续帧中能继续检测到这个目标，分类去需要能够在很多位置上都能进行详尽的评估，同时在滑动的过程中都会提供一个新的图像Patch来帮助提升模型的性能。

在上述任务中，我们将感兴趣的对象——即待追踪的目标称作正样本（positive samples），将目标所在的环境或者背景称作负样本。一幅图像中包含的负样本数量几乎是无限的，在达到合并尽可能多的样本和保持较低的计算量之间的平衡上，前人作出了很多尝试。KCF算法中，引入了轮转矩阵（circulant matrices）这一工具来合并大量的样本，实现了一个基于“核岭回归”的跟踪器（tracer），该跟踪器可以认为是一个核化版的线性相关滤波器（Linear Correlation Filter ）。

Correlation Filter（CF）源于信号处理领域，其用于tracking方面的想法是：相关是衡量两个信号相似值的度量，如果两个信号越相似，那么其相关值就越高。在tracking的应用里，设计一个滤波模板，利用该模板与目标候选区域做相关运算，最大输出响应的位置即为当前帧的目标位置。

用数学语言来表述，xxx表示一幅图像，w\mathcal{w}w表示相应的滤波模板，yyy表示模板与相关候选区域进行相关运算后的相应输出，则有如下关系：y=x⊗wy=x\otimes\mathcal{w}y=x⊗w，为减小计算量，可将x,y,wx,y,\mathcal{w}x,y,w分别转换到其傅里叶空间后计算点积，即y^=x^⋅w^∗\hat{y} = \hat{x}· \hat{\mathcal{w}}^*y^=x^⋅w^∗。相关滤波的任务则是寻找最优的滤波模板w\mathcal{w}w。

目下以相关滤波类的方法和深度学习方法为代表的判别类模型方法效果普遍好于生成模型，其显著特点是分类器的训练过程中运用到了前景与背景信息，分类器的任务专注于区分前景与背景。Kernelized Correlation Filter(KCF)是相关滤波方法中的典型算法，该方法的一般流程是：在跟踪的过程中训练一个目标检测器，使用目标检测器去检测下一帧预测位置是否是目标，然后再使用新检测结果去更新训练集进而更新目标检测器。

2. HOG特征

KCF算法中使用了图像的HOG特征替代了传统跟踪器所用的灰度特征。

HOG特征全称为Histogram of Oriented Gradients，即方向梯度直方图。作为一种图像的特征描述子（图像的一种表示，通过提取有用的信息并扔掉多余的信息来简化图像），HOG特征将一张大小为width×height×channelswidth \times height \times channelswidth×height×channels的图像转化为一个长度为nnn的特征向量，例如输入图像大小为64×128×364 \times 128 \times 364×128×3，HOG特征的输出为长度378037803780的向量。HOG特征的计算过程如下：

(1)图像预处理

最早提出的HOG特征通过在一张64×12864 \times 12864×128的图像上计算得到，预处理的操作要求图片保持1:2的横纵比，原文中提及的γ\gammaγ矫正已知其增益效果较小，不再考虑。

(2)计算梯度图

首先计算图像的水平方向梯度gxg_xgx和竖直方向的梯度gyg_ygy，以如下公式来计算梯度的强度值ggg和梯度方向θ\thetaθ：
g=gx2+gy2θ=arctangxgyg=\sqrt{g_x^2+g_y^2}\\ \theta = arctan \frac{g_x}{g_y} g=gx2+gy2θ=arctangygx
同时梯度方向有如下的性质：θ∈[0,π]\theta \in [0,\pi]θ∈[0,π]。

(3)计算梯度直方图

这一步中，首先将图像分割成多个8×88 \times 88×8的cell，在这些cell中计算梯度直方图。(2)中计算的每个像素点的位置处包含了2个值：g和θg和\thetag和θ，一个cell就保存了128个值，单个像素的梯度信息往往包含了噪声，采用8×88 \times 88×8的图片块表示后能够是的直方图对噪声不敏感。

由(2)的性质，将0-180度分成9个区间：0,20,40,…160，之后统计每个像素点所在的区间——将这个区间命名为bin，采取的原则是对每个像素点处的ggg值，按θ\thetaθ的比例将ggg分配给相邻的bin，如下图所示。

最终统计得到如下图所示的直方图：

(4)对16×1616 \times 1616×16大小的Block进行标准化

标准化（Normalization）也称归一化，即将每个向量的分量除以向量的模长。Block选取示意图如下：

(5)得到HOG向量

将(4)中计算的单个Block的向量连接起来得到整个图片块的最终的特征向量，其中可以由这样的认识：每个16×16块由36×1向量表示。

OpenCV中封装了提取HOG特征的类HOGDescriptor，实现过程中只需构造一个hog对象来实现，见于HOG.py文件中的构造方法：

    def __init__(self, winSize):self.winSize = winSizeself.blockSize = (8, 8)self.blockStride = (4, 4)self.cellSize = (4, 4)self.nbins = 9self.hog = cv2.HOGDescriptor(winSize,self.blockSize, self.blockStride,self.cellSize, self.nbins)

3. 傅里叶变换

3.1 傅里叶变换概念简述

傅里叶变换（Fourier transform，FT）是一种线性积分变换，用于函数（应用上称作“信号”）在时域和频域之间的变换，在物理学和工程学中有许多应用，其作用是将函数分解为不同特征的正弦函数的和，如同化学分析来分析一个化合物的元素成分。对于一个函数，也可对其进行分析，来确定组成它的基本（正弦函数）成分。

经过傅里叶变换生成的函数f^\hat {f}f^称作原函数fff的傅里叶变换，应用意义上称作频谱。在特定情况下，傅里叶变换是可逆的，即将f^\hat{f}f^通过逆变换可以得到其原函数fff。通常情况下，fff是一个实函数，而f^\hat{f}f^ 则是一个复数值函数，其函数值作为复数可同时表示振幅和相位。高斯函数是傅里叶变换的本征函数（也称固有函数，即对已定义的函数空间中任意一个非零函数fff 进行变换仍然是函数 fff或者其标量倍数的函数。更加精确的描述就是Af=λf\mathcal{A}f = \lambda fAf=λf其中λ\lambdaλ是标量，它是对应的特征值），即若一个函数(或分布)在真实空间中是一个高斯函数，则在傅里叶空间中它还是一个高斯函数，如下图所示：

维基百科中对连续可积函数fff的傅里叶变换给出的定义形式如下：
f^(ξ)=∫−∞∞f(x)e−2πixξdx，ξ为定义在频域上的任意实数\hat{f}(\xi)=\int _{-\infty}^{\infty}f(x)e^{-2 \pi ix\xi}dx， \xi为定义在频域上的任意实数 f^(ξ)=∫−∞∞f(x)e−2πixξdx，ξ为定义在频域上的任意实数
同时在适当的条件下，f^\hat{f}f^的逆傅里叶变换可得到fff：
f(x)=∫−∞∞f^(ξ)e−2πixξdξ，x为定义在时域上的任意实数f(x)=\int _{-\infty}^{\infty}\hat{f}(\xi)e^{-2 \pi ix\xi}d\xi， x为定义在时域上的任意实数 f(x)=∫−∞∞f^(ξ)e−2πixξdξ，x为定义在时域上的任意实数

3.2 傅氏空间与离散傅里叶矩阵

傅里叶变换的奇妙之处就在于，将真实空间中复杂的波形通过“分解”之后可以得到若干简单的正弦波，这个“分解”操作是为傅里叶变换，从而“傅里叶空间”也可理解为“频率空间”，如下图所示：

离散傅里叶变换矩阵则是将离散傅里叶变换以矩阵乘法来表示的一种表达形式，其一般形式为：
F=1n[11...111w...wn−2wn−11w2...w2(n−2)w2(n−1)...............1wn−1...w(n−1)(n−2)w(n−1)2]F=\frac{1}{\sqrt{n}}\begin{bmatrix} 1&1&...&1&1\\ 1&\mathcal{w}&...&\mathcal{w}^{n-2}&\mathcal{w}^{n-1}\\ 1&\mathcal{w}^2&...&\mathcal{w}^{2(n-2)}&\mathcal{w}^{2(n-1)}\\ {...}&{...}&{...}&{...}&{...}&\\ \\ 1&\mathcal{w}^{n-1}&...&\mathcal{w}^{(n-1)(n-2)}&\mathcal{w}^{(n-1)^2}\end{bmatrix} F=n1⎣⎢⎢⎢⎢⎢⎢⎡111...11ww2...wn−1...............1wn−2w2(n−2)...w(n−1)(n−2)1wn−1w2(n−1)...w(n−1)2⎦⎥⎥⎥⎥⎥⎥⎤

4.快速核相关计算

在part1中提及，实现快速相关计算的方法为：将模板、图像以及相应输出转换到傅里叶空间做点积运算。同时，KCF中创新性地使用了循环转换模型(cyclic shift model)，即前文提及的循环矩阵来解决了计算量的问题，将计算复杂度由O(n2)\mathcal{O}(n^2)O(n2)降为了O(nlongn)\mathcal{O}(nlong\mathcal {n})O(nlongn)。

具体来说，由如下定理：循环矩阵都能够在傅氏空间中使用离散傅里叶矩阵进行对角化，即X=Fdiag(x^)FHX = Fdiag(\hat{x})F^HX=Fdiag(x^)FH。

对于多项式K(X,X′)=(XTX′+a)b\mathcal{K}(X,X')=(X^TX'+a)^bK(X,X′)=(XTX′+a)b，可有如下的表示形式：KXX′=(F−1(x^∗⊙x′^)+a)b\mathcal{K}^{XX'} = (\mathcal{F^{-1}}(\hat{x}^* \odot\hat{x'})+a)^bKXX′=(F−1(x^∗⊙x′^)+a)b，相关的证明过程参见原文，其中的F−1\mathcal{F^{-1}}F−1代表离散傅里叶逆变换(IDFT)，x^∗与x′^\hat{x}^* 与\hat{x'}x^∗与x′^代表各自在傅里叶空间中的映射。

径向基函数可以认为是∣∣xi−xj∣∣2||x_i-x_j||^2∣∣xi−xj∣∣2的函数，对于一些核函数h，径向基核如有如下的形式：
k(X,X′)=h(∣∣xi−xj∣∣2)\mathcal{k}(X,X')=h(||x_i-x_j||^2) k(X,X′)=h(∣∣xi−xj∣∣2)
于是在高斯核中，径向基函数的形式为：Kxx′=exp(−1σ2(∣∣xi∣∣2+∣∣xj∣∣2)−2(F−1(x^∗⊙x′^))T)K^{xx'}=exp(-\frac{1}{\sigma^2}(||x_i||^2+||x_j||^2)-2(\mathcal{F}^{-1}(\hat{x}^* \odot\hat{x'}))^T)Kxx′=exp(−σ21(∣∣xi∣∣2+∣∣xj∣∣2)−2(F−1(x^∗⊙x′^))T)，其中关于循环矩阵对角化及其推导的内容参见原文。这一部分计算的实现过程如下：

    def kernel_correlation(self, x1, x2, sigma):"""核化的相关滤波操作:param x1::param x2::param sigma:   :return:"""# 转换到傅里叶空间fx1 = fft2(x1)fx2 = fft2(x2)# 相乘并返回共轭复数tmp = conj(fx1) * fx2# 离散傅里叶逆变换转换回真实空间idft_rbf = ifft2(np.sum(tmp, axis=0))# 将零频率分量移到频谱中心。idft_rbf = fftshift(idft_rbf)d = np.sum(x1 ** 2) + np.sum(x2 ** 2) - 2.0 * idft_rbf# 径向基函数k = np.exp(-1 / sigma ** 2 * np.abs(d) / d.size)return k

5. 响应值的计算与分类器训练

训练的过程即是在核空间中做岭回归。具体来说，线性不可分的样本，经过非线性映射函数ϕ(x)\phi (x)ϕ(x)后可变为线性可分的样本，此时所处的新空间即时核空间，此时的任务是：在核空间中使用岭回归寻找一个分类器f(xi)=wTϕ(xi)f(x_i)=\mathcal{w}^T \phi(x_i)f(xi)=wTϕ(xi)，目标在于寻找最优的w\mathcal{w}w，根据梯度取零可以表示出w\mathcal{w}w的向量表达式，详细推导过程参见原文。将w\mathcal{w}w表示为如下形式后，可得到其对偶问题w=∑iαiϕ(xi)\mathcal{w}=\sum_i\alpha_i\phi(x_i)w=∑iαiϕ(xi)，对偶表达为
α=minα∣∣ϕ(X)ϕ(X)Tα−y∣∣2+λ∣∣ϕ(X)Tα∣∣\alpha = min_{\alpha}||\phi(X)\phi(X)^T\alpha-y||^2+\lambda||\phi(X)^T\alpha|| α=minα∣∣ϕ(X)ϕ(X)Tα−y∣∣2+λ∣∣ϕ(X)Tα∣∣
同样令其梯度为零可求得其表达式：α∗=(ϕ(X)ϕ(x)T+λI)−1y\alpha ^* = (\phi(X)\phi(x)^T+\lambda I)^{-1}yα∗=(ϕ(X)ϕ(x)T+λI)−1y，由于并不清楚非线性映射函数ϕ(X)\phi(X)ϕ(X)的具体形式，并且只关心刻画在核空间的核矩阵(ϕ(X)ϕ(x)T(\phi(X)\phi(x)^T(ϕ(X)ϕ(x)T，故令KKK表示这个核矩阵。

KCF算法始终未离开一个工具——循环矩阵，借助循环矩阵的傅里叶对角化来简化计算，核矩阵虽然得到了表示但是仍旧不知道其具体数值、形式，接下来的任务就变成了将KKK对角化——希望找到一个核函数能使得对应的核矩阵是循环矩阵。原文中给出了如下的定理：

Given circulant data C (x), the corresponding kernel matrix K is circulant if the kernel function satisfifies

K(x,x′)=K(Mx,Mx′)K(x,x') = K(Mx, Mx')K(x,x′)=K(Mx,Mx′), for any permutation matrixM。

即：第一个样本和第二个样本都是由生成样本循环移位产生的，可以不是同一个样本。从而，在train函数中核矩阵的计算，以及α\alphaα的计算方式有如下表示：

    def train(self, x, y, sigma, lambdar):"""原文所给参考train函数:param x::param y::param sigma::param lambdar::return:"""k = self.kernel_correlation(x, x, sigma)return fft2(y) / (fft2(k) + lambdar)

在计算得到α\alphaα之后，即可借助α\alphaα来计算响应值：

    def detect(self, x, z, sigma):"""原文所给参考detect函数:param x::param z: x的HOG特征:param sigma: 高斯参数:return:"""k = self.kernel_correlation(x, z, sigma)# 傅里叶逆变换的实部，对应响应的值return np.real(ifft2(self.alphaf * fft2(k)))

之后即可根据响应值来更新检测的目标，详细代码参见update()函数。

6. 运行结果

详见于result文件夹。

Reference

High-Speed Tracking with Kernelized Correlation Filters： https://arxiv.org/pdf/1404.7584
Histogram of Oriented Gradients explained using OpenCV (learnopencv.com) : https://learnopencv.com/histogram-of-oriented-gradients/
傅里叶变换 - 维基百科：https://zh.wikipedia.org/wiki/傅里叶变换
从真实空间到傅立叶空间 - 知乎： https://zhuanlan.zhihu.com/p/37061414
KCF目标跟踪方法分析与总结 - 一只有恒心的小菜鸟 - 博客园 : https://www.cnblogs.com/YiXiaoZhou/p/5925019.html
GitHub代码仓库： https://github.com/chuanqi305/KCF

Appendix：库函数说明以及原文相关说明

1. HOG.py

项目1	对应内容
函数名	cv2.HOGDescriptor.compute（）
功能描述	计算给定图像的HOG特征描述子
参数	img: 待计算的图像Patch winStride: 窗口滑动的步长 padding : Padding
返回值	拼接好的n*1维HOG特征向量

2. KCF.py

项目1	对应内容
函数名	numpy.conj()
功能描述	帮助用户对任何复数进行共轭
参数	x[array_like]: 输入值
返回值	numpy.ndarray: x的复共轭

3. runKCF.py

项目1	对应内容
函数名	cv2.VideoCapture(video_path or device index)
功能描述	实例化一个VideoCapture对象
参数	video_path：以.mp4/.avi结束的string类型，代表视频所在路径 or device index：指定camare的序号，值为0或-1
返回值	VideoCapture对象

项目2	对应内容
函数名	cv2.VideoCapture.read()
功能描述	逐帧读取视频，返回值分两部分：是否正确读取以及每一帧的图像
参数	无
返回值	bool ret：是否正确读取&&是否读取到视频末尾 np.ndarray frame：每一帧的图像，是一个三维矩阵

项目3	对应内容
函数名	cv2.selectROI(windowName, img, showCrosshair=None, fromCenter=None)
功能描述	在一幅图像中对特定图像区域以矩形框的形式进行截取
参数	windowName：选择的区域被显示在的窗口的名字 img：要在什么图片上选择ROI showCrosshair：是否在矩形框里画十字线. fromCenter：是否是从矩形框的中心开始画
返回值	tuple (min_x, min_y, w, h)，依次表示：矩形框中最小的x值矩形框中最小的y值矩形框的宽矩形框的高

项目4	对应内容
函数名	cv2.rectangle(image, start_point, end_point, color, thickness)
功能描述	在一幅图像上绘制一个矩形
参数	image：待绘制矩阵的图像 start_point：矩形的起始坐标，tuple类型 end_point：矩阵的结束坐标，tuple类型 color：绘制边框线的颜色 thickness：边框线的粗细
返回值	np.ndarray：image

项目5	对应内容
函数名	cv2.VideoCapture.isOpen()
功能描述	判断摄像头是否初始化成功
参数	无
返回值	bool res：True or False

项目6	对应内容
函数名	argparse.ArgumentParser()
功能描述	创建一个命令行参数解释器
参数	无(本例中)
返回值	一个ArgumentParser对象

项目7	对应内容
函数名	argparse.ArgumentParser.parse_args()
功能描述	给ArgumentParser对象中的参数赋值
参数	无(本例中)
返回值	argparse.Namespace 参数名等

4. 原文参考源码与超参数

参考源码：

单目标跟踪--KCF算法(核化相关滤波算法)Python实现(超详细)相关推荐

KCF算法（相关滤波算法）跟踪目标
这是笔者目前碰到的性能最好的传统目标跟踪算法,所以这里也简单介绍下这个算法的由来. 一.算法介绍 KCF全称为Kernel Correlation Filter 核相关滤波算法.是在2014年由Joa ...
单目标跟踪小结（2）：算法
内容主要提取自 edX 平台上Chalmers的micromaster项目:Emerging Automotive Technologies: Sensor Fusion and Non-linear ...
单目标跟踪（一）相关滤波器方法综述
<Handcrafted and Deep Trackers: Recent Visual Object Tracking Approachesand Trends> MOSSE: Min ...
python 粒子滤波目标追踪_Python实现基于相关滤波的单目标跟踪算法
最近在阅读一些单目标跟踪的算法论文,主要看了一系列基于相关滤波的算法,尝试着用python实现了其中一些比较经典的算法,地址在https://github.com/wwdguu/pyCFTracker ...
人脸跟踪：KCF核相关滤波算法
一直以来没有很想写这个,以为这个东西比较简单,还算是比较容易理解的一个算法,但是在知乎上回答过一个问题之后就有朋友私信我一些关于细节的东西,我一直以为关于细节的东西大家可以自己去理解,大家都是想快速了 ...
kcf 跟随_KCF相关滤波跟踪算法
写作理由:看了网上各种版本的kcf原理,都不是特别理解(主要是不能特别说服我,哈哈).详细琢磨了下,谈谈自己的理解,本人水平有限,难免有理解不到位的地方,请大家批评指正本文主要会从下面几个方面讲:目 ...
基于嵌入式设备的单目标跟踪算法
基于嵌入式设备的单目标跟踪实现最近基于嵌入式设备(Khadas Vim3)做了一套单目标跟踪算法,跟踪效率可以做到每秒25帧左右. 算法运行耗时记录 time is:37.6241 ms time ...
Video Target Tracking Based on Online Learning—TLD单目标跟踪算法详解
视频目标跟踪问题分析视频跟踪技术的主要目的是从复杂多变的的背景环境中准确提取相关的目标特征,准确地识别出跟踪目标,并且对目标的位置和姿态等信息精确地定位,为后续目标物体行为分析提供足 ...
单目标跟踪——个人笔记
单目标跟踪--个人笔记以<Handcrafted and Deep Trackers: A Review of Recent ObjectTracking Approaches>为主线看 ...

单目标跟踪--KCF算法(核化相关滤波算法)Python实现(超详细)