基于信息熵的流線質量評估和布局算法

張 沙，解利軍，桂立業(yè)，鄭 耀

浙江大學 航空航天學院，杭州 310027

1 引言

矢量場可視化是科學可視化的一個經(jīng)典分支，在氣象分析、航空航天等眾多領域具有無可替代的作用。矢量場可視化的目的是將矢量場仿真數(shù)據(jù)以一種視覺可見的方式正確表達出來，同時這樣表達出來的信息必須便于人們理解和分析矢量場的特征。流線可視化，能夠直觀地反映矢量場的特征，同時計算相對簡便快捷，一直是矢量場可視化分析中最重要的方法之一。

流線是從某個種子點出發(fā)，沿著該點流向的正向和反向，通過矢量場積分計算得到的曲線。種子點的選取對流線的質量具有非常關鍵的作用。流線種子點選取不當可能會引起視覺混亂、特征不明顯等一系列問題。

近年來，矢量場種子點布局算法呈現(xiàn)百家爭鳴的局面。2D流場中，Jobard等提出等距離種子點放置方法[1]，并經(jīng)過多次改進[2-4]。等距離種子點放置方法的目的是在矢量場中等距離地放置流線以避免視覺聚類。但這種方法不對矢量場結構進行分析，可能會造成矢量場中重要特征不明顯的問題。最遠種子點放置算法[5]在距現(xiàn)有流線最遠的地方放置種子點，保證矢量場中不會出現(xiàn)大塊空白的區(qū)域。這個算法的缺點與等距離種子點的缺點類似，不對矢量場結構進行分析可能錯過某些重要特征。基于特征的可視化方法主要關注矢量場中的主要特征提取，Helman等人[6]對2D 矢量場的拓撲進行分析，對臨界點進行檢測和分類，并以此作為基礎執(zhí)行可視化算法。Verma等人提出了基于拓撲的種子點放置方法[7]，這個算法首先分析矢量場得到關鍵點，然后利用預先定義的模板在關鍵點周圍進行布點，這樣可以保證重要特征的突出顯示。本文后面初始布點時用到的模板就與Verma等人的模板相同。3D流場中，Ye等人在[8]中將Verma等的方法擴展到三維流場中，并定義了三維流場的模板。Li和Shen提出了基于圖像的三維矢量場種子點布置方法[9]，有效地避免了流線雜亂和視覺的混亂。Xu等人提出了矢量場信息化的信息框架，以獨特的視角詮釋了流線生成的信息意義。國內(nèi)的王少榮等人[10]提出了一種基于出入度的種子點布局算法。

在已經(jīng)有了如此多的種子點選擇算法后，對種子點的質量評估算法卻寥寥無幾。本文提出了一種基于信息論的矢量場種子點評估和布局算法，并利用CUDA（Compute Unified Device Architecture，計算統(tǒng)一設備架構）對其做了不錯的性能加速[11]。

本文的核心思想是：矢量場中矢量變化越大的區(qū)域，其局部熵值越大，擁有重要特征的可能性越大（詳見2.1節(jié)）。初始布點，選擇擁有極大局部熵值的網(wǎng)格點，對其周圍進行模板布點，突出顯示矢量場的重要特征。然后利用現(xiàn)有流線擴散出一個中間矢量場，求解中間矢量場與初始矢量場的條件熵，來衡量初始矢量場中除中間矢量場中的共有信息外的未知信息的熵值，確定現(xiàn)有流線是否已經(jīng)能夠展示出初始矢量場中的足夠信息，從而確定是否要繼續(xù)在矢量場中加入新的種子點。新的種子點的放置采取基于重要性取樣的方法。本文提出的算法最終目的是以較少的流線展示出矢量場中大部分信息。

2 基于信息熵的流線評估和布局算法

2.1 程序流程

