基于近似最近鄰搜索的并行光流計算

楊昕欣，姜精萍

北京航空航天大學 電子信息工程學院，北京 100191

1 引言

1.1 光流計算簡介

光流（Optical Flow）的概念是Gibson于1950年首先提出的，它是由相機和場景相對運動產生的，表示的是圖像中各像素瞬時相對運動的速度情況。光流場是指圖像灰度模式的表面運動[1]。光流計算是當今計算機視覺的熱點之一，光流算法被廣泛應用于運動目標檢測與跟蹤[2]、機器人導航[3]、三維重建[4]等領域。20世紀80年代初期，Horn和Schunck[5]提出了建立在光流平滑性假設基礎上的稠密光流算法，為光流計算的發展起了奠基性作用，之后學者們相繼提出了多種計算光流的算法，Barron[6]等人根據數學方法和相應的理論基礎，將計算光流的方法分為基于梯度的方法、基于塊匹配的方法、基于能量的方法和基于相位的方法。

與其他計算光流的方法相比，基于塊匹配的光流法具有運算速度快、易于硬件實現的優點。塊匹配是圖像處理中很常見的一種問題，尤其在視頻編碼中的運動估計中應用廣泛，許多年來，學者們已經提出了許多塊匹配的算法，如雙十字搜索法[7]、六邊形搜索法[8]、十字-菱形搜索法[9]等。為了進一步提高運算速度，學者們提出了近似塊匹配算法，如局部搜索[10]、降低維度[11]等，雖然這些方法得到的結果只是精確匹配的近似，但由于其大大減少了運算量而被廣泛應用于各種高層圖像處理中。最近，Barnes[12]等提出了一種近似最近鄰域塊匹配算法，該算法利用了圖像的局部相關性，基于“相鄰塊的最近鄰域也是相鄰的”這一假設，將前面已經匹配好的結果傳遞給相鄰塊，與其他近似塊匹配算法相比，該算法利用了塊匹配的特性，極大地減少了運算量，本文將其應用于光流計算中。

1.2 Barnes近似塊匹配算法

Barnes近似最近鄰搜索算法利用了塊匹配結果的特點，即相鄰的查詢塊的匹配塊也有很大的可能性是相鄰的。同時Barnes算法還基于另一個假設：大范圍隨機的初始坐標也可能包含有某些塊的正確匹配塊的坐標。

該算法首先定義了一個最近鄰域場NNF（Nearest-Neighbor Field）。假設查詢圖像為A，參考圖像為B，每個圖像塊的塊坐標用其中心點的坐標表示。NNF是方程 f:A??2，定義域為A中所有塊可能的坐標，值域為匹配塊相對于查詢塊的坐標偏移量，如果查詢圖像A中坐標為a的查詢塊的最近鄰域塊為在參考圖像B中的坐標為b，則有：

給定某一特定距離計算函數E，E(f(a))表示查詢圖像中坐標為a的塊與參考圖像中坐標為a+f(a)的塊的誤差，一般用歐式距離。Barnes算法包括初始化和迭代修正兩個階段。

初始化是將NNF初始化為隨機的值或是使用先驗知識。例如，采用參考圖像B中全部區域的獨立均勻采樣偏移值，作為NNF的初始值。

迭代修正包括傳播和隨機搜索，這兩個步驟交替進行，對NNF進行修正。傳播：假設當前要計算的是點p(x,y)，可以利用 p的相鄰像素來改善 p的結果，因為對于相鄰的像素來說，其偏移量是差不多的。對于p來說，f(x-1,y)和 f(x,y-1)是已經計算好了的結果，可以用它們來修正當前塊 p的匹配結果 f(x,y)，即令

這樣做的意義在于，如果(x,y)的映射是正確的，但其相鄰區域?的映射不正確，這樣做之后，?區域右方和下方的點都能得到正確的映射值。進一步，在奇數次迭代中采取右下方向的傳播，在偶數次迭代中采取左上方向的傳播，這樣就解決了修正過程中的方向問題。隨機搜索過程是以當前查詢塊的最近鄰域值為中心，在一個窗口大小指數遞減的搜索窗內隨機搜索候選塊，利用這些候選塊來改善當前查詢塊的匹配信息。候選塊產生方式如下：

