張繼凱，張然，趙君，顧蘭君

(內(nèi)蒙古科技大學 信息工程學院，內(nèi)蒙古 包頭 014010)

近年來，隨著高分衛(wèi)星的應用和遙感技術的飛速發(fā)展，地理信息相關的遙感數(shù)據(jù)量呈爆發(fā)式增長.因遙感數(shù)據(jù)的地理跨度較大，地表弧度帶來的影響無法忽略不計，導致傳統(tǒng)體可視化中的立方體法無法適用；另外，傳統(tǒng)的全球范圍數(shù)據(jù)可視化方式多采用經(jīng)緯度網(wǎng)格進行數(shù)據(jù)存儲，這就造成數(shù)據(jù)分布在全球各處密度不同，不僅會造成各緯度的可視化細節(jié)表現(xiàn)差異，更是會帶來嚴重的極點問題和數(shù)據(jù)存儲冗余現(xiàn)象.國內(nèi)外研究人員對此提出過多種不同的球面網(wǎng)格剖分方式，其中六邊形離散網(wǎng)格因其良好的數(shù)據(jù)分布特性，成為地理信息領域和數(shù)據(jù)可視化領域近年來的一個研究熱點.

本文在對球面進行體可視化的操作過程中將六邊形離散網(wǎng)格引入其中，使得球面各處得數(shù)據(jù)密度可以保持一致.同時，為了簡化光線投射算法的插值過程，設計一種基于三角形單元的編碼方式，以紋理作為其存儲載體，使光線投射過程的數(shù)據(jù)插值點更容易查找，更適合硬件加速實現(xiàn).

1 相關工作

隨著硬件計算能力地提升和圖形顯示手段地發(fā)展，對海量數(shù)據(jù)進行可視化成為數(shù)據(jù)分析和知識發(fā)現(xiàn)的有效手段之一，而其中三維遙感數(shù)據(jù)的體可視化是信息含量最大難度最高的可視化形式.

直接體可視化技術經(jīng)過多年的發(fā)展，已經(jīng)產(chǎn)生諸如拋雪球算法[1，2]、頻域算法[3]、錯切-變形算法[4]以及光線投射算法[5]等多類算法.其中，因算法原理相對簡單，可視化效果相對更符合實際數(shù)據(jù)特征等優(yōu)勢，光線投射算法成為當前體可視化中應用最為廣泛的一類算法.

在球面離散網(wǎng)格的研究上，趙學勝等[6，7]于2009年提出一種基于菱形塊四叉樹節(jié)點的鄰近搜索算法，使用三角形單元作為頂層節(jié)點，而其余各層使用四邊形單元且層次越低越接近于矩形，但因各緯度帶的剖分差異性造成復雜的拓撲結構.針對六邊形離散網(wǎng)格，賁進、童曉沖和張永生等人進行了一系列探索[8-11]，他們通過分析各種六邊形剖分方法，實現(xiàn)一種新的全球六邊形離散網(wǎng)格剖分及編碼算法，并證明了其在計算幾何領域有較好地應用潛力[12,13].但因該方法的編碼規(guī)則比較復雜且涉及大量的歸屬性判斷，使得很難將其應用于數(shù)據(jù)可視化尤其是三維數(shù)據(jù)場可視化中.

XIE J等[15]于2013年提出了基于六邊形網(wǎng)格的光線投射算法，才使得球面離散網(wǎng)絡應用于體可視化領域的探索有所突破.該算法通過采用六邊形網(wǎng)格實現(xiàn)了直接體繪制，但是由于需提前按照六邊形網(wǎng)格相互之間的對應關系將其頂點、邊以及網(wǎng)格單元全都連接起來，然后實時的依據(jù)該相互之間的對應關系檢索相鄰的網(wǎng)格單元，從而實現(xiàn)光線投射，但該算法復雜程度高.為了提高渲染速度，他們之后還研究使用多GPU并行處理的方式對全球體數(shù)據(jù)進行體可視化[16].

