MOT中的Data Association（二）：最小代價流

4 最小代價流

4.1 算法形式

在了解最小代價流之前，我們需要先鋪墊一下幾個常見圖模型，以幫助我們理解，比如最短路、最大流、最小費用最大流，最小割（閉嘴，我暫時沒看懂）。下圖是一個很常見的圖網絡：

我們可以看到，圖上有很多節點和邊，這兩個元素是組成圖模型的核心。其次，每條邊上都會有對應的數值，比如最短路問題中的相鄰兩節點的距離，最大流中的邊容量，最小費用最大流問題中的邊容量和費用。那么我們來看看幾個問題的具體定義：

最短路問題

最短路問題一般特指單源單匯最短路問題，即給定起點和終點，從各種路徑中選擇最短的路徑。

上面公式中如果結合圖模型來思考會簡單很多，即中間節點無論會不會通過，其流入邊和流出邊一定有且只有0或1個，不可能經過這個點而不經過與之相鄰的邊。不過對於起點和終點則允許有一條流出邊或一條流入邊。

最大流問題

最大流問題就是選擇從起點到終點的最大流量分配，與最短路的最大區別在於 最短路問題中每個節點只能選擇一條流出邊和一條流入邊，而最大流問題則只要滿足邊容量限制，則可任意選擇流入流出邊數量。

 

上面公式的意思是，即中間節點無論會不會通過，其流入邊流量之和=流出邊流量之和，從起點流出的總流量=流入終點的總流量，每條邊的流量有上限。

最小費用流問題

最小費用流的約束條件和最大流的一樣，只不過為了更好描述目標函數我改寫成了類似於最短路問題的形式。 其目標是選擇費用最短的流，當然，它跟最大流問題不同，這里需要設置起始點的流出流量，而且，如果在最大流限制下求解最小費用流，那么就是最小費用最大流問題了。

 

 聯系到多目標跟蹤任務，其數據關聯任務從短期來看就是一個二分圖匹配問題，從長期來看就是一個圖網絡模型。

作者：黃飄
鏈接：https://zhuanlan.zhihu.com/p/111397247
來源：知乎
著作權歸作者所有。商業轉載請聯系作者獲得授權，非商業轉載請注明出處。

1）如果我們要用最短路模型來描述數據關聯問題，其節點就是跟蹤對象id，邊代表跟蹤軌跡和檢測之間的代價。那單源單匯最短路模型就遠遠不足以描述，因為跟蹤軌跡和檢測數量是大於1的，所以從形式上來講是多源多匯最短路，但是最短路沒有限制中間節點的可重復性，所以這個問題應該用路由問題中的K –最短不相交路線（K shortest disjoint paths）來描述；

（2）如果我們要用最大流模型來描述該問題，那么不同於最短路，邊容量代表跟蹤軌跡和檢測之間的連接可能性，所以只可能是0和1。最終要求的就是最大流量，由於邊容量的限制，所以不可能重復，也就是最多可能軌跡。而且這么看來，匈牙利算法很像是最大流模型的特例；

（3）如果我們要用最小費用流模型來描述該問題，那么就跟第一部分中的最短路問題一樣了，只不過多源多匯問題變成了給定初始流量的情形，距離變成了費用。最大的區別在於需要合理設定初始流量（代表了最終有多少條軌跡），還要設定邊容量，不然容易所有流都流向同一條邊；

（4）結合（2）（3）來看，最大流模型只需要設定跟蹤軌跡和檢測的連接可能性，但是缺乏了相對性。而最小費用流則只需要設定代價值，但是需要設定初始流量。這里的初始流量代表了軌跡數量，所以先用最大流模型求出最大流，即可作為初始的軌跡數量，然后再求最小費用流即可，也就是最小費用最大流。不過兩個任務都有着相同的任務，那就是尋找目標軌跡，所以這樣來說時間效率會較低。一般來說我們會直接使用最大流模型/最小割模型，或者直接使用最小費用流+搜索算法。

