單目標跟蹤之相關濾波 MOSSE
論文名稱：Visual Object Tracking using Adaptive Correlation Filters
文獻地址：https://www.cs.colostate.edu/~vision/publications/bolme_cvpr10.pdf
原始碼地址：https://github.com/xingqing45678/Mosse_CF
由於基於CNNs深度學習在單目標跟蹤方法的引數量和計算量都較大，難以與目標檢測演算法一起移植到嵌入式中。同時，受CVPR2020 AutoTrack
本文主要介紹相關濾波系列演算法的開篇——MOSSE基本原理及其python程式碼實現流程。
由於論文中涉及到大量傅立葉變換的公式，可能顯得MOSSE晦澀難懂，本篇儘量淡化傅立葉變換的公式內容，突出MOSSE在跟蹤過程中所做的工作。
相關濾波
相關濾波（CF）源於訊號處理領域，有這麼一句話"兩個訊號越相似，其相關值越高。在跟蹤，就是找到與跟蹤目標響應最大的項" 。通過後續程式碼的解讀，會體會到相關值越高在MOSSE中的作用。
相關操作
假設有兩個訊號f和g，則兩個訊號的相關性（correlation）為：
其中，f∗表示f的複共軛...（說好的淡化公式的... 看不明白不要緊，只要明白卷積就可以往下看）
對於影象而言，相關性可以理解為相關核在影象上進行點積操作，如下圖所示。影象的相關公式可以表達為g = f ⓧ h，具體而言，對應元素相乘在求和：
從上圖的結果可以看出，輸出正好是相關核旋轉了180°。
具體步驟：
a. 相關核，中心分別位於輸入影象的f(i, j)畫素上進行滑動；（會補零）
b. 利用上式進行點積操作，並求和，得到輸出影象對應的畫素值g(i, j);
c. 在輸出影象上進行滑動，重複b的操作，求出輸出影象上的所有畫素值。
你可能會覺得相關操作和卷積操作一樣，但二者實際上是有區別的。
卷積操作
影象卷積操作的公式：g = f ★ h，具體而言，對應到每個畫素的表示式如下所示，可以看出，卷積操作滿足交換律。
In Convolution operation, the kernel is first flipped by an angle of 180 degrees and is then applied to the image.
也就是說， 卷積操作與相關操作的差異在於卷積將核kernel旋轉了180度（順逆都一樣）；相關可以反應兩個訊號的相似度，卷積不可以；卷積滿足交換律，相關不可以；卷積可以直接通過卷積定理（時域上的卷積等於頻域上的乘積）來加速運算，相關不可以。
參考文章：https://towardsdatascience.com/convolution-vs-correlation-af868b6b4fb5
MOSSE
基本思想
MOSSE的基本思想：以上一幀目標位置為bbox，在當前幀影象上擷取目標，使用相關核找到當前擷取影象上的最大響應，這裡的最大響應就是當前目標的中心，然後更新當前幀目標的中心，並通過第一幀給出的bbox的寬高為寬高，生成此時物體的框體。（這裡有一個假設，就是相鄰幀物體不會發生過大的位移，即使用上一幀的bbox擷取當前幀的影象，仍可以獲得到目標的中心） --- 可以看出，MOSSE對於物體大小的變化，以及物體的消失沒有絲毫抵抗能力。
具體操作流程
通過基本思想可以看出，需要首先獲得一個可以在目標中心獲得最大響應的相關核h。即下述公式操作後，g中最大值是物體的中心。
由於上式在計算機中計算量較大，作者對上式採用快速傅立葉變換(FFT):
便得到論文中給出的
因此，跟蹤的任務便是找到相應的H*
在實際跟蹤中，需要考慮目標外觀的變化等因素的影響，需要進行一個優化的求解，即輸出的響應與期望的響應越接近越好：
該優化的求解可以根據偏導數進行求取，求取結果可以直接表達為：
以上的操作均是在跟蹤第一幀完成的。也就是說，相關核H的確定是通過第一幀標出的目標線上得出的。而上式中的i，對應到程式碼中即為若干次影象的仿射變換。而對應的輸出響應的期望G*是高斯核，後續會結合程式碼具體解讀。
為了提升演算法的魯棒性，應對目標在跟蹤過程中的外觀變換，需要對核進行線上的更新，更新策略如下所示：
值得一提的是，F、G、H的尺寸都是一致的。並且，在工程實踐中確實將H分解為A，B兩部分，並進行迭代更新。
程式碼解讀
初始化
初始化應用於第一幀的功能實現，包括初始bbox資訊的獲取（中心點、寬、高、特徵），以及核的生成。其中，包含影象的預處理和準備工作：
- 灰度化（輸入特徵為灰度，特徵表達能力有限）
- 歸一化
- 餘弦窗的生成
- 期望輸出響應G的生成（高斯核，同時讓圖片的左邊緣與右邊緣連線，上邊緣與下邊緣連線）
之後，便是最終要的核的生成。程式碼如下：
1 def init(self,first_frame,bbox): 2 if len(first_frame.shape)!=2: 3 assert first_frame.shape[2] == 3 4 first_frame=cv2.cvtColor(first_frame, cv2.COLOR_BGR2GRAY) 5 first_frame=first_frame.astype(np.float32)/255 # 歸一化 6 x,y,w,h=tuple(bbox) # x,y為bbox的左下角點，和寬高 7 self._center=(x+w/2, y+h/2) # 中心點座標 8 self.w,self.h=w,h 9 w,h=int(round(w)), int(round(h)) # 取整 10 self.cos_window=cos_window((w,h)) # 生成寬高的一維hanning窗 後再進行矩陣乘法 形成h*w的hanning矩陣 11 self._fi=cv2.getRectSubPix(first_frame,(w,h), self._center) # 根據寬高和中心點的位置擷取影象/特徵 12 self._G=np.fft.fft2(gaussian2d_labels((w,h), self.sigma)) # 高斯相應圖 13 self.crop_size=(w,h) # 裁剪尺寸 14 self._Ai=np.zeros_like(self._G) 15 self._Bi=np.zeros_like(self._G) 16 for _ in range(8): 17 fi=self._rand_warp(self._fi) # 對裁剪的影象進行8次仿射變化(旋轉)，從而求得H* 18 Fi=np.fft.fft2(self._preprocessing(fi,self.cos_window)) # 進行餘弦窗處理 19 self._Ai += self._G*np.conj(Fi) # conj共軛負數 20 self._Bi += Fi*np.conj(Fi) # 求解H*所需
其中，程式碼第4-5行：影象轉化為灰度圖；
程式碼第6-9行：獲取第一幀影象的資訊，將左下角點，寬高資訊轉化為中心點座標。（bbox是通過cv2.selectROI()獲得到目標的左下角點座標和寬高資訊。）
程式碼第10行：獲取餘弦窗；
程式碼第11行：通過中心點座標和寬高在第一幀影象上擷取目標影象；
程式碼第12行：獲得高斯響應作為期望輸出G。
程式碼第16-20行：通過8次仿射變換，獲取相關核H。
餘弦窗的生成
餘弦窗通過漢寧窗生成。
1 def cos_window(sz): 2 cos_window = np.hanning(int(sz[1]))[:, np.newaxis].dot(np.hanning(int(sz[0]))[np.newaxis, :]) 3 return cos_window
漢寧窗的模樣
通過下圖餘弦矩陣的顏色分佈可以看出，餘弦窗具有兩邊小中間大的特性，並且四周的值趨近於零，有利於突出靠近中心的目標。
期望響應輸出G的生成
G相當於是真值，但此處不需要人為的標註。而是通過高斯函式生成即可。高斯函式產生的最大值在影象的中心。生成後，需要進行傅立葉變換轉化到頻域。
1 def gaussian2d_labels(sz,sigma): 2 w, h=sz 3 xs, ys = np.meshgrid(np.arange(w), np.arange(h)) # 生成兩個矩陣 一個矩陣各行是0-w-1 一個矩陣各列是0-h-1 4 center_x, center_y = w / 2, h / 2 5 dist = ((xs - center_x) ** 2 + (ys - center_y) ** 2) / (sigma**2) 6 labels = np.exp(-0.5*dist) 7 return labels
為什麼可以用高斯函式產生的值來作為期望的響應輸出呢？這是由於高斯函式自帶中間大兩邊小的特性，十分符合物體中心是最大值的要求。並且在第一幀框定目標時，也會標出十分貼合的bbox，bbox的中心可以近似為物體的中心。
隨機仿射變換
通過隨機數的生成控制仿射變換的力度。
1 def _rand_warp(self,img): # 輸入是裁剪下的圖形 2 h, w = img.shape[:2] 3 C = .1 4 ang = np.random.uniform(-C, C) 5 c, s = np.cos(ang), np.sin(ang) 6 W = np.array([[c + np.random.uniform(-C, C), -s + np.random.uniform(-C, C), 0], 7 [s + np.random.uniform(-C, C), c + np.random.uniform(-C, C), 0]]) 8 center_warp = np.array([[w / 2], [h / 2]]) 9 tmp = np.sum(W[:, :2], axis=1).reshape((2, 1)) 10 W[:, 2:] = center_warp - center_warp * tmp # 仿射矩陣 11 warped = cv2.warpAffine(img, W, (w, h), cv2.BORDER_REFLECT) # 進行仿射變換 W是仿射矩陣 12 return warped
線上更新
線上更新核心操作在於相關核的更新（下述程式碼第6行）以及影響相關核的兩個因素的更新（下述程式碼22-23行）。
從程式碼14-18行以及返回值可以看出，線上更新的過程只有中心點位置的更新，寬高不更新，如果物體的大小發生明顯的變化，很難自適應。
1 def update(self,current_frame,vis=False): 2 if len(current_frame.shape)!=2: 3 assert current_frame.shape[2]==3 4 current_frame=cv2.cvtColor(current_frame,cv2.COLOR_BGR2GRAY) 5 current_frame=current_frame.astype(np.float32)/255 6 Hi=self._Ai/self._Bi # Ht kernel 7 fi=cv2.getRectSubPix(current_frame,(int(round(self.w)),int(round(self.h))),self._center) # 裁剪 8 fi=self._preprocessing(fi,self.cos_window) 9 Gi=Hi*np.fft.fft2(fi) # fft 10 gi=np.real(np.fft.ifft2(Gi)) # 返回負數型別引數的實部， 以及二維傅立葉反變換 11 if vis is True: 12 self.score=gi 13 curr=np.unravel_index(np.argmax(gi, axis=None),gi.shape) # 找到最大那個相應的位置 14 dy,dx=curr[0]-(self.h/2),curr[1]-(self.w/2) # 中心點位置移動量 15 x_c,y_c=self._center 16 x_c+=dx 17 y_c+=dy # 中心點更新 18 self._center=(x_c,y_c) 19 fi=cv2.getRectSubPix(current_frame,(int(round(self.w)),int(round(self.h))),self._center) 20 fi=self._preprocessing(fi,self.cos_window) 21 Fi=np.fft.fft2(fi) 22 self._Ai=self.interp_factor*(self._G*np.conj(Fi))+(1-self.interp_factor)*self._Ai 23 self._Bi=self.interp_factor*(Fi*np.conj(Fi))+(1-self.interp_factor)*self._Bi 24 return [self._center[0]-self.w/2,self._center[1]-self.h/2,self.w,self.h]
缺點
- 輸入的特徵為單通道灰度影象，特徵表達能力有限
- 沒有尺度更新，對於尺度變化的跟蹤目標不敏感