2 球面六邊形網(wǎng)格數(shù)據(jù)預處理

2.1 數(shù)據(jù)組織方式

通過合理的方式構造和儲存六邊形網(wǎng)格信息，從而應用于光線投射算法.圖1所示為六邊形網(wǎng)格的構造組成，六邊形網(wǎng)格單元的每個頂點使用灰色數(shù)字進行編碼，六邊形的中心則用黑色圓點表示，并用黑色數(shù)字對其中心進行編碼，由圖可知每個六邊形網(wǎng)格的中心都與其6個頂點相對應，而每個六邊形的頂點又與其周圍相鄰的3個六邊形中心相對應，網(wǎng)格數(shù)據(jù)根據(jù)該結構進行構造.

為了滿足不同分辨率數(shù)據(jù)下的可視化需求，依據(jù)文獻[15]提出的球面離散網(wǎng)格生成算法對不同分層的六邊形網(wǎng)格數(shù)據(jù)進行構造生成.對于每個六邊形網(wǎng)格單元，使用XML的文件形式構造保存網(wǎng)格數(shù)據(jù)，主要包含單元網(wǎng)格編碼、六邊形網(wǎng)格各個頂點及其中心點的經(jīng)度坐標數(shù)據(jù)與緯度坐標數(shù)據(jù).表1展示了不同分層的網(wǎng)格構造過程，其中在當前分辨率下，每個六邊形網(wǎng)格單元的面積用表示，地球的表面積使用表示，網(wǎng)格分層層數(shù)用N表示.

表1 網(wǎng)格剖分說明表

2.2 數(shù)據(jù)項處理

在對體數(shù)據(jù)進行可視化時，若通過具有相同分辨率的第一層網(wǎng)格中心為原點按不同步長向徑向進行縮放，用同樣的層次數(shù)且分辨率也一致的六邊形網(wǎng)格作為數(shù)據(jù)載體，以充分體現(xiàn)體數(shù)據(jù)的多層次特性，即可得到具有不同高度但分辨率相同的多層次網(wǎng)格.

六邊形單元都是由6個頂點和1個中心點組成，將體數(shù)據(jù)與六邊形網(wǎng)格互相關聯(lián)時，可以在頂點和中心點間進行選擇.若選擇前者，則圖2所示的每個黑色網(wǎng)格代表一個“數(shù)據(jù)單元”，任一層六邊形單元在位置上都與其他層單元在徑向上呈對應關系，且相鄰2層的“數(shù)據(jù)單元”可以形成1個六棱柱結構，如圖3(a)所示；若選擇中心點綁定，每3個相鄰六邊形的中心點相連得到1個三角形“數(shù)據(jù)單元”，圖2的紅色4虛線所示，相鄰兩層的三角形“數(shù)據(jù)單元”之間可構成三棱柱，如圖3(b)示.

基于六邊形網(wǎng)格的光線投射插值算法會涉及2步耗時量的工作：一是確定采樣點所在的棱柱，二是計算采樣點的數(shù)據(jù)值.前者需要獲取采樣點的經(jīng)緯度以及高度信息，后者需要計算采樣點在棱柱內(nèi)部所處的位置，并根據(jù)棱柱頂點數(shù)據(jù)進行插值計算.鑒于中心綁定的計算量比頂點綁定要少，因此選擇六邊形中心綁定數(shù)據(jù).在減少計算時間的同時，為進一步減少存儲空間，在縱向上利用棱柱體形成的各層數(shù)據(jù)的一一映射關系，在橫向上利用每個六邊形中心與一系列處于同一經(jīng)緯度的數(shù)據(jù)形成的關聯(lián)關系，可得到其他層的數(shù)據(jù)，所以只需存儲最底層網(wǎng)格數(shù)據(jù).

