基于GPU的二維矢量場(chǎng)LIC算法研究

陳 丁，萬(wàn) 剛，王 龍，趙宸立

（1．信息工程大學(xué)，河南 鄭州450052；2．總參第一測(cè)繪導(dǎo)航基地大連測(cè)繪大隊(duì)，遼寧 大連116000）

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

目前，矢量場(chǎng)的可視化方法主要包括4類(lèi)［2］：直接映射可視化、基于幾何的可視化、基于紋理的可視化和基于特征的可視化?；诩y理的可視化將矢量場(chǎng)的方向變化等細(xì)節(jié)信息利用紋理表達(dá)，綜合了直接映射和基于幾何的可視化方法的優(yōu)點(diǎn)，在矢量場(chǎng)表達(dá)研究中占據(jù)著重要地位。另一方面，隨著科學(xué)計(jì)算產(chǎn)生的數(shù)據(jù)量越來(lái)越大，科研人員希望快速準(zhǔn)確地可視化大規(guī)模數(shù)據(jù)，并實(shí)現(xiàn)針對(duì)可視化結(jié)果的實(shí)時(shí)交互與分析，對(duì)矢量場(chǎng)可視化提出了更高的要求。線積分卷積方法是基于紋理的可視化方法中最常用的一種。其成像效果好但帶來(lái)大量的計(jì)算，本文針對(duì)LIC在二維矢量場(chǎng)表達(dá)中的應(yīng)用進(jìn)行分析并提出改進(jìn)措施，考慮矢量場(chǎng)的相關(guān)性并對(duì)流線跟蹤計(jì)算進(jìn)行優(yōu)化。

1 LIC算法及相關(guān)研究

1．1 LIC算法原理

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

圖1 LIC算法原理示意圖

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

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

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

LIC算法能夠表達(dá)矢量場(chǎng)的整體和局部連續(xù)信息，當(dāng)卷積函數(shù)為周期函數(shù)時(shí)，通過(guò)變換相位還可實(shí)現(xiàn)矢量場(chǎng)的動(dòng)態(tài)效果。

1．2 提高LIC效率的研究

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

本文提出利用GPU實(shí)現(xiàn)LIC的加速計(jì)算，從以下幾個(gè)方面進(jìn)行改進(jìn)：首先，為提高矢量場(chǎng)的表達(dá)效果，矢量場(chǎng)數(shù)據(jù)的插值采用反距離加權(quán)方法；其次，流線計(jì)算在整個(gè)算法中耗費(fèi)時(shí)間最多，通過(guò)用戶(hù)的自主選擇或矢量場(chǎng)的雷諾數(shù)計(jì)算，從而實(shí)現(xiàn)流線的分級(jí)生成；最后，依據(jù)矢量場(chǎng)的變化率來(lái)決定積分步長(zhǎng)，建立自適應(yīng)步長(zhǎng)積分。

2 基于GPU的加速LIC改進(jìn)

LIC算法中每個(gè)像素進(jìn)行流線追蹤和積分卷積計(jì)算是相互獨(dú)立的，最終輸出結(jié)果只影響當(dāng)前像素，這說(shuō)明LIC算法的并行實(shí)現(xiàn)十分可行。通過(guò)GPU進(jìn)行實(shí)現(xiàn)，將數(shù)據(jù)存儲(chǔ)在顯存大大減少數(shù)據(jù)的I／O操作，并且不需要對(duì)圖像像素進(jìn)行顯式遍歷，片元處理器的操作對(duì)象即是逐像素的。

如圖2所示，將離散矢量數(shù)據(jù)場(chǎng)進(jìn)行顏色映射生成連續(xù)的矢量紋理，將矢量紋理和噪聲紋理同時(shí)傳入GPU，利用矢量場(chǎng)對(duì)噪聲紋理進(jìn)行積分卷積計(jì)算。GPU加速部分通過(guò)Open GL實(shí)現(xiàn)，對(duì)輸出紋理進(jìn)行邊緣檢測(cè)，使矢量場(chǎng)的整體特征更加明顯。

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

2．1 數(shù)據(jù)插值

在LIC算法中，流線的計(jì)算是十分重要的一個(gè)過(guò)程，在求解流線時(shí)，必須要確定采樣點(diǎn)的位置，通過(guò)插值計(jì)算得到采樣點(diǎn)的矢量信息。GPU自帶的顏色插值功能采用雙線性插值方法，將單元格內(nèi)矢量信息的變化看作是線性的。其優(yōu)點(diǎn)是方法簡(jiǎn)單，運(yùn)算速度快，但沒(méi)有考慮矢量場(chǎng)的相關(guān)性，當(dāng)矢量信息變化較大時(shí)誤差不能忽略?？紤]矢量場(chǎng)的局部強(qiáng)相關(guān)性，本文采用反距離加權(quán)法進(jìn)行矢量場(chǎng)的插值計(jì)算替代GPU的固有算法，既保證矢量場(chǎng)的局部相關(guān)性，運(yùn)算量也適中［8－9］。反距離加權(quán)算法的實(shí)現(xiàn)原理如下：設(shè)矢量場(chǎng)中的一點(diǎn)P，其矢量值與周?chē)?個(gè)采樣點(diǎn)的相關(guān)性由P到各點(diǎn)的距離決定，各點(diǎn)的權(quán)重值與到P點(diǎn)距離成反比。P點(diǎn)的矢量值由式（3）給出。