總的來說，最大流模型的優點是參數量少，但是確定跟匈牙利算法一樣，我們無法對於0.9和0.5相似度的邊進行相對選擇，因為都是1。最小費用流模型的優點是保留了相似度，但是初始流量這一超參數不好設定。K-最短不相交路模型跟最小費用流一樣，都需要設定軌跡數量。所以我們會選擇用搜索算法使用最小費用流，通過搜索閾值使用最大流模型.

4.2 基於最大化后驗概率模型的網絡流建圖

對於最小費用流而言，最難的地方在於設定初始流量和邊容量，使得跟蹤軌跡不交叉，而且跟蹤軌跡盡可能多而合理。最重要的是，我們不知道在網絡模型中哪個節點是軌跡的起點或者終點，這些都需要我們去建模。再加上我們的目標是使得代價最小，極可能最終出現每條軌跡只有一個節點的情形。

下面我們要設定幾個代價值，由於每個點都有屬於軌跡起點和終點的可能性，所以網絡會非常大，為了更好地借鑒已有的最小費用流模型，我們可以轉化為單一起點和終點的網絡圖：

我們可以看到，簡單的最小代價流結構存在幾個問題：

• 一開始我們就要選定哪些節點有可能成為起點，哪些節點有可能成為終點，這無疑增加了參數量；
• 類似於最大流模型，我們通過設定邊容量為1可以保證每條邊最多被選擇一次。但是，我們無法確保最多只有一個節點可以連接到目標節點，這就不能保證跟蹤軌跡的不重疊。

針對以上問題，我們可以引入過渡節點的概念，同時也就引入了過渡邊，每個節點連接一條過渡邊，這樣通過設定過渡邊容量，可以限制每個節點的流出流量，相應地就可以限制最多只有一個節點可以連接到目標節點。而且，我們讓每個節點都連接起點，每個節點的過渡節點連接終點，這樣就保證上面兩個問題都解決了。具體網絡結構如下：

我們可以看到每個節點u都連接了起點，每個節點v都連接了終點，每個節點u都連接了過渡節點v。正如我們之前說的，每條邊容量都是1，可以有效防止軌跡重疊。另外我在圖中注明了每條邊費用的取值范圍，我們定義每條包含起點和終點的邊的費用都是比較大的正數，節點與過渡節點的邊的費用是負數，這樣可以避免過早的終止軌跡，過渡節點與節點之間的邊的費用就是跟蹤軌跡和檢測的代價值，取正數，不然每條軌跡都會在最后一幀終止。所以這里的參數有：初始流量的大小（軌跡數量）、節點屬於軌跡起點/終點的概率、節點到過渡節點的補償（選擇這個節點的補償）、過渡節點到節點的概率（匹配代價）。

下面我們聯系多目標跟蹤模型的形式來為這些參數賦予特殊的含義，首先給出后驗概率形式，T表示已有軌跡，Z表示觀測值：

 

 如果我們不考慮目標之間的聯系，即假設目標相互獨立，將聯系歸於代價值之中。那么上式就可以轉化為：

這里我們就需要了解三個概率：

 

 

另外，我們還需要補充節點的概念來完善數據關聯模型，因為上面的幾個概念中我們還沒有加入軌跡的終點的概率，所以明確一下聯合概率數據關聯模型：

上圖中虛線部分代表非當前時刻的節點，這樣我們就將虛警和雜波利用起點和終點消除了，由於我們可以跨幀連接，所以虛擬目標就可以近似忽略。

接下來我們開始分析概率模型，我們可以利用對數似然概率來描述：

 

其中 就是當前跟蹤軌跡和檢測之間的相似度， 就是當前跟蹤軌跡存在的概率， ​表示該觀測存在的概率，因為觀測有可能是雜波或者虛警，我們可以用於過渡邊的代價描述。最后就只剩下 和 ​兩個概率值，這個我們可以當做是一個參數進行試驗。就這樣，我們最小代價流模型中的每個節點和每條邊都賦予了有意義的概念。

4.3 在線和離線跟蹤分析