關于空間某點所屬三棱柱的確定，首先根據(jù)三棱柱的投影三角形與坐標點的經(jīng)度、維度進行關聯(lián)，得到該點底部三角形的位置，而后通過點的高度信息向上正方向映射從而確定該點所在網(wǎng)格單元的層級，得到該點所屬三棱柱的信息.假設在各種條件都滿足的條件下，三角形的網(wǎng)格單元與經(jīng)度信息和緯度信息理應具有線性轉換關系，也就是說即使經(jīng)度值與緯度值具有非常高的精度，也都可以通過相關的轉換或檢索計算，最終獲得相應的三角形網(wǎng)格單元.利用體數(shù)據(jù)分布的離散性，構造一張精度極高的“檢索表”，該表可以對經(jīng)緯度到三角形單元的映射進行詳細記錄，進而有效地減少巨大的計算量，如對空間中某點所屬的三角形單元進行檢索時，通過“檢索表”快速縮小查找范圍，再經(jīng)少量計算即可完成轉換.該“檢索表”的原理如圖4(a)所示，其具有個大小均勻的單元，為最大經(jīng)度，最小經(jīng)度，為最大緯度，為最小緯度，每個單元都一個有獨一無二的經(jīng)緯度范圍，并且固定經(jīng)緯度范圍的單元中存放了所屬的三角形單元編號.在查找某一經(jīng)緯度所屬三角形單元時，只需查找其所在單元格即可.

如圖4(a)所示，圖中所屬三角形的邊緣處用藍色網(wǎng)格單元表示，當劃分該單元格所屬的三角形時，由于其或許橫跨多個三角形，所以存儲于該藍色網(wǎng)格單元中的編碼信息也許存在稍許誤差.對此，如圖4(b)所示，將單元格的密度增大，減小誤差，從而也使得處于三角形邊緣的單元格區(qū)域大大減小.然而，如果單純的增大劃分單元格的密度，則會占用“檢索表”大量的存儲空間，由此需要找到一個合適劃分的密度，既不會讓誤差增大，也無需占用大量存儲空間.通過多次對比實驗表明，當將“檢索表”的單元格精度設置為(單位：度)時，最終可視化呈現(xiàn)的效果最佳，的計算公式如下：

(1)

式中：dc為每2個相鄰的六邊形中心的距離；β為緯度值.式(1)顯示緯度與劃分的網(wǎng)格精度成反比.

“檢索表”初始化之前必須對其精度進行確認，而每張“檢索表”的單元格的標準均設置為其中心點的經(jīng)緯度坐標位置，因為已知每個三角形單元的所有頂點的坐標，由此可運用重心法計算得到各個三角形和每個“檢索表”單元的從屬關系，每個單元存儲屬于其三角形的編碼.為了在可視化進行中可以快速讀取，通過二進制或文本的結構形式將該表按照“行優(yōu)先”的存儲方式置于外存中，該表的生成均在預處理階段完成.

重心法的基本原理是:已知三角形的所有頂點均處于同一個平面，若將該三角形的任一頂點設為原點坐標，那么由該點指向其他2點的方向即為2個坐標軸的方向.如圖5所示，將點A定為原點，則AB方向為坐標軸V的正坐標方向A方向為坐標軸U的正坐標方向.

因此，對于同一平面內(nèi)的任意一點，均可通過式(2)表示得出：

p=A+(C-A)·x+(B-A)·y.

(2)

可見，如果x或y為負數(shù)，則點是向U或V軸的負方向移動，就處于三角形外部，因此為了保證點處于三角形內(nèi)部，需要同時滿足如下條件：x≥0，y≥0且x+y≤1.其中，當x=0且y=0時，p點與A點重合；同理，當x=0且y=1時，p點B與點重合；當x=1且y=0時，p點與C點重合.

為了求取已知直角坐標的點在三角形中的UV坐標，對式(2)做如下變換：

p-A=(C-A)·x+(B-A)·y.

(3)

令向量v0=C-A，v1=B-A，v2=p-A，則v2=v0·x+v1·y，有

