曾 理 李宗劍 劉長(zhǎng)江

①(重慶大學(xué)光電技術(shù)及系統(tǒng)教育部重點(diǎn)實(shí)驗(yàn)室ICT研究中心 重慶 400044)②(重慶大學(xué)數(shù)理學(xué)院 重慶 400044)

③(重慶大學(xué)光電工程學(xué)院 重慶 400044)

④(四川理工學(xué)院數(shù)學(xué)系 自貢 643000)

1 引言

工業(yè)CT利用射線源在無損狀態(tài)下對(duì)被檢測(cè)材料進(jìn)行射線掃描檢測(cè)和圖像重建，具有直觀、可靠、靈敏度和分辨率高、影像不重疊等優(yōu)點(diǎn)，被國際無損檢測(cè)界譽(yù)為最佳的無損檢測(cè)手段[1]。

工業(yè)CT體數(shù)據(jù)的主要特征由點(diǎn)、線、面、管狀物和絲狀物刻畫[2]。其中，面特征常常對(duì)應(yīng)目標(biāo)對(duì)象的輪廓或邊界面。抽取出工業(yè)CT體數(shù)據(jù)中的面特征對(duì)于描述和解譯工業(yè)CT體數(shù)據(jù)本身是很重要的。

對(duì)體數(shù)據(jù)面特征的提取策略包括兩大類：一類是直接對(duì)體數(shù)據(jù)進(jìn)行面特征的提取；另一類是將體數(shù)據(jù)分解成一系列2D圖像，對(duì)每幅2D圖像進(jìn)行線特征的提取，再將各幅2D圖像的線特征組合成面特征。最經(jīng)典的處理方法是邊緣差分算子[3]，其它常見的方法如主動(dòng)輪廓模型及其改進(jìn)模型[4]、Facet模型[5]、基于小波的方法[6]等等。上述方法提取得到的面特征，其結(jié)果并不是直接描述面特征的方程，而是由一個(gè)個(gè)點(diǎn)特征組成，因此，只能從主觀視覺效果上說得到了面特征。

多尺度幾何分析是近十多年來，圖像稀疏分解發(fā)展的一個(gè)分支，其非常適合檢測(cè)、表示、處理某些類型的高維數(shù)據(jù)。Wedgelet變換作為多尺度幾何分析方法家族中的一員[7]，由斯坦福大學(xué)的Donoho于1999年提出。Wedgelet以區(qū)域基的方式分析數(shù)據(jù)，對(duì)以不同灰度區(qū)域形式存在的邊緣線和邊緣面有很好的分析能力[8,9]。

2 離散3D-Wedgelet變換

Wedgelet被定義為N×N(N=2J)圖像中被一條Beamlet分成的兩塊區(qū)域[10]，同理，3D-Wedgelet被定義為N×N×N(N=2J)數(shù)據(jù)中被一張Planelet分成的兩塊區(qū)域[11]，如圖1所示。

2.1 3D-Wedgelet分解

記I為一N×N×N(N=2J)的體數(shù)據(jù)，Cj,k=表示尺度為j，位置為k=(k1, k2, k3)的立方體塊，且0≤j≤J，0≤k1, k2, k3＜2j，j, k1, k2, k3∈Z。記Cj,k的邊界點(diǎn)集為V(如圖2黑點(diǎn)所示，其中，1小立方體塊代表1體素)。任取不在同一邊上的3點(diǎn)(v1, v2, v3)∈V，則平面將Cj,k劃分成兩塊區(qū)域，兩塊區(qū)域的體素均值分別記為ma和mb，則立方體塊Cj,k的3D-Wedgelet分解集W[ I(Cj,k)]為

圖1 Wedgelet(a)和3D-Wedgelet(b)

圖2 立方體塊中的邊界點(diǎn)

2.2 快速3D-Wedgelet分解

對(duì)于N×N×N(N=2J)的立方體塊，其邊界點(diǎn)的個(gè)數(shù)為12N?4，理論上說，有張Planelet。但其中存在重復(fù)的情形，如圖1(b)所示的Planelet包含了4個(gè)邊界點(diǎn)，若按每3個(gè)邊界點(diǎn)計(jì)算一張Planelet，則該P(yáng)lanelet將被計(jì)算4次。因此，若對(duì)立方體塊中的每張Planelet直接進(jìn)行計(jì)算，將存在大量的重復(fù)計(jì)算。鑒于此，考慮通過找出單尺度下不同Planelet間的內(nèi)在聯(lián)系，以實(shí)現(xiàn)減少3D-Wedgelet分解的計(jì)算量。