圖1給出了本文算法的程序流程圖。算法包括三個比較大的模塊，分別為基于信息熵的初始種子點布置、流線質量評估以及基于重要性種子點取樣?；谛畔㈧氐某跏挤N子點布置模塊對矢量場特征進行分析，找到其主要特征并生成流線。流線質量評估模塊對現(xiàn)有流線進行質量評估，以確定現(xiàn)有流線是否已經(jīng)展示出矢量場中足夠的信息，是否需要繼續(xù)添加流線?；谥匾苑N子點取樣的目的是根據(jù)條件熵向矢量場中添加新的種子點。

圖1 基于信息熵種子點評估和布局算法流程圖

后續(xù)的章節(jié)中，將對各模塊的具體算法實現(xiàn)進行介紹。

2.2 基于信息熵的初始布點

在本文中，將矢量及其所有鄰近網(wǎng)格點（本文網(wǎng)格點的鄰近點規(guī)模：2D數(shù)據(jù)為132個，3D數(shù)據(jù)為133個）作為一個隨機系統(tǒng)并求解該系統(tǒng)的熵值。矢量變化越大即矢量越混亂的系統(tǒng)，其熵值越大，含有重要特征的可能越大。利用這個方法選取的種子點，可以很好地抓住矢量場的主要特征。

2.2.1 概率分布函數(shù)求解

要求解網(wǎng)格點的局部熵值，首先要統(tǒng)計網(wǎng)格點及其所有鄰近點的方向區(qū)間直方圖，得到該矢量場的概率分布函數(shù)p(xi)。

對所有矢量，先按方向進行量化，將其分為若干個方向區(qū)間。

對于2D矢量，先將矢量表示為角度θ∈[0,360)，角度空間被劃分為60個區(qū)間，每個區(qū)間6°。對于矢量為(a,b)的網(wǎng)格點，如果a，b均不為0，則該矢量角度區(qū)間求法如式（2）：

1948年，信息論之父C.E.Shannon在論文[12]中提出了信息熵的概念，用來解決信息的量化度量問題。信息熵描述了隨機系統(tǒng)的不確定程度，越有序的系統(tǒng)，其熵值越小，不確定的信息也越少；反之越混亂的系統(tǒng)，熵值越大，不確定的信息越多。

對于一個離散變量X，有n個可能的值(x1,x2,…,xn)每個值出現(xiàn)的概率為(p1,p2,…,pn)，則X的離散熵為：

對于3D矢量，為了簡化計算，本文將z軸單獨劃分為6個區(qū)間，x軸和y軸的劃分與2D矢量相同。z軸區(qū)間值的計算方法為z×3/r+3其中，3D矢量的角度區(qū)間為z軸的區(qū)間值乘以xOy平面的區(qū)間值。3D數(shù)據(jù)的角度區(qū)間共有360個。

統(tǒng)計網(wǎng)格點及其所有鄰近點矢量的方向區(qū)間。假如落在某區(qū)間Ci內(nèi)的矢量個數(shù)為ni，鄰近點的總個數(shù)為N，則該角度區(qū)間出現(xiàn)的概率為pi=ni/N。

2.2.2 基于模板的種子點布置

按照2.1節(jié)所示，求取每一個網(wǎng)格點的熵值，所有點的熵值構成EM（熵值矩陣），選取所有局部熵值大于0.7倍熵最大值的極值點，作為初始布點的模板中心。本文采用的模板是由Verma等人在文[7]中首先提出來的所有模板的綜合體。對于2D數(shù)據(jù)，采用了菱形模板，共9個點，分別為4個頂點，4條邊的中點及1個中心點。對于3D數(shù)據(jù)采用了八面體模板，共27個點：1個中心點，6個頂點，8個面的中心點，12條邊的中點。

圖2（a）是Isabel的某個時間片的切片矢量圖，圖2（b）展示了根據(jù)本文的算法對該數(shù)據(jù)求解熵矩陣，并對熵極值進行模板布點后的流線。

圖2 Isabel數(shù)據(jù)切片的矢量圖及其對應的熵圖1）

2.3 流線質量評估

這一模塊主要包括兩個步驟：首先根據(jù)現(xiàn)有流線生成中間矢量場，然后根據(jù)初始矢量場與中間矢量場的條件熵量化評估流線對矢量場信息的反映程度。

2.3.1 條件熵

