高藝　裘杭萍　羅健欣　吳波

地形渲染是三維地形仿真研究最重要的內(nèi)容之一.數(shù)字高程模型(Digital elevation model,DEM)[1]作為三維地形可視化的一個重要數(shù)據(jù)來源,是實(shí)現(xiàn)地形渲染的基礎(chǔ)支撐.

當(dāng)前DEM數(shù)據(jù)主要有兩種表示方法:不規(guī)則三角網(wǎng)(TIN)和格網(wǎng)(GRID).TIN能用較少的三角形表示地形的高細(xì)節(jié),但由于三角網(wǎng)的不規(guī)則性帶來的存儲困難,使其難以表達(dá)多分辨率地形.

格網(wǎng)DEM是以二維數(shù)組的形式存儲數(shù)據(jù),每個數(shù)組包含一個高程值.因其數(shù)據(jù)結(jié)構(gòu)簡單,易于更新等優(yōu)點(diǎn)成為了當(dāng)前比較常用的一種數(shù)據(jù)格式.

隨著地理數(shù)據(jù)獲取技術(shù)的發(fā)展,一方面,DEM自身的數(shù)據(jù)量日漸龐大,給計(jì)算機(jī)實(shí)時渲染三維地形場景帶來了極大的困難[2];另一方面,在地形渲染時,存儲空間大小隨著分辨率的增加呈指數(shù)級增長.為了減少計(jì)算機(jī)圖形硬件的數(shù)據(jù)處理量,需要尋找一種能夠減少計(jì)算機(jī)實(shí)時處理數(shù)據(jù)量的方法,這對于大規(guī)模地形可視化系統(tǒng)有著重要的意義.

傳統(tǒng)的靜態(tài)渲染算法已經(jīng)不足以處理地形的實(shí)時渲染問題.目前,相關(guān)的研究主要集中在地形的LOD(Level of Detail).Lindstrom提出了連續(xù)LOD(Continues LOD,CLOD)技術(shù),將地形分塊并按精度需求調(diào)用不同的塊來渲染整個場景,減少了實(shí)時渲染的地形網(wǎng)格數(shù)量,提高了實(shí)時渲染的效率[3].

Duchaineau在CLOD算法基礎(chǔ)上提出了一種基于二叉樹模型的ROAM(Real-time Optimally Adapting Meshes)算法,該算法以自頂向下的方式構(gòu)建二叉樹,隨視點(diǎn)的移動,動態(tài)地更新各區(qū)域的LOD級別[4].Zhao等基于ROAM算法結(jié)合四叉樹模型提出了一種修正的LOD技術(shù)[5].

Gobbetti等提出的P-BDAM[6](Planet sized batched dynamic adaptive meshes)和C-BDAM[7](Compressed batched dynamic adaptive meshes)算法處理LOD渲染并采用了小波變換的數(shù)據(jù)壓縮.

Marc Treib提出一種稀疏小波變換算法[8]能達(dá)到較高的壓縮比.小波變換算法不使用傳統(tǒng)的幾何處理的方式去改變地形網(wǎng)格,誤差率低,但是算法涉及復(fù)雜的編解碼運(yùn)算.

本文在基于四叉樹的LOD算法基礎(chǔ)上,提出了一種高效的基于表面近似的渲染方法.這一方法是用近似輪廓取代格網(wǎng)來表示地形,再通過q-Bernstein多項(xiàng)式將離散近似輪廓拼接成G1-連續(xù)曲面.相較于傳統(tǒng)的DEM渲染方法,本文算法不需要復(fù)雜的編解碼過程,直接用參數(shù)映射生成的連續(xù)平滑表面優(yōu)化了DEM數(shù)據(jù),節(jié)省了系統(tǒng)存儲空間;同時,本文算法渲染時不需要加載全部的DEM數(shù)據(jù),極大地提高了渲染效率.

1　近似輪廓

