翟 銳，呂 科，代雙鳳，潘衛(wèi)國

(中國科學院大學工程管理與信息技術學院，北京，100049)

基于地形高度域的數據壓縮算法研究

翟 銳，呂 科，代雙鳳，潘衛(wèi)國

(中國科學院大學工程管理與信息技術學院，北京，100049)

隨著遙感技術的發(fā)展，地形數據規(guī)模越來越大，遠遠超過了內存處理的范圍，成為急需解決的問題.通過數據壓縮提高系統(tǒng)吞吐量是常用技術之一，隨著GPU技術的快速發(fā)展，傳統(tǒng)的壓縮算法無法充分利用GPU的能力.鑒于此，本文提出了一種基于GPU的地形數據壓縮方法，實現(xiàn)了高度域和位置信息的壓縮.不同于其他的算法僅對高度或位置進行壓縮，本文的主要貢獻在于將地形的位置和高度同時進行處理，當前頂點的所有信息都可以根據當前分段計算得到.算法對地形的高度域進行貝塞爾曲線的近似，保存每個頂點的差值，實現(xiàn)有損和無損的相結合的高比率的壓縮.通過與傳統(tǒng)方法的比較，實驗結果表明，能夠取得很好的壓縮效果.

數據壓縮；地形渲染；圖形處理器

1 引言

大規(guī)模地形渲染是計算機圖形學的主要研究內容之一，廣泛應用于虛擬現(xiàn)實、地理信息系統(tǒng)、飛行模擬和游戲等領域.隨著遙感技術的發(fā)展，數字地形數據的分辨率日益增高，數據規(guī)模越來越大.近幾年，GPU計算能力得到了飛速提升，處理速度比從內存?zhèn)鬏斨翀D形顯卡的速度更快，目前的渲染算法中已從早期處理單個三角形變?yōu)樘幚砣切螇K，因此數據傳輸成為地形渲染的瓶頸之一，而通過對地形數據采用壓縮技術增大系統(tǒng)的吞吐量可以有效的解決這一問題，所以壓縮技術的研究已成為國內外研究的熱點.

在地形渲染中，數據壓縮可用于多種數據：包括紋理、高度域或位置等.根據壓縮方法的不同，可以分為有損壓縮、無損壓縮和兩種相結合的方法.針對于不同的應用場景，需求也有一定的差異.例如，對于無交互的應用，解碼速度要求相對較低，可以使用無損的壓縮方法.而對渲染速度要求較高的場合，可能需要采用有損的壓縮方法.結合目前的GPU技術，本文提出了一種無損和有損相結合的地形位置信息和高度域信息的壓縮算法，可以實現(xiàn)高比率的數據壓縮，提高系統(tǒng)的吞吐量.

2 國內外研究現(xiàn)狀

根據數據結構的不同，地形渲染算法可以分為兩大類：基于規(guī)則網格的算法如實時優(yōu)化自適應網格算法翟 銳：基于地形高度域的數據壓縮算法研究(Real-Time Optimally Adapting Mesh，ROAM)[1]和基于不規(guī)則三角網的算法(Triangulated Irregular Network，TIN)，如漸進網格(Progressive Mesh)[2]等.與不規(guī)則網格相比，規(guī)則結構更適合基于GPU的并行的環(huán)境，統(tǒng)一的結構更易于編碼與實現(xiàn)，在壓縮領域也是如此.本文的研究內容是基于約束四叉樹的位置信息和高度域信息的壓縮.

根據壓縮策略的不同，可以分為有損壓縮、無損壓縮或兩者相結合的方法.有損壓縮常被用于實時渲染領域，對渲染速度要求較高的系統(tǒng).例如Gerstner[3]處理原始地形數據的子集來壓縮高度域，通過線性的插值計算被刪除的頂點.Kim[4]等將地形數據進行小波變換來構建近似的三角化網格.但是這些方法編碼速度較慢.除了實時渲染系統(tǒng)，有損壓縮還用于分布的網絡傳輸中[5].