其中v0=f(x,y)，Ri為內的一個均勻采樣的隨機點，w為搜索窗大小，為圖像的寬度和高度中的最大值，α為搜索窗半徑衰減率，i=1,2,…直到wαi的值小于一個像素。

2 基于Barnes算法的塊匹配光流法

本文參考facebook在surround 360中的源碼，將Barnes算法應用于計算稠密光流時，有以下三個調整。

2.1 基于梯度的匹配準則

Barnes算法中的誤差計算函數D常用的是絕對誤差和準則（SAD），SAD在塊匹配中應用的最多，因為其不需要進行乘法計算，因此節省時間容易實現。SAD公式表示為：

其中塊大小為M×N，塊的運動矢量即為SAD最小時的 (u,v)，即

SAD雖然常用，但當場景中存在光照變化時，得到的光流場會有明顯的離散化，因此本文使用另一種匹配準則，基于梯度守恒的準則[13]。

圖像序列可以看成是圖像灰度值與時間、空間位置的函數，可以將圖像中像素點的灰度值表示為函數的形式：

使用泰勒公式將 f(x+u,y+v,t+Δt)展開，忽略高階無窮小項并認為時間無限趨近于0，得到：

將式（7）代入式（6），得到：

將式（8）代入式（5）中，得：

要求上式的最小值，可以用最小二乘法得到關于(u,v)T的非齊次方程組，然后可以得到新的匹配準則

其中m(i,j)代表窗口權重，表示了塊中每個像素點對匹配結果的影響，在實際計算中，一般用以ρ為半徑的高斯函數代替，并計算高斯函數與圖像序列梯度差值之間的卷積。可以看到，在新的匹配規則中，匹配塊與待匹配塊之間的差異只由圖像的梯度決定，而不是灰度。而光照會對圖像的灰度產生影響但不會對圖像的梯度產生影響，因此，采用基于梯度的匹配準則可以很好地避免光照帶來的不利因素，使得到的光流場更準確。

本文在計算稠密光流時，使用的塊大小為1，即m(i,j)=1，用G0x,G0y,G1x,G1y分別表示圖像G0,G1在x方向和y方向的梯度值矩陣。則式（10）可以簡化表示為：

為了使E(u,v)在定義域內連續可導，將式（11）等價地寫為下面的形式：

式（12）就是本文使用的匹配誤差計算公式。其中u,v表示該像素點的位移，為保證計算的精確度，將其類型定義為浮點數，而G1x,G1y中只有整數坐標有值，因此計算G1x(i+u,j+v)、G1y(i+u,j+v)時需要進行線性插值。例如，在G1x中，假設點(i+u,j+v)周圍最近的四個整數坐標分別為(m,n),(m+1,n),(m,n+1),(m+1,n+1)，令：

則根據雙線性插值原理，得到：

于是由式（12）～（14）便可以求解出每個點的匹配誤差。

2.2 多分辨率迭代

將圖像按金字塔式分層，是圖像處理中的一種常見的策略，經典的Lucas-Kanade光流算法后來就應用了這一策略[14]。采用由粗到精的金字塔算法，不僅可以減少運算量，還可以避免陷入局部點，增強魯棒性。在光流算法中，采用金字塔策略，可以有效獲得大位移光流。

具體方法為設圖像的原始尺寸為(W0,H0)，將其設為圖像金字塔的最底層（即第0層），定義縮放因子λ（小于1的常數，本文設為0.9），則第一層圖像尺寸為(λW0,λH0),依次類推，第 L層的圖像尺寸(WL,HL)=(λLW0,λLH0)，當圖像的寬度或高度小于某閾值MinSize時，停止分層。

設圖像金字塔共Lm+1層，光流的計算從最頂層（即第Lm層）圖像開始。每層的計算方法是利用上一層的光流場作為初始值，然后逐點進行傳播和修正。其中，最頂層圖像的光流場初始值為(0,1)之間的隨機均勻采樣值。以第L層圖像的光流計算為例：