2.2.1 Planelet的形狀 根據(jù)Planelet的邊數(shù)，可將Planelet分為三邊形、四邊形、五邊形和六邊形4種，如圖3所示。

圖3 多邊形Planelet

2.2.2 單尺度下Planelet間的聯(lián)系 為討論方便，建立如圖4的坐標(biāo)系，原點(diǎn)為立方體塊的前左下頂點(diǎn)。

圖4 坐標(biāo)系

(1)三邊形Planelet 三邊形Planelet的3頂點(diǎn)相對(duì)位置關(guān)系只可能為圖5所示的8種情形。

圖5 三邊形Planelet

設(shè)f(x, y, z)=0是圖5(a)的一個(gè)三邊形，則f(x, y, N?z)=0，f(N?x, y, z)=0，f(N?x, y, N?z)=0，f(x, N?y, z)=0，f(x, N?y, N?z)=0，f(N?x, N?y, z)=0，f(N?x, N?y, N?z)=0分別是圖5(b)-5(h)的一個(gè)三邊形。圖5(a)中，三邊形的3頂點(diǎn)坐標(biāo)分別為(x1,0,0)，(0,y2,0)，(0,0,z3)。

(2)四邊形Planelet 四邊形Planelet包含兩大類情況：一類是兩組對(duì)邊均平行；另一類是只有一組對(duì)邊平行。構(gòu)成四邊形Planelet的立方體面可以是兩組對(duì)面，也可以是相鄰4個(gè)面，如圖6所示，其中圖6(a)-6(c)由兩組對(duì)面形成，圖6(d)-6(o)由相鄰4個(gè)面形成。

對(duì)于由兩組對(duì)面形成的四邊形Planelet，設(shè)f(x, y, z)=0是圖6(a)的一個(gè)四邊形，則f(N?z,y, x)=0、f(x, N?z, y)=0分別是圖6(b)，6(c)的一個(gè)四邊形。圖6(a)中，頂點(diǎn)A、B、C的坐標(biāo)分別為(0,0,z1)，(N,0,z2)，(0,N, z3)。

對(duì)于由相鄰4個(gè)面形成的四邊形Planelet，設(shè)f(x, y, z)=0是圖6(d)的一個(gè)四邊形，則f(N?y,x, z)=0，f(N?x, N?y, z)=0，f(y, N?x, z)=0，f(x, N?z, y)=0，f(z, x, y)=0，f(N?x, z, y)=0，f(N?z, N?x, y)=0，f(x, N?y, N?z)=0，f(y, x,N?z)=0，f(N?x, y, N?z)=0，f(N?y, N?x,N?z)=0分別是圖6(d)-6(o)的一個(gè)四邊形。圖6(d)中，頂點(diǎn)A，B，C的坐標(biāo)分別為(0,0,z1)，(x2,0,N)，(0,N, z3)。

(3)五邊形Planelet 五邊形Planelet與立方體塊的5個(gè)面相交，換句話說，只與1個(gè)面不相交，據(jù)此可粗分為6種情形。而對(duì)于每一種情形，又可細(xì)分為4種。設(shè)f(x, y, z)=0是圖7的一個(gè)五邊形，則f(N?x, y, z)=0，f(z, y, N?x)=0，f(z, y, x)=0，f(x, N?y, z)=0，f(N?x, N?y, z)=0，f(z,N?y, N?x)=0，f(z, N?y, x)=0，f(y, x, z)=0，f(y, N?x, z)=0，f(y, N?z, N?x)=0，f(y, N?z, x)=0，f(N?y, x, z)=0，f(N?y, N?x, z)=0，f(N?y, N?z, N?x)=0，f(N?y, N?z, x)=0，f(x, z, N?y)=0，f(N?x, z, N?y)=0，f(N?z, x, N?y)=0，f(N?z, N?x, N?y)=0，f(x, z, y)=0，f(N?x, z, y)=0，f(N?z, x, y)=0，f(N?z, N?x, y)=0分別是立方體中的一個(gè)不重復(fù)五邊形。圖7中，頂點(diǎn)A，B，C的坐標(biāo)分別為(0,N,z1)，(x2,0,N)，(x3,0,0)。且x2≠x3，因?yàn)椋魓2=x3，則為四邊形Planelet。

圖6 四邊形Planelet

