多目标跟踪Sort算法及绘制多目标跟踪的运动轨迹

参考博文:

https://www.freesion.com/article/9284823172/

https://blog.csdn.net/App_12062011/article/details/51758989

目标跟踪就是在检测出某视频序列初始帧的目标位置后,预测后续帧中该目标的大小与位置。

背景:

       可能有的人会问为什么不每帧检测,但前提是检测算法得足够鲁棒,对每一帧图像中的目标都能检测出来,否则就会出现不连续的情况,即明明存在目标,却没有检测出来,这时候就需要用到跟踪算法了。这里用Sort(SIMPLE ONLINE AND REALTIME TRACKING)多目标跟踪算法。

原始论文:http://arxiv.org/pdf/1602.00763.pdf

原始论文中的开源代码:https://github.com/abewley/sort

Sort算法原理:

       Sort多目标跟踪算法是一种解决跟踪的BBox与检测的BBox相关联的算法,关联两个BBox主要的思想是:通过计算两个BBox之间的IOU + 匈牙利算法选择最优关联结果。Sort算法是在卡尔曼滤波和匈牙利算法的基础上将卡尔曼滤波预测的BBox与目标检测的BBox进行匹配,对于匹配上的,选择目标检测的BBox更新下一帧的目标预测BBox而已。所以理解卡尔曼滤波预测和匈牙利关联这两个算法也就理解Sort算法了。

       卡尔曼滤波:可以在任何含有不确定信息的动态系统中使用卡尔曼滤波,对系统的下一步走向做出有根据的预测,即使伴随着各种干扰,卡尔曼滤波总是能指出真实发生的情况。比如:假设物体当前时刻在BBox1的位置,卡尔曼滤波可以预测出物体下一时刻在BBox2的位置(但其实是个很粗糙的预测)。一般包括三个过程:(1)初始化:卡尔曼滤波器的状态变量和观测输入;(2)跟踪器列表中目标框的预测;(3)对于匹配上的目标框用检测框更新(并非直接替换)。

(1)初始化卡尔曼滤波器的状态变量(7维)和观测输入(4维):

首先状态变量x(self.kf.x)是一个七维向量:

x=[u,v,s,r,\dot{u},\dot{v},\dot{s}]

前四维分别表示目标框中心位置的x,y坐标,面积s和当前目标框的长宽比(即为观测输入,通过实参传入),最后三维则是横向,纵向,面积的变化速率,其中变化速率初始化为0。

(2)跟踪器列表中目标框的预测

首先插入一个形象的例子,这样后面理解起来会容易点:

       假如现在有一辆在路上做直线运动的小车,该小车在 t 时刻的状态可以用一个向量来表示\binom{p_{t}}{v^{_{t}}},其中pt 表示当前的位置,vt表示该车当前的速度。当然,司机还可以踩油门或者刹车来给车一个加速度ut,ut相当于是一个对车的控制量。显然,如果司机既没有踩油门也没有踩刹车,那么ut就等于0,此时车就会做匀速直线运动。

       如果已知上一时刻 t-1时小车的状态,现在来考虑当前时刻t 小车的状态,显然有:

v_{t} =v_{t-1}+u_{t}*\Delta t

p_{t}=p_{t-1}+v_{t-1}*\Delta t+1/2*u_{t}*\Delta t^{2}

从以上公式可以看出:输出变量pt和vt都是输入变量的线性组合,这也就是卡尔曼滤波器被称为线性滤波器的原因。既然上述公式表征了一种线性关系,那么就可以用一个矩阵来表示:

其中的分别就是状态转移矩阵F(表示如何从上一状态来推测当前状态)和控制矩阵B(表示控制量ut(加速度)如何作用于当前状态),至此就是卡尔曼滤波器的状态预测公式。

在卡尔曼滤波器中状态转移矩阵F根据运动学公式可以确定为:

F=(1000100

        0100010

        0010001

        0001000

        0000100

        0000010

        0000001)

在卡尔曼滤波器中运动形式(加速度)和状态转移矩阵的确定都是基于匀速运动模型,所以当前帧的状态可以如下表示:

当前帧的状态变量 = F*上一帧的状态变量

这时的状态是根据上一状态推测而来的,是一个估计值,应该对其进行修正以得到最优估计,那么就理应考虑到噪声的影响,若要引入噪声的影响,只要考虑噪声的方差σ即可,维度提高后,为了综合考虑各个维度偏离其均值的程度,需要引入协方差矩阵P。噪声同目标的状态一样也会进行传递,且预测模型本身并不是绝对准确的,用协方差矩阵 Q 来表示预测模型本身的噪声,那么所有噪声的传递可以表示如下(表示噪声在各个时刻间的传递关系):

