基于ＧＰＵ的等值面提取與繪制

（國防科學技術大學 多媒體研發中心， 長沙 410073）



摘要：利用圖形硬件的并行性將六面體網格數據映射為紋理，從每個六面體中提取等值面片，并將其繪制到紋理而得到最終等值面。基于Cg著色器編程語言實現三維電磁環境表現的實驗結果表明，該方法有效地減輕了CPU負擔，提高了等值面提取速度，適合實時應用。

關鍵詞：等值面；移動立方體；硬件加速

中圖分類號：TP391文獻標志碼：A

文章編號：1001-3695(2008)11-3468-03



Isosurface extraction and rendition using GPU



WU Ling-da，YANG Chao，CHEN Peng

（Multimedia RD Center， National University of Defense Technology， Changsha 410073， China）

Abstract:The development of GPU gave a new platform for general purpose computations. First the input cube-dataset was mapped to textures. Then based on the parallelity of graphic hardware， a new algorithm for isosurface extraction from cubes was designed.At last the result isosurface was rendered to output texture. Experiment result of representation for 3D electromagnetic environment based on Cg shader programming language indicates that the method effectively lightens CPU’s burden， accelerates isosurface extraction process and can be used in real-time application.

Key words:isosurface; marching cubes; hardware accelerated



0引言

自從1987年Lorensen等人^[1]提出marching cubes（MC）算法用于三維醫學數據等值面可視化以來，國內外學者對等值面提取進行了深入廣泛的研究，取得了很大進展。基于CPU實現的等值面提取算法，按照所使用數據結構和搜索算法的不同，基本上可以分為三類：a)基于數據值區間劃分，其典型代表是1991年Gallagher^[2]提出的span-filtering方法。b)基于空間坐標劃分，典型代表是1992年Wilhelms等人^[3]提出的八叉樹分解方法。這兩種方法均能有效地減少與等值面相交單元的搜索時間，從而提高了等值面提取的速度。c）基于曲面的一致性，如Hung和Itoh等人利用曲面的鄰接關系從一個很小的種子集開始不斷增長生成最終等值面^[4，5]。

隨著圖形硬件的發展，GPU(graphics processing units，圖形處理器)的原始計算能力超出CPU峰值性能至少一個數量級，可編程能力和浮點支持能力在不斷提高。GPU不再僅僅局限于處理特定的圖形繪制任務，利用其高度并行計算和多數據流處理能力，逐步成為能夠輔助CPU計算的通用計算單元。實時繪制語言的出現無疑更是加快了GPU在通用計算上的應用，如Cg、HLSL以及GLSL或其他更高級的面向通用計算的繪制語言^[6]的出現，使得人們更容易編程實現這些運算。采用圖形硬件作為通用計算的動力來自這些新硬件所具有的一定的并行性、高密集的運算、減少了GPU 與CPU 的數據通信^[7]等優勢。

過去的幾年中，很多人利用GPU提供的并行性將它作為流處理器來做一些通用計算方面的工作，甚至用來求解有限差分方程組。例如Matsumura等人^[8]和Thomas等人^[9]分別在2003和2004年嘗試利用GPU提取和繪制等值面，但是當時由于顯存以及GPU功能的限制，等值面提取和繪制的過程并不能完全在圖形硬件上實現。本文將基于六面體的等值面提取和繪制過程完全移到GPU上執行，從而降低了CPU計算負荷，同時減少了內存與GPU之間的數據傳輸量。

1等值面提取算法分析

MC 算法處理的基本單位是邏輯立方體，基本思想是逐一處理數據場中的每個立方體，找出與等值面相交的所有立方體，采用線性插值計算出等值面與立方體邊的交點。根據立方體每一個頂點與等值面的相對位置，將等值面與立方體邊的交點按一定方式連接生成等值面，作為整個等值面在該立方體內的一個逼近表示，將所有的單元連接在一起就形成了最終的等值面。

雖然MC算法中立方體與等值面相交的情況共有28=256種之多，但就等值面提取本質而言，可以分為兩個階段，即搜索與查詢閾值相交的單元和生成等值面幾何圖元。針對第一階段，前述三類基于CPU的等值面提取算法分別利用空間坐標進行數據劃分、利用值范圍進行數據劃分、利用曲面一致性來提高等值面提取的速度，減少所占用的內存空間。對于第二個階段，通常利用查找表的方法來提高處理速度，如果僅從軟件的角度考慮，其可以提升的空間很小，因此很少有人從第二個階段去研究如何提高效率。然而在等值面提取的第二個階段，幾何圖元索引生成、坐標點插值以及法向量插值等關鍵運算目前基本都是在CPU和內存中進行。隨著圖形硬件的發展，利用其高度的并行性和高密集的運算能力來做通用計算方面的工作已不再新奇，而且各種等值面提取算法都需要對大量體元進行獨立的處理，這是非常適合于應用高度并行運算的。本文正是從這一點出發，充分利用GPU高密集和高并行的運算能力，從硬件的角度針對第二個階段來加速等值面提取和繪制，減少CPU的負擔，提高整體效率。

