周大鑫，周茂林，鄒 文，魯 才

(1.電子科技大學(xué) 計(jì)算機(jī)科學(xué)與工程學(xué)院，四川成都 611731;2.川慶鉆探工程有限公司 地球物理勘探公司，四川成都 610213;3.電子科技大學(xué)寬帶光纖傳感與通信教育部重點(diǎn)實(shí)驗(yàn)室，成都 611731)

基于立方體元的Shear－warp體繪制加速算法

周大鑫1，周茂林2，鄒 文2，魯 才3

(1.電子科技大學(xué) 計(jì)算機(jī)科學(xué)與工程學(xué)院，四川成都 611731;2.川慶鉆探工程有限公司 地球物理勘探公司，四川成都 610213;3.電子科技大學(xué)寬帶光纖傳感與通信教育部重點(diǎn)實(shí)驗(yàn)室，成都 611731)

大規(guī)模三維地震數(shù)據(jù)的直接體繪制，是一個(gè)計(jì)算和數(shù)據(jù)雙重密集型問題。為了提高三維圖像的重建效率，提出了一種基于立方體元的Shear－warp地震數(shù)據(jù)體繪制算法。該算法通過構(gòu)建立方體元在相鄰的體素點(diǎn)之間建立聯(lián)系;然后根據(jù)地震數(shù)據(jù)的特征對(duì)體元進(jìn)行分類，在繪制過程中通過二叉樹索引，快速定位分類結(jié)果。通過索引結(jié)果，避免了對(duì)等值體元的插值計(jì)算和跳過透明體元，提高了繪制速度。經(jīng)實(shí)驗(yàn)結(jié)果表明，Shear-warp算法得到了有效加速。

可視化;體繪制;地震數(shù)據(jù);二叉樹;Shear-warp

0 前言

直接體繪制是近年來迅速發(fā)展起來的，對(duì)三維數(shù)據(jù)場(chǎng)的一種體可視化繪制方法。該方法可以非常清晰地反應(yīng)出三維數(shù)據(jù)場(chǎng)的細(xì)節(jié)，并且能夠?qū)⒚胬L制不能顯示的三維數(shù)據(jù)內(nèi)部構(gòu)造顯示出來。由于直接體繪制能夠?qū)碜杂诘刭|(zhì)勘探、核磁共振、計(jì)算機(jī)斷層掃描、計(jì)算機(jī)流體力學(xué)等的三維體數(shù)據(jù)進(jìn)行有效的繪制，所以體繪制技術(shù)在石油勘探、生物醫(yī)學(xué)工程、地質(zhì)科學(xué)、氣象、有限元分析后處理、化學(xué)等領(lǐng)域得到了廣泛地應(yīng)用和重視［1］。近年來，體繪制技術(shù)又被應(yīng)用于考古和文物保護(hù)、宇航、地震預(yù)測(cè)、自然現(xiàn)象仿真等領(lǐng)域。

三維地震是目前石油勘探的主要手段，由于三維地震數(shù)據(jù)的數(shù)據(jù)量很大(例如一個(gè)普通三維工區(qū)有1024×1024個(gè)道，每道有5000個(gè)點(diǎn)，每個(gè)點(diǎn)為一個(gè)浮點(diǎn)數(shù)存儲(chǔ)，其大小為20G，相對(duì)現(xiàn)有的普通微機(jī)內(nèi)存是很難將其全部載入)，因此其體繪制是一個(gè)計(jì)算和數(shù)據(jù)雙重密集型問題。如何在不損失最終繪制結(jié)果質(zhì)量的前提下，通過提高體繪制的成像速度來達(dá)到交互式體繪制，對(duì)于推動(dòng)直接體繪制技術(shù)在石油勘探領(lǐng)域進(jìn)一步應(yīng)用具有重要的實(shí)際意義，并已成為國(guó)內(nèi)、外石油工業(yè)界和學(xué)術(shù)界共同關(guān)心的前沿問題。

傳統(tǒng)的體繪制方法有光線投影算法(Ray-Casting)［2］，足跡表法 (Splatting)［3］和錯(cuò)切變換(Shear-warp)［4］等。近年來，國(guó)內(nèi)、外學(xué)者也提出了很多加速體繪制的方法，如包圍盒空間八叉樹(Octree)［5］等數(shù)據(jù)結(jié)構(gòu)和提前不透明度截至［6］等策略來加速體繪制算法，也有通過基于三維紋理硬件［7］、基于可編程圖形處理器(GPU)［8］等硬件方式加速體繪制的方法。但由于計(jì)算量大，計(jì)算復(fù)雜度高，很難達(dá)到交互式體繪制的要求。基于硬件的體繪制算法，由于硬件價(jià)格昂貴和紋理顯存有限等原因也難以普及。相對(duì)而言，目前Shear-warp算法的速度最快［9］。