(3)对于匹配上的目标框用检测框更新(并非直接替换)

从目标的真实状态到我们看到的观测状态之间还有一个变换关系,这个变换关系叫做观测函数H,在卡尔曼滤波器中观测函数H可以确定为:

H =(1000000

         0100000

         0010000

         0001000)

那么真实状态就可以表示为:

Z = H*x+ a 

Z表示真实状态,x表示观测状态,a表示误差

现在对前面得出的状态估计值进行修正(这里就是对匹配上的检测框与预测框进行处理的具体操作:即通过检测框Z去更新预测框):

修正值 = 估计值 + 卡尔曼系数*(Z-H*估计值)

还有对噪声的传递状态进行修正:

Sort算法流程:

前提:在跟踪之前,对所有目标已经完成检测。

1、第一帧进来时,以检测到的目标初始化,对检测到的每个目标都分别创建新的卡尔曼滤波器对象,并将每个创建的对象都加入到跟踪器列表中;

2、后面帧进来时,先到卡尔曼滤波器中得到由前面帧BBox产生的状态预测和协方差预测,计算跟踪器列表中所有目标的状态预测与当前帧中检测到的BBox的IOU,再通过匈牙利算法分别得到所有目标IOU最大的唯一匹配(数据关联部分),最后再去掉匹配值小于iou_threshold的匹配对;

3、针对匹配成功的目标,用当前帧中匹配到的目标检测BBox去更新卡尔曼滤波器的状态预测BBox。对于当前帧中没有匹配到的目标(两种情况:当前帧只有检测框,没有跟踪框;当前帧有跟踪框,但没有检测框),分两种情况进行处理,具体处理见问题解答和代码解析部分。

问题解答:

当时我看完原始文献和开源代码之后,还存在几个方面的问题,花了几天也没弄明白,现在大概弄清楚了,这里做个小结:

(1)没有匹配到的目标在当前帧没有跟踪框吗?(即当前帧没有检测到目标,在当前帧会有跟踪框吗?)【这个问题现在看起来很傻,但在当时还是花了一点时间的,比如跟踪框到底跟踪几帧视频)

       在原始论文的代码中如果连续两帧都没有检测框的话,即没有匹配到目标,那么只会跟踪一帧,即只有一帧是有跟踪框的,也就出现了在连续两帧的第二帧中没有了跟踪框。如果要改跟踪帧数的话,可以按如下修改代码:

if (self.time_since_update > 1):#如果连续两帧都没有检测框的话,那么跟踪2帧
    self.hit_streak = 0

(2)当前帧检测到目标,但为什么没有跟踪框?

       这种情况的本质原因是self.min_hits设置的大了,虽然进行了更新,self.hit_streak+1,但是还没有达到self.min_hits,这样更新的目标框不会加入到跟踪结果列表中,也就不会显示跟踪框了,比如一般把self.min_hits设置为1,这样就是当第一次检测到目标的时候只有检测框,没有跟踪框,从第二帧开始就会有跟踪框了(前提:第二帧还存在检测框),如果self.min_hits设置为0,那就是从第一次检测到目标开始,就会出现跟踪框,这样的弊端就是如果检测的是错的话,跟踪的就一定是错的,没有反应的时间,所以一般都把self.min_hits设置为1,即连续两帧都检测及匹配到的话,就在第二帧开始绘制跟踪框。

       这样就能理解为什么当前帧检测到目标,但为什么没有跟踪框了,因为当前帧(不包括第一帧,第一帧是会显示跟踪框的)是第一次出现目标(前提是self.min_hits=1)。

(3)预测和跟踪怎么理解,预测了就会有跟踪框吗?

       不是的,只有匹配上更新过的才会有跟踪框(前提self.min_hits=1),如果没有匹配上做更新,永远也不会触发

trk.hit_streak >= self.min_hits

也就不会有跟踪框了,但在跟踪器列表中的框都会做预测,除非连续预测的次数已经大于max_age,这时框会从跟踪器列表中删除,认为目标已经彻底的从画面中消失了。

代码解析:

"""
    SORT: A Simple, Online and Realtime Tracker
    Copyright (C) 2016-2020 Alex Bewley alex@bewley.ai

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
"""
from __future__ import print_function

import os
import numpy as np
import matplotlib

matplotlib.use('TkAgg')
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from skimage import io

import glob
import time
import argparse
from filterpy.kalman import KalmanFilter