設第L+1層的光流場為 fL+1,其尺寸為(WL+1,HL+1),由于第L層圖像的尺寸(WL+1/λ,HL+1/λ)，因此首先需要對 fL+1進行上采樣，為使結果盡量平滑，上采樣過程中使用雙三次插值[15]。設經上采樣后的光流場為那么第L層圖像光流場的初始值為：

乘上1/λ是因為圖像放大后，相應的像素點的運動位移也會變大。接下來對每個像素點進行傳播修正，用偽代碼描述如下：

其中，匹配誤差的計算公式見式（12），梯度下降法修正光流將在2.3節論述。至此，第L層圖像的光流場計算完畢。每一層都如此計算，當計算完第0層時，就得到了最終的光流場。

2.3 光流修正

Barnes算法的迭代修正過程中，傳播和隨機搜索交替進行，隨機搜索是在大范圍內設置搜索窗口來改善匹配結果，而視頻序列中連續的兩幀圖像的時間間隔一般很小，大部分像素的運動矢量也較小，這種情況下，對每個塊在大范圍隨機搜索來改善匹配結果的效率就不高。同時，由于采取了多分辨率迭代策略，已經可以有效獲得大位移光流，因此，在傳播后，可以不再用隨機搜索，而是使用梯度下降法來修正光流。

對于某個點來說，其匹配誤差可以看成是其運動矢量的函數（見式（12）），其真實的光流矢量就是使誤差函數最小時的運動矢量。

在式（12）中，G0x和 G0y與 u,v無關，而由式（13）、（14）可知，G1x和G1y是關于u,v的連續可導的函數，于是誤差函數E(u,v)也是關于u,v的連續可導的函數，由微積分原理可知：

注意上式僅在逼近點(u,v)的鄰域內成立，在實際計算中，用下式近似計算：

式（18）中，du和 dv取0.001，根據前面的分析，每次傳播后，修正光流矢量為：

由式（16）可知，修正后誤差總是會減小。式（19）就是機器學習中常用的梯度下降法。其中，α為步長，本文采用動態步長的策略，因為在接近圖像金字塔頂層時，誤差一般較大，可以將步長設置大一些，當接近金字塔底層時，步長應該小一些。將步長α可以設置為金字塔層數的函數，第L層的步長為：

綜上所述，整個算法流程用偽代碼描述如下：

//開始

讀入兩幀圖像G0,G1；轉換為灰度圖并做歸一化；

初始化第Lm層的光流場 fLm為（0，1）之間的隨機均勻采樣值；

對光流場 fL進行一次右下方向的傳播修正；

對光流場 fL進行一次左上方向的傳播修正；

對光流場 fL進行中值濾波；

if L＞0:

對 fL進行上采樣，作為下一層的初始值；

}

//結束

現對整個算法的復雜度進行簡要分析。對于n×n的原始圖像，總共處理的像素點個數為：

對于每層圖像的每個像素點，其傳播過程僅與相鄰的4個像素點有關，所用時間為常數；其光流修正過程與周圍像素點無關，所用時間也為常數；因此傳播和修正算法的時間復雜度僅為O(c?1)，其中c為常數。因此整個算法的時間復雜度為，即等效于平方復雜度O(n2)。

3 并行加速方法

圖形處理器GPU是由成千上萬的更小、更高效的核心組成的大規模并行計算架構，專為同時處理多任務而設計。如今，GPU的運算能力已經遠遠超過了CPU。GPU相比CPU在浮點運算能力，內存帶寬，以及價格上都體現了極大的優勢。

關于圖像的算法由于要處理數量龐大的像素點，通常比較耗時。例如，本文在計算每個像素點的光流時，要計算相鄰4個點的匹配誤差，用梯度下降法修正光流時，還要計算一次匹配誤差，而每次計算匹配誤差都要進行雙線性插值運算，因此每個像素點總的運算量也不小，隨著圖像尺寸的增大，總的運行時間會以快速增加。因此，對算法用GPU進行并行加速就非常有意義。