Shear-warp算法通過矩陣變換的方式，將繪制過程分解成三維數(shù)據(jù)場(chǎng)shear和二維圖像的warp兩個(gè)過程，以達(dá)到將三維數(shù)據(jù)場(chǎng)的重采樣和插值過程轉(zhuǎn)變?yōu)閷?duì)二維平面的重采樣和插值過程，從而簡(jiǎn)化了計(jì)算復(fù)雜度和計(jì)算量。

在Shear-warp算法實(shí)現(xiàn)過程中，可以采用多種加速方法。如文獻(xiàn)［5］中的八叉樹編碼體數(shù)據(jù)，實(shí)現(xiàn)對(duì)體數(shù)據(jù)的快速編碼，還有游程編碼(RLE)［10］數(shù)據(jù)結(jié)構(gòu)加速的方式，都能跳過大規(guī)?！巴该鳌钡臒o效體素，從而達(dá)到加速算法的目的。但是，當(dāng)體數(shù)據(jù)的透明度發(fā)生變化之后，八叉樹編碼和RLE都要重新對(duì)體數(shù)據(jù)進(jìn)行編碼，需要額外的預(yù)處理時(shí)間，這樣對(duì)體繪制的交互性造成了影響。而且這兩種方法都只考慮了體數(shù)據(jù)中的透明無效體素，沒有考慮到在體數(shù)據(jù)中存在的相鄰數(shù)據(jù)值相同的情況。同時(shí)RLE方法需要在每個(gè)主方向存儲(chǔ)一套切片，這對(duì)于大規(guī)模的三維地震數(shù)據(jù)而言，增大了內(nèi)存消耗，從而增加了數(shù)據(jù)讀取時(shí)間。

針對(duì)上述問題，作者在本文中提出了基于立方體元的Shear-warp地震數(shù)據(jù)體繪制算法。該算法通過立方體元的方式，使相鄰的體素點(diǎn)之間建立聯(lián)系，從而在計(jì)算的過程中，能夠有效地利用地震數(shù)據(jù)空間關(guān)系對(duì)體繪制進(jìn)行加速。體元作為存儲(chǔ)和繪制的基本單位，當(dāng)繪制的主方向發(fā)生變化時(shí)，通過改變Shear-warp算法掃描線的方向，對(duì)一份數(shù)據(jù)進(jìn)行操作，相對(duì)RLE的編碼方式節(jié)省了存儲(chǔ)空間。同時(shí)，可根據(jù)地震數(shù)據(jù)的特點(diǎn)對(duì)立方體元進(jìn)行分類，然后采用二叉樹對(duì)分類進(jìn)行快速索引。通過索引結(jié)果，在繪制過程中不僅跳過了“透明”的無效體素，而且還避免了對(duì)相同地質(zhì)體中重采樣點(diǎn)的插值計(jì)算，從而提高了Shear-warp算法的繪制速度。

1 利用立方體元實(shí)現(xiàn)Shear-warp加速算法

Shear-warp地震數(shù)據(jù)體會(huì)制算法的整個(gè)流程，可以劃分為兩個(gè)階段:

(1)數(shù)據(jù)預(yù)處理階段。

與大型出版企業(yè)相比，中小出版企業(yè)在資源壟斷、出版許可、網(wǎng)點(diǎn)布局、文化沉淀、社會(huì)責(zé)任等方面差距較大。因此，中小出版企業(yè)轉(zhuǎn)型升級(jí)不能效仿大企業(yè)的“高、大、上、全”。

(2)繪制階段。

如圖1所示，預(yù)處理階段包括立方體元數(shù)據(jù)結(jié)構(gòu)的壓縮和存儲(chǔ)，二叉樹索引結(jié)構(gòu)的建立。繪制階段包括透明立方體元的跳躍，對(duì)不同體元的插值處理方式，以及顏色和不透明度的映射。

1.1 數(shù)據(jù)預(yù)處理

1.1.1 地震數(shù)據(jù)壓縮和道的立方體元轉(zhuǎn)換