# try:
#     from numba import jit
# except:
#     def jit(func):
#         return func

np.random.seed(0)


def linear_assignment(cost_matrix):
    try:
        import lap
        _, x, y = lap.lapjv(cost_matrix, extend_cost=True)
        return np.array([[y[i], i] for i in x if i >= 0])  #
    except ImportError:
        from scipy.optimize import linear_sum_assignment
        x, y = linear_sum_assignment(cost_matrix)
        return np.array(list(zip(x, y)))


# @jit
def iou(bb_test, bb_gt):
    """
    Computes IUO between two bboxes in the form [x1,y1,x2,y2]
    """
    xx1 = np.maximum(bb_test[0], bb_gt[0])
    yy1 = np.maximum(bb_test[1], bb_gt[1])
    xx2 = np.minimum(bb_test[2], bb_gt[2])
    yy2 = np.minimum(bb_test[3], bb_gt[3])
    w = np.maximum(0., xx2 - xx1)
    h = np.maximum(0., yy2 - yy1)
    wh = w * h
    o = wh / ((bb_test[2] - bb_test[0]) * (bb_test[3] - bb_test[1])
              + (bb_gt[2] - bb_gt[0]) * (bb_gt[3] - bb_gt[1]) - wh)
    return (o)


def convert_bbox_to_z(bbox):
    """
    Takes a bounding box in the form [x1,y1,x2,y2] and returns z in the form
      [x,y,s,r] where x,y is the centre of the box and s is the scale/area and r is
      the aspect ratio
    """
    w = bbox[2] - bbox[0]
    h = bbox[3] - bbox[1]
    x = bbox[0] + w / 2.
    y = bbox[1] + h / 2.
    s = w * h  # scale is just area
    r = w / float(h)
    return np.array([x, y, s, r]).reshape((4, 1))


def convert_x_to_bbox(x, score=None):
    """
    Takes a bounding box in the centre form [x,y,s,r] and returns it in the form
      [x1,y1,x2,y2] where x1,y1 is the top left and x2,y2 is the bottom right
    """
    w = np.sqrt(x[2] * x[3])
    h = x[2] / w
    if (score == None):
        return np.array([x[0] - w / 2., x[1] - h / 2., x[0] + w / 2., x[1] + h / 2.]).reshape((1, 4))
    else:
        return np.array([x[0] - w / 2., x[1] - h / 2., x[0] + w / 2., x[1] + h / 2., score]).reshape((1, 5))


class KalmanBoxTracker(object): #卡尔曼滤波器的python实现
    """
    This class represents the internal state of individual tracked objects observed as bbox.
    """
    count = 0

    def __init__(self, bbox):
        """
        Initialises a tracker using initial bounding box.
        """
        # define constant velocity model
        self.kf = KalmanFilter(dim_x=7, dim_z=4)
        self.kf.F = np.array(
            [[1, 0, 0, 0, 1, 0, 0], [0, 1, 0, 0, 0, 1, 0], [0, 0, 1, 0, 0, 0, 1], [0, 0, 0, 1, 0, 0, 0],
             [0, 0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 0, 0, 1]]) #状态变量模型
        self.kf.H = np.array(
            [[1, 0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0]])#观测函数

        self.kf.R[2:, 2:] *= 10. #测量噪声矩阵
        self.kf.P[4:, 4:] *= 1000.  # give high uncertainty to the unobservable initial velocities 协方差矩阵
        self.kf.P *= 10.
        self.kf.Q[-1, -1] *= 0.01 #过程噪声矩阵
        self.kf.Q[4:, 4:] *= 0.01

        self.kf.x[:4] = convert_bbox_to_z(bbox)
        self.time_since_update = 0
        self.id = KalmanBoxTracker.count
        KalmanBoxTracker.count += 1
        self.history = []
        self.hits = 0
        self.hit_streak = 0
        self.age = 0

    def update(self, bbox):#用检测框替换跟踪器self.trackers列表中对应的跟踪框
        """
        Updates the state vector with observed bbox.
        """
        self.time_since_update = 0
        self.history = []
        self.hits += 1
        self.hit_streak += 1
        self.kf.update(convert_bbox_to_z(bbox))

    def predict(self):#卡尔曼滤波对跟踪器列表中的目标进行下一时刻位置的预测
        """
        Advances the state vector and returns the predicted bounding box estimate.
        """
        if ((self.kf.x[6] + self.kf.x[2]) <= 0):
            self.kf.x[6] *= 0.0
        self.kf.predict()
        self.age += 1
        if (self.time_since_update > 0):#这里更改跟踪的帧数,默认只跟踪一帧
            self.hit_streak = 0
        self.time_since_update += 1
        self.history.append(convert_x_to_bbox(self.kf.x))
        return self.history[-1]

    def get_state(self):#坐标转换
        """
        Returns the current bounding box estimate.
        """
        return convert_x_to_bbox(self.kf.x)