無損壓縮方法中，根據相鄰頂點的信息預測頂點的高度，使用通用的無損壓縮方法對預測誤差進行編碼.早期的一些算法中，使用各種方法對編碼進行預測，如霍夫曼編碼、拉格朗日乘數優(yōu)化的線性方法或算術編碼[6,7]等.這些方法的壓縮效率比直接使用數字圖像編碼的壓縮效率要高，但處理數據的方式是順序的，不適用于并行的壓縮.還有一些算法將高度域映射到更高階的表面[8]，實現(xiàn)基于點的并行，不過這些方法提出的時間較早，在可編程的GPU之前出現(xiàn)，并沒有考慮到目前的并行的框架.也有一些算法將圖像的壓縮算法應用于地形的壓縮方法中.與圖像的灰度壓縮不同，在一些地形的壓縮算法中，地形的壓縮算法需要考慮到地形的結構，而在灰度圖的壓縮中通常沒有考慮到這一點.

最近的一些壓縮方法在GPU上處理數據.Dick等[9]使用與文獻[3]中類似的編碼，在GPU上進行了實現(xiàn)，不過該方法只對位置進行了壓縮，并沒有針對處理高度域進行處理.Lindstrom[10]等使用線性預測方法，不過這種方法中每一塊的解壓實際上是順序的過程，因為待解壓的頂點的高度是根據之前壓縮的頂點獲得，無法獲取任意頂點的值.此外，Stookey[11]采用了有損壓縮的方法并在網絡上實現(xiàn)了分布傳輸.Olanda[12]使用小波分塊的金字塔用于在線的三維可視化系統(tǒng).Li[13]等人采用類似于文獻[11]中的方法進行5D的數據壓縮，并使用CUDA解方程.但是上述方法需要迭代過程，本質上仍然是順序的方法.Dore[14,15]等人最近提出了一種基于高度域壓縮的算法，將一塊地形的高度值近似到貝塞爾曲面，只需保存控制點的信息，運行時動態(tài)的計算頂點的高度.但這種方法中的分塊只能是一個層次，而且只對頂點進行了壓縮，并沒有結合位置信息.Niski[16]將高度域作為紋理進程處理，利用圖像壓縮的算法對高度域進行壓縮.但該方法沒有考慮到位置的壓縮，因此不夠靈活.隨著通用圖形處理器(General Purpose GPU,GPGPU)的發(fā)展，研究者使用并行技術實現(xiàn)經典的編碼，如基于CUDA實現(xiàn)的霍夫曼編碼[17]，JPEG2000編碼[18]等.不過這些編碼靈活性有限，而且將其用于地形壓縮時的效果還需要進一步的考慮.

根據上述的研究基礎，本文提出了一種對位置信息和高度信息同時進行壓縮的壓縮方法.與之前的算法不同，應用本文的算法，在獲取任意點的位置和高度時，只需要局部的參考信息，并且結合GPU并行環(huán)境，實現(xiàn)了完全的并行數據處理.實驗結果表明，本文所提算法在數據壓縮比和運行效率上都得到了很大的提高.

3 算法

在地形渲染領域，約束四叉樹是一種常見的數據結構，與不規(guī)則的三角形網格相比，更加適用于基于GPU的環(huán)境.本文所提方法采用約束四叉樹的數據結構，并對地形數據的位置信息采用無損壓縮方法，對地形數據的高度信息可以采用有損或無損的壓縮方法.

地形數據的預處理流程如圖1所示，包含如下幾個步驟：首先是將顯示所需的地形數據以約束四叉樹的結構表示并三角化.然后對地形的位置信息和高度信息分別進行壓縮.

3.1 位置信息壓縮

對于位置信息的壓縮時，借鑒了Gerstiner提出約束四叉樹的序列化方法[3].主要思想如下：對于每個網格，初始化階段時，選擇兩個頂點并記錄其位置信息，然后根據前兩個頂點和標識符計算第三個頂點的信息，重復該步驟直至訪問所有的頂點.