(4)六邊形Planelet 六邊形Planelet與立方體塊的6個(gè)面均相交，且3組對(duì)邊兩兩平行，可以將其分為4種情形。設(shè)f(x, y, z)=0是圖8的一個(gè)六邊形，則f(N?y, x, z)=0，f(N?x, N?y, z)=0，f(y, N?x, z)=0分別是立方體中的一個(gè)不重復(fù)六邊形。圖8中，頂點(diǎn)A，B，C的坐標(biāo)分別為(0,y1, N)，(x2,0,0)，(N, N, z3)。

圖7 五邊形Planelet代表

圖8 六邊形Planelet代表

2.2.3 單尺度下3D-Wedgelet的分解 記I0=被一張Planelet分開的兩塊3D-Wedgelet區(qū)域分別為w1, w2，各自的體素均值分別為m1, m2，包含的體素個(gè)數(shù)分別為N1, N2，N1+ N2=N3。

第1步 使用2.2.2節(jié)的方法，選定頂點(diǎn)A，B，C，計(jì)算平面Planelet的方程f(x, y, z)=0。

第2步 計(jì)算被平面f(x, y, z)=0分成的兩塊3D-We d g e l e t區(qū)域的體素均值m1, m2，m1=/N2。

第3步 利用2.2.2節(jié)討論的坐標(biāo)變換得到與該3D-Wedgelet區(qū)域相關(guān)的其它3D-Wedgelet區(qū)域，并計(jì)算出各自的體素均值。

重復(fù)上述步驟，依次計(jì)算完立方體內(nèi)的三邊形、四邊形、五邊形和六邊形Planelet。

3 工業(yè)CT體數(shù)據(jù)面特征提取

對(duì)于立方體中的兩互補(bǔ)3D-Wedgelet，按圖4所示的x，y，z方向分別進(jìn)行切片劃分，均會(huì)得到一序列組的Wedgelet。因此，對(duì)于工業(yè)CT體數(shù)據(jù)面特征的提取，給出如下的兩種思路：(1)基于3D-Wedgelet分解的提取方法；(2)先對(duì)體數(shù)據(jù)在x，y，z方向分別進(jìn)行切片劃分，然后對(duì)每張切片作基于Wedgelet分解的線特征提取，最后融合3個(gè)方向的線特征從而得到體數(shù)據(jù)的面特征。

3.1 單尺度下的特征提取

對(duì)于體數(shù)據(jù)而言，根據(jù)所選尺度j和位置k=(k1, k2, k3)，求出立方體塊Cj,k內(nèi)滿足條件(t為閾值)的其中，m1, m2為由Planelet分成的兩互補(bǔ)3DWedgelet的體素均值。對(duì)于切片圖像而言，處理的思想一樣。

3.2 多尺度下的特征提取

多尺度與單尺度相比，具有更大的靈活性。多尺度Wedgelet或3D-Wedgelet分解由不同尺度下的單尺度Wedgelet或3D-Wedgelet分解組成。經(jīng)多尺度Wedgelet或3D-Wedgelet分解后的數(shù)據(jù)形成一種四叉樹或八叉樹結(jié)構(gòu)。對(duì)四叉樹或八叉樹進(jìn)行由下至上地最優(yōu)修剪，就能得到線特征或面特征的最優(yōu)表示。

由文獻(xiàn)[11]知道：大尺度下的1條Beamlet，由小尺度下的1，2或3條Beamlet組成；大尺度下的1張Planelet，由小尺度下的1，2，3，4，5或6張Planelet組成。記大尺度下的正方形塊為S，立方體塊為C，小尺度下的正方形塊分別為S1～S4，立方體塊分別為C1～C8。沿用3.1節(jié)單尺度下的特征提取思想，給出多尺度下特征提取的判斷條件。

對(duì)體數(shù)據(jù)而言，當(dāng)

同時(shí)成立時(shí)，將C分解為C1～C8，其它情況則保留C。

對(duì)切片圖像而言，當(dāng)

同時(shí)成立時(shí)，將S分解為S1～S4，其它情況則保留S。

條件式(2)，式(4)說明子立方體塊或圖像塊中存在比父立方體塊或圖像塊中更明顯的面或線特征。條件式(3)，式(5)保證了不會(huì)將1張大的面特征劃分為多張小的面特征或1條長(zhǎng)的線特征劃分為多條短的線特征。

3.3 線特征的融合