條件熵反映的是兩個隨機變量之間的不確定性。在本文中初始矢量場和中間矢量場的條件熵描述的是：在完全已知中間矢量場Y(x)的情況下，初始矢量場X(x)中除了與Y(x)共享的信息之外的信息的熵值。也就是，Y(x)中剩余信息的不確定性。

條件熵有兩種計算方法，如公式（3）和公式（5）：

本文中，全局條件熵的計算直接放置CPU中進行，沒有空間限制，采用公式（5）只需對所有網(wǎng)格進行一次遍歷，能獲得較高的時間性能。而在基于重要性布點過程中求解局部條件熵時，算法需要求解每一個網(wǎng)格點及其所有鄰居點的條件熵。如果采用串行算法，將面臨高昂的時間代價，所以將這一部分移動到GPU中利用CUDA進行加速計算。實驗表明，此時采用公式（3）將取得更好的時間性能。

2.3.2 求解中間矢量場

本文通過最小化能量方程（7）的方法來求解中間矢量場Y(x)。這種方法能夠保證整個矢量場被恰當覆蓋，同時，還能保證中間矢量場經(jīng)過光滑處理，并且光滑后流線上的點的中間值與真實值之間的差異最小。

本文應用的能量方程與參考文獻[13]中提出的外部力相同。具體的求解過程參見文獻[13]。上述能量方程的求解結果為：

其中μ=0.1（經(jīng)驗值），?f(x)=?(x)，dt自己設定，但須滿足CFL約束條件[14]。

圖3展示了經(jīng)過若干次的迭代后求解出的中間矢量（藍色）與初始矢量（紅色）的比較結果。只有藍色矢量的網(wǎng)格點表示初始矢量與中間矢量相同，在顯示時被覆蓋了。圖3（d）中，中間矢量場中絕大部分的網(wǎng)格點或者因為被流線覆蓋或者因為有效的矢量恢復，得到了與原始矢量近似或者相同的結果。

對于大規(guī)模的數(shù)據(jù)，此模塊耗費的時間非常長。為了減少此模塊的時間，本文將這一模塊放到了GPU上進行并行計算。

圖3 各步驟流線恢復出的中間矢量場

圖3中圖（a）顯示判斷網(wǎng)格點是否在流線上，其中紅色的網(wǎng)格點表示在流線上，在建立近似矢量場時，這些網(wǎng)格點的矢量值將等于初始值，其他藍色點均為0。圖（b）展示了根據(jù)（a）中流線生成的中間矢量場與初始矢量場的比較。圖（d）為迭代終止后的中間矢量場。圖（b）（c）（d）三幅圖反映了中間矢量場逐漸與原始矢量場重合的過程。

2.4 基于重要性取樣

初始布點完成以后，矢量場中大部分的重要特征已經(jīng)顯示出來。后續(xù)迭代的重要性取樣是為了更加強調(diào)這種特征顯示同時填充矢量場中的空白區(qū)域。這個迭代過程的結束條件是矢量場中所有的信息可見，即：矢量場的全局條件熵趨于收斂。

根據(jù)求解出的中間矢量場Y(x)和初始矢量場X(x)求解每個網(wǎng)格點的條件熵，更新CEM。對于某個網(wǎng)格點來說，該處的條件熵越大，則成為新一輪種子點的可能性也越大。新一輪的布點過程中，網(wǎng)格點(x,y)被選中的概率為：

其中，h(x,y)為點(x,y)處的條件熵。

所有網(wǎng)格點的布點概率構成了概率分布函數(shù)（PDF），然后利用反變換的方法選取新一輪的網(wǎng)格點。

2.5 終止條件

重要性取樣迭代進行，在每次循環(huán)中，將一定數(shù)量的種子點布置到矢量場中，直到初始矢量場與中間矢量場的條件熵值趨于收斂。對于不同的數(shù)據(jù)，每次循環(huán)放置的種子點個數(shù)也不相同，通常情況下，每次循環(huán)放置的種子點數(shù)為網(wǎng)格點數(shù)的平方根。

3 結果和討論

