一種固體發動機光線投射體繪制加速方法

馬 強，朱 敏，盧洪義，于光輝，李 朋

（海軍航空工程學院 a.訓練部；b.七系；c.飛行器工程系；d.研究生管理大隊，山東 煙臺 264001）

0 引言

三維規則體數據可視化通常根據繪制過程中數據描述方法的不同而分為2類[1]：一類是表面繪制方法(Surface Rendering)；另一類是體繪制方法(Volume Rendering)。其中體繪制方法能不同程度表現物體中的細微結構和細小變化，可以同時將物體各組成結構的質地屬性、形狀特征及相互之間的層次關系表現出來。因此，在基于ICT 斷層圖像的固體發動機三維可視化故障診斷中用此方法實現三維重建。

目前，體繪制算法中繪制質量最高、使用最廣泛的是光線投射法[2]。但是，傳統的光線投射法重建速度比較慢，且一旦觀察方向發生變化，需要重新進行重采樣及圖像合成計算，效率低下。而固體發動機ICT 斷層圖像構造的三維規則體數據場數據量大，更顯現出繪制速度慢這一問題。

近些年來，國內外研究人員提出了許多光線投射法的改進算法，歸納起來主要分為2類：一類是利用圖像空間的相關性減少射線的數目，如趙建[3]、M.Levoy[4]并不從圖像平面所有像素投射光線，而是利用相鄰像素的相關性，有選擇地投射光線，但這種方法的顯示效果不是很好，不是當前研究的主要方向；另一類是利用對象空間的相關性減少不必要的采樣點數目，如洪歧[5]、R.Yagel[6]、J.Willtelms[7]、M.Agate[8]、K.R.Subramanian[9]等人采用分層的體數據結構重新規劃三維規則體數據場，跳過其中對繪制效果沒有影響的空體素，C.Suneup[10]對空體素賦一反映其與鄰近非空體素距離的值，在繪制時直接跳過空體素，這些方法不同程度的減少了計算量，在不降低圖像質量的前提下，提高了繪制速度。

上述改進方法都是針對空體素比例較高的醫學三維體數據場，如對試驗用的人體頭部CT數據進行了統計分析，發現其中73%的體素均為空體素[5]。而對于固體發動機，體數據場規模大，空體素比例有限，上述的加速算法效果不很明顯。為此，本文在作者申請的發明專利[11]中闡述的分割結果的基礎上，借鑒文獻[5-9]中的分層體數據結構加速方法，提出一種已分割規則體數據的光線投射加速方法。

1 光線投射法及計算量分析

光線投射法將三維規則體數據場按照一定的規則轉換為像平面上二維離散信號，即生成每個像素點顏色的R、G、B值，其實質是重新采樣和圖像合成過程[12]，基本原理如圖1所示。

圖1 光線投射法算法

圖1中，像平面的長和寬分別為W和H，采樣點間距相等。重新采樣過程中，首先需要將三維規則體數據場從物體空間坐標轉換為圖像空間坐標，然后從像平面的每一個像素根據用戶的觀察方向發出一條射線，這條射線穿過三維規則體數據場，沿著這條射線選擇N個等間距的采樣點，并由距離某一采樣點最近的8個數據點的顏色值和不透明度作三線性插值[13]，求出該采樣點的顏色值和不透明度。圖像合成過程中，將每條射線上各采樣點的顏色值及不透明度值用Over 算子[14]合成得到像平面上該像素點的最終顏色值。

重采樣過程中，在圖像坐標系中需要向三維規則體數據場發出W×H根射線，要運用矩陣變換式(1)計算每一個采樣點Pi在物體空間的坐標 'iP：

式中，M?1是一四階矩陣，每次變換需要16次乘法和12次加法操作。

另外，需要在物體空間中對每一采樣點作三線性插值，以得到Pi'點的灰度值，采用三線性插值需要7次乘法，14次加法。以本文實驗用到的固體發動機 ICT數據為例，三維規則體數據場大小為1024×1024×120，則一次光線投射繪制總的運算量為：(16+7)×1024×1024×120≈2.894×109次 乘法，(12+14)×1024×1024×120 ≈ 3.272×109次加法，可見光線投射法計算量相當大。

2 加速方法描述

本文加速方法的設計思想是：在分割預處理操作的基礎上，構造出包含分割結果信息的三維規則體數據場，再結合分層體數據結構加速方法，設計出一種通過用戶自定義三維重建的組成部分，而將不需重建的部分置為空體素的光線投射加速方法。

2.1 構造分割信息體數據場

由高能ICT檢測設備可以得到固體發動機的二維序列斷層圖像，把它們組織在一起構成三維規則體數據場，一般可假設其規模為L×M×N，如圖2所示。