對(duì)于三維數(shù)據(jù)場(chǎng)D的數(shù)據(jù)分布，通常可以通過某個(gè)函數(shù)f3D(x)來定義，為了達(dá)到壓縮三維數(shù)據(jù)場(chǎng)的目的，我們通過門檻分類法把數(shù)據(jù)分為256個(gè)區(qū)間。如n+1個(gè)門檻值0=S0＜S1＜S2、…、＜Sn－1＜Sn=256，把數(shù)據(jù)場(chǎng)分為n部份，對(duì)數(shù)據(jù)場(chǎng) D=D0∪ D1∪ 、…、∪ Dn－1。其中Di={x|Si＜=f3D(x)＜Si+1}，i=0、1…、n－1。函數(shù)f3D(x)就是將每個(gè)采樣點(diǎn)x映射到特定的區(qū)間［Si，Si+1)。映射函數(shù)f3D(x)如下:

f3D(x)=

圖1 繪制流程Fig． 1 Rendering process

式中 dmin、dmax為原始數(shù)據(jù)的最小值與最大值;x為原始數(shù)據(jù);f3D(x)為壓縮后的數(shù)據(jù)值。

三維地震數(shù)據(jù)的數(shù)據(jù)值分別代表地質(zhì)體的不同屬性，如沙層、巖層等之間不同的密度值。而三維地震數(shù)據(jù)中有大量相鄰的相同地質(zhì)體存在(如下頁圖2所示)，合理利用地震數(shù)據(jù)中相同地質(zhì)體的空間關(guān)系，可以對(duì)算法進(jìn)行有效加速。原始的三維地震數(shù)據(jù)通常是以點(diǎn)為單位，然后按道的方式進(jìn)行存儲(chǔ)(見下頁圖3(a))，由于點(diǎn)的結(jié)構(gòu)難以反映三維地震數(shù)據(jù)的空間關(guān)系。因此作者采用立方體元，為基本單位的結(jié)構(gòu)進(jìn)行存儲(chǔ)，將以點(diǎn)為單位的地震道轉(zhuǎn)換為以體元為單位的數(shù)據(jù)道(如圖3所示)。使每個(gè)體數(shù)據(jù)點(diǎn)和周圍七個(gè)點(diǎn)通過立方體的方式發(fā)生聯(lián)系，同時(shí)按體元的方式對(duì)體數(shù)據(jù)進(jìn)行存儲(chǔ)和處理。

圖2 三維地震數(shù)據(jù)剖面Fig． 2 3D seismic data profile

圖3 地震道到體元道的轉(zhuǎn)換Fig． 3 Transition of seismic trace to volume element

因?yàn)橄嗤牡刭|(zhì)體具有相同的體數(shù)據(jù)值，而相同的體數(shù)據(jù)在通過傳遞函數(shù)映射之后，必然得到相同的顏色值和透明度。Shear-warp算法對(duì)重采樣點(diǎn)的值進(jìn)行計(jì)算的時(shí)候，采用的是線性插值算法，若C為AB連線上的一點(diǎn)，公式(2)為一次線性插值算法，其中 va和 vb表示點(diǎn) a和 b的值，α=CB/AB。若 va=vb，那么 vc=va=vb。

因此如果一個(gè)立方體元的八個(gè)頂點(diǎn)值相同或值相近，對(duì)在這個(gè)立方體元六個(gè)面上的重采樣點(diǎn)的值，就與頂點(diǎn)的值相同。在此我們認(rèn)為，一個(gè)立方體元的八個(gè)頂點(diǎn)的差值如果小于定義的閘值，則這個(gè)立方體元為一個(gè)等值體元。因此，通過對(duì)等值體元進(jìn)行提取分類和有效索引，可以減少插值次數(shù)和跳過透明體元。

如圖4所示，作者對(duì)地震數(shù)據(jù)體元道采用了二叉樹索引結(jié)構(gòu)。二叉樹的節(jié)點(diǎn)中包括了指向數(shù)據(jù)起始體元的指針p和等值體元個(gè)數(shù)size，以及標(biāo)識(shí)體元是否透明的flag。若一個(gè)立方體元道內(nèi)的所有體元都是等值的，那么二叉樹只有一個(gè)根節(jié)點(diǎn)。否則將體元道進(jìn)行上下等分，使二叉樹的子節(jié)點(diǎn)只有父節(jié)點(diǎn)的n/2個(gè)體元。通過遞歸分類之后，將體元道所有的體元都進(jìn)行分類。如果顏色和不透明度發(fā)生改變，只需要改變傳遞函數(shù)的映射規(guī)則，將二叉樹節(jié)點(diǎn)中的透明標(biāo)志改為不透明，然后在第一次繪制的過程中，改變透明體元的二叉樹索引節(jié)點(diǎn)中的透明標(biāo)志，不需要重新構(gòu)建二叉樹索引結(jié)構(gòu)(八叉樹編碼和游程編碼都要重新構(gòu)建索引結(jié)構(gòu))。二叉樹索引結(jié)構(gòu)中的非葉子節(jié)點(diǎn)的size值，代表當(dāng)前節(jié)點(diǎn)子孫節(jié)點(diǎn)的體元個(gè)數(shù)，葉子節(jié)點(diǎn)的size值，代表當(dāng)前體元(包括當(dāng)前體元)之后等值體元的個(gè)數(shù)。在繪制過程中，可以通過二分查找在二叉樹索引中快速查找到對(duì)應(yīng)體元的索引節(jié)點(diǎn)。