這里我們所說的在線和離線模型的意思是一幀一幀使用min cost flow或者多幀一起優化。對於在線的優化，我們就不需要考慮有多個節點連接同一個節點的特殊情況了，也就是說我們可以直接忽略過渡邊，這樣就跟KM算法一模一樣了，所以我們可以認為匈牙利算法是最大流模型的特殊情況，KM算法是最小費用流的特殊情況。

接下來我們分別對在線和離線的最小代價流模型進行對比實驗，其中在線最小代價流模型我們還是采用Kalman+馬氏距離的方式構建代價矩陣。而離線的方式下我們則直接使用IOU和HSV直方圖作為構建代價矩陣的指標。而對於觀測量的概率，即決定過渡邊權的指標，我們采用檢測的置信度(ln(a*confidence+b))作為指標，而對於起點和終點的判定，我們將其作為超參數，連同代價閾值和特征衰減變量作為超參數。

其中特征衰減變量是對軌跡短暫消失的懲罰：

 

 

作者：黃飄
鏈接：https://zhuanlan.zhihu.com/p/111397247
來源：知乎
著作權歸作者所有。商業轉載請聯系作者獲得授權，非商業轉載請注明出處。

另外，無論是在線跟蹤還是離線跟蹤，MinCostFlow這個任務本身都需要設定初始流量，也就是跟蹤軌跡數量，這個值我們都知道是最少是1，最多是總id數。那么我們就需要用搜索算法來解決，為了保證求解效率，我們簡單假設這個問題是一維凸優化問題，采用二分搜索或者斐波那契搜索來進行。

其中二分搜索很簡單，對於斐波那契搜索，我們知道斐波那契數列{0,1,1,2,3…}，即f(n)=f(n-1)+f(n-2)。對於這個通項公式，我們可以看到對於長度為f(n)的搜索空間，可以將其分為f(n-1)和f(n-2)兩個部分，這樣就實現了搜索空間的縮減。下面給出具體的算法：

可以看到我上面用了哈希表來存儲搜索過程中的結果，避免重復運算，保證O(1)的查詢效率。另外，對於求解結果，一般返回的是匹配點對，我們如何將其變成軌跡呢？這就是一個經典的“朋友圈”問題，可以采用並查集來求解，只需要O(n)的時間復雜度和O(n)的空間復雜度。不過我們這個問題簡單一點，不存在一個點對應多個點的情況，所以可以簡單利用數組或哈希表建立多叉樹求解。

4.4 代碼實現

為了方便，我們直接用IOU和HSV顏色直方圖作為特征進行試驗。由於代碼太長，我這里這放一部分，其中的Fibonacci搜索過程代碼如下：

def fibonacci(self, n):
    """Use Fibonacci Search to speed up Searching
    there can exist u~v flows(id), so we need to find the min cost flows
​
    Parameters:
    -----------
    n: int
​
    Returns:
    -----------
    fn: int
        the n th fibonacci number
    """
    assert n > -1, "n must be non-negative number"
​
    if n in self.fib:
        return self.fib[n]
    else:
        return self.fib.setdefault(n, self.fibonacci(n - 1) + self.fibonacci(n - 2))
​
def fibonacci_search(self):
    """Run Fibonacci Searching to find the min cost flow
​
    Returns
    -------
    trajectories: List[List]
        List of trajectories
    min_cost: float
        cost of assignments
    """
    k = 0
    r = max(0, self.max_flow - self.min_flow)
    s = self.min_flow
    cost = {}
    trajectories = []
​
    # find the nearest pos of fibonacci
    while r > self.fibonacci(k):
        k = k + 1
​
    while k > 1:
        u = min(self.max_flow, s + self.fibonacci(k - 1))
        v = min(self.max_flow, s + self.fibonacci(k - 2))
​
        if u not in cost:
            self.graph.SetNodeSupply(0, u)
            self.graph.SetNodeSupply(1, -u)
            if self.graph.Solve() == self.graph.OPTIMAL:
                cost[u] = self.graph.OptimalCost()