圖2 三維規則體數據場

數據場中定義每個體素坐標為(i,j,k)，其中，i=0,1,???,L?1，j=0,1,???,M?1，k=0,1,???,N?1。此時每個體素記錄的值為經ICT檢測后生成的每個體素的灰度值。

對于分割處理后的三維規則體數據場，每個體素都有其歸屬類別，如文獻[11]中分割的固體發動機三維規則體數據為例，其分為背景、空氣環偽影、殼體、推進劑、星孔和內部缺陷6種類別。不失一般性，假設數據場中有K類組成部分。我們再構造一個體數據場，其空間坐標分布和圖2所示完全一致，只是每個體素記錄的值不是灰度值，而是該體素的類別，設為C (i,j,k)，C (i,j,k)=1,2,???,K。

2.2 構建分層體數據結構

2.2.1 分層體數據結構

假設L×M×N的三維規則體數據場中，N值最大，且N=2M+1(若N 不滿足，用空體素填補滿足)，M為一正整數。再將L×M×N的體數據場拓展為N×N×N的體數據場，多出部分用空體素填充。在此前提下，我們可以將N×N×N的體數據場分解為M+1個分層二值化體數據場，每一個分層二值化體數據場可用層編號m 檢索，其中 m=0,1,???,M，第m層體數據場表示為Vm。體數據場中，8個相鄰體素所圍成的立方體區域稱為體元。體數據場 V0為原始狀態，每條邊有N?1個體元，體數據場 V1每條邊有(N?1)/2個體元，以此類推，體數據場 VM中只有1個體元。圖3所示的M=2的分層體數據結構中，V0每條邊有4個體元，V1每條邊有2個體元，V2中僅有1個體元。

圖3 M=2的分層體數據結構

第m層體數據場中的體元用編號m和當前層中體素坐標 i=(i,j,k)檢索，表示為Vm(i)。第m層體數據場中體元的大小為原始體素間隔（體數據場 V0中體素間隔）的2m倍，體數據場 Vm中的體元要比體數據場Vm?1中的體元大8 倍，體數據場 Vm中體元 Vm(i)和其在體數據場Vm?1中的8個八叉元[15]子體元互為父子結構關系。圖3所示，體素坐標(0,0,0)定義為體數據場的右前下角，則體數據場 V0中體元 V0(0,0,0)為以體素(0,0,0)和(1,1,1)為對角頂點的立方體空間，體元 V0(0,0,0)是體元V1(0,0,0)的子體元。

2.2.2 基于分割信息的分層體數據結構構建

包含分割信息體數據場中，并不是所有的體數據都有意義，如固體發動機體數據場中，重建出背景和空氣環偽影就毫無用處。不失一般性，不妨設重建第a類和第b類組成部分，即取C (i)=a,C (i)=b，a,b=1,???,K，且a≠b，那么其他類別組成部分置為空體素。

我們按以下規則構造各層體數據場：對于初始體數據場 V0，如果其中體元 V0(i)8個頂點的體素類別不為a和b，則體元值為0，否則為1。除體數據場 V0外，如果第m層體數據場 Vm中體元 Vm(i)在體數據場Vm?1中的8個八叉元子體元均為0，則Vm(i)=0，否則Vm(i)=1。用公式定義可表示為下述2式。

第0層體數據場 V0中體元值為：

第m(m=1,???,M)層體數據場 Vm中體元值為：

式(2)中，{1，2，???，N?1}3表示由1，2，???，N?1任意組合構成的三維向量的集合。同理可知{0,1}3和所表示含義。

2.3 加速方法實現過程

結合上節所構建分層結構，將光線投射法作如下改進：光線投射從只包含1個體元的第M層體數據場VM開始，進入體元判斷體元的值，如果為0，穿越此體元不作采樣；如果為1，降低到下一層體數據場，繼續判斷體元的值。若穿越的相鄰體元不為同一個父體元，則跳到上一層體數據場，判斷體元的值，如此反復，遍歷整個三維規則體數據場。這樣，當體數據場中含空體素比例越大，節省的計算時間越多。當光線投射到第0層初始體數據場V0中，可知體元的8個體素頂點至少有1個點屬于a類或b類組成部分，在此體元中沿光線方向采樣。

更形象地，用圖4所示的流程框圖將上述過程表述如下。

圖4 加速方法實現流程框圖

2.4 加速方法實例