測試算法在64位Windows 7系統(tǒng)，Lenovo ideapad Y470，Intel core i5，3 GB內(nèi)存機器上運行。機器配備的顯卡為Nvidia GeForce GT525。

3.1 颶風Isabel矢量場切片

圖4展示了兩種不同的算法對Isabel數(shù)據(jù)集切片的可視化結果。圖4（a）是本文算法執(zhí)行基于信息熵的初始布點的結果。由圖中可以看出即使此時流場中只存在較少數(shù)量的流線，流場中的特征區(qū)域已經(jīng)被全部顯示出來。圖 4（a）中III部分存在顯著源或匯，結合圖 2（a）點圖標矢量方向，可知此處存在源，II和V部分存在比較明顯鞍面，II、IV和VI部分存在需要用戶重視的特征區(qū)域。這證明了本文的算法在初始布點后能夠很好地捕獲流場中主要特征區(qū)域。圖4（b）算法經(jīng)過兩次迭代后的中間結果。

圖4（c）展示算法最終結果圖。圖中，各個特征區(qū)域周圍的流線數(shù)目更多，特征區(qū)域得到加強顯示；而非特征區(qū)域流線數(shù)目相對較少但不影響用戶對流場的理解。圖4（d）是最遠種子點算法生成的流線。圖中幾個特征區(qū)域并不明顯且短流線較多，尤其對應于圖4（a）I區(qū)域位置的特征區(qū)域很難識別出來。與最遠種子點算法相比，本文提出的算法能夠更好的捕獲流場中的特征區(qū)域，并予以突出顯示。

圖（e）顯示了在算法迭代過程中，中間矢量場與初始矢量場的條件熵逐漸趨于收斂，矢量場中未被顯示的信息逐漸減少。

3.2 Rayleigh-Bénard數(shù)據(jù)集

圖6展示的是本文算法應用于Rayleigh-Bénard數(shù)據(jù)集的實驗結果。Rayleigh-Bénard對流現(xiàn)象（RBC）是指當流體的Rayleigh數(shù)達到某一特定數(shù)值時產(chǎn)生對流的現(xiàn)象[15]。RBC現(xiàn)象是一種在自然界和工業(yè)中普遍存在的熱流動現(xiàn)象。本文獲取矢量場數(shù)據(jù)的方法是利用CFD模擬計算兩無限平行板間的RBC對流。圖5是模擬采用的幾何模型示意圖。模擬網(wǎng)格采用127×34×40的結構網(wǎng)格。

圖5 Rayleigh-Bénard對流模擬的模型

圖6（a）中紅色線條是通過李明等人提出的算法提取的渦核[16]。

圖6（b）展示的是對矢量場經(jīng)過初始布點以后的流線圖。圖中可以看出，流線大多纏繞渦核，很好地完成了抓取特征區(qū)域的任務。

圖6（c）是流線最終放置結果。最終的結果圖中，在初始布點生成流線的基礎上增加了流線，特征區(qū)域得到進一步加強顯示，非特征區(qū)域也適當增加了較少流線，具有全局性。

圖6（d）顯示隨著流線數(shù)目增多，初始矢量場與中間矢量場的條件熵逐步減少直至收斂，最終完成整個評估。

3.3 算法性能分析

上面兩個算例中，Isabel單次評估和生成流線需要9.26 s，Rayleigh-Bénard數(shù)據(jù)需要46.88 s，算法的速度當前還達不到交互式可視化的要求。

本文算法主要分為3個模塊，分別為基于信息熵的初始布點、流線質量評估以及基于重要性取樣。以RB對流數(shù)據(jù)為例，算法各模塊的時間效率如表1。其中，更新LCM時，需要對每個網(wǎng)格點掃描鄰域，且計算過程中需建立一個360×360的矩陣，這遠遠超出了GPU 16K的共享內(nèi)存存儲范圍，需要采取以時間換取空間的手段，加大了計算量，耗費了大量的時間。

圖6 對流數(shù)據(jù)的可視化結果

4 結論