圖4 二叉樹索引結(jié)構(gòu)建立過程Fig． 4 Process of creating binary tree index

1.2 基于立方體元的Shear-warp體繪制

1.2.1 基于存儲(chǔ)方向的體掃描線

傳統(tǒng)的Shear-warp算法采用游程編碼(RLE)的數(shù)據(jù)結(jié)構(gòu)加速數(shù)據(jù)場(chǎng)的遍歷。RLE需要建立三份垂直于主觀察方向的編碼數(shù)據(jù)，每組編碼數(shù)據(jù)重采樣時(shí)，都以掃描線為遍歷的步長(zhǎng)單位。對(duì)于本來就規(guī)模龐大的三維地震數(shù)據(jù)，三份編碼數(shù)據(jù)增大內(nèi)存消耗的同時(shí)也提高了數(shù)據(jù)訪問的時(shí)間。

作者采用基于體元道存儲(chǔ)順序的體元掃描線方式(如下頁圖5所示)。該方式在任意的主觀察方向上，都采用和體元的存儲(chǔ)方向一致的體元掃描方式。按照存儲(chǔ)方向訪問數(shù)據(jù)，提高了高速緩沖存儲(chǔ)器的命中率，只存儲(chǔ)一份數(shù)據(jù)，也減少了內(nèi)存的消耗。在從前向后融合的過程中，把每個(gè)體元道中對(duì)采樣點(diǎn)的融合的中間結(jié)果記錄在一個(gè)二維數(shù)組中，該二維數(shù)組記錄了Shear－warp錯(cuò)切空間中二維平面像素融合的值。在對(duì)下一個(gè)體元道進(jìn)行計(jì)算的時(shí)候，再讀取該二維數(shù)組的值帶入融合計(jì)算。

圖5 基于存儲(chǔ)方向的體掃描線Fig． 5 The volume scanning line based on storage direction

1.2.2 繪制流程

體繪制的繪制流程就是對(duì)重采樣點(diǎn)進(jìn)行融合，形成最終圖片的過程，而減少對(duì)重采樣點(diǎn)值的計(jì)算量是對(duì)體繪制加速的主要環(huán)節(jié)。作者提出的算法通過對(duì)二叉樹索引的快速查詢，避免或減少了對(duì)重采樣點(diǎn)的計(jì)算。每個(gè)重采樣點(diǎn)的繪制流程可以分以下步驟:

(1)遍歷每一個(gè)重采樣點(diǎn)(x，y，z)，計(jì)算其所在的體元編號(hào)，并在二叉樹索引中進(jìn)行二分查找對(duì)應(yīng)體元索引節(jié)點(diǎn)N。

(2)判斷節(jié)點(diǎn)N的透明標(biāo)識(shí)flag，如果flag為真，則表示該體元全透明，返回步驟(1)處理下一個(gè)重采樣點(diǎn);否則，則進(jìn)入步驟(3)。

(3)判斷該體元是否為等值體體元，如果為等值體，則將該重采樣點(diǎn)的值用體元值替換，并通過采用Shear-warp算法采樣步長(zhǎng)單位等距的特點(diǎn)，快速處理下一個(gè)重采樣點(diǎn)而無需再次查找該采樣點(diǎn)的體元索引，轉(zhuǎn)到步驟(1)處理下一個(gè)采樣點(diǎn);如果該節(jié)點(diǎn)是非等值體，則進(jìn)入步驟(4)。

(4)對(duì)于非等值體節(jié)點(diǎn)采用二維線性插值獲得采樣點(diǎn)的值。

在繪制過程中遇到等值體時(shí)的時(shí)候，可以通過Shear－warp算法采樣步長(zhǎng)單位等距的特點(diǎn)，快速處理下一個(gè)重采樣點(diǎn)而無需再次查找該采樣點(diǎn)的體元索引。同時(shí)，當(dāng)Shear-warp錯(cuò)切空間中像素融合結(jié)果值的透明度大于設(shè)定閘值的時(shí)候，提前終止該像素點(diǎn)從前向后的融合過程。

2 實(shí)驗(yàn)結(jié)果與分析

