基于GPU的二維矢量場LIC算法研究

陳 丁，萬 剛，王 龍，趙宸立

（1．信息工程大學，河南 鄭州450052；2．總參第一測繪導航基地大連測繪大隊，遼寧 大連116000）

矢量場可視化技術在科學計算和工程分析中有著重要意義，其目的就是采用計算機圖形（圖像）來同時表示場中每點的方向和大小信息［1］。通過可視化，把大量數據轉變成直觀的圖形圖像信息，并提供交互控制，以輔助信息挖掘與理解，進行有效分析。

目前，矢量場的可視化方法主要包括4類［2］：直接映射可視化、基于幾何的可視化、基于紋理的可視化和基于特征的可視化。基于紋理的可視化將矢量場的方向變化等細節信息利用紋理表達，綜合了直接映射和基于幾何的可視化方法的優點，在矢量場表達研究中占據著重要地位。另一方面，隨著科學計算產生的數據量越來越大，科研人員希望快速準確地可視化大規模數據，并實現針對可視化結果的實時交互與分析，對矢量場可視化提出了更高的要求。線積分卷積方法是基于紋理的可視化方法中最常用的一種。其成像效果好但帶來大量的計算，本文針對LIC在二維矢量場表達中的應用進行分析并提出改進措施，考慮矢量場的相關性并對流線跟蹤計算進行優化。

1 LIC算法及相關研究

1．1 LIC算法原理

LIC通過紋理內紋素間的相關性來表達矢量場的方向信息，它由Brian Cabral和Leit h Leedo m［3］在SIGGRAPH’93上提出。其基本思想是通過矢量場中任一點的流線，沿流線方向按給定的卷積核函數進行積分，積分結果作為最終的輸出，從而能夠表示矢量場的數據特征（方向等），如圖1所示。

圖1 LIC算法原理示意圖

LIC算法把矢量場數據集定義成笛卡爾網格下的規則結構，將矢量場數據集和相同分辨率的白噪聲紋理作為輸入。針對矢量場中任一點（x，y），其局部特征由以該點為中心的流線段描述。運用1D濾波卷積核，沿正反兩個方向進行積分卷積得到每個像素的輸出紋理。采用白噪聲作為輸入保證了輸入像素之間的不相關，通過沿流線方向進行卷積，從而合成圖像具有了矢量場方向相關性。

首先，進行流線的跟蹤計算。從點（x，y）出發，沿正、反兩方向分別追蹤生成長度為L（預先定義的流線段長度）的流線。正反兩方向的長度要一致，當遇特殊情況（流線跟蹤到區域邊界或矢量大小為0）需要提前結束時，兩方向仍要截取相同長度，以保證矢量場特征不變形。然后，進行積分卷積計算。LIC常采用盒式卷積核函數作為卷積核，將正、反兩方向流線上的輸入噪聲紋理按卷積核函數進行卷積運算。由此可得到點（x，y）對應的輸出紋理值Fout（x，y）為

式中：Fin（p）為點（xi，yi）處的輸入像素值；si為第i步的流線長度；Δsi為第i步流線的積分步長；l和l′為正、反方向的流線長度；k（ω）為盒式卷積函數。

LIC算法能夠表達矢量場的整體和局部連續信息，當卷積函數為周期函數時，通過變換相位還可實現矢量場的動態效果。

1．2 提高LIC效率的研究

利用LIC算法能夠很好地描述矢量場，但它還存在很多不足。LIC的流線計算和卷積計算過程要針對紋理中的每個像素進行計算，需要消耗大量的計算時間。矢量場的表達效果嚴重依賴輸出紋理圖像的分辨率，算法的效率很低，遠遠不能滿足實時可視化的需求。許多基于LIC的改進算法相繼提出，致力于提高LIC的效率和面向三維矢量場的擴展。D．Stalling和 H．C．Hege［4］對LIC算法做了很大的改進，提出快速卷積法，利用流向上像素間共用流線，減少重復計算，提高LIC算法的速度和精度。在此基礎上，Zockler［5］引入并行算法，充分利用幀間紋理的時間相關性，實現矢量場的動畫效果。秦勃［6］等利用異構算法，解決計算的耗時問題，但數據的預處理與繪制顯示仍然存在不平衡問題。