近似輪廓的概念在文獻(xiàn)[9]中被提出,用于存儲體素數(shù)據(jù).這個輪廓為幾何模型定義了一系列的平行平面,如果輪廓能夠很好地近似原始幾何,即用平面來取代原始幾何.這一方法對具有較多平坦表面的場景有很好的壓縮性能.受這一思想啟發(fā),我們將近似輪廓用于地形渲染中,用近似輪廓取代格網(wǎng)來表示地形.

如圖1所示,在一個四叉樹結(jié)構(gòu)中,首先判斷地形變化的總趨勢,用整個地形DEM數(shù)據(jù)計(jì)算最佳匹配平面,然后沿著匹配平面的法線方向計(jì)算最大平面和最小平面,這兩個平面即為近似輪廓的上界和下界.在地形渲染時,當(dāng)近似輪廓可以很好地近似原始幾何時,如圖中黑色線框區(qū)域,無論分辨率如何增加都只用近似輪廓來代替.

圖1　使用近似輪廓表達(dá)DEM數(shù)據(jù)

利用特征值最小二乘平面法來構(gòu)造DEM的最佳匹配平面.DEM數(shù)據(jù)中的2D點(diǎn)像素坐標(biāo)為用齊次坐標(biāo)表示為P2=R3?(0,0,0)是二維映射空間.一個3D平面用齊次坐標(biāo)表示為:

圖2　平面示意圖

給定DEM中一組n個點(diǎn)P={P1,P2,P3,···,Pn},(xi,yi,zi)≡Pi.計(jì)算矩心則定義一組點(diǎn)Q為相對于矩心center的DEM點(diǎn),要獲得最佳匹配平面,在條件a2+b2+c2=1約束下,滿足最小.利用拉格朗日乘數(shù)法求函數(shù)極值構(gòu)成特征值方程:

由文獻(xiàn)[10]的推導(dǎo)可得:

令M=即為e的最小特征值,最小特征值對應(yīng)的特征向量即為平面方程的參數(shù)a,b,c.

給出點(diǎn)到匹配平面的最大最小距離:

則,原點(diǎn)到最大最小平面的距離為:

求得近似輪廓平面為:

實(shí)際情況中,DEM數(shù)據(jù)在形成和傳輸過程中會引入多種噪聲,噪聲往往表現(xiàn)為或大或小的極值.如果不考慮噪聲的影響,將會導(dǎo)致獲得的匹配平面產(chǎn)生偏差,算法將不具有魯棒性.因此,在生成匹配平面前需進(jìn)行DEM的去噪操作.

2　近似輪廓的曲面連續(xù)

上述離散的近似輪廓在渲染時會出現(xiàn)不同程度的地形裂縫問題.為解決這一問題,本文引入q-Bernstein多項(xiàng)式[11]將離散近似輪廓拼接成G1-連續(xù)曲面.

為每一個匹配平面定義一個離散函數(shù)f(向量形式F):

其中,V為平面頂點(diǎn)集,定義一個連續(xù)函數(shù)f′,如果f′(V)=f(V)則f′為f的插值.給出一個權(quán)值因子是一個q-Bernstein多項(xiàng)式,與標(biāo)準(zhǔn)的Bernstein多項(xiàng)式相比,表達(dá)式中含有一個可調(diào)控形狀的參數(shù)值q.其定義如下:

其中,參數(shù)q為正實(shí)數(shù),

當(dāng)n≥i≥1且r=0時,值為1,否則值為0.當(dāng)q=1時,即簡化為一般的二項(xiàng)式系數(shù).

根據(jù)文獻(xiàn)[12-13],定義含有兩個張量積形式的n×m次q-B′ezier曲面為:

其中,Vij為控制頂點(diǎn),和是參數(shù)為q1和q2的q-Bernstein多項(xiàng)式.

定理.P1(u,v)和P2(u,v)為n×m次張量積q-B′ezier曲面,{pi,j}和{ri,j}分別為其q-Bernstein多項(xiàng)式控制點(diǎn)矩陣,P1(u,v)和P2(u,v)拼接達(dá)到G1-連續(xù)的充要條件為滿足式(3)和式(4):