構建鏈表的步驟如下：首先根據觀察者與該塊的距離等信息決定該塊的細分層次，隨后進行三角化.三角化后，選擇其中一個三角形的兩個頂點作為初始頂點，根據該三角形的最后一個頂點的位置和與前兩個頂點的相對位置決定三角形的類型.完成當前三角形的構建后，根據該頂點和其中一個起始頂點(由方向決定)作為下一個三角形的初始頂點，重復上述構建過程直至訪問所有的頂點.

三角形包括三種類型：直角邊到直角邊，直角邊到斜邊，斜邊到直角邊.方向為兩種：順時針和逆時針.圖2表示了這幾種情況，壓縮時，已知頂點為V0和V1，根據V2的位置決定三角形的類型.解壓時，對于直角三角形的三個頂點V0、V1、V2，已知頂點V0和V1，可以根據兩個頂點的位置和三角形的類型計算得到另一個頂點V2的位置信息.

除了頂點的位置信息進行編碼外，還需要處理頂點的高度信息.本文采用了一種分層的方法對高度域進行壓縮，將高度域分為多個層次，然后分別對每一層進行處理.第一層是對整個高度域的粗糙近似，第二層在第一層的基礎上增加了細節(jié)，最后一層是與真實地形的最終的差值，其中第二層和第三層的值可能都為0，因此在實際中可能不需要保存這些數據.

一塊數據可能包含多個子帶，每個子帶包含頭信息和數據信息.一塊數據可能包含多個頭信息和多個頂點，之所以這樣設計是可以靈活的添加和刪除頂點，除非結構有大的變化，無需對結構進行大的調整.

表1為本文所采用的數據結構，可以分為頭信息和每個頂點的信息.頭信息包含當前子帶中的頂點的共有信息，包括頂點個數、子帶中第一個頂點的高度和最后一個頂點的高度.頂點數表明每一個貝塞爾曲線覆蓋多少個頂點.對于地形變化較大的區(qū)域，如果包含的頂點數過多，可能導致差值1和差值2較大，反而不能得到好的壓縮效率.相比于粗糙的區(qū)域，平坦區(qū)域中的貝塞爾曲線可能可以包含更多的頂點也能有較好壓縮比率.除了第一個和最后一個頂點，其他的頂點使用兩個字段保存差值信息，實際應用中，根據實際需求，可以決定是否對兩個差值進行無損壓縮.運行時，根據頭信息中的貝塞爾參數和每個頂點自帶的信息計算相應頂點的高度.

表1 數據結構

3.2 高度域壓縮

上一節(jié)中對數據的編碼和結構進行了介紹，本節(jié)中主要介紹如何使用貝塞爾曲線對高度域進行模擬.由于地形高度域具有規(guī)則性，可以使用一些專門的曲線對其進行近似.貝塞爾曲線是一種常用的構建地形的工具，只需要幾個控制點信息即可計算曲線上頂點的位置.

在大多數的場合中，原始的高度域和根據貝塞爾曲線構建的層次1的差值的范圍集中在0的周圍.這樣的分布在基于預測的方法中是非常常見的.大多數的差值可以使用一個相對較小的位數保存.若還不滿足需求，將最終的值保存在層次3中.此外，層次2和層次3還可以進一步的壓縮.

實際中，可以使用二階或更高階的貝塞爾進行模擬，階數越高，模擬結果越精確.但是如果階數越高，導致計算開銷過大，因此，在實際應用中，通常使用二階的貝塞爾曲線進行模擬，通過修改每個貝塞爾曲線覆蓋的頂點個數修正壓縮率.二階貝塞爾曲線包含兩個頂點V0，Vn和一個控制點C0.其中V0和Vn高度為H0和Hn，C0是貝塞爾參數的控制點的值.

已知頂點V0，Vn的信息，只需計算控制點C0的位置，可以考慮使用最小二乘法對貝塞爾曲線進行擬合.

解壓時，中間頂點的高度通過下述的公式計算得到.每個頂點保存了一個參數t，相應頂點的高度H(t)的計算方法如式(1)所示：

H(t)=(1-t)2H0+2t(1-t)H1+t2H2,t∈[0,1]

(1)

