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

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

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

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

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

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

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

可視化;體繪制;地震數據;二叉樹;Shear-warp

0 前言

直接體繪制是近年來迅速發展起來的，對三維數據場的一種體可視化繪制方法。該方法可以非常清晰地反應出三維數據場的細節，并且能夠將面繪制不能顯示的三維數據內部構造顯示出來。由于直接體繪制能夠對來自于地質勘探、核磁共振、計算機斷層掃描、計算機流體力學等的三維體數據進行有效的繪制，所以體繪制技術在石油勘探、生物醫學工程、地質科學、氣象、有限元分析后處理、化學等領域得到了廣泛地應用和重視［1］。近年來，體繪制技術又被應用于考古和文物保護、宇航、地震預測、自然現象仿真等領域。

三維地震是目前石油勘探的主要手段，由于三維地震數據的數據量很大(例如一個普通三維工區有1024×1024個道，每道有5000個點，每個點為一個浮點數存儲，其大小為20G，相對現有的普通微機內存是很難將其全部載入)，因此其體繪制是一個計算和數據雙重密集型問題。如何在不損失最終繪制結果質量的前提下，通過提高體繪制的成像速度來達到交互式體繪制，對于推動直接體繪制技術在石油勘探領域進一步應用具有重要的實際意義，并已成為國內、外石油工業界和學術界共同關心的前沿問題。

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

Shear-warp算法通過矩陣變換的方式，將繪制過程分解成三維數據場shear和二維圖像的warp兩個過程，以達到將三維數據場的重采樣和插值過程轉變為對二維平面的重采樣和插值過程，從而簡化了計算復雜度和計算量。

在Shear-warp算法實現過程中，可以采用多種加速方法。如文獻［5］中的八叉樹編碼體數據，實現對體數據的快速編碼，還有游程編碼(RLE)［10］數據結構加速的方式，都能跳過大規模“透明”的無效體素，從而達到加速算法的目的。但是，當體數據的透明度發生變化之后，八叉樹編碼和RLE都要重新對體數據進行編碼，需要額外的預處理時間，這樣對體繪制的交互性造成了影響。而且這兩種方法都只考慮了體數據中的透明無效體素，沒有考慮到在體數據中存在的相鄰數據值相同的情況。同時RLE方法需要在每個主方向存儲一套切片，這對于大規模的三維地震數據而言，增大了內存消耗，從而增加了數據讀取時間。

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

1 利用立方體元實現Shear-warp加速算法

Shear-warp地震數據體會制算法的整個流程，可以劃分為兩個階段:

(1)數據預處理階段。

與大型出版企業相比，中小出版企業在資源壟斷、出版許可、網點布局、文化沉淀、社會責任等方面差距較大。因此，中小出版企業轉型升級不能效仿大企業的“高、大、上、全”。

(2)繪制階段。

如圖1所示，預處理階段包括立方體元數據結構的壓縮和存儲，二叉樹索引結構的建立。繪制階段包括透明立方體元的跳躍，對不同體元的插值處理方式，以及顏色和不透明度的映射。

1.1 數據預處理

1.1.1 地震數據壓縮和道的立方體元轉換

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

f3D(x)=

圖1 繪制流程Fig． 1 Rendering process

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

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

圖2 三維地震數據剖面Fig． 2 3D seismic data profile

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

因為相同的地質體具有相同的體數據值，而相同的體數據在通過傳遞函數映射之后，必然得到相同的顏色值和透明度。Shear-warp算法對重采樣點的值進行計算的時候，采用的是線性插值算法，若C為AB連線上的一點，公式(2)為一次線性插值算法，其中 va和 vb表示點 a和 b的值，α=CB/AB。若 va=vb，那么 vc=va=vb。

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

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

圖4 二叉樹索引結構建立過程Fig． 4 Process of creating binary tree index

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

1.2.1 基于存儲方向的體掃描線

傳統的Shear-warp算法采用游程編碼(RLE)的數據結構加速數據場的遍歷。RLE需要建立三份垂直于主觀察方向的編碼數據，每組編碼數據重采樣時，都以掃描線為遍歷的步長單位。對于本來就規模龐大的三維地震數據，三份編碼數據增大內存消耗的同時也提高了數據訪問的時間。

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

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

1.2.2 繪制流程

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

(1)遍歷每一個重采樣點(x，y，z)，計算其所在的體元編號，并在二叉樹索引中進行二分查找對應體元索引節點N。

(2)判斷節點N的透明標識flag，如果flag為真，則表示該體元全透明，返回步驟(1)處理下一個重采樣點;否則，則進入步驟(3)。

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

(4)對于非等值體節點采用二維線性插值獲得采樣點的值。

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

2 實驗結果與分析

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

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

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

表1 繪制時間比較Tab． 1 Time-consumes

3 結束語

作者在本文中根據三維地震數據自身的特點，用立方體的方式模擬了地震數據的空間關系。通過對相同地質體的分類，該算法在繪制時只對周圍點的值都不相同的重采樣點進行插值計算，忽略了透明體元和相同地質體上重采樣點的插值計算。當地震數據中有大規模的相同地質體時，本文中算法的加速效果將更加明顯。

通過對地震數據的繪制結果證明，本算法相對游程編碼的Shear-warp算法繪制速度有了顯著的提高。

［1］ 唐澤圣.三維數據場可視化［M］.北京:清華大學出版社，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］ 宋濤，歐宗瑛，王瑜，等.八叉樹編碼體數據的快速體繪制算法［J］.計算機輔助設計與圖形學學報，2005，9:1991.

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

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

［8］ 蘇超軾，趙明昌，張向文.GPU加速的八叉樹體繪制算法［J］.計算機應用，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

國家自然科學基金重點項目資助(40839905);教育部“新世紀優秀人才支持計劃”資助(NCET-07-0148)

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

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