​
            else:
                cost[u] = np.inf
​
        if v not in cost:
            self.graph.SetNodeSupply(0, v)
            self.graph.SetNodeSupply(1, -v)
            if self.graph.Solve() == self.graph.OPTIMAL:
                cost[v] = self.graph.OptimalCost()
​
            else:
                cost[v] = np.inf
​
        if cost[v] == cost[u]:
            s = v
            k = k - 1
        elif cost[v] < cost[u]:
            k = k - 1
        else:
            s = u
            k = k - 2
​
    self.graph.SetNodeSupply(0, s)
    self.graph.SetNodeSupply(1, -s)
​
    if self.graph.Solve() == self.graph.OPTIMAL:
        min_cost =  self.graph.OptimalCost() / multi_factor
        hashlist = {0: []}
        # create disjoint set
        for arc in range(self.graph.NumArcs()):
            if self.graph.Flow(arc) > 0:
                if self.graph.Tail(arc) == 0:
                    hashlist[0].append(self.graph.Head(arc))
                else:
                    hashlist[self.graph.Tail(arc)] = self.graph.Head(arc)
        for entry in hashlist[0]:
            tracklet = [(
                        self.node[entry]['frame_idx'],
                        self.node[entry]['box_idx'],
                        self.node[entry]['box']
                         )]
            point = hashlist[entry]
            while point != 1:
                if self.node[point]['type'] == 'object':
                    tracklet.append((
                        self.node[point]['frame_idx'],
                        self.node[point]['box_idx'],
                        self.node[point]['box']
                         ))
                if point in hashlist:
                    point = hashlist[point]
                else:
                    break
            trajectories.append(tracklet)
​
    else:
        min_cost = inf_cost
​
    return trajectories, min_cost

跟蹤部分代碼：

def process(self, boxes, scores, image = None, features = None, **kwargs):
    """Process one frame of detections.
    Parameters
    ----------
    boxes : ndarray
        An Nx4 dimensional array of bounding boxes in
        format (top-left-x, top-left-y, width, height).
    scores : ndarray
        An array of N associated detector confidence scores.
    image : Optional[ndarray]
        Optionally, a BGR color image;
    features : Optional[ndarray]
        Optionally, an NxL dimensional array of N feature vectors
        corresponding to the given boxes. If None given, bgr_image must not
        be None and the tracker must be given a feature model for feature
        extraction on construction.
    **kwargs : other parameters that model needed
​
    Returns
    -------
    trajectories: List[List[Tuple[int, int, ndarray]]]
        Returns [] if the tracker operates in offline mode. Otherwise,
        returns the set of object trajectories at the current time step.
    entire_trajectories: List[List[Tuple[int, int, ndarray]]]
        entire time steps trajectories
    """
    # save the first node id in current frame
    first_node_id = deepcopy(self.node_idx)
​
    # initialize graph in every time step when online
    if self.mode == "online" and self.current_frame_idx > 1:
        self.graph = pywrapgraph.SimpleMinCostFlow()
        self.trajectories = []
        if self.powersave:
             self.node = {key: self.node[key] for key in self.node \
                          if key not in range(2, self.last_frame_id)}
        for i in range(self.last_frame_id, self.node_idx):
            self.graph.AddArcWithCapacityAndUnitCost(0, int(i), 1, \
                                                     int(multi_factor * self.entry_exit_cost))
​
    # Compute features if necessary.
    parameters = {'image': image, 'boxes': boxes, 'scores': scores,
                  'miss_rate': self.miss_rate, 'batch_size': self.batch_size}
    parameters.update(kwargs)
    if features is None:
        assert self.feature_model is not None, "No feature model given"
        features = self.feature_model(**parameters)
​
​
    # Add nodes to graph for detections observed at this time step.
    observation_costs = (self.observation_model(**parameters)
        if len(scores) > 0 else np.zeros((0,)))
    node_ids = []
    for i, cost in enumerate(observation_costs):
        self.node.update({self.node_idx:
            {
                "type": 'object',
                "box": boxes[i],
                "feature": features[i],
                "frame_idx": self.current_frame_idx,
                "box_idx": i,
                'cost': cost
            }
        } )