v2·v0=(v0·x+v1·y)·v0

v2·v1=(v0·x+v1·y)·v1，

(4)

對其求解可得：

(5)

所以，將三角形網(wǎng)格的單元中心點與其頂點的位置坐標代入式(5)，若和值是滿足條件，那么可以得出三角形與單元格之間滿足相對應的從屬關系.

3 基于六邊形網(wǎng)格的體繪制

本文使用光線投射算法對體數(shù)據(jù)進行可視化的操作過程中，各光線的積分、插值等運算沒有任何融合等現(xiàn)象存在.由于GPU有著強大的并行處理能力，利用該方法進行計算，并且每條光線獨立積分運算.算法的大致流程圖如圖6所示.

球面光線投射算法的具體步驟如下：

3.1 球面光線插值

傳統(tǒng)的光線投射算法將數(shù)據(jù)分別存放于六面體的6個面中，且所有數(shù)據(jù)均通過長、寬、高3個維度相互對應.當有屏幕像素隨著視線穿過六面體時，六面體的表面就留下了2個交點，分別為輸入點與輸出點，如圖7(a)所示，射線與六面體分別相交于和.光線投射算法的基本原理為，均勻的在和之間取若干個數(shù)據(jù)，通過顏色疊加和變化等運算獲得屏幕上光線的顏色值.

在全球數(shù)據(jù)體可視化中，體數(shù)據(jù)范圍一般分布在一定高度的球殼中，坐標系也從直角坐標系轉換成球面坐標系，所以數(shù)據(jù)維度也就不存在上述所說的一一對應的關系，數(shù)據(jù)被壓縮到地球表面的空心球殼中，光線在該球殼體中的投射過程如圖7(b)所示.本文采用六邊形網(wǎng)格作為數(shù)據(jù)載體，由1.2節(jié)獲取光線中的采樣點數(shù)據(jù)可知，求得各采樣點的經(jīng)緯高坐標才是關鍵，由于光線只在直角坐標系椽筆，所以均勻采樣也要在該坐標系進行，因此球面光線投射中采樣點經(jīng)緯高坐標的計算步驟如下：

Step1 利用幾何著色器和渲染到紋理技術將代理幾何體前后面的坐標信息保存成紋理形式，三維位置信息中的XYZ分別保存至紋理的RGB通道中.

Step2 在Step1中的2個紋理中，按對應坐標查找光線的2個交點，將其坐標分別設為和.

Step3 按照設定步長從前到后選擇采樣點，假設采樣步長為，對于第個采樣點，其空間坐標計算式(6)所示，其中l(wèi)ength()為求取向量長度的函數(shù).

(6)

Step4 利用式(7)將采樣點坐標從直角坐標(X，Y，Z)轉化為地球球面經(jīng)緯高坐標(lon，lat，alt).

lon=arctan(Y/X)

lat=arctan(U/V)

alt=p/cos(lat)-N.

(7)

式中：P，N，U，V的計算公式如下：

(8)

式中：re，rp分別為地球的赤道半徑和2極半徑，ds，theta的計算公式如下：

ds=(re2-rp2)/rp2

theta=arctan(Z·re/(p·rp)) .

(9)

2.2 采樣點數(shù)據(jù)計算

通過在數(shù)據(jù)預處理階段生成的三角網(wǎng)格編碼的“檢索表”，可以從采樣點的經(jīng)緯度直接得到相應的三角網(wǎng)格編碼TCn，然后，根據(jù)1.1中的XML文件，查詢到相應的3個六邊形編碼HCn1，HCn2和HCn3，并找到它們對應的數(shù)據(jù)鏈表，為了確定對應的三棱柱6個頂點的相關數(shù)據(jù),應根據(jù)式(10)計算采樣點的網(wǎng)格層次，即采樣點相對于地面的高度h，體數(shù)據(jù)的空間高度hv，體數(shù)據(jù)的層數(shù)l，以及體數(shù)據(jù)所在棱鏡的上表面索引index.

