MeanShift：
如果仅仅使用特征检测来跟踪物体会丢失对应物体的信息。为了解决对应性问题，我们可以使用之前学过的特征匹配和光流法，或者，可以使用meanshift track方法。
meanshift track是一个追踪任意物体的简单但是有效的方法。meanshift的思想就是考虑在一个小的我们感兴趣的区间（ROI）里的从潜在的概率密度函数采样的像素来表示一个目标。
其中，灰色的点代表从某个概率分布里采样的点。点的距离越近他们就越相似。直觉上，meanshif就是在图里找到密度最大的区域然后画一个圈。这个算法会从不稠密的区域画圈，然后移动圆圈到稠密的地方然后固定在哪里。如果这个场景不是点，而是对应彩色直方图，可以用Meanshift找到最接近目标直方图的对象。
meanshift非常适合做目标跟踪。在opencv中使用cv2.meanShift来调用，但是需要一个预过程。大致过程如下：
- 固定一个包围数据点的窗口：可以是一个ROI(range of interest)的bounding box
- 计算窗口内数据的平均值，一般用像素值的直方图，通常会转化为HSV的色彩空间
- 反复移动窗口到平均值直到收敛。这个使用cv2.meanShift完成的。我们可以控制迭代的长度和精度。
CamShift
Meanshift有个不好的地方时，它框出的区域不会随着目标的扩大（或缩小）而扩大（或缩小），CamShift则不会有这样的缺点。
Saliency map
在介绍实例之前，我们介绍一下特征图。使用傅里叶分析可以得到我们对自然图片数据的一般性理解，帮助我们建立一般图片背景的模型。通过比较背景模型与指定图片帧的不同，我们可以得到去掉背景的子区域，这样我们的注意力就会在这些子区域上。这个技术叫做Visual saliency（视觉显著性）。
传统模型试图通过特征匹配或视觉变换来检测目标，这需要手动的标记与训练，可是，如果特征和对象数量都是未知的情况呢？
这里的方法是应用视觉显著性的方法，立即定位到抓住我们注意力的区域（那些脱离常规的数据），这样这个算法就可以跟踪任意数量的物体。
就像我们的大脑会被某些图片上的区域抓住注意力一样，视觉显著性试图去描述对象们的视觉质量，越高则越能够引起重视，同时忽视掉较低的不重要的部分。这在信息丰富的环境中无疑是一个创造性的策略。就像下面两张图中的两根小棍，瞬间就引起你的注意：
但是如果把这些小棍杂糅在一起那么红色的近乎垂直的小棍就很难找到：
所以，要想找到唯一的脱颖而出的目标是比较困难的，那么，应该怎么做呢，就是要让计算机知道自己要把注意力放在什么样的目标上。
- 傅里叶谱图
要想找到视觉显著的区域，我们需要查看它的频率光谱图。一般我们都在空间域查看图片，分析像素或者图片在各个子通道的强度。然而，图片也可以在频率域（frequency domain）中进行表示：通过分析像素的频率或者学习像素出现在图片的周期性。
我们可以通过傅里叶变换将图片从空间域转换到频率域。在频率域中，我们不在考虑图片的坐标，而是将焦点放在图片的谱图上。傅里叶的根本思想就是这个问题：是否能将信号或者图片转化为一系列圆弧路径（或者说：谐波），就像将太阳光里不同的光谱展示出来。对比彩虹中的频率：电磁频率，图片中的是空间频率——像素值的空间周期性。比如一所监狱的图片，空间频率就好比两间邻接监狱的距离。
傅里叶谱图有两个部分，一个是量级（magnitude），一个是相位（phase）。量级描述的是一张图片里不同频率的数量，相位说的是这些频率的空间位置。
右图是左图的傅里叶谱图，右图告诉我们在左图的灰度版本中哪个频率组件是最突出的（最亮）。谱图是调整过的，所以图片的中心对应于x和y的频率0。所以，越往外，频率越高。这张图告诉我们左图里有很多低频率的组件（因为亮度都集中在中间）。
在Opencv中，这个转换可以通过Discrete Fourier Transform(DFT)，包含在Saliency类的plot_magnitude方法中。具体过程如下：
- 转换图片成灰度图：
def plot_magnitude(self):
if len(self.frame_orig.shape)>2:
frame = cv2.cvtColor(self.frame_orig,cv2.COLOR_BGR2GRAY)
else:
frame = self.frame_orig
- 扩展图片到一个最佳的大小：图片大小是2的倍数时DFT的转换速度是最快的，因此一般都用0来填充图片。
rows, cols = self.frame_orig.shape[:2]
nrows = cv2.getOptimalDFTSize(rows)
ncols = cv2.getOptimalDFTSize(cols)
frame = cv2.copyMakeBorder(frame, 0, ncols-cols, 0, nrowsrows, cv2.BORDER_CONSTANT, value = 0)
- 应用DFT：使用Numpy中的fft2函数，返回一个2维复数矩阵。
img_dft = np.fft.fft2(frame) - 将实数值和复数值转换为量级：对复数取绝对值，就是取模。
magn = np.abs(img_dft) - 转换到对数量度（logarithmic scale）：傅里叶系数通常都大到无法在屏幕上显示。一些小的和大的改变值无法观察到。因此，高的值会成为白点，小的值为黑点。为了实现灰度值的可视化，将线性量度转换为对数量度：
log_magn = np.log10(magn) - 平移：将谱图放在中间，方便观察。
spectrum = np.fft.fftshift(log_magn) - 返回结果：
return spectrum/np.max(spectrum)*255
- 自然统计规则
自然世界有很多统计规则，最普遍知道的大概是1/f规则。它陈述了自然图片的集成的振幅遵从1/f 分布，也被成为尺度不变形（scale invariance维基百科）。
一张2维图片的一维功率谱图可以用Saliency类里的plot_power_spectrum函数来视觉化查看。我们可以用与之前量级谱图相似的方式，但我们要保证正确的将2维谱图降到单轴上。
- 转换图片成灰度图：
def plot_power_spectrum(self):
if len(self.frame_orig.shape)>2:
frame = cv2.cvtColor(self.frame_orig,cv2.COLOR_BGR2GRAY)
else:
frame = self.frame_orig
- 扩大图片到最优大小：
rows, cols = self.frame_orig.shape[:2]
nrows = cv2.getOptimalDFTSize(rows)
ncols = cv2.getOptimalDFTSize(cols)
frame = cv2.copyMakeBorder(frame, 0, ncols-cols, 0, nrowsrows, cv2.BORDER_CONSTANT, value = 0)
- 应用DFT得到对数谱图：这里可以选择numpy的傅里叶方法或者opencv的傅里叶方法。
if self.use_numpy_fft:
img_dft = np.fft.fft2(frame)
spectrum = np.log10(np.real(np.abs(img_dft))**2)
else:
img_dft = cv2.dft(np.float32(frame), flags=cv2.DFT_COMPLEX_OUTPUT)
spectrum = np.log10(img_dft[:, :, 0] ** 2 + img_dft[:, :, 1] ** 2)
- 径向平均（radial averaging）：简单的将二维spectrum在x或y方向上取平均数是错误的。这被称为径向平均功率谱图（radially averaged power spectrum ）。
L = max(frame.shape)
freqs = np.fft.fftfreq(L)[:L/2]
dists = np.sqrt(np.fft.fftfreq(frame.shape[0])[:,np.newaxis]**2 + np.fft.fftfreq(frame.shape[1])**2)
dcount = np.histogram(dists.ravel(), bins=freqs)[0]histo, bins = np.histogram(dists.ravel(), bins=freqs, weights=spectrum.ravel())
- 画出结果：记得用bin的值规范化上一步累加的值
centers = (bins[:-1] + bins[1:]) / 2
plt.plot(centers, histo, dcount)
plt.xlabel('frequency')
plt.ylabel('log-spectrum')
plt.show()
结果和频率成反比例的。如果要确认1/f特性，你可以将x的值进行np.log10操作来查看曲线是否是大致上线性递减的。这里仅仅是对y的值取对数操作，结果如下：
这个特性告诉我们如果我们将所有自然图片的所有谱图平均化，也会得到上图的图片。那么，如何利用这个特性来告诉我们的算法把焦点放在Limmat river图的水上的船只而不是旁边的树木呢？
- 利用光谱残留（spectral residual）生成特征图
知道1/f 规则后，我们可以想到，图片上的什么内容会引起我们的注意呢，那就是那些不遵从1/f 规则的内容，那些异常的数据，被叫做光谱残留，对应的就是潜在的引起我们兴趣的对象。将这些异常数据以白点的形式显示的图片就被称为特征图（saliency map）。
光谱残留的方法是在这里提出的：Xiaodi Hou and Liqing Zhang (2007). Saliency
Detection: A Spectral Residual Approach. IEEE Transactions on Computer
Vision and Pattern Recognition (CVPR), p.1-8. doi:
10.1109/CVPR.2007.383267.
生成特征图需要单独处理图片的每个通道。
图片的每个通道可以通过私有方法Saliency.getchannel_sal_magn来生成特征图：
def getchannel_sal_magn(self, channel):
# 计算图片傅里叶谱图的量级和相位
if self.use_numpy_fft:
img_dft = np.fft.fft2(channel)
magnitude, angle = cv2.cartToPolar(np.real(img_dft), np.imag(img_dft))
else:
img_dft = cv2.dft(np.float32(channel), flags=cv2.DFT_COMPLEX_OUTPUT)
magnitude, angle = cv2.cartToPolar(img_dft[:, :, 0], img_dft[:, :, 1])
# 计算傅里叶谱图的振幅的log值，量级的下边界限制在1e-9，防止计算log值时除数为0
log_ampl = np.log10(magnitude.clip(min=1e-9))
# 近似计算典型自然图片的平均谱图
log_ampl_blur = cv2.blur(log_amlp, (3, 3))
# 计算光谱残留，大体上包括了场景里不平凡的部分。
magn = np.exp(log_amlp – log_ampl_blur)
# 通过反向的傅里叶变换计算特征图
if self.use_numpy_fft:
real_part, imag_part = cv2.polarToCart( magn, angle)
img_combined = np.fft.ifft2(real_part + 1j*imag_part)
magnitude, = cv2.cartToPolar(np.real(imgcombined), np.imag(img_combined))
else:
img_dft[:, :, 0], img_dft[:, :, 1] = cv2.polarToCart(residual, angle)
img_combined = cv2.idft(img_dft)
magnitude, = cv2.cartToPolar(imgcombined[:, :, 0], img_combined[:, :, 1])
return magnitude
得到的单通道特征图然后返回给Saliency.get_saliency_map函数，输入图片的每个通道都要经过这样的过程，如果是灰度图就比较简单，只处理一个通道：
def get_saliency_map(self):
if self.need_saliency_map:
# 这一帧的特征图还没有计算
num_channels = 1
if len(self.frame_orig.shape)==2:
# 单个通道
sal = self.getchannel_sal_magn(self.frame_small)
else:
# 考虑每一个通道
sal = np.zeros_like(self.frame_small).astype(np.float32)
for c in xrange(self.frame_small.shape[2]):
sal[:, :, c] = self.getchannel_sal_magn(self.frame_small[:, :, c])
# 多个通道取平均
sal = np.mean(sal, 2)
# 可选的后处理，例如使结果更平滑
if self.gauss_kernel is not None:
sal = cv2.GaussianBlur(sal, self.gauss_kernel,sigmaX=8, sigmaY=0)
# 将sal的值平方，这样可以突出高特征的区域，并将它还原到原始的分辨率并且归一化
sal = sal**2
sal = np.float32(sal)/np.max(sal)
sal = cv2.resize(sal, self.frame_orig.shape[1::-1])
# 为了避免下一次进行这么密集的计算，将当前结果保存，并更改flag
self.saliency_map = sal
self.need_saliency_map = False
return self.saliency_map
下面是得出的特征图：
- 检测对象原型
在一定程度上，saliency map 已经拥有了对象原型的信息，我们要做的就是通过设置阈值获取对象原型。
这里的代码可以让用户选择3倍于平均特征的阈值，或者使用Otsu阈值：
def get_proto_objects_map(self, use_otsu=True):
saliency = self.get_saliency_map()
if use_otsu:
# saliency的范围是0到1所以要乘以255，并且转为uint8型
_, img_objects = cv2.threshold(np.uint8(saliency*255), 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
else:
thresh = np.mean(saliency)*255
_, img_objects = cv2.threshold(np.uint8(saliency*255), thresh, 255, cv2.THRESH_BINARY)
return img_objects
结果如下图：
实例：自动跟踪足球场上所有的选手
将特征检测得到的结果，即特征图作为meanshift的目标输入，视频来自Alfheim数据库，可以从http://home.ifi.uio.no/paalh/dataset/alfheim/下载。
实例主要包括两个过程：
- 将某一帧中的所有原型的bounding boxes标出。特征检测器会在当前帧图片中进行操作，同时，meanshift tracker会在当前帧中寻找上一帧里存在的原型目标。
- 只保留两个算法检测到的目标bounding boxes的交集，就是指两个算法都认为这些bounding boxes框出的是正确的目标。
难点就是多目标的跟踪，这里介绍一下它的实现，具体实现在MultiObjectTracker 类中：
新的一帧获取后调用advance_frame函数，它有一个frame帧参数，另外接受一个目标图proto_objects_map作为参数，然后对frame进行一个深拷贝：
def advance_frame(self, frame, proto_objects_map):
self.tracker = copy.deepcopy(frame)
然后该方法会建立多个bounding boxes作为候选，这些候选框中包括特征图里的也包括从上一帧图片到这一帧的meanshift tracking的结果：
box_all = []
# 添加从当前目标原型图里获得的bouding boxes
box_all = self.appendboxes_from_saliency(proto_objects_map, box_all)
# 找到meanshift tracking在上一帧中的所有bounding boxes
box_all = self.appendboxes_from_meanshift(frame, box_all)
然后要将所有的bounding boxes合并在一起，并且去掉重复的项。通过cv2.groupRectangles方法来完成，如果有group_thresh+1个或者更多的bounding boxes重复了，该方法会返回一个唯一的bounding box，这样一来，如果仅有单独的一个bounding boxes，它就会排除，即只保留交集：
if len(self.object_roi) == 0:
group_thresh = 0 # 没有前一帧，全部bounding boxes来自于特征图
else:
group_thresh = 1 #前一帧加上特征图
box_grouped, _ = cv2.groupRectangles(box_all, group_thresh, 0.1)
要想使meanshift正常工作，需要记录存留的boxes：
# 更新留存的boxes记录
self.updatemean_shift_bookkeeping(frame, box_grouped)
然后将这些没有重复的boxes画出来，并将图片返回：
for (x, y, w, h) in box_grouped:
cv2.rectangle(self.tracker, (x, y), (x+w, y+h), (0, 255, 0), 2)
return self.tracker
-
获取原型目标的bounding boxes
将一个原型目标图以及一个bounding boxes的列表作为输入，首先检测原型目标图的轮廓：
def appendboxes_from_saliency(self, proto_objects_map, box_all):
box_all = []
cnt_sal, = cv2.findContours(proto_objects_map, 1, 2)
去掉比阈值小的块：
for cnt in cnt_sal:
if cv2.contourArea(cnt) < self.min_cnt_area:
continue
将结果存入box_all:
box = cv2.boundingRect(cnt)
box_all.append(box)
return box_all
-
为meanshift tracking建立必须的记录
这个方法需要两个参数，分别是输入的图片和一系列的bounding boxes：
def updatemean_shift_bookkeeping(self, frame, box_grouped):
Bookkeeping主要保存了每个bounding box的HSV彩色值，因此需要将输入图片转换为HSV颜色空间：
hsv = cv2.cvtColor(frame, cv2.COLOR_BGR2HSV)
我们要保存bounding box的位置和大小，也包括HSV彩色值得彩色直方图：
self.object_roi = [] # 彩色直方图
self.object_box = [] # bounding box位置和大小
从box列表里抽取出box的大小和位置，从HSV图片里获取ROI：
for box in box_grouped:
(x, y, w, h) = box
hsv_roi = hsv[y : y+h, x : x+w]
然后计算ROI中H值的直方图，设置一个mask滤掉昏暗的区域，并规范化直方图：
mask = cv2.inRange(hsv_roi, np.array((0., 60., 32.)),np.array((180., 255., 255.)))
roi_hist = cv2.calcHist([hsv_roi], [0], mask, [180], [0, 180])
cv2.normalize(roi_hist, roi_hist, 0, 255, cv2.NORM_MINMAX)
然后存储这个信息到对应的私有变量中，在循环的下一帧可以被使用，同时在下一帧中会使用meanshift算法来找到ROI的区域：
self.object_roi.append(roi_hist)
self.object_box.append(box)
-
使用meanshift算法进行目标跟踪
最后，通过记录的前一帧的bookkeeping信息来跟踪原型目标。与appendboxes_from_meanshift类似，建立一个bounding boxes的列表，它有两个参数，分别是输入的图片和用来存储bounding boxes的列表：
def appendboxes_from_meanshift(self, frame, box_all):
hsv = cv2.cvtColor(frame, cv2.COLOR_BGR2HSV)
然后这个方法解析之前存储的原型目标：
for i in xrange(len(self.object_roi)):
roi_hist = copy.deepcopy(self.object_roi[i])
box_old = copy.deepcopy(self.object_box[i])
为了获取新的ROI的位置，我们将ROI的反向投影作为meanshift算法的参数。终止条件（self.term_crit）保证足够的迭代次数（100）并且均值平移至少一个像素：
dst = cv2.calcBackProject([hsv], [0], roi_hist, [0, 180], 1)
ret, box_new = cv2.meanShift(dst, tuple(box_old),self.term_crit)
在添加新的检测到的平移的bounding box到列表之前，我们希望确认一下，是否是正确的目标。那些不动的目标通常是错误的正例，例如线标或者其他的与任务无关的特征块。
为了丢弃掉无关的跟踪结果，我们比较了新旧box的所在位置：
(xo, yo, wo, ho) = box_old
(xn, yn, wn, hn) = box_new
如果它们的中心没有移动至少sqrt(self.min_shift2)像素，我们就将他们排除：
co = [xo + wo/2, yo + ho/2]
cn = [xn + wn/2, yn + hn/2]
if (co[0] - cn[0])**2 + (co[1] - cn[1])**2 >= self.min_shift2:
box_all.append(box_new)
然后将结果返回：
return box_all
-
将上面所说的组合在一起
这个实例的关键点就在于，通过特征图和meanshift跟踪的相结合，排除掉了那些不动的错误正例。
