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

（國防科學(xué)技術(shù)大學(xué) 多媒體研發(fā)中心， 長沙 410073）



摘要：利用圖形硬件的并行性將六面體網(wǎng)格數(shù)據(jù)映射為紋理，從每個六面體中提取等值面片，并將其繪制到紋理而得到最終等值面?；贑g著色器編程語言實現(xiàn)三維電磁環(huán)境表現(xiàn)的實驗結(jié)果表明，該方法有效地減輕了CPU負擔(dān)，提高了等值面提取速度，適合實時應(yīng)用。

關(guān)鍵詞：等值面；移動立方體；硬件加速

中圖分類號：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）算法用于三維醫(yī)學(xué)數(shù)據(jù)等值面可視化以來，國內(nèi)外學(xué)者對等值面提取進行了深入廣泛的研究，取得了很大進展。基于CPU實現(xiàn)的等值面提取算法，按照所使用數(shù)據(jù)結(jié)構(gòu)和搜索算法的不同，基本上可以分為三類：a)基于數(shù)據(jù)值區(qū)間劃分，其典型代表是1991年Gallagher^[2]提出的span-filtering方法。b)基于空間坐標劃分，典型代表是1992年Wilhelms等人^[3]提出的八叉樹分解方法。這兩種方法均能有效地減少與等值面相交單元的搜索時間，從而提高了等值面提取的速度。c）基于曲面的一致性，如Hung和Itoh等人利用曲面的鄰接關(guān)系從一個很小的種子集開始不斷增長生成最終等值面^[4，5]。

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

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

1等值面提取算法分析

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

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

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

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

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

21輸入輸出數(shù)據(jù)到紋理的映射

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

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

22等值面提取

在將輸入數(shù)據(jù)映射到紋理并設(shè)計好輸出數(shù)據(jù)結(jié)構(gòu)后，就可以在片段著色程序中實現(xiàn)等值面提取。從前面設(shè)計的輸入紋理以及輸出結(jié)構(gòu)可知，每個六面體對應(yīng)12個輸出像素。假設(shè)輸入紋理大小為ITX×ITY×ITZ，則輸出紋理大小為12×（ITX-1）×（ITY-1）×（ITZ-1）。對于紋理坐標為(x，y，z)的當前輸出像素，對應(yīng)的六面體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個頂點對應(yīng)的損失值Si（i=0…7）；再按照閾值c可以將六面體8個頂點進行分類。通過式（1）來計算掩碼值β以記錄分類結(jié)果：

βi=0Si＜c

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

其中:βi表示β的第i位。根據(jù)這里計算的比特掩碼就可以決定六面體與等值面相交的邊。在β所有可能的256種狀態(tài)中，等值面與六面體相交形成的三角形面片最多只有4個，當少于4個時認為其中包含了重復(fù)頂點。對于每一種情況，按照MC算法，三角形面片的位置是確定的，因此頂點位置也是固定的，于是引入三角形頂點索引π=0，1，…，11來分別標志這4個三角形的12個頂點（可能包含重復(fù)頂點）。通過該索引就可以定位每個三角形三條邊的端點所在的六面體棱邊，進而可以確定棱邊端點在輸入紋理中的坐標。圖3展示了狀態(tài)為01001001的情況，若當前處理的是第二個三角形，即π=3，4，5，則其3個頂點所在的棱邊為(0~1)、(1~3)和(1~5)。根據(jù)其所在的六面體就可以得到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種關(guān)系列表（表1）并以索引紋理的形式保存。采用8位浮點RGB格式，每一種六面體頂點函數(shù)值分布狀態(tài)對應(yīng)24個紋理單元；這24個紋理單元分別對應(yīng)六面體與等值面的4個相交三角形的12個頂點所在六面體的棱邊。一個紋理單元存儲一個三角形頂點所在棱邊的一個端點在六面體中的局部頂點序號0~7。這里需要強調(diào)的是，本文假定每個六面體和等值面有4個相交三角形，少于4個的情況則認為存在有三角形頂點重合的情形。

當?shù)玫搅讼嘟蝗切?個頂點所在六面體棱邊的6個端點在輸入紋理中的坐標后，就可以根據(jù)紋理坐標計算相應(yīng)的坐標值Vi（i=0，…，5）并在紋理中查找到相應(yīng)函數(shù)值Si（i=0，…，5），那么輸出三角形3個頂點的位置就可以通過式（2）線性插值得到：

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

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

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

23程序?qū)崿F(xiàn)

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

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

3實驗結(jié)果

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

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

4結(jié)束語

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

參考文獻：

［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]吳恩華，柳有權(quán).基于圖形處理器(GPU)的通用計算[J].計算機輔助設(shè)計與圖形學(xué)學(xué)報，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].吳秉瑋，趙楊，劉元林，等譯.南京：機械電子工業(yè)部第十四研究所，1990.