1 前言
跟蹤是很多視覺系統中的一個核心模塊，有很多算法都需要使用到跟蹤的信息。比如在基於視頻的行為識別，我們就需要獲得視頻中每個個體的行為片段。在我們項目的pipeline中，跟蹤采用的是DeepSORT算法，而DeepSORT的基礎是SORT算法，所以本文主要先介紹SORT算法，后面另開一篇介紹DeepSORT算法。
2 SORT
2.1 SORT是什么
SORT是論文《Simple Online and Realtime Tracking》的縮寫，它是一個解決多目標跟蹤(Multiple Object Tracking: MOT)問題的算法，該算法基於“tracking-by-detection”框架，且是一個在線跟蹤器(Online Tracker)。而所謂Online Tracker，就是跟蹤器只能利用當前和之前幀的檢測結果去實現跟蹤算法。
SORT算法在設計時的建模有以下特點：
- 不考慮遮擋，無論是短時的還是長時的
- 未使用外觀特征(appearance feature)，在運動估計和數據關聯時只利用了檢測框的位置(postiion)和大小(size)
- 沒有過多考慮跟蹤中的一些corner case以及檢測錯誤，因此算法對detection error的魯棒性可能不是那么好，或者說跟蹤效果的好壞很大程度上受到檢測的影響
2.2 SORT原理
SORT算法主要包括4個模塊：1）檢測模塊；2）運動估計模塊；3）數據關聯模塊；4）被跟蹤物體的建立與銷毀模塊。
檢測模塊
其中檢測模塊采用的是Faster RCNN，這個在實際項目中可以被其它檢測算法替換，比如我們項目中使用的就是YOLO算法。
運動估計模塊
每個物體的狀態定義為\(\mathbf{x}=[u, v, s, r, \dot{u}, \dot{v}, \dot{s}]^{T}\)。假如當前幀檢測出3個物體，運動估計模塊利用Kalman Filter，得到下一幀(或下幾幀)這3個物體的狀態。預測過程使用的是勻速模型，如果下一幀的某個檢測結果和其中某個物體關聯起來了，那么用檢測結果作為觀測值來更新該物體的狀態。
數據關聯模塊
數據關聯就是回答當前檢測得到的物體是現有的哪個物體(assigning detections to existing targets)。假設上一幀檢測出3個物體，當前幀檢測出4個物體，那么我們可以得到一個\(3 \times 4\)的矩陣，矩陣中每個元素表示預測的bbox和當前檢測的bbox的IoU距離(Intersection-over-union distance)。基於這個矩陣，為了使總的IoU距離最小，可以利用匈牙利算法進行匹配/指配，從而完成數據關聯。
被跟蹤物體的建立與銷毀模塊
- 如果某個檢測得到的物體和所有現有的trakcer的IoU都很小(\(< IoU_{min}\))，那么會基於它建立一個新的trakcer，這個物體的速度初始化為0，速度相關的協方差分量會初始化為一個很大的值；
- 如果一個被跟蹤的物體(tracker)如果在\(T_{Lost}\)幀未被檢測關聯上，那么該tracker會被銷毀。
為了跟細致地理解運動估計模塊和數據關聯模塊，下面簡單介紹一下卡爾曼濾波算法和匈牙利算法。
2.2.1 卡爾曼濾波(Kalman Filter)算法
卡爾曼濾波是一種最優估計算法，它是融合/結合各種帶有不確定信息的一種有力工具。卡爾曼濾波算法對內存的占用很少，因為它只需要存儲上一個狀態的信息，同時它的實時性非常好，這使得卡爾曼濾波成為一種得到廣泛應用的算法，1969年的阿波羅登月計划中就使用到了該算法。
下面以運動估計為例，大致看一下Kalman Filter，詳細的卡爾曼濾波及其推導會另開一篇文章說明。
參考《How a Kalman filter works, in pictures》。
假設我們關注的狀態為一個物體的位置和速度:\(\vec{x}=\left[\begin{array}{l} \vec{p} \\ \vec{v} \end{array}\right]\)，且它們間的協方差矩陣(co-variance matrix)為\(\mathbf{P}_{k}=\left[\begin{array}{ll} \Sigma_{p p} & \Sigma_{p v} \\ \Sigma_{v p} & \Sigma_{v v} \end{array}\right]\)。
-
預測
根據運動模型(比如在很短的一小段時間\(\Delta t\)內,可以把物體視作勻速運動),可以得到下一狀態的預測值: \(\begin{aligned} \hat{\mathbf{x}}_{k} &=\left[\begin{array}{cc} 1 & \Delta t \\ 0 & 1 \end{array}\right] \hat{\mathbf{x}}_{k-1} \\ &=\mathbf{F}_{k} \hat{\mathbf{x}}_{k-1} \end{aligned}\)。
上面的預測是基於一個非常理想的模型，但在實際場景中，我們還需要考慮一些已知的外部影響以及環境中的一些額外的不確定性。綜合考慮，可以得到下一時刻的預測：
\(\begin{aligned} \hat{\mathbf{x}}_{k} &=\mathbf{F}_{k} \hat{\mathbf{x}}_{k-1}+\mathbf{B}_{k} \overrightarrow{\mathbf{u}_{k}} \\ \mathbf{P}_{k} &=\mathbf{F}_{\mathbf{k}} \mathbf{P}_{k-1} \mathbf{F}_{k}^{T}+\mathbf{Q}_{k} \end{aligned}\)
其中\(\mathbf{B}_{k}\)表示已知的外部影響，\(\mathbf{Q}_{k}\)表示環境中額外的不確定性。 -
測量更新
通過傳感器，我們可以獲得對狀態的一組觀測值。如果要獲得一組靠譜的預測值，需要調和通過運動模型得到的預測值和實際測量的觀測值。如果上述兩個值分別用一個多維高斯分布進行表示，那么更好的預測應該是這兩個多維高斯分布的乘積。
在預測階段，我們以運動模型為基礎，給定\(k-1\)時刻的狀態，預測了\(k\)時刻的狀態，但是我們的測量結果並不在狀態空間，通過\(H_k\)矩陣，可以將狀態空間轉換到測量空間：
\(\begin{aligned} &\vec{\mu}_{\text {expected }}=\mathbf{H}_{k} \hat{\mathbf{x}}_{k} \\ &\mathbf{\Sigma}_{\text {expected }}=\mathbf{H}_{k} \mathbf{P}_{k} \mathbf{H}_{k}^{T} \end{aligned}\) -
兩正太分布相乘仍為正態分布
這里直接給出公式，其中\(K\)矩陣稱為卡爾曼增益:
- 完整的Kalman FIlter
結合上面對卡爾曼濾波的簡單介紹，可以看看源碼對它的實現:
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):
"""
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
# dim_x: 狀態值是一個7維向量[u, v, s, r, u', v', s']
# dim_z: 測量值是一個4維向量[u, v, s, r]
# 代碼中的F、H、R、P、Q都可以和上面的公式對應
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]])
# F和H基於模型進行設置，但是R、P、Q的設置好像是根據經驗進行設置？
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):
"""
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)
2.2.2 匈牙利(Hungarian)算法
匈牙利算法在有些資料中也稱作KM算法。SORT利用匈牙利算法解決了二分圖中的最小權重匹配問題，我們可以用以下例子來理解這個問題。假設二分圖的第一部分有N個節點，表示N個工人，第二部分有M個節點，表示M件工作。矩陣C中的元素\(C[i, j]\)表示將地j個工作分配給第i個工人需要的花費。試問如何分配工作，才能使得花費最小？
SORT中數據關聯的源碼如下，其中linear_sum_assignment()函數是sklearn中匈牙利算法的實現。
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)))
def iou_batch(bb_test, bb_gt):
"""
From SORT: Computes IOU between two bboxes in the form [x1,y1,x2,y2]
"""
bb_gt = np.expand_dims(bb_gt, 0)
bb_test = np.expand_dims(bb_test, 1)
# 求兩個bbox的iou，橫軸表示x，右為正，縱軸表示y，下為正
# (x1, y1)表示左上角點，(x2, y2)表示右下角點
# (xx1, yy1)和(xx2, yy2)分別表示兩個bbox交集的左上角點和右下角點
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 associate_detections_to_trackers(detections,trackers,iou_threshold = 0.3):
"""
Assigns detections to tracked object (both represented as bounding boxes)
Returns 3 lists of matches, unmatched_detections and unmatched_trackers
"""
if(len(trackers)==0):
return np.empty((0,2),dtype=int), np.arange(len(detections)), np.empty((0,5),dtype=int)
iou_matrix = iou_batch(detections, trackers)
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)
else:
matched_indices = np.empty(shape=(0,2))
# 沒有匹配到現有tracker的檢測物體
unmatched_detections = []
for d, det in enumerate(detections):
if(d not in matched_indices[:,0]):
unmatched_detections.append(d)
# 沒有匹配到檢測物體的tracker
unmatched_trackers = []
for t, trk in enumerate(trackers):
if(t not in matched_indices[:,1]):
unmatched_trackers.append(t)
#filter out matched with low IOU
matches = []
for m in matched_indices:
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)
return matches, np.array(unmatched_detections), np.array(unmatched_trackers)