且拼接邊界上u的偏導(dǎo)數(shù)對所有的v都有

即

α為大于0的常數(shù).

證明.將式(2)轉(zhuǎn)換為標(biāo)準(zhǔn)B′ezier曲面P′(u,v),為其標(biāo)準(zhǔn)Bernstein多項(xiàng)式控制點(diǎn)矩陣.式(2)可寫成

由文獻(xiàn)[13]可知,q-B′ezier曲面和標(biāo)準(zhǔn)B′ezier曲面存在如下的矩陣變換關(guān)系:

矩陣M有元素:

mi,j=s(n,j)為生成函數(shù)的系數(shù),滿足:(1+t)(1+[2]t)···(1+[n]t)=詳細(xì)的公式推導(dǎo)過程請參見文獻(xiàn)[13].

由式(5)推出式(4)有如下的轉(zhuǎn)化公式:

由B′ezier曲面的仿射不變性可知,式(2)和式(3)成立.

因此,本文通過q-B′ezier曲面可以實(shí)現(xiàn)將離散近似輪廓拼接成G1-連續(xù)曲面,從而消除地形裂縫.

3　算法結(jié)果

本文用C++與Cg語言實(shí)現(xiàn)了本文算法,并在微軟Window XP系統(tǒng)上運(yùn)行.機(jī)器CPU為Intel Core TMi7 950(3.07GHz),內(nèi)存3GB.顯卡為NVGF GTX 480,顯存1536MB.實(shí)驗(yàn)數(shù)據(jù)使用一為Puget Sound數(shù)據(jù)集.二為整個亞洲地區(qū)的Google 19層紋理數(shù)據(jù)集及S高程數(shù)據(jù),數(shù)據(jù)塊的調(diào)度使用文獻(xiàn)[14]的屏幕誤差公式:

式中,h為視口的像素高度,θ為垂直視野角,τ為屏幕誤差容忍,Pz為相機(jī)深度,ε為世界空間誤差容忍.在本文實(shí)驗(yàn)中,屏幕誤差設(shè)置為2/3像素.

3.1　存儲及渲染效率對比分析

將Puget Sound數(shù)據(jù)集比較不同的分辨率下與文獻(xiàn)[15]提出的地形壓縮存儲算法進(jìn)行存儲量的對比(每個高程值用一個2字節(jié)的無符號整型數(shù)來存放.從表1可以看出,本文算法在低分辨率時比文獻(xiàn)[15]的存儲量大,因?yàn)樗牟鏄涞钠史謺a(chǎn)生更多的節(jié)點(diǎn),但隨著分辨率的不斷增加,本文算法只存儲近似輪廓即可,四叉樹將不再進(jìn)行遞歸剖分,分辨率越大,存儲優(yōu)勢越明顯.

表1　存儲量對比結(jié)果

圖3顯示了本文算法的渲染結(jié)果,可以看出算法可以獲得較高質(zhì)量的渲染平面地形.

圖3　場景截圖

將本文算法與經(jīng)典的Chunked LoD[16]進(jìn)行渲染效率的對比.平面分辨率為4096×4096的Puget Sound數(shù)據(jù)集,窗口大小為1024×768.幀速率公式[17]為:

式中,dt為幀時間,w?h為屏幕分辨率.如圖4所示,圖中橫坐標(biāo)為渲染過程中的時間記錄,縱坐標(biāo)為記錄的幀速率.Chunked LoD方法在渲染過程中為了避免裂縫的產(chǎn)生,會強(qiáng)制添加過渡三角形來填補(bǔ)塊間裂縫,造成網(wǎng)格數(shù)量的增加,從而影響渲染效率.本文算法的平均幀速率大約為102.334fps,而Chunked LoD的平均幀速率大約為86.349fps,由此可見本文算法具有更高的渲染效率.

3.2　裂縫消除效果

圖5顯示了用本文方法不使用q-B′ezier曲面渲染和使用q-B′ezier曲面渲染的不同效果.不使用q-B′ezier曲面時,能從圖5(a)中看到明顯的裂縫,而使用q-B′ezier曲面后,這樣的問題得到了明顯的改善(圖5(b)).