def associate_detections_to_trackers(detections, trackers, iou_threshold=0.1):
    #将物体检测的BBox与卡尔曼滤波器预测的跟踪BBox匹配
    """
    Assigns detections to tracked object (both represented as bounding boxes)

    Returns 3 lists of matches, unmatched_detections and unmatched_trackers
    """
    if (len(trackers) == 0):#第一帧没有跟踪框,只有检测框,所以返回3个值:(1)匹配到的[d,t](空的);(2)没有匹配到的检测框;(3)没有匹配到的跟踪框(空的)
        return np.empty((0, 2), dtype=int), np.arange(len(detections)), np.empty((0, 5), dtype=int)
    iou_matrix = np.zeros((len(detections), len(trackers)), dtype=np.float32)

    for d, det in enumerate(detections):#遍历物体检测BBox集合,每个BBox标识为d
        for t, trk in enumerate(trackers):#遍历卡尔曼滤波器预测的BBox集合,每个BBox标识为t
            iou_matrix[d, t] = iou(det, trk)

    if min(iou_matrix.shape) > 0:
        a = (iou_matrix > iou_threshold).astype(np.int32)
        if a.sum(1).max() == 1 and a.sum(0).max() == 1:
            matched_indices = np.stack(np.where(a), axis=1)
        else:
            matched_indices = linear_assignment(-iou_matrix) #通过匈牙利算法匹配卡尔曼滤波器预测的BBox与物体检测BBox以[[d,t],,,]的二维矩阵保存到matched_indices
    else:
        matched_indices = np.empty(shape=(0, 2))

    unmatched_detections = []
    for d, det in enumerate(detections):#没有匹配上的物体检测BBox放入unmatched_detections列表,表示有新的物体进入画面了,后面要新增跟踪器追踪新物体
        if (d not in matched_indices[:, 0]):
            unmatched_detections.append(d)
    unmatched_trackers = []
    for t, trk in enumerate(trackers):#没有匹配上的卡尔曼滤波器预测的BBox放入unmatched_trackers列表,表示之前跟踪的物体离开画面了,后面要删除对应的跟踪器
        if (t not in matched_indices[:, 1]):
            unmatched_trackers.append(t)

    # filter out matched with low IOU
    matches = []
    for m in matched_indices:#遍历matched_indices矩阵,将IOU值小于iou_threshold的匹配结果分别放入unmatched_detections,unmatched_trackers列表中
        if (iou_matrix[m[0], m[1]] < iou_threshold):
            unmatched_detections.append(m[0])
            unmatched_trackers.append(m[1])
        else:
            matches.append(m.reshape(1, 2))
    if (len(matches) == 0):
        matches = np.empty((0, 2), dtype=int)
    else:
        matches = np.concatenate(matches, axis=0)#匹配上的卡尔曼滤波器预测的BBox与物体检测的BBox以[[d,t],,,]的形式放入matches矩阵

    return matches, np.array(unmatched_detections), np.array(unmatched_trackers) #返回跟踪成功的物体矩阵,新增物体的矩阵,离开画面的物体矩阵