但從上面Barnes的算法描述可知，Barnes算法并不易于在GPU上并行實現，原因在于其傳播過程是以逐塊掃描式的，即當前塊(x,y)的計算依賴于(x-1,y)和(x,y-1)兩個塊的計算，這樣的串行性使得算法難以并行實現。

為此，本文對傳播過程稍作修改。每次傳播時，將當前塊的計算僅依賴于相鄰的一個塊，而不是之前的兩個塊，同時將好的結果也只傳播給相鄰的一個塊，這樣就可使得每一行或每一列的傳播變得獨立，因此很容易在GPU上并行實現。

具體來說，對每層圖像的光流場，由之前的傳播兩次變為傳播4次，依次是從左到右、從上到下、從右到左、從下到上。每次傳播時，每行或每列的計算是并行的。

該方法的可行性分析如下，假設其中的某個塊p獲得了最好的匹配結果，其傳播的過程如圖1所示。

圖1 傳播過程示意圖

圖1 中深色的塊表示當前傳播到的塊。可以看出，經過4次傳播，其好的結果能傳遍圖像中所有的點。

4 實驗分析

實驗將本文算法及本文并行加速的算法與OpenCV實現的稠密光流Farneback[16]算法和Dual TV L1[17-18]算法進行對比。

Farneback算法的基本思想是用一個二次多項式來對圖像進行近似建模：，其中A是一個2×2對稱矩陣，b是2×1的矩陣。對于每幀圖像中的每個像素點周圍設定一個鄰域(2n+1)×(2n+1)，利用鄰域內的共(2n+1)2個像素點作為最小二乘法的樣本點，擬合得到中心像素點的系數。該過程中使用了分離卷積的方式，有效地降低了運算量。對于某個像素點，假設其原灰度值為 f(x)，增加全局位移后的值為f(x+d)，由灰度一致性假設可知 f(x)=f(x+d)，然后根據多項式相應的系數相等，可以求解出該像素點的全局位移d。對每個像素點均計算后，就得到了圖像的稠密光流場。

Dual TV L1算法是在TV L1光流模型的基礎上提出的。TV L1光流模型是指基于L1范數的整體變分（Total Variation-TV）光流模型，它是用L1范數代替文獻[5]中HS光流模型的L2范數，因為基于L1范數的方程具有更優的間斷處理能力。但由于TV L1方程的數據項和平滑項沒有完全、連續的可微性，因此求解較困難。Zach[17]等為此提出了一種對偶的方法來求解TV L1光流模型，它將原模型數據項中的原光流變量(u,v)用對偶變量(u′,v′)替換，從而構造了TV L1的對偶模型（Duality TV L1）。同時將模型分為兩部分凸替求解(u,v)和 (u′,v′)。這樣求得的結果相比傳統的TV L1光流模型具有更好的全局最優性。

4.1 評估指標

光流場的評估理論和方法一直是學者們研究的一個重要主題。2011年，Baker[19]等對光流場的評估方法做了一個詳細的總結。本次實驗的評估指標采用常用的平均角誤差（Average Angle Error，AAE）和平均點誤差（Average Endpoint Error，AEE）。

AAE的計算公式為：其中，N表示圖像中像素點的總數，(ui,vi)表示計算得到的圖像某像素點的光流矢量表示圖像某像素點標準的光流失量。AAE反映了計算得到的光流矢量整體偏移標準光流矢量的角度，是評估光流算法準確性的重要指標。

平均點誤差(AEE)的計算公式為：

AEE反映的計算得到的光流場矢量端點偏離標準光流場矢量端點的距離，這是個絕對的量，有一定局限性，比如當標準光流場矢量很小時，這個值一般較小，而當標準光流場矢量較大時，這個值一般較大，為方便不同圖片的比較，本文還使用Mccane[20]定義的歸一化點誤差（Normalized Endpoint EEROR，NEE），其計算公式為：

其中c表示真實光流，e表示計算光流，T為顯著性閾值，設置T是為了避免小光流帶來的不利影響。本文使用歸一化化點誤差的平均值（ANEE）作為統計指標：