應當指出的是，每個頂點的高度可以并行獲得.已知V1和Vn的高度，其他頂點的高度V2,V3,…,Vn-1可以獨立計算.差值1可以以稀疏矩陣方式存儲在層次2.差值2存儲在第3層.如果為了直接和快速的訪問每個最終的殘留，可以不對差值進行壓縮.

3.3 基于GPU的實現(xiàn)

目前GPU不僅用于繪制，由于其具有很強的浮點數計算能力，被廣泛應用在通用計算領域[19～21].

在基于CUDA的壓縮中，壓縮在不同的核中執(zhí)行，其中一個核計算近似的貝塞爾曲線，另外兩個核計算層次2和層次3，將中間結果保存在訪問速度較快的共享存儲器中.計算近似的貝塞爾曲線時，分別對每個高度域頂點分配一個線程，使用最小二乘法實現(xiàn)對曲線的近似，對于大小為32×32的線程塊，可以計算128個分段的信息.構建層次2時，將高度值和上一個步驟的近似做差值計算，然后將最終的殘留保存在層次3中.

目前的圖形硬件支持可編程的著色器，在渲染每一幀時，將可見的未緩存的塊進行壓縮，并傳輸到顯存中，在GPU中解壓.在對某一塊數據進行解碼時，首先讀取帶頭信息，加載到頂點緩沖區(qū)中.根據帶頭信息得到第一個頂點和最后一個頂點的高度值，結合子帶中兩個頂點V0、V1和下一個頂點的類型得到第三個頂點的位置.然后按照上文中的步驟構建其他的三角形.通過貝塞爾曲線得到近似高度，根據實際需要與兩個差值相加，得到最終的頂點高度.

壓縮和解壓步驟分別可以并行的執(zhí)行，壓縮過程中同時對多個帶進行計算.在解壓過程中，可以分別計算每個段中的頂點的詳細位置，而無需考慮其他的信息.整個系統(tǒng)可以并行的執(zhí)行，壓縮和解壓過程中的最小單位即貝塞爾曲線包含的頂點個數，例如一個貝塞爾曲線包含10個頂點，實現(xiàn)了細粒度的并行.

4 實驗與結果分析

本節(jié)中，我們通過實驗對本文算法進行驗證和分析.數據對象是結構較為復雜的Pudget Sound和結構較為平坦的夏威夷的數據，如圖5所示.實驗運行環(huán)境為：微軟Win7 32位系統(tǒng)，Intel i7-3610，2.3GHz內存，3.14GB內存，顯卡為NVIDIA GeForce GT630M.

在未經壓縮的結構中，每個頂點占用32bit，其中高度域12bit，位置x和y分別是10bit，共32bit.每個三角形占用96bit.經過本文方法的壓縮后，頭文件中，頂點數為6bit，高度1和高度2共為24bit，兩個貝塞爾參數共為16bit.共計42bit.對于單個頂點而言，第一個頂點只需要保存位置信息20bit，第二個頂點20bit和差值2～7bit，共計22～27.頂點3至(n-1)為類型3bit和插值2～7bit，共計5～10bit.最后一個頂點3bit.這些頂點所占用的bit數共計：[5(n+6),10(n+2)]bit，共計可組成n-2個三角形.如果n=8，8個頂點能夠組成6個三角形，占用的空間共計為70～100bit，與之前每個三角形占用96bit相比，壓縮效率為：5.76～8.2.在實際應用中，由于差值可能為0，壓縮效率可能更高.

表2和表3分別表示了Puget Sound和夏威夷的數據使用不同的貝塞爾個數的壓縮效率，第一行表示貝塞爾參數包含的頂點個數，第二行表示其壓縮效率.圖6和圖7是與表2和表3相對應的柱狀圖的顯示.橫軸表示貝塞爾曲線包含不同的頂點時的壓縮效率，縱軸表示相應的壓縮率.

表2 Puget Sound壓縮結果

表3 夏威夷壓縮結果