(10)

因此，可以從3個數(shù)據(jù)鏈接表中獲得index和index-1層的6個三角形頂點a，b，c，a′，b′和c′，及其對應的6個數(shù)據(jù)值，并且它們與采樣點的關系如圖8所示.其中，O點是地球的中心點，p為采樣點，O和P的連線在上下表面的交點分別是p′和p″.首先，采用重心插值法計算P′和P″的數(shù)據(jù)值，然后采用線性插值則可獲得P的數(shù)據(jù)值.

數(shù)據(jù)值和顏色值之間的轉換需要使用傳遞函數(shù).在文中，線性映射函數(shù)被用作體積可視化的傳遞函數(shù)，即圖9所示的顏色映射表的藍色末端.對應于數(shù)據(jù)的最小值，紅色端對應于數(shù)據(jù)的最大值,并根據(jù)數(shù)據(jù)大小通過線性插值得到其他數(shù)據(jù)的相應位置.相應的屏幕像素的最終顏色值可以通過累積光線上的采樣點的顏色來獲得.

2.3 繪制模式

由于體繪制過程中需要進行大量復雜的計算，如果對圖像的每一幀都進行計算，將占用很大的計算資源.根據(jù)光線投射原理，當視覺參數(shù)不變，即視點、場景、窗口等參數(shù)一致時，光線投射產(chǎn)生完全相同的結果，即無需重復計算；除此之外，視點的視覺辨別力在移動的過程中會下降，從而對場景渲染精度的要求也會隨之減小，所以，在視點移動的過程中，可以通過適當?shù)慕档凸饩€的計算精度，從而使繪制速度有效提升.

為了有效利用繪制資源，避免額外的開銷，提升預測的效率，設置了3種繪制方式：方式Ⅰ，零繪制，即將上一幀的繪制結果直接用于下一幀的繪制；方式Ⅱ，粗繪制，即光線投射計算時設置較大的步長進行采樣；方式Ⅲ，精細繪制，即實施完整的光線投射計算，同時在計算過程中按規(guī)定步長采樣.這3種方式的選擇和調(diào)度過程如圖10所示，其中閾值f是指當視覺參數(shù)保持f幀靜止時設置精細計算,設置閾值的目的是為了預防交互過程中出現(xiàn)假靜止，避免不必要的繪制.

4 實驗結果分析

為了證明該算法的有效性，在windows7操作系統(tǒng)下，分別使用OSG三維圖形引擎、GPU編程語言CUDA以及GLSL著色語言對算法進行了實驗；實驗所使用的硬件環(huán)境是Intel i7-4770k 3.5 GHz CPU，NVIDIA GTX780TI GPU，以及16G RAM；測試數(shù)據(jù)分為單層數(shù)據(jù)和多層體數(shù)據(jù)，其中前者為全球海洋某一時刻的溫度場數(shù)據(jù)，后者數(shù)據(jù)是來自太平洋西部海域某時刻的海洋鹽度數(shù)據(jù)，經(jīng)度在(99.0°E-150.0°E)和(0.5°N-52°N)范圍之內(nèi)，分辨率為 409×497×51.接下來將從以下方面對實驗的結果進行論述.

4.1 六邊形網(wǎng)格的全球分布及其可視化

圖11為單層六邊形網(wǎng)格在全球的分布效果示意圖，其中圖11(a)，(b)，(c)分別為地球中部低緯度地區(qū)和兩極地區(qū)的可視化結果，可見在地球球面的任意位置，六邊形網(wǎng)格的分布都呈現(xiàn)比較均勻的特性.將單層六邊形網(wǎng)格用于全球海洋溫度場可視化的效果如圖12所示.

4.2 渲染效果對比