ANEE反映了點誤差占標準光流矢量大小的平均百分比。

4.2 實驗過程

實驗的硬件環境為：CPU為Core i7 3.4 GHz，RAM 32 GB，顯卡為NVIDIA GeForce GTX 1080。

實驗的軟件環境為：Windows 10，Visual Studio 2015，OpenCV3.1.0，cuda 8.0。

實驗的測試圖片采用Middlebury[21]網站提供的包含真實光流場的訓練組序列。訓練組序列包含自然場景序列和合成場景序列，自然場景序列包括Dimetrodon、Hydrangea、RubberWhale，其中包含了大位移運動，非剛性運動以及真實而多樣性的光照效果；合成場景序列包括Grove類序列以及Urban類序列，其中Grove類序列是由多個程序合成的石頭和樹組成，地面紋理和表面位移是隨機生成的，含有明顯的遮擋，攝像機還有自旋和3D運動。Urban類序列是模擬的城市環境，含有許多的運動間斷點，較大的位移，還有獨立的運動物體。目前所有的合成序列明顯超越了先前的Yosemite序列。

實驗時記錄各個算法在每組序列中的各項誤差指標，同時記錄各算法的運行時間，對于各算法的運行時間均采用多次實驗（10次）取平均的方法。另外，由于本文算法和本文并行加速的算法是隨機初始化的，故每次實驗各項誤差數據也會有些微不同，因此也采用了多次實驗（10次）取平均的方法記錄誤差數據。對于各算法得到的光流場，采用文獻[18]中使用的彩色圖表示法以方便直觀的對比，彩色圖中不同的顏色代表不同方向的光流矢量，顏色的深淺表示該點光流矢量的大小；各算法計算的光流場與標準的光流場的誤差用灰度圖表示，圖中黑色的點表示沒有誤差，灰度值越大表示該點的誤差越大。

4.3 實驗結果及分析

以自然場景序列Dimetrodon為例，實驗結果如圖2。

圖2 Dimetrodon序列測試結果

表1 不同算法的各項指標對比

圖2中第二行彩色圖為各算法所計算的光流場經可視化后的圖像，第三行灰度圖為各算法計算得到的光流場的誤差。詳細的測試數據如表1所示。

結合表1和表2可以看出，Farneback算法各項誤差最大，Dual TV L1算法的對于自然場景序列效果較好，而對于有大位移的合成序列，算法誤差較大，比如對于序列Urban2，其平均點誤差超過了3個像素長度。而本文基于Barnes算法實現的塊匹配光流法對于所有序列均表現良好，從表2的平均性能來看，其平均角度誤差不超過5°，平均點誤差不超過半個像素長度，歸一化點誤差約為10%。而經并行加速的算法的準確度與原算法基本一致，均明顯好于OpenCV實現的兩種稠密光流算法。另外，本文原算法采用的是基于梯度的匹配準則，對于光照變化的場景有良好的魯棒性，但對于有旋轉、縮放的場景，梯度守恒的假設便不再適用，因此本文算法對Grove3序列的準確性次于其他序列。

表2 不同算法在所有測試圖片中的平均誤差

從時間性能來看，Farneback算法最快，本文原算法次之，Dual TV L1算法最慢，而本文的并行加速算法相比原算法有大約2.5倍的加速比，速度與最快的Farneback算法十分接近。

綜合準確度和速度來看，本文的并行加速光流算法性能最佳。

5 結束語

本文在Barnes算法基礎上實現的稠密光流算法，使用了梯度作為匹配準則，采用了金字塔分層策略，使用了梯度下降法修正光流，每個像素點的時間復雜度僅為O(1)，并且其準確性明顯優于OpenCV實現的兩種稠密光流算法。本文對原算法的傳播過程進行改進，從原來的兩次傳播變為4次傳播，使算法易于在GPU上實現并行加速，由此得到的并行加速算法相比本文的原始算法，而速度快了約2.5倍，幾乎與OpenCV中最快的Farneback算法相當，而其各項誤差指標明顯優于Farneback算法。