class Sort(object):
    def __init__(self, max_age=3, min_hits=1):
        """
        Sets key parameters for SORT
        """
        self.max_age = max_age
        self.min_hits = min_hits
        self.trackers = []
        self.frame_count = 0

    def update(self, dets=np.empty((0, 5))):
        """
        Params:
          dets - a numpy array of detections in the format [[x1,y1,x2,y2,score],[x1,y1,x2,y2,score],...]
        Requires: this method must be called once for each frame even with empty detections (use np.empty((0, 5)) for frames without detections).
        Returns the a similar array, where the last column is the object ID.

        NOTE: The number of objects returned may differ from the number of detections provided.
        """
        self.frame_count += 1
        # get predicted locations from existing trackers.
        trks = np.zeros((len(self.trackers), 5))#根据当前所有的卡尔曼跟踪器的个数(等于上一帧中被跟踪的物体个数)创建二维矩阵,trks行号为卡尔曼跟踪器标识,列向量为跟踪BBox与物体跟踪ID
        to_del = []
        ret = []
        for t, trk in enumerate(trks):#循环遍历卡尔曼跟踪器列表
            pos = self.trackers[t].predict()[0]#用卡尔曼跟踪器t产生对应物体在当前帧中预测的BBox
            trk[:] = [pos[0], pos[1], pos[2], pos[3], 0]
            if np.any(np.isnan(pos)):
                to_del.append(t)
        trks = np.ma.compress_rows(np.ma.masked_invalid(trks))#trks中存放了上一帧中被跟踪的所有物体在当前帧中预测的BBox
        for t in reversed(to_del):
            self.trackers.pop(t)
        matched, unmatched_dets, unmatched_trks = associate_detections_to_trackers(dets, trks)#将物体检测的BBox与卡尔曼滤波器预测的跟踪BBox匹配,获得跟踪成功的物体矩阵,新增物体的矩阵,离开画面的物体矩阵

        # update matched trackers with assigned detections
        for m in matched:
            self.trackers[m[1]].update(dets[m[0], :])#跟踪成功的物体BBox信息更新到对应的卡尔曼滤波器

        # create and initialise new trackers for unmatched detections
        for i in unmatched_dets:#新增的物体要创建新的卡尔曼滤波器对象用于跟踪
            trk = KalmanBoxTracker(dets[i, :])
            self.trackers.append(trk)
        i = len(self.trackers)
        for trk in reversed(self.trackers):#跟踪成功的物体BBox与ID放入ret列表中
            d = trk.get_state()[0]
            #3表示连续预测的次数,self.min_hits不设置为0是因为第一次检测到的目标不用跟踪,不能设大,一般就是1
            if (trk.time_since_update < 3) and (trk.hit_streak >= self.min_hits or self.frame_count <= self.min_hits):#self.min_hits为最小更新的次数
                ret.append(np.concatenate((d, [trk.id + 1])).reshape(1, -1))  # +1 as MOT benchmark requires positive
            i -= 1
            # remove dead tracklet
            if (trk.time_since_update > self.max_age):#self.max_age为连续预测的最大次数
                self.trackers.pop(i)
        if (len(ret) > 0):
            return np.concatenate(ret)#返回当前画面中所有被跟踪物体的BBox与ID,二维矩阵[[x1,y1,x2,y2,id1],,,[,,,]],在main函数的帧循环中将当前帧的跟踪结果绘制和显示到屏幕上
        return np.empty((0, 5))

if __name__ == '__main__':
    mot_tracker = sort.Sort(max_age=5, min_hits=1)
    total_frames = 0
    for image_name in images:
        im = cv.imread(image_name)
        total_frames += 1
        trackers = mot_tracker.update(detect_result)
        trail_dic = dict()
        for t in trackers:
            start_X = int(d[0])
            start_Y = int(d[1])
            end_X = int(d[2])
            end_Y = int(d[3])
            id = int(d[4])
            center = cal_center(t)
            trail_dic.setdefault(id, []).append(center.tolist())
            cv.rectangle(im, (start_X, start_Y), (end_X, end_Y), (0, 0, 255), 2)#绘制跟踪框
            if total_frames > 1:
                for key, value in trail_dic.items():
                    for i in range(len(value)-1):
                        index_start = i
                        index_end = index_start + 1
                        cv.line(im, tuple(value[index_start]), tuple(value[index_end]), (0, 255, 0), thickness=2, lineType=8)#绘制运动轨迹
    cv.imshow("result", im)
    cv.waitKey(100)

关键函数/参数解释:

1、time_since_update:当前连续预测的次数,只要调用KalmanBoxTracker类中update函数,time_since_update就会清零
2、hit_streak:判断当前是否做了更新,大于等于1的说明做了更新,只要连续帧中没有做连续更新,hit_streak就会清零
3、max_age:连续预测的最大次数,就是放在self.trackers跟踪器列表中的框用卡尔曼滤波器连续预测位置的最大次数 
4、min_hits:最小更新的次数,就是放在self.trackers跟踪器列表中的框与检测框匹配上,然后调用卡尔曼滤波器类中的update函数的最小次数,min_hits不设置为0是因为第一次检测到的目标不用跟踪,只需要加入到跟踪器列表中,不会显示,这个值不能设大,一般就是1,表示如果连续两帧都检测到目标,
5、ret:跟踪结果列表,返回形式:[[x1,y1,x2,y2,id1],,,[,,,]],id表示多个目标的编号,以进行多个目标的区分

可以修改的地方以修改跟踪的效果:

1、if (self.time_since_update > 0):#更改跟踪的帧数
    self.hit_streak = 0
2、def __init__(self, max_age=3, min_hits=1):
3、trk.time_since_update < 3

 

评论 11
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值