標籤：

Farneback 光流演算法詳解與 calcOpticalFlowFarneback 源碼分析

04-17

真實的三維空間中，描述物體運動狀態的物理概念是運動場。三維空間中的每一個點，經過某段時間的運動之後會到達一個新的位置，而這個位移過程可以用運動場來描述。

而在計算機視覺的空間中，計算機所接收到的信號往往是二維圖片信息。由於缺少了一個維度的信息，所以其不再適用以運動場描述。光流場（optical flow）就是用於描述三維空間中的運動物體表現到二維圖像中，所反映出的像素點的運動向量場。

  • 當描述部分像素時，稱為：稀疏光流
  • 當描述全部像素時，稱為：稠密光流

光流法是利用圖像序列中的像素在時間域上的變化、相鄰幀之間的相關性來找到的上一幀跟當前幀間存在的對應關係，計算出相鄰幀之間物體的運動信息的一種方法。光流法理解的關鍵點有：

  • 核心問題：同一個空間中的點，在下一幀即將出現的位置
  • 重要假設：光流的變化（向量場）幾乎是光滑
  • 角點處的光流能夠通過角點的鄰域完全確定下來，因此角點處的運動信息最為可靠；其次是邊界的信息

光流法有著各種各樣的分支，本文介紹的則是一種被廣泛使用的經典稠密光流演算法：Farneback 光流演算法

一、OpenCV 主函數源碼解讀：

源代碼中為了解決孔徑問題，主函數中引入了圖像金字塔。

孔徑問題（Aperture Problem）：

  • blog.csdn.net/hankai102
  • 形象理解：從小孔中觀察一塊移動的黑色幕布觀察不到任何變化。但實際情況是幕布一直在移動中
  • 解決方案：從不同尺度（圖像金字塔）上對圖像進行觀察，由高到低逐層利用上一層已求得信息來計算下一層信息

主函數 calcOpticalFlowFarneback 中需要的變數參數包括：

  1. _prev0：輸入前一幀圖像
  2. _next0：輸入後一幀圖像
  3. _flow0：輸出的光流
  4. pyr_scale：金字塔上下兩層之間的尺度關係
  5. levels：金字塔層數
  6. winsize：均值窗口大小，越大越能 denoise 並且能夠檢測快速移動目標，但會引起模糊運動區域
  7. iterations：迭代次數
  8. poly_n：像素鄰域範圍大小，一般為 5、7 等
  9. poly_sigma：高斯標準差，一般為 1～1.5（函數處理中需要高斯分布權重）
  10. flags：計算方法，包括 OPTFLOW_USE_INITIAL_FLOW 和 OPTFLOW_FARNEBACK_GAUSSIAN

實際的 Farneback 演算法在每一層金字塔上的處理過程為：

輸入的圖像默認為灰度圖片，默認只有亮度值

圖像輸入與輸出時進行的高斯模糊操作（Gaussian Blur）操作都是使得光流場結果平滑，從而滿足假設：光流的變化幾乎是光滑的

所以需要關注的子函數有 4 個：

  1. FarnebackPolyExp
  2. FarnebackUpdateMatrices
  3. FarnebackUpdateFlow_GaussianBlur
  4. FarnebackUpdateFlow_Blur

二、OpenCV 子函數 FarnebackPolyExp：

2.1. 理論基礎

2.1.1 圖像建模

將圖像視為二維信號的函數（輸出圖像是灰度圖像），因變數是二維坐標位置,

並且利用二次多項式對於圖像進行近似建模的話，會得到:

其中，A是一個2×2的矩陣，b是一個2×1的矩陣

因此係數化之後，以上公式等號右側可以寫為：

2.1.2 求解空間轉換

如果將原有（笛卡爾坐標系）圖像的二維信號空間，轉換到以 (1,x,y,x^2,y^2,xy)作為基函數的空間，則表示圖像需要一個六維向量作為係數，代入不同像素點的位置 x,y 求出不同像素點的灰度值。

Farneback 演算法對於每幀圖像中的每個像素點周圍設定一個鄰域(2n+1)×(2n+1)，利用鄰域內的共(2n+1)^2個像素點作為最小二乘法的樣本點，擬合得到中心像素點的六維繫數。因此對於圖像中的每個像素點，都能得到一個六維向量。

在一個鄰域內灰度值的 (2n+1)×(2n+1) 矩陣中，將矩陣按列優先次序拆分組合成 (2n+1)^2×1 的向量f，同時已知 (1,x,y,x^2,y^2,xy)作為基函數的轉換矩陣 B 維度為 (2n+1)^2×6（也可以視為 6 個列向量 bi 共同組成的矩陣），鄰域內共有的係數向量r 維度為 6×1，則有：

2.1.3 權重分配

利用最小二乘法求解時，並非是鄰域內每個像素點樣本誤差都對中心點有著同樣的影響力，函數中利用二維高斯分布將影響力賦予權重。 在一個鄰域內二維高斯分布的 (2n+1)×(2n+1) 矩陣中，將矩陣按列優先次序拆分組合成 (2n+1)^2×1 的向量 a。因此原本的基函數的轉換矩陣 B 將變為：

2.1.4 對偶轉換（並不清楚轉換的具體過程，只能參考論文中已有的公式）

為了「進一步加快求解」得到係數矩陣r，博士論文 4.3 小節中提出使用對偶的方式再次轉換基函數矩陣B，此時的對偶轉換矩陣為 G，經過轉換後的基函數矩陣的列向量為 bi~。博士論文中 G 的計算方式為：

