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
來源：知乎