本文提出了一種種子點布局及流線評估方法。矢量場中熵值越大的區(qū)域，其含有特征點的可能性越大。初始布點中，選取了矢量場中若干個極值點進行基于模板的布點，最終生成的流線能很好地抓住矢量場的主要特征。利用初始矢量場與流線生成的中間矢量場的條件熵對現(xiàn)有流線進行質量評估，條件熵值越小說明矢量場中未被現(xiàn)有流線挖掘出的信息越少，流線質量越好，需要添加的新流線也越少。加入新的種子點采用基于重要性取樣的方法，保證了矢量場特征區(qū)域的顯著加強顯示。當評估結果收斂時終止迭代過程，此時初始矢量場中大部分的信息已經(jīng)被現(xiàn)有流線顯示出來。

實驗結果表明，基于信息熵的初始布點能夠有效地捕獲流場中的主要特征。當需要流線數(shù)量少的時，僅通過初始布點也能得到很好的效果。初始布點結束時，矢量場主要特征雖然被捕獲，但矢量場中仍有大量信息被隱藏，此時條件熵值也是比較大的。隨著新的流線的不斷加入，流場中隱藏的信息逐漸減少，條件熵值也逐步減小直到收斂于一個較小的值。整個評估過程有效地量化評估了現(xiàn)有流線對矢量場的反應程度，即量化評估了流線的質量。

表1 算法各模塊的運行時間

目前，本文算法只能支持結構化網(wǎng)格，網(wǎng)格種類比較單一，算法性能仍需進一步提升。未來將研究本算法應用于非結構網(wǎng)格的方法。

[1]Jobard B，Lefer W.Creating evenly-based streamlines of arbitrary density[C]//Proceedings of Eighth Eurographics Workshop on Visualization in Scientific Computing，1997：45-55.

[2]Liu Z，Moorhead R，Groner J.An advanced evenly-spaced streamline placement algorithm[J].IEEE Transactions on Visualization and Compute Graphics，2007，13（3）：630-640.

[3]Spencer B，Laramee R S，Chen G，et al.Evenly spaced streamlines forsurfaces：An image-based approach[J].Computer Graphics Forum，2009，28（6）：1618-1631.

[4]Wu K，Kao D，Pang A.Strategy for seeding 3D streamlines[C]//Vis’05：Processing of IEEE Visualization 2005，2005：471-478.

[5]Mebarki A，Alliez P，Devillers O.Farthest point seeding forefficientplacementofstreamlines[C]//Vis’05：Proceedings of the IEEE Visualization 2005，2005：479-486.

[6]Helman J L，Hesselink L.Visualizing vector field topology in fluid flows[J].IEEE Computer Graphics and Applications，1991，11（3）：36-46.

[7]Verma V，Kao D T，Pang A.A flow-guided streamline seeding strategy[C]//Vis’00：ProceedingsoftheIEEE Visualization 2000，2000：163-170.

[8]Ye X，Kao D T，Pang A.Strategy for seeding 3D streamlines[C]//Vis’05：Proceedings of the IEEE Visualization 2005，2005：471-478.

[9]Li L，Shen H W.Image-based streamline generation and rendering[J].IEEE Transactions on Visualization and Computer Graphics，2007，13（3）：630-640.

[10]王少榮，吳迪，汪國平.一種流線放置方法[J].軟件學報，2012，23（S（2））：42-52.

[11]Nvida.CUDA programming guide[M/OL].（2013）.http：//docs.nvidia.com/cuda/pdf/CUDA_C_Programming_Guide.pdf.

[12]Shannon C E.A mathmatical theory of communication[J].Mobile Computing and Communication Reviews，5（1）.

[13]Xu C，Prince J L.Snakes，shapes，and gradient vector flow[J].IEEE Transactions on Image Processing，1998，7（3）：359-369.

[14]Ames W F.Numerical methods for partial differential equations[M].3rd ed.New York：Academic，1992.

[15]Lohse D，Xia K.Small-scale properties of turbulent Rayleigh-Bénard convection[J].AnnualReview of Fluid Mechanics，2010，42：335-364.

[16]李明.非定常三維流場渦結構分析方法研究[D].杭州：浙江大學，2011.