圖6和圖7中的數據表明，對于單個地形而言，貝塞爾曲線包含的頂點增多，壓縮率會逐步提高，但是到達某臨界值之后，壓縮效率會有所下降，需要額外空間保存差值2和差值3的值.對于不同的地形而言，壓縮效率也有所不同，在地形較為平坦的區(qū)域，壓縮的效率更高.而且出現(xiàn)拐點的貝塞爾曲線包含的頂點的值也不同，平坦區(qū)域的拐點的值大于粗糙的區(qū)域的拐點的值.總體上而言，平坦的數據的壓縮效果要好于結構復雜的數據的壓縮效果.

使用經典的算法ZIP對數據進行壓縮，其壓縮比率大概為2.5∶1，而本文所提算法相對于結構復雜的數據的平均壓縮比率為6.3∶1，相對于平坦的數據的平均壓縮比率為7.06∶1.可以看出，本文所提算法具有更好的壓縮效果.由于壓縮過程中包含復雜的計算，因此解壓速度要快于壓縮速度.不過貝塞爾曲線包含的頂點個數對壓縮和解壓的速度不會有太大的影響，這是因為實際的操作基本是基于每個頂點執(zhí)行.實際應用中，每秒可以解壓108的頂點，滿足了實時的需求.如果提高壓縮的速度，可以更進一步的提高吞吐量.

5 結語

本文提出了一種新的高度域快速壓縮和解壓方法，可以充分利用GPU的并行功能.實現(xiàn)了有損和無損的壓縮.我們的方法在兩個無關的數據中進行了試驗，實驗結果表明，能夠實現(xiàn)非常高效的壓縮效果.與傳統(tǒng)的壓縮方法相比，我們的壓縮率要遠遠高于其壓縮效率，隨著數據規(guī)模的增大，可能實現(xiàn)進一步的提升.

本文的方法可以用于大規(guī)模的地形數據處理的二次處理中.解壓速度和經典的基于GPU的解壓速度相當.此外，系統(tǒng)還可以與其他壓縮方法相結合.下一步的工作中，我們將考慮更高階的貝塞爾曲線能否取得更好的壓縮效率，以及如何快速的計算高階貝塞爾曲線的控制點.

[1]Duchaineau M,Wolinsky M,Sigeti D E,et al.ROAMing terrain:real-time optimally adapting meshes[A].Visualization'97[C].Phoenix,AZ,USA:IEEE,1997.81-88.

[2]Hoppe H.View-dependent refinement of progressive meshes[A].ProcSiggraph[C].Los Angeles,USA:ACM,1997.189-198.

[3]Gerstner T.Multi-resolution visualization and compression of global topographic data[J].Geoinformatica,2003,7(1):7-32.

[4]Kim J K,Ra J B.A real-time terrain visualization algorithm using wavelet-based compression[J].Visual Computer,2004,20(2):67-85.

[5]Xie Z,Marcus A,Randolph F,Barbara C,et al.Progressive transmission of lossily compressed terrain[A].Conferencialatinoamericana de informática (CLEI 2008)[C].Santa Fe,Argentina:2008.8-12.

[6]Kidner B,Derek S.Advances in the data compression of digital elevation models[J].Computers & Geosciences,2003,29(8):985-1002.

[7]Inanc M.Compressing terrain elevation datasets[D].NY,USA:Rensselaer Polytechnic Institute Troy,2008.

[8]Kidner B,Derek S.Storageefficient techniques for representing digital terrain models[A].Innovations in GIS 4[C].London,England:Taylor & Francis,1997.25-41.

[9]Dick C,Schneider J,Westermann R.Efficient geometry compression for GPU-based decoding in realtime terrain rendering[J].Computer Graphics Forum,2009,28(1):67-83.

[10]Lindstrom P,Cohen D.On-the-fly decompression and rendering of multiresolution terrain[A].Proceedings of the 2010 ACM SIGGRAPH Symposium on Interactive 3D Graphics and Games[C].Bethesda,MD,USA:ACM,2010.65-73.

[11]Jared S,Randolph F,Xie Z,et al.Parallel ODETLAP for terrain compression and reconstruction[A].Proceedings of the 16th ACM SIGSPATIAL International Conference on Advances in Geographic Information Systems[C].Irvine,California,USA:ACM,2008.17.