本文提出利用GPU實現LIC的加速計算，從以下幾個方面進行改進：首先，為提高矢量場的表達效果，矢量場數據的插值采用反距離加權方法；其次，流線計算在整個算法中耗費時間最多，通過用戶的自主選擇或矢量場的雷諾數計算，從而實現流線的分級生成；最后，依據矢量場的變化率來決定積分步長，建立自適應步長積分。

2 基于GPU的加速LIC改進

LIC算法中每個像素進行流線追蹤和積分卷積計算是相互獨立的，最終輸出結果只影響當前像素，這說明LIC算法的并行實現十分可行。通過GPU進行實現，將數據存儲在顯存大大減少數據的I／O操作，并且不需要對圖像像素進行顯式遍歷，片元處理器的操作對象即是逐像素的。

如圖2所示，將離散矢量數據場進行顏色映射生成連續的矢量紋理，將矢量紋理和噪聲紋理同時傳入GPU，利用矢量場對噪聲紋理進行積分卷積計算。GPU加速部分通過Open GL實現，對輸出紋理進行邊緣檢測，使矢量場的整體特征更加明顯。

圖2 基于GPU的LIC算法可視化流程

2．1 數據插值

在LIC算法中，流線的計算是十分重要的一個過程，在求解流線時，必須要確定采樣點的位置，通過插值計算得到采樣點的矢量信息。GPU自帶的顏色插值功能采用雙線性插值方法，將單元格內矢量信息的變化看作是線性的。其優點是方法簡單，運算速度快，但沒有考慮矢量場的相關性，當矢量信息變化較大時誤差不能忽略。考慮矢量場的局部強相關性，本文采用反距離加權法進行矢量場的插值計算替代GPU的固有算法，既保證矢量場的局部相關性，運算量也適中［8－9］。反距離加權算法的實現原理如下：設矢量場中的一點P，其矢量值與周圍4個采樣點的相關性由P到各點的距離決定，各點的權重值與到P點距離成反比。P點的矢量值由式（3）給出。

其中，

2．2 流線生成

流線的構造一般采用數值積分的形式給出，常用的有一階歐拉法、二階Runge－Kuta和四階Runge－Kuta積分法。通常情況下，數值積分的階數越高計算的精度也越高，但對應的計算量也成倍增加。采用精度高的積分方法，使得可視化的效果較好，但可視化的實時性將會大打折扣。在積分方法的選取上，必須折中考慮可視化的精度和效率。

為了得到較好的效果，本文采用三種積分方法的融合方式，通過矢量場的雷諾數來進行積分公式的選擇。雷諾數是反映流場形態的一個無量綱參數，當雷諾數較小時，流場的變化緩慢，采用一階歐拉法得到的流線即可滿足精度要求；當雷諾數較大時，流場的變化紊亂，需要采用高階數值積分公式以獲取更好的可視化效果。可按式（5）進行積分方法的選擇。

其中：n表示選擇積分公式的階數；Re表示雷諾數。

另外，通過數值積分生成的流線并不是光滑的，它由分段折線段構成。傳統的LIC算法采用固定的積分步長來生成流線，對于變化劇烈的地方，流線積分無法精確表達；而對于變化平緩的地方，又會造成大量計算的浪費。針對這一問題，Buning［7］提出自適應步長積分方法，認為積分步長受網格單元和速度大小影響，但該方法計算較為復雜，在GPU實現中存在一定限制。因此，本文采用流線上兩相鄰采樣點的速度方向變化來確定積分步長，算式為

依據夾角的大小，設置積分步長的選擇算式為

式中：Δs為預設的積分步長，θ1，θ2為夾角閾值。

3 算法實現及分析

基于GPU的LIC算法實現中，將矢量場數據作為輸入紋理，本文采用實驗數據作為中國東南部海域海流數據和臺風數據，需要對矢量場數據建立顏色映射關系，如式（8）所示。