圖13是基于經(jīng)緯度網(wǎng)格的體可視化效果，圖14是基于六邊形網(wǎng)格的體可視化效果，使用了不同的顏色代表的了該位置此刻的鹽度強度，越接近藍色說明鹽度越低，相反越接近紅色說明鹽度越高.由圖13(b)和圖14(b)進行對比可以看出，基于經(jīng)緯度網(wǎng)格的體可視化過程中會出現(xiàn)顆粒感，特別是在一些起伏較大的區(qū)域，在使用文中的算法之可以使可視化后的效果更加平滑，整體效果表現(xiàn)良好.

4.3 渲染幀率對比

文中在場景渲染效率方面基于六邊形網(wǎng)格的體可視化方法與傳統(tǒng)方法進行了比較.如圖15所示，步長越短，采樣的點就越多，積分的速度就越慢，導致光線投射的時間延長，但是最終繪制的可視化效果也會因此變佳，需要在最終呈現(xiàn)的可視化效果和性能權衡中取1個平衡值，將其用于實際應用中.

從圖15中的4條曲線可以看出，基于經(jīng)緯度的方法比基于六邊形網(wǎng)格渲染幀率高，前者在渲染幀率上增加了19.7%，造成該現(xiàn)象的原因首先是傳統(tǒng)經(jīng)緯網(wǎng)格只需1次三線性插值即可獲取數(shù)據(jù)值，而基于六邊形的方法在平頭截錐體中進行多次插值才能獲得；其次另一個原因是六邊形網(wǎng)格在光線投射過程中要進行幾次的空間坐標與經(jīng)緯度坐標的轉換，這同時就會增加運算的時間.但在實際的應用中兩種方法在繪制幀率上的差別幾乎難以發(fā)覺.

4.4 數(shù)據(jù)量對比

六邊形網(wǎng)格方法不僅能夠有效提升繪制的效果，而且還能夠減少體數(shù)據(jù)量，以本文為例，經(jīng)緯度網(wǎng)格的分辨率是409×497×51，總共有10 366 923個標量值.如果按照六邊形網(wǎng)格進行存儲，僅需9 246 934個標量值就可以達到低緯數(shù)據(jù)密度，數(shù)據(jù)量減少了10.8%，如在全球范圍內(nèi)使用六邊形網(wǎng)格存儲將比經(jīng)緯度存儲減少33.7%的數(shù)據(jù)存儲量.由于硬件的發(fā)展不均衡，目前顯存價格是遠遠要比同等容量的內(nèi)存要昂貴，所以縮小體數(shù)據(jù)的容量以減少顯存占用明顯具有重要的意義.

4.5 算法復雜度

文獻[14]中提出的基于六邊形網(wǎng)格的光線投射算法，使用了4張二維表將各個六邊形單元的中心點、頂點還有邊的對應關系進行了存儲，用于緊鄰搜索，因為要遍歷網(wǎng)格單元，所以算法的復雜度是O(n).而本文算法實在數(shù)據(jù)預處理階段就生成了一張經(jīng)緯度-三角形網(wǎng)格的檢索表，所以將算法的復雜度降為了O(1)，有效的減少了計算量.

5 結論

將球面六邊形網(wǎng)格引入球面體繪制算法中，有效地降低了數(shù)據(jù)存儲空間的占用率，在渲染效果方面有顯著提升；通過建立經(jīng)緯度-三角形檢索表以及設計三棱柱采樣點插值算法，降低了算法的復雜度，同時經(jīng)過優(yōu)化渲染模式，在對結果影響很小的情況下提升了渲染的速率.

由于全球空間有著大范圍、大數(shù)據(jù)量的特性，將這些完整的數(shù)據(jù)在普通計算機上進行可視化，傳統(tǒng)的算法是難以實現(xiàn)的，考慮到三維空間中人眼對于較近位置的場景觀察力要比較遠處高，所以下一步將考慮把光線投射中的自適應采樣技術引入到六邊形網(wǎng)格體繪制中，進一步提高數(shù)據(jù)規(guī)模較大時的計算速度.