圖4　光線投射渲染效率對比

圖5　圖像比較

4　結(jié)論

本文提出了利用近似輪廓取代格網(wǎng)來表示地形的算法,并將離散近似輪廓通過q-Bernstein多項(xiàng)式生成了G1-連續(xù)的q-B′ezier曲面,避免了離散近似輪廓在渲染時帶來的裂縫問題.實(shí)驗(yàn)結(jié)果驗(yàn)證了本文算法.下一步研究將該算法移植到GPU中,使其分擔(dān)CPU中數(shù)值計(jì)算的任務(wù),從而提高渲染效率.

1 TAO Y,YANG X,TANG G,et al.Primary studies of local shape correction of coarse DEM based on 3D terrain features[C]//International Conference on Geoinformatics.IEEE,2010:1?5.

2 ZHANG R H.Real-time optimization technology and its application in terrain rendering[C]//International Congress on Image and Signal Processing.IEEE,2011:1349?1352.

3 LINDSTROM P,KOLLER D,RIBARSKY W,et al.Real-time,continuous level of detail rendering of height fields[C//Proceedings of the 23rd annual conference on Computer graphics and interactive techniques.ACM,1996:109?118.

4 DUCHAINEAU M,WOLINSKY M,SIGETI D E,et al.ROAMing terrain:real-time optimally adapting meshes[C]//Proceedings of the 8th Conference on Visualization’97.IEEE Computer Society Press,1997:81?88.

5 ZHAO Y,MA Y.A modifieLOD Terrain Model Based on QuadTree Algorithm[C]//InternationalJointConferenceonComputationalSciences and Optimization.IEEE,2009:259?263.

6 CIGNONI P,GANOVELLI F,GOBBETTI E,et al.Planet-sized batched dynamic adaptive meshes(P?BDAM)[C]//IEEE Visualization Conference,2003:20.

7 GOBBETTI E,MARTON F,CIGNONI P,et al.C-BDAM–compressed batched dynamic adaptive meshes for terrain rendering[J].Computer Graphics Forum,2006,25(3):333?342.

8 TREIB M,REICHL F,AUER S,et al.Interactive Editing of GigaSample Terrain Fields[C]//Computer Graphics Forum,2012:383?392.

9 LAINE S,KARRAS T.Efficient sparse voxel octrees[J].IEEE Transactions on Visualization and Computer Graphics,2011,17(8):1048?1059.

10官云蘭,程效軍,施貴剛.一種穩(wěn)健的點(diǎn)云數(shù)據(jù)平面擬合方法[J].同濟(jì)大學(xué)學(xué)報(自然科學(xué)版),2008,36(7):981?984.

11 ORUC?H,PHILLIPS G M.q-Bernstein polynomials and B′ezier curves[J].Journal of Computational and Applied Mathematics,2003,151(1):1?12.

12 ZHOU B,SONG L,LI J.Quasi rational B′ezier curves with shape parameter[J].Journal of Residuals Science&Technology,2016,13(7):143.1?143.6.

13 HAN L W,CHU Y,QIU Z Y.Generalized B′ezier curves and surfaces based on Lupa q-analogue of Bernstein operator[J].Journal of Computational and Applied Mathematics,2014,261:352?363.

14 COHEN J D.Appearance-preserving simplificatioof polygonal models[M].Chapel Hill:The University of North Carolina at Chapel Hill,1998.

15白建軍,閆超德.基于二叉樹的多分辨率地形數(shù)據(jù)壓縮存儲[J].西北大學(xué)學(xué)報(自然科學(xué)版),2010,40(6):1088?1092.

16 ULRICH T.Rendering massive terrains using chunked level of detail control[C]//Course Notes of Computer Graphics Annual Conference Series,SanAntonio,ACM,2002.

17 WINGERDEN T V.Real-time ray tracing and editing of large voxel scenes[D].Utrecht:Utrecht University,2015.