以文獻[11]中分割的固體發動機結果為例，在二維空間中解釋在分割基礎上的分層體數據結構構建和光線投射加速方法實現過程，見圖5。圖中，a)為固體發動機ICT 原始斷層圖像，可分割成背景、空氣環偽影、殼體、推進劑、星孔和缺陷6種組成部分，分別用1、2、3、4、5、6表示。取圖a)中部分區域，大小為16像素×16像素，用分割類別數字形式表現為圖b)。在圖b)基礎上繪制感興趣的固體發動機內部缺陷區域，則可將其余部分置為空體素，二維分層體數據結構光線投射過程見圖c)，圖中空白區域為置為0的空體素，深色區域為缺陷部分。

圖5 已分割固體發動機光線投射加速方法

3 實驗結果

采用C++和OpenGL 編程，計算機配置為酷睿2 雙核2.20GHzCPU、2GB 內存、GeForce 8400GS 256M 顯存顯卡。數據源為海軍無損檢測中心提供的固體發動機 ICT 斷層序列圖像，規模為1024×1024×120，經Z軸8 倍插值生成體數據場規模為1024×1024×960，用本文提出的光線投射加速方法進行重建，結果如圖6所示。

圖6 三維重建結果

圖6中：a)為在原始的未經分割的ICT 序列圖像上重建的結果，可見其外圍有一圈環狀噪聲；b)為經分割去除空氣環偽影的重建結果；c)為分割后選擇殼體、缺陷和星孔重建的結果；d)為分割后僅重建缺陷和星孔的結果。由a)到d)體數據場中空體素越來越多，本文提出的算法的加速效果也就越明顯，表1中給出傳統光線投射算法和本文算法的運算時間比較。

表1 兩種方法重建效率比較

4 結束語

固體發動機的三維可視化故障診斷中，感興趣的是內部缺陷的空間位置、大小等結構信息。由ICT斷層圖像構成的體數據場中，并不是所有體素都對缺陷的檢測有意義。本文提出的在分割基礎上的光線投射加速方法，其實質是將分割后不需重建的組成部分置為空體素，建立基于分割信息的分層體數據結構，動態調整采樣步長，減少采樣點。實驗表明，在固體發動機三維可視化故障診斷中，使用本文的加速方法，合理組合重建出固體發動機已分割的組成部分，在保證缺陷檢測準確度和精確度同時，體繪制速度得到明顯提升。

[1]沈海戈,柯有安.醫學體數據三維可視化方法的分類與評價[J].中國圖像圖形學報,2000,85(7)∶545-550.

[2]GRIMM S,BRUEKNER S,KANITSAR A.Memory efficient acceleration structures and techniques for CPU— based volume ray casting of large data[C]//Proceedings of IEEE Symposium on Volume Visualization and Graphics 2000.2004∶1-8.

[3]趙建.基于體繪制的醫學數據可視化方法的研究[D].大連∶大連理工大學,2006∶21-24.

[4]LEVOY M.Volume rendering by adaptive refinement[J].The Visual Computer,1990,6(1)∶2-7.

[5]洪歧,張樹生,劉雪梅.基于二值體金字塔數據結構的自適應快速體繪制方法[J].計算機工程與應用,2007,43(9)∶99-100.

[6]YAGEL R,SHI Z.Accelerating volume animation by space-leaping[C]//Proc.of IEEE Visualization '93.1993∶62-69.

[7]WILLTELMS J,GELDER A V.Octrees for faster isosurface generation[J].ACM Transactions on Graphics,1992,11(3)∶201-227.

[8]AGATE M,GRIMSDELE R L,LISTER P F.The HERO algorithm for ray tracing octrees[J].Advance in Computer Graphics Hardware,1997,12(3)∶61-73.

[9]SUBRAMANIAN K R,FUSSELL D S.Applying space subdivision techniques to volume rendering[C]//Proceedings of Visualization '90.San Francisco,1990∶150-159.

[10]SUNEUP C,HNGEONGDO K.Efficient space-leaping method for volume rendering[C]//SPIE Conference on Visual Data Exploration and Analysis.1999∶263-270.

[11]盧洪義,朱敏.一種固體發動機CT圖像的分割方法∶中國,200710301441.5[P].2007-12-28.

[12]唐澤圣.三維數據場可視化[M].北京∶清華大學出版社,1999∶8-15.

[13]WITTENBRINK C,MALZBENDER T,GOSS M.Opacity-weighted color interpolation for volume sampling[C]//Symposium on Volume Visualization.1998∶135-142.

[14]CLINE H E,LORENSEN W E,LUDKE S.Two algorithms for the 3d reconstruction of tomograms[J].Medical Physics,1988,15(3)∶320-327.

[15]LEVOY M.Efficient ray tracing of volume data[J].ACM Transactions on Graphics,1990,9(3)∶245-261.