馬 強(qiáng)，朱 敏，盧洪義，于光輝，李 朋

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

0 引言

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

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

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

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

1 光線投射法及計(jì)算量分析

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

圖1 光線投射法算法

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

重采樣過程中，在圖像坐標(biāo)系中需要向三維規(guī)則體數(shù)據(jù)場發(fā)出W×H根射線，要運(yùn)用矩陣變換式(1)計(jì)算每一個(gè)采樣點(diǎn)Pi在物體空間的坐標(biāo) 'iP：

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

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

2 加速方法描述

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

2.1 構(gòu)造分割信息體數(shù)據(jù)場

由高能ICT檢測設(shè)備可以得到固體發(fā)動(dòng)機(jī)的二維序列斷層圖像，把它們組織在一起構(gòu)成三維規(guī)則體數(shù)據(jù)場，一般可假設(shè)其規(guī)模為L×M×N，如圖2所示。

圖2 三維規(guī)則體數(shù)據(jù)場

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

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

2.2 構(gòu)建分層體數(shù)據(jù)結(jié)構(gòu)

2.2.1 分層體數(shù)據(jù)結(jié)構(gòu)

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

圖3 M=2的分層體數(shù)據(jù)結(jié)構(gòu)

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

2.2.2 基于分割信息的分層體數(shù)據(jù)結(jié)構(gòu)構(gòu)建

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

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

第0層體數(shù)據(jù)場 V0中體元值為：

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

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

2.3 加速方法實(shí)現(xiàn)過程

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

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

圖4 加速方法實(shí)現(xiàn)流程框圖

2.4 加速方法實(shí)例

以文獻(xiàn)[11]中分割的固體發(fā)動(dòng)機(jī)結(jié)果為例，在二維空間中解釋在分割基礎(chǔ)上的分層體數(shù)據(jù)結(jié)構(gòu)構(gòu)建和光線投射加速方法實(shí)現(xiàn)過程，見圖5。圖中，a)為固體發(fā)動(dòng)機(jī)ICT 原始斷層圖像，可分割成背景、空氣環(huán)偽影、殼體、推進(jìn)劑、星孔和缺陷6種組成部分，分別用1、2、3、4、5、6表示。取圖a)中部分區(qū)域，大小為16像素×16像素，用分割類別數(shù)字形式表現(xiàn)為圖b)。在圖b)基礎(chǔ)上繪制感興趣的固體發(fā)動(dòng)機(jī)內(nèi)部缺陷區(qū)域，則可將其余部分置為空體素，二維分層體數(shù)據(jù)結(jié)構(gòu)光線投射過程見圖c)，圖中空白區(qū)域?yàn)橹脼?的空體素，深色區(qū)域?yàn)槿毕莶糠帧?/p>

圖5 已分割固體發(fā)動(dòng)機(jī)光線投射加速方法

3 實(shí)驗(yàn)結(jié)果

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

圖6 三維重建結(jié)果

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

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

4 結(jié)束語

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

[1]沈海戈,柯有安.醫(yī)學(xué)體數(shù)據(jù)三維可視化方法的分類與評價(jià)[J].中國圖像圖形學(xué)報(bào),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]趙建.基于體繪制的醫(yī)學(xué)數(shù)據(jù)可視化方法的研究[D].大連∶大連理工大學(xué),2006∶21-24.

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

[5]洪歧,張樹生,劉雪梅.基于二值體金字塔數(shù)據(jù)結(jié)構(gòu)的自適應(yīng)快速體繪制方法[J].計(jì)算機(jī)工程與應(yīng)用,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]盧洪義,朱敏.一種固體發(fā)動(dòng)機(jī)CT圖像的分割方法∶中國,200710301441.5[P].2007-12-28.

[12]唐澤圣.三維數(shù)據(jù)場可視化[M].北京∶清華大學(xué)出版社,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.