我們使用C++和OpenGL實(shí)現(xiàn)了基于游程編碼(RLE)的Shear－warp算法和作者在文中所描述的算法，所用的機(jī)器配置為IntelDual-core2.5 GHz雙核，2G內(nèi)存，ATIHD3600顯卡，256MB顯存，windows7系統(tǒng)。實(shí)驗(yàn)數(shù)據(jù)為200×400×1100的地質(zhì)斷層數(shù)據(jù)，實(shí)驗(yàn)結(jié)果如下頁表1。圖像如圖5所示。表1中對(duì)等值體元定義的閘值大小分別對(duì)應(yīng)圖5中的繪制效果。

圖6 不同閘值在相同顏色和透明度下的繪制結(jié)果Fig． 6 The drawing result of different brake value in thesame color and transparency

從表1的結(jié)果可以看出，新算法的加速比受等值體元進(jìn)行定義的閘值大小影響。若閘值增大，等值體元增多，插值計(jì)算量減小，算法速度提高。但是通過圖6(見上頁)的繪制結(jié)果可以看出，閘值越大，被近似當(dāng)做透明和等值的重采樣點(diǎn)增多，使重采樣點(diǎn)誤差增大，從而降低了體繪制圖像的質(zhì)量。因此如何定義閘值的大小，要通過實(shí)際應(yīng)用中的需求調(diào)整，才能達(dá)到不影響質(zhì)量的情況下最大限度地提高體繪制速度。

表1 繪制時(shí)間比較Tab． 1 Time-consumes

3 結(jié)束語

作者在本文中根據(jù)三維地震數(shù)據(jù)自身的特點(diǎn)，用立方體的方式模擬了地震數(shù)據(jù)的空間關(guān)系。通過對(duì)相同地質(zhì)體的分類，該算法在繪制時(shí)只對(duì)周圍點(diǎn)的值都不相同的重采樣點(diǎn)進(jìn)行插值計(jì)算，忽略了透明體元和相同地質(zhì)體上重采樣點(diǎn)的插值計(jì)算。當(dāng)?shù)卣饠?shù)據(jù)中有大規(guī)模的相同地質(zhì)體時(shí)，本文中算法的加速效果將更加明顯。

通過對(duì)地震數(shù)據(jù)的繪制結(jié)果證明，本算法相對(duì)游程編碼的Shear-warp算法繪制速度有了顯著的提高。

［1］ 唐澤圣.三維數(shù)據(jù)場(chǎng)可視化［M］.北京:清華大學(xué)出版社，1999.

［2］ LEVOYM.Displayofsurfacesfromvolumedata［J］.IEEEComputerGraphics＆Applications，1988，8(3):29.

［3］ WESTOVERL.Footprintevaluationforvolumerendering［J］.ComputerGraphics，1990，24(4):367.

［4］ LACROUTEP，LEVOYM.FastVolumeRenderingUsing Shear－WarpFactorizationoftheViewingTransformation［J］.ComputerGraphicsproceedin － gs.July.1994:451.

［5］ 宋濤，歐宗瑛，王瑜，等.八叉樹編碼體數(shù)據(jù)的快速體繪制算法［J］.計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2005，9:1991.

［6］ MATSUIM，INOF，HAGIHARAK.Parallelvolumerenderingwithearlyrayterminationforvisualizinglargescaledatasets［A］.ISPA2004［C］.2004.

［7］ 許慶功，劉慶偉，張永勝.基于硬件加速的紋理映射體繪制［J］.計(jì)算機(jī)工程與應(yīng)用，2009(26):190.

［8］ 蘇超軾，趙明昌，張向文.GPU加速的八叉樹體繪制算法［J］.計(jì)算機(jī)應(yīng)用，2008，9:1232.

［9］ SORENG，STEFANB，ARMINK，etal.MemoryefficientaccelerationstructuresandtechniquesforCPU-basedvolumeraycastingoflargedata［A］.IEEESymposiumonVolumeVisualizationandGraphics2004［C］.Austin，Texas，USA:2004.

［10］ LACROUTEP.FastVolumeRenderingUsingaShear-WarpFactorizationoftheViewingTransformation［D］.(Ph.D.dissertation).America.StanfordUniversity，1995.

TP391.41

A

1001—1749(2011)05—0563—05

國(guó)家自然科學(xué)基金重點(diǎn)項(xiàng)目資助(40839905);教育部“新世紀(jì)優(yōu)秀人才支持計(jì)劃”資助(NCET-07-0148)

2011－07－04 改回日期:2011－07－14

周大鑫(1986－)，男，四川木里人，碩士，研究方向:計(jì)算機(jī)圖形。