而通過對偶轉換之後，計算係數向量r 的方式就簡單了很多，其中 * 表示兩個信號的互相關（實質上類似於兩個函數的卷積 zh.wikipedia.org/wiki/互相關 ）

過程：

2.1.5 函數輸出

子函數 FarnebackPolyExp 輸出得到的是單張圖像中每個像素點的係數向量 r（不包括常數項係數 r1，因為之後的計算光流過程中沒有用到）

2.1.6 Separable Normalized Convolution

博士論文 4.4 小節中提出的 Separable Normalized Convolution 計算方法提出將卷積操作由一維的直接計算拆分成兩個維度的分別計算，可以降低計算複雜度，拆分的依據源自：

2.2. 源碼解讀

2.2.1【基於鄰域】產生二維高斯分布的基礎是一維高斯分布，一維高斯分布存儲於數組 g 中，並且進行了求和後歸一化：

2.2.2【基於鄰域】基於 Separable Normalized Convolution 方法，分別求解(1,x,x^2)上的權重分布為：

2.2.3【基於鄰域】求解對偶轉換矩陣 G，注意此時的二維高斯分布權重已經被按照列向量拉伸成為了一個(2n+1)^2×1的向量。而根據一維高斯分布數組 g 生成二維高斯分布權重也很簡單，嵌套兩層循環即可。循環時可以利用矩陣的對稱性減少計算量。並且因為鄰域中心點為 (0,0)(0,0) 點，所以循環求和時G(0,1)+=g[y]*g[x]*x這類計算在x和-x上相互抵消，結果必然為 0 無需計算：

注意：實際計算中 G 的特殊結構會使得 G^-1 的特殊結構，所以只需要保存逆矩陣中幾個特殊位置的元素即可。證明過程見博士論文的附錄 A

2.2.4 分別在 vertical 和 horizontal 兩個方向進行卷積計算，卷積結果分別存在 row 數組和 b1~b6 中，最終的係數向量存在於 drow 數組中（不需要保存常量係數 r1，因為不涉及後面子函數的計算過程）。計算時注意 x 和 -x 在計算時導致的正負性差異：

三、OpenCV 子函數 FarnebackUpdateMatrices

子函數 FarnebackUpdateMatrices 中需要的變數參數包括：

  1. _R0：輸入前一幀圖像（編號：0）
  2. _R1：輸入後一幀圖像（編號：1）
  3. _flow：已知的前一幀圖像光流場（用於 Blur 子函數的迭代，以及主函數中圖像金字塔的不同層數間迭代）
  4. _M：儲存中間變數，不斷更新
  5. _y0：圖像求解光流的起始行
  6. _y1：圖像求解光流的終止行

3.1. 理論基礎

利用 FarnebackPolyExp 函數分別得到了視頻前後幀中每個像素點的係數向量之後，則需要利用係數向量計算得到光流場。這個計算過程在博士論文 7.6 小節中有描述。

首先，每個像素點都有著初始位移（最開始設置為全 0 變數），將上一幀的初始位移增加到第一幀圖像上的像素點位置 x 上，得到此像素點在下一幀圖像上的大致位置 x~：

x~ 可能是浮點位置，則無法估計此位置的係數向量值。對待此問題的方法有兩種：

1）將初始位移四捨五入得到整數值（博士論文中採用）

2）利用二次線性插值（OpenCV 中的默認做法）

由此得到用於計算的中間變數A(x),Δb(x)：

如果涉及到尺度變換，例如圖像金字塔技術（子函數 FarnebackUpdateMatrices 中再次使用 multi-scale 的思路，但需要和主函數中的圖像金字塔區分開來）。還需要另外與尺度縮放矩陣 S(x) 有關。總之，子函數 FarnebackUpdateMatrices 的目的是為了得到中間變數 G 與 h：

3.2. 源碼解讀

3.2.1 二次插值得到新一幀圖像位置中的係數向量：

3.2.2 處理不同尺度上的縮放，這裡借鑒的是 multi-scale 思路，詳見博士論文的 7.7 小節，目的是為了提高本演算法的魯棒性：

3.2.3 計算中間變數 G 與 h，存儲到矩陣 M 中：

四、OpenCV 子函數 FarnebackUpdateFlow_GaussianBlur、FarnebackUpdateFlow_Blur

4.1. 理論基礎

根據光流的基本假設：光流的變化（向量場）幾乎是光滑的。因此利用中間變數 G與 h 求解光流場 dout 前，需要進行一次局部模糊化處理（主函數中 winsize 輸入變數控制），可選均值模糊 (FarnebackUpdateFlow_Blur)、高斯模糊 (FarnebackUpdateFlow_GaussianBlur)。對於模糊後的中間變數，可以直接求解光流場：

4.2. 源碼解讀

4.2.1 根據中間變數的元素，計算得到光流場存儲於 flow 數組。flow 數組也是主函數最後的輸出量：

4.2.2 因為模糊化操作可能會執行多次，需要調用子函數 FarnebackUpdateMatrices 更新中間變數矩陣 M：

附：函數使用舉例

以上代碼的運行效果截圖為：

參考資料：

  1. blog.csdn.net/zouxy09/a
  2. blog.csdn.net/carson200
  3. Farneback G. Two-frame motion estimation based on polynomial expansion[C]. scandinavian conference on image analysis, 2003: 363-370.
  4. Farneback G. Polynomial Expansion for Orientation and Motion Estimation[J]. Link?ping University Sweden, 2002.

推薦閱讀：

TAG:三維 |