[12]Olanda R,Pérez M,et al.Terrain data compression using wavelet-tiled pyramids for online 3D terrain visualization[J].International Journal of Geographical Information Science,2013,28(2):407-425.

[13]Li Y.CUDA-accelerated HD-ODETLAP:A high dimensional geospatial data compression framework[D].NY,USA:Rensselaer Polytechnic Institute Troy,2011.

[16]Niski K,Purnomo B,Cohen J.Multi-grained level of detail using a hierarchical seamless texture atlas[A].Proceedings of the 2007 Symposium on Interactive 3D Graphics and Games[C].Seattle,WA,USA:ACM,2007.153-160.

[17]Rahmani H,CihanT,Cuneyt A.A parallel huffman coder on the CUDA architecture[A].Visual Communications and Image Processing Conference,2014 IEEE[C].Valletta,Malta:IEEE,2014.311-314.

[18]Lee J,Kim B,Yoon K.CUDA-based JPEG2000 encoding scheme[A].International Conference on Advanced Communication Technology[C].Pyeongchang,Korea:IEEE,2014.671-674.

[19]趙星,胡晶晶,潘曉川,張朋.一種新的基于GPU實現(xiàn)的錐束CT正投影算法[J].電子學報,2009,37(6):1165-1169. Zhao Xing,Hu Jingjing,Pan Xiaochuan,Zhang Peng.A novel GPU based cone beam CT forward projection method[J].Acta Electronica Sinica,2009,37(6):1165-1169.(in Chinese)

[20]楊正龍,金林,李蔚清.基于GPU的圖形電磁計算加速算法[J].電子學報,2007,35(6):1056-1060. Yang Zhenglong,Jin Lin,Li Weiqing.Accelerated GRECO based on GPU[J].Acta Electronica Sinica,2007,35(6):1056-1060.(in Chinese)

[21]王綱,季振洲,張澤旭.大規(guī)模真實感雪景實時渲染[J].電子學報,2012,40(9):1746-1751. Wang Gang,Ji Zhenzhou,Zhang Zexu.Large scale realistic snow scene real-time rendering[J].Acta Electronica Sinica,2012,40(9):1746-1751.(in Chinese)

翟 銳 男，1981年生于河南焦作,中國科學院大學在讀博士生，主要研究方向為數字圖像處理、計算機圖形學、三維可視化.

E-mail:zhairui11b@mails.ucas.ac.cn

呂 科 男，1971出生于寧夏西吉，教授，博士生導師，主要研究方向為數字圖像處理、計算機圖形學、智能信息處理技術.

E-mail:luk@ucas.ac.cn

Research on Terrain Height Field Compression Algorithm

ZHAI Rui,Lü Ke,DAI Shuang-feng,PAN Wei-guo

(CollegeofEngineeringandInformationTechnology,UniversityofChineseAcademyofSciences,Beijing100049,China)

With the development of remote sensing technology,the size of terrain is growing rapidly,and far beyond the scope of main memory,has become an urgent problem.Data compression is a popular technology to increase system throughput.With the rapid development of GPU(Graphics Processing Unit) technology,the traditional compression algorithms cannot take full advantage of the ability of the current GPU.In this paper,we propose a GPU-based terrain data compression method,and achieve a high rate compression of terrain height field and location.Comparing to other algorithms,the main contribution of our algorithm is that the compression of terrain height filed and position is executed in the same time,and all the information of a node can be calculated according to the presentstrip.For terrain height domain,we firstly make Bezier curve approximation,then save the difference.After the steps above,we can achieve high compression ratio.By comparison with traditional methods,we get reasonable experimental results.

data compression;terrain rendering;graphics processing unit (GPU)

2015-04-09;

2015-05-05;責任編輯：梅志強

國家自然科學基金(No.61271435);北京市自然科學基金重點項目(No.4141003)

TP302.7

A

0372-2112 (2016)12-2894-06

??學報URL:http://www.ejournal.org.cn

10.3969/j.issn.0372-2112.2016.12.012