其中，

2．2 流線生成

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

為了得到較好的效果，本文采用三種積分方法的融合方式，通過(guò)矢量場(chǎng)的雷諾數(shù)來(lái)進(jìn)行積分公式的選擇。雷諾數(shù)是反映流場(chǎng)形態(tài)的一個(gè)無(wú)量綱參數(shù)，當(dāng)雷諾數(shù)較小時(shí)，流場(chǎng)的變化緩慢，采用一階歐拉法得到的流線即可滿(mǎn)足精度要求；當(dāng)雷諾數(shù)較大時(shí)，流場(chǎng)的變化紊亂，需要采用高階數(shù)值積分公式以獲取更好的可視化效果?？砂词剑?）進(jìn)行積分方法的選擇。

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

另外，通過(guò)數(shù)值積分生成的流線并不是光滑的，它由分段折線段構(gòu)成。傳統(tǒng)的LIC算法采用固定的積分步長(zhǎng)來(lái)生成流線，對(duì)于變化劇烈的地方，流線積分無(wú)法精確表達(dá)；而對(duì)于變化平緩的地方，又會(huì)造成大量計(jì)算的浪費(fèi)。針對(duì)這一問(wèn)題，Buning［7］提出自適應(yīng)步長(zhǎng)積分方法，認(rèn)為積分步長(zhǎng)受網(wǎng)格單元和速度大小影響，但該方法計(jì)算較為復(fù)雜，在GPU實(shí)現(xiàn)中存在一定限制。因此，本文采用流線上兩相鄰采樣點(diǎn)的速度方向變化來(lái)確定積分步長(zhǎng)，算式為

依據(jù)夾角的大小，設(shè)置積分步長(zhǎng)的選擇算式為

式中：Δs為預(yù)設(shè)的積分步長(zhǎng)，θ1，θ2為夾角閾值。

3 算法實(shí)現(xiàn)及分析

基于GPU的LIC算法實(shí)現(xiàn)中，將矢量場(chǎng)數(shù)據(jù)作為輸入紋理，本文采用實(shí)驗(yàn)數(shù)據(jù)作為中國(guó)東南部海域海流數(shù)據(jù)和臺(tái)風(fēng)數(shù)據(jù)，需要對(duì)矢量場(chǎng)數(shù)據(jù)建立顏色映射關(guān)系，如式（8）所示。

式中，紋理的RGB顏色通道分別和速度矢量的方向分量uxuy及速率v相對(duì)應(yīng)，vminvmax為速度最大最小值。因矢量方向分量存在負(fù)值，將其轉(zhuǎn)換到（0，1）區(qū)間存儲(chǔ)在紋理空間，顏色映射完成后的紋理數(shù)據(jù)如圖3所示。

圖3 矢量數(shù)據(jù)分區(qū)映射紋理

基于GPU加速主要通過(guò)可編程片元著色器實(shí)現(xiàn)，主要工作包括數(shù)據(jù)插值、流線計(jì)算和卷積積分。

通過(guò)Open Gl＋GLSL來(lái)實(shí)現(xiàn)，算法實(shí)現(xiàn)代碼如下：

Sampler2D FlowI mage；／／矢量場(chǎng)紋理數(shù)據(jù)

Sampler2D NoiseI mage；／／噪聲紋理

Unif or m Lengt h；／／積分長(zhǎng)度

vec2 f Pos＝tex Coor d；／／點(diǎn)定位

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）；／／流線計(jì)算

Angle Cal（f Pre，f Pos）；／／夾角計(jì)算

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為實(shí)驗(yàn)效果圖，圖4（a）為海流數(shù)據(jù)的可視化效果，通過(guò)對(duì)LIC紋理進(jìn)行邊緣檢測(cè)操作，使矢量場(chǎng)的平流特征更加明顯。圖4（b）為臺(tái)風(fēng)數(shù)據(jù)場(chǎng)可視化效果，依托已有平臺(tái)，將生成LIC紋理圖像映射到地球曲面顯示，并添加光照條件的渲染結(jié)果。

圖4 LIC算法可視化效果圖

通過(guò)實(shí)驗(yàn)可以得出：基于GPU加速的LIC算法能夠取得較好的渲染效果，通過(guò)并行性提煉，保證繪制速度；進(jìn)一步的算法優(yōu)化使可視化效果基本滿(mǎn)足實(shí)時(shí)繪制的需求。

4 結(jié)束語(yǔ)

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

［1］ 唐澤勝．三維數(shù)據(jù)場(chǎng)可視化［M］．北京：淸華大學(xué)出版社，1999．

［2］ 李小梅，黃朝輝．科學(xué)計(jì)算可視化導(dǎo)論［M］．北京：國(guó)防科技大學(xué)出版社，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］．測(cè)繪工程，2013，23（5）：17．

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