​
        # save object node id to this time step
        node_ids.append(self.node_idx)
​
        if self.mode == 'online':
            if self.current_frame_idx == 0:
                self.graph.AddArcWithCapacityAndUnitCost(0, int(self.node_idx), 1, \
                                                     int(multi_factor*self.entry_exit_cost))
            else:
                self.graph.AddArcWithCapacityAndUnitCost(int(self.node_idx), 1, 1, \
                                                         int(multi_factor * self.entry_exit_cost))
            self.node_idx += 1
​
        else:
            self.node.update({self.node_idx + 1:
                {
                    "type": 'transition',
                }
            })
            self.graph.AddArcWithCapacityAndUnitCost(0, int(self.node_idx), 1, \
                                                     int(multi_factor * self.entry_exit_cost))
            self.graph.AddArcWithCapacityAndUnitCost(int(self.node_idx), int(self.node_idx + 1), \
                                                     1, int(multi_factor * cost))
            self.graph.AddArcWithCapacityAndUnitCost(int(self.node_idx + 1), 1, 1, \
                                                     int(multi_factor * self.entry_exit_cost))
            self.node_idx += 2
​
    # Link detections to candidate predecessors.
    predecessor_time_slices = (
        self.nodes_in_timestep[-(1 + self.max_num_misses):])
    for k, predecessor_node_ids in enumerate(predecessor_time_slices):
        if len(predecessor_node_ids) == 0 or len(node_ids) == 0:
            continue
        predecessors = [self.node[x] for x in predecessor_node_ids]
        predecessor_boxes = np.asarray(
            [node["box"] for node in predecessors])
        if isinstance(features,np.ndarray):
            predecessor_features = np.asarray(
                [node["feature"] for node in predecessors])
        else:
            predecessor_features = torch.cat(
                [node["feature"].unsqueeze(0) for node in predecessors])
​
        time_gap = len(predecessor_time_slices) - k
​
        transition_costs = self.transition_model(
            miss_rate = self.miss_rate,
            time_gap = time_gap, predecessor_boxes = predecessor_boxes,
            predecessor_features = predecessor_features,
            boxes = boxes, features = features, **kwargs)
​
        for i, costs in enumerate(transition_costs):
            for j, cost in enumerate(costs):
                if cost > self.cost_threshold:
                    continue
                if self.mode == 'online':
                    last_id = int(predecessor_node_ids[i])
                else:
                    last_id = int(predecessor_node_ids[i] + 1)
                self.graph.AddArcWithCapacityAndUnitCost(last_id, int(node_ids[j]), 1,
                                                         int(multi_factor * cost))
    self.nodes_in_timestep.append(node_ids)
​
    # Compute trajectories if in online mode
    if self.mode == 'online':
        if self.current_frame_idx > 0:
            min_cost, n_flow = self.binary_search(high = min(len(predecessor_time_slices[0]), len(node_ids)))
            if n_flow > 0:
                self.trajectories = self.get_trajectory()
            else:
                self.trajectories = self.node2trajectory(first_node_id, self.node_idx)
            self.entire_trajectories = self.merge_trajectories(self.trajectories, self.entire_trajectories)
        else:
            self.trajectories = self.node2trajectory(2, self.node_idx)
            self.entire_trajectories = deepcopy(self.trajectories)
​
    self.current_frame_idx += 1
    self.last_frame_id = first_node_id
​
    return self.trajectories, self.entire_trajectories

完整版代碼：

https://github.com/nightmaredimple/libmot

在線MOTA=0.676，離線的MOTA=0.594，離線的特征關聯方法很簡單，而在線的用的是Kalman+馬氏距離。以上都是我自己根據自己理解寫的，可能理解有誤，也有可能代碼實現有問題。

作者：黃飄
鏈接：https://zhuanlan.zhihu.com/p/111397247
來源：知乎