2基于GPU的等值面提取算法

經典的MC算法中，頂點值分布情況一共有28=256種，基于CPU的算法通常使用查找表的方法來加快提取速度。本章說明如何在GPU上的片斷著色程序中實現基于六面體的等值面提取。圖1給出了從一個六面體通過硬件加速提取等值面幾何圖元并輸出其中一個插值點的過程。

首先將頂點數據、六面體網格數據映射為紋理輸入以供GPU讀取，然后將輸出等值面的各幾何圖元（四邊形）映射到輸出紋理，以真正地在GPU上執行等值面提取與繪制。

21輸入輸出數據到紋理的映射

對于輸入的雷達波損失數據，首先將其映射為三維紋理。采用32位浮點的RGBA格式，每個數據點對應一個紋理單元，用來存儲該位置的坐標和雷達波損失值。

根據等值面與六面體的相交情況，對于立方體頂點函數值分布的256種情況，根據互補和旋轉對稱性，可以減少為15種，如圖2所示。從這15種狀態中可以發現，無論等值面與六面體如何相交，最終生成的三角形等值面片最多只有4個。同時，當在圖形管線中繪制一個三角形時，如果輸入三角形的三個頂點相同，在柵格化過程中將去掉該圖元。因此，可以用一種一致的方式來存儲每個六面體輸出的等值面片幾何圖元——對每個六面體都使用12個紋理單元（32位RGB格式）來保存輸出數據，即所生成的4個三角形等值面片的12個頂點，每個紋理單元對應一個三角形頂點。

22等值面提取

在將輸入數據映射到紋理并設計好輸出數據結構后，就可以在片段著色程序中實現等值面提取。從前面設計的輸入紋理以及輸出結構可知，每個六面體對應12個輸出像素。假設輸入紋理大小為ITX×ITY×ITZ，則輸出紋理大小為12×（ITX-1）×（ITY-1）×（ITZ-1）。對于紋理坐標為(x，y，z)的當前輸出像素，對應的六面體8個頂點的紋理坐標分別為(t，y，z)、(t+1，y，z)、(t，y+1，z)、(t+1，y+1，z)、(t，y，z+1)、(t+1，y，z+1)、(t，y+1，z+1)和(t+1，y+1，z+1)。其中t=x mod 12。由此就可以從輸入紋理中查找到這8個頂點對應的損失值Si（i=0…7）；再按照閾值c可以將六面體8個頂點進行分類。通過式（1）來計算掩碼值β以記錄分類結果：

βi=0Si＜c

1Si≥c（i=0，…，7）（1）

其中:βi表示β的第i位。根據這里計算的比特掩碼就可以決定六面體與等值面相交的邊。在β所有可能的256種狀態中，等值面與六面體相交形成的三角形面片最多只有4個，當少于4個時認為其中包含了重復頂點。對于每一種情況，按照MC算法，三角形面片的位置是確定的，因此頂點位置也是固定的，于是引入三角形頂點索引π=0，1，…，11來分別標志這4個三角形的12個頂點（可能包含重復頂點）。通過該索引就可以定位每個三角形三條邊的端點所在的六面體棱邊，進而可以確定棱邊端點在輸入紋理中的坐標。圖3展示了狀態為01001001的情況，若當前處理的是第二個三角形，即π=3，4，5，則其3個頂點所在的棱邊為(0~1)、(1~3)和(1~5)。根據其所在的六面體就可以得到6個端點在輸入紋理中的坐標分別是（t，y，z）、（t+1，y，z）；（t+1，y，z）、（t+1，y+1，z）；（t+1，y，z）、（t+1，y，z+1）。其中t=x mod 12。

將所有15×4×3種關系列表（表1）并以索引紋理的形式保存。采用8位浮點RGB格式，每一種六面體頂點函數值分布狀態對應24個紋理單元；這24個紋理單元分別對應六面體與等值面的4個相交三角形的12個頂點所在六面體的棱邊。一個紋理單元存儲一個三角形頂點所在棱邊的一個端點在六面體中的局部頂點序號0~7。這里需要強調的是，本文假定每個六面體和等值面有4個相交三角形，少于4個的情況則認為存在有三角形頂點重合的情形。