考慮圖9(a)中的體數(shù)據(jù)模型，背景為白色，目標(biāo)對(duì)象為一實(shí)心立方體(為了視覺效果，用黑色標(biāo)記出12條邊)。記上，下，左，右，前，后表面分別為平面1，2，3，4，5，6。分別按x，y，z方向劃分切片，切片圖像均如圖9(b)(含目標(biāo))和9(c)(不含目標(biāo))所示。對(duì)于這個(gè)體數(shù)據(jù)模型，我們感興趣的面特征，實(shí)際上為立方體的6個(gè)面。

圖9 體數(shù)據(jù)模型

若沿x，y，z方向分別劃分切片，并進(jìn)行基于Wedgelet的線特征提取，則提取結(jié)果分別是平面1，2，5，6；平面1，2，3，4；平面3，4，5，6。

從上面的分析可以看出，若按一個(gè)方向?qū)w數(shù)據(jù)劃分切片，最終得到的面特征提取結(jié)果將比實(shí)際感興趣的面特征要少。因此，本文考慮將3個(gè)方向的提取結(jié)果進(jìn)行數(shù)據(jù)融合。

記x，y，z方向各自的線特征提取結(jié)果分別為out_x，out_y，out_z。從上面的分析可以知道，若體素點(diǎn)v位于感興趣的面特征上，則滿足下式

即將3個(gè)方向上各自切片序列組的線特征提取結(jié)果進(jìn)行合并，從而得到整個(gè)體數(shù)據(jù)的面特征。

但式(6)只是體素點(diǎn)位于感興趣面特征上的必要條件，不是充分條件。對(duì)于實(shí)際的體數(shù)據(jù)而言，由于其結(jié)構(gòu)組成更加復(fù)雜，利用式(6)得到的面特征中除了真正的面特征外還有一部分線特征，保留這部分線特征并不影響后續(xù)對(duì)體數(shù)據(jù)的分析。

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

實(shí)驗(yàn)的實(shí)際數(shù)據(jù)采用Internet上公布的工業(yè)CT體數(shù)據(jù)engine[12]，尺寸為256×256×128，其一張切片及3維顯示分別如圖10(a)，10(b)所示。實(shí)驗(yàn)用計(jì)算機(jī)配置如下：AMD3100+，1 G內(nèi)存，VC++6.0，Windows XP。基于3D-Wedgelet和Wedgelet方法的分解尺度均為尺度4到尺度6，面特征提取結(jié)果的3維顯示分別如圖10(c)、10(d)所示，計(jì)算耗時(shí)分別為5507.320 s和1180.308 s。3維顯示由3DMed軟件[13]實(shí)現(xiàn)。

從圖10(c)、10(d)中可以看出，本文的兩種策略都是有效的。這是因?yàn)樵诠I(yè)CT體數(shù)據(jù)中，大部分的線特征和面特征并不是以獨(dú)立的形式存在，而是表現(xiàn)為不同灰度區(qū)域的邊界線或邊界面。而Wedgelet或3D-Wedgelet基的表現(xiàn)形式為線(Beamlet)或面(Planelet)分開的兩塊區(qū)域，根據(jù)這兩塊區(qū)域的像素或體素均值，就能判斷出該線或面是否為線特征或面特征。

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

從提取結(jié)果的主觀效果上看，基于Wedgelet比基于3D-Wedgelet方法提取的面特征更平滑，尤其體現(xiàn)在曲面部分。出現(xiàn)這種結(jié)果有兩方面的原因。一方面，基于3D-Wedgelet方法得到的面特征是由一系列的Planelet組成，而Planelet本身是平面，對(duì)體數(shù)據(jù)中曲面特征的描述實(shí)際上是一個(gè)平面逼近曲面的過程。另一方面，基于Wedgelet的方法是將面特征轉(zhuǎn)化為3個(gè)相互垂直方向上的線特征進(jìn)行提取。在轉(zhuǎn)化的過程中同一張曲面經(jīng)不同方向的劃分后，有可能出現(xiàn)曲線，有可能出現(xiàn)直線。Beamlet本身是直線，對(duì)切片圖像中曲線特征的描述實(shí)際上是一個(gè)直線逼近曲線的過程。但是融合3個(gè)相互垂直方向的線特征之后，會(huì)讓這種逼近的痕跡減弱甚至消失。