式中，紋理的RGB顏色通道分別和速度矢量的方向分量uxuy及速率v相對應，vminvmax為速度最大最小值。因矢量方向分量存在負值，將其轉換到（0，1）區間存儲在紋理空間，顏色映射完成后的紋理數據如圖3所示。

圖3 矢量數據分區映射紋理

基于GPU加速主要通過可編程片元著色器實現，主要工作包括數據插值、流線計算和卷積積分。

通過Open Gl＋GLSL來實現，算法實現代碼如下：

Sampler2D FlowI mage；／／矢量場紋理數據

Sampler2D NoiseI mage；／／噪聲紋理

Unif or m Lengt h；／／積分長度

vec2 f Pos＝tex Coor d；／／點定位

vec2b Pos＝tex Coor d；

vec2f Pre＝f Pos；

vec2b Pre＝b Pos；

vec2 f Point List［lengt h］，b Point List［length］；／／采樣列表

int n＝ Re Cal（tex Coor d）；／／積分方法選擇

f or（int k＝0，k＜lengt h，k＋＋）

｛

Strea mLine Cal（n，f Pos）；／／流線計算

Angle Cal（f Pre，f Pos）；／／夾角計算

f wd［k］＝f Pos；

f Pre＝f Pos；

Strea mLine Cal（n，b Pos）；

Angle Cal（b Pre，b Pos）；

b wd［k］＝b Pos；

b Pre＝b Pos；

｝

圖4為實驗效果圖，圖4（a）為海流數據的可視化效果，通過對LIC紋理進行邊緣檢測操作，使矢量場的平流特征更加明顯。圖4（b）為臺風數據場可視化效果，依托已有平臺，將生成LIC紋理圖像映射到地球曲面顯示，并添加光照條件的渲染結果。

圖4 LIC算法可視化效果圖

通過實驗可以得出：基于GPU加速的LIC算法能夠取得較好的渲染效果，通過并行性提煉，保證繪制速度；進一步的算法優化使可視化效果基本滿足實時繪制的需求。

4 結束語

LIC算法作為紋理可視化方法的經典之一，在矢量場的可視化技術中占據著重要地位。本文所實現的基于GPU的算法降低了對硬件的要求，具有一定通用性，初步達到矢量場實時交互可視化的可能性，但還存在一定問題。下一步將研究LIC的優化算法：數據的壓縮、傳遞和基于GPU集群的并行實現，以適應更多的需求。

［1］ 唐澤勝．三維數據場可視化［M］．北京：淸華大學出版社，1999．

［2］ 李小梅，黃朝輝．科學計算可視化導論［M］．北京：國防科技大學出版社，1996．

［3］ CABRAL B，LEEDOM L C．Imaging vector fields using line integral convolution［C］／／Proceedings of the 20th annual conference on Computer graphics and interactive techniques．ACM，1993：263－270．

［4］ STALLING D，HEGE H C．Fast and resolution independent line integral convolution［C］／／Proceedings of the 22nd annual conference on Co mputer graphics and interactive techniques．ACM 1995 249－256．

［5］ Z?CKLER M，STALLING D，HEGE H C．Parallel line integral convolution［J］．Parallel Co mputing，1997，23（7）：975－989．

［6］ QIN B，WU Z，SU F，et al．GPU－based parallelization algorith m for 2D line integral convolution［M］／／Advances in Swar m Intelligence．Springer Berlin Heidelberg，2010：397－404．

［7］ POST F H，VROLIJK B，HAUSER H，et al．The state of the art in flow visualisation：Feature extraction and tracking［C］／／Computer Graphics Forum．Black well Publishing，Inc，2003，22（4）：775－792．

［8］ 焦永清，李斌，張坤．基于空間可視化的氣溫插值方法比較［J］．測繪工程，2013，23（5）：17．

［9］ 侯溯源，安曉亞，許劍，等．地理信息可視化新技術綜述與分析［J］．測繪與空間地理信息，2014，37（1）：30－32．