當得到了相交三角形3個頂點所在六面體棱邊的6個端點在輸入紋理中的坐標后，就可以根據紋理坐標計算相應的坐標值Vi（i=0，…，5）并在紋理中查找到相應函數值Si（i=0，…，5），那么輸出三角形3個頂點的位置就可以通過式（2）線性插值得到：

Vi=tiV2i+(1-ti)V2i+1(2)

其中ti=（c-S2i）/（S2i+1-S2i）；0≤i≤2。

對于等值面法線的計算，先在預處理階段計算每個頂點的梯度值，得到等值面坐標后插值計算等值面各頂點的法向。這個過程與頂點坐標插值完全類似，只需要在片段處理程序中簡單地加入法向插值代碼即可以實現。這里不再給出詳細說明。

23程序實現

在MS WindowsXP平臺基于C++和OpenGL實現該算法。使用的顯卡芯片是ATI Radeon X1900 GT，核心頻率為575 MHz，顯存頻率為1.2 GHz，顯存為256 MB。這雖不是頂級顯卡，但具有36個像素渲染單元、8個頂點渲染單元、12條像素渲染管線，完全支持微軟DirectX 9.0和OpenGL 2.0。

具體運用Cg語言編寫的像素著色器程序代碼如下：

3實驗結果

本文將該方法在電磁環境三維可視化應用系統中進行了實驗和應用。體數據來源于自由空間下電磁波傳播方程計算生成的三維空間電磁波功率密度值^[10]。其中三維環境中隨機布置了20個電磁設備，每個電磁設備采用全向天線。具體參數如表2所示。針對體數據分辨率分別為100×100×100和300×100×100的兩種情況，對在CPU和通過GPU加速的運行結果進行了比較，如表3所示。

由于所使用的硬件和文獻[9]的不同，無法進行比較。但可以肯定的是，本文的GPU要先進一些，效率肯定有所提高。實驗最終生成的三維電磁環境效果如圖4、5所示。

4結束語

不同于以往單純地從軟件的角度加速等值面的提取方法，本文充分利用當前圖形硬件在通用計算方面的能力，提出了基于圖形處理器的等值面提取算法，適合于任意六面體網格形式的體數據。利用Cg語言實現三維電磁環境可視化表現的實驗證明，該算法在顯存允許的條件下，可以用于實時等值面提取與繪制，對體數據實時可視化研究具有一定的意義。但是由于顯卡內存有限，對于超過顯卡內存容量的體數據，如何利用硬件進行加速提取和繪制仍然有待進一步的研究。

參考文獻：

［1］

LORENSEN W E，CLINE H E.Marching cubes:a high resolution 3D surfaces construction algorithm[J].Computer Graphics，1987，21(4):163-169.

[2]GALLAGHER R S.Span filtering: an optimization scheme for volume visualization of large finite element models[C]//Proc of the 2nd Conference on Visualization.Los Alamitos:IEEE Computer Society Press，1991.

[3]WILHELMS J，GELDER A V.Octrees for faster isosurface generation[J].ACM Trans on Graphics，1992，11(3):201-227.

[4]HUNG C H，YANG Chuan-kai.A simple and novel seed-set finding approach for isosurface extraction[C]//Proc of Eurographics/IEEE-VGTC Symposium on Visualization.2005:125-132.

[5]ITOH T，KOYAMADA K.Automatic isosurface propagation using an extrema graph and sorted boundary cell lists[J].IEEE Trans on Visualization and Computer Graphics，1995，1(4): 319-327.

[6]BUCK I.Data parallel computing on graphics hardware[EB/OL].(2007).http://graphics.stanford.edu/projects/brookgpu/. 

[7]吳恩華，柳有權.基于圖形處理器(GPU)的通用計算[J].計算機輔助設計與圖形學學報，2004，16(5):601-612.

[8]MATSUMURA M，ANJYO K.Accelerated isosurface polygonization for dynamic volume data using programmable graphics hardware in visua-lization and data analysis[C]//Proc of Conference on Visualization and Data Analysis.2003:145-152.

[9]THOMAS K，SIMON S，THOMAS E.Hardware-accelerated reconstruction of polygonal isosurface representations on unstructured grids[C]//Proc of the 12th Pacific Conference on Computer Graphics and Applications.Washington DC:IEEE Computer Society，2004:185-196.

[10]BLAKE L V.雷達距離性能分析[M].吳秉瑋，趙楊，劉元林，等譯.南京：機械電子工業部第十四研究所，1990.