從描述面特征的簡(jiǎn)潔性上看，基于3DWedgelet的方法得到的面特征是由一系列的平面(Planelet)方程組成，而基于Wedgelet的方法得到的面特征是由一系列的直線(Beamlet)方程組成，因此，前者對(duì)體數(shù)據(jù)面特征的描述會(huì)更簡(jiǎn)潔。

從計(jì)算時(shí)間上看，基于3D-Wedgelet方法的計(jì)算時(shí)間大致是基于Wedgelet方法的4.667倍。

5 結(jié)束語

本文先分析了單尺度下Planelet的組成及內(nèi)在聯(lián)系，得到了一種快速3D-Wedgelet分解方法。以此為基礎(chǔ)，提出了分別基于3D-Wedgelet，Wedgelet的兩種提取工業(yè)CT體數(shù)據(jù)中面特征的方法。實(shí)驗(yàn)結(jié)果表明，這兩種方法都是有效的。如何更好地描述工業(yè)CT體數(shù)據(jù)中的曲面特征并減少計(jì)算時(shí)間，是下一步研究的重點(diǎn)。

[1] 曾理, 安貝貝, 馬睿. 脊波在工業(yè)CT圖像裂紋邊緣檢測(cè)中的應(yīng)用. 儀器儀表學(xué)報(bào), 2007, 28(6): 981-986.Zeng Li, An Bei-bei, and Ma Rui. Application of ridgelet in crack edge detection of ICT image. Chinese Journal of Scientific Instrument, 2007, 28(6): 981-986.

[2] 焦李成, 譚山. 圖像的多尺度幾何分析: 回顧和展望. 電子學(xué)報(bào), 2003, 31(12A): 1975-1981.Jiao Li-cheng and Tan Shan. Development and prospect of image multiscale geometric analysis. Acta Electronica Sinica,2003, 31(12A): 1975-1981.

[3] 張丹丹. 基于邊緣算子的結(jié)構(gòu)損傷診斷方法研究. [碩士論文],湖南大學(xué), 2008.Zhang Dan-dan. The studies on structural damage identification method based on edge operator. [MA.dissertation], Hunan University, 2008.

[4] Vasilache S and Najarian K. A unified method based on wavelet filtering and Active Contour Models for segmentation of Pelvic CT images. International Conference on Complex Medical Engineering, Tempe, Arizona, USA, April 9-11, 2009:1-5.

[5] 馬睿, 曾理, 盧艷平. 改進(jìn)的基于Facet模型的亞像素邊緣檢測(cè). 應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報(bào), 2009, 17(2): 296-302.Ma Rui, Zeng Li, and Lu Yan-ping. Improved sub-pixel edge detection method based on Facet model. Journal of Basic Science and Engineering, 2009, 17(2): 296-302.

[6] 胡正磊, 孫進(jìn)平, 袁運(yùn)能等. 基于小波邊緣提取和脊線跟蹤技術(shù)的SAR圖像河流檢測(cè)算法. 電子與信息學(xué)報(bào), 2007, 29(3):524-527.Hu Zheng-lei, Sun Jin-ping, and Yuan Yun-neng, et al.. River detection in SAR images based on edge extraction in wavelet domain and ridge tracing technique. Journal of Electronics &Information Technology, 2007, 29(3): 524-527.

[7] Donoho D L. Wedgelet: nearly-minimax estimation of edges.Annals of Statistics, 1999, 27(3): 859-897.

[8] Liu Ze-yi, Sun Zi-qiang, and Xu Ling, et al.. Contour representation based on wedgelet. Journal of Systems Engineering and Electronics, 2006, 17(2): 251-257.

[9] Ren Chao, Wu Si-liang, and Jiao Li-cheng. Edge detection algorithm of SAR images with wedgelet filter. Journal of Beijing Institute of Technology, 2008, 17(3): 346-350.

[10] Donoho D L and Huo Xiao-ming. Beamlets and multiscale image analysis. Lecture Notes in Computational Science and Engineering, Berlin, Germany: Springer Press, 2001, 20:149-196.

[11] Ofer L. Multiscale geometric analysis of three-dimensional data. [Ph.D. dissertation], Scientific Computing and Computational Mathematics of Stanford University, 2005.

[12] General Electric. Engine. http://www.gris.uni-tuebingen.de/areas/scivis/volren/datasets/datasets.html.

[13] 中國科學(xué)院自動(dòng)化研究所醫(yī)學(xué)影像研究室. 3DMed. http://www.3dmed.net.Institute of Automation Chinese Academy of Science, Medical Image Processing Group. 3DMed. http://www. 3dmed.net.