3D-Wedgelet分解與工業(yè)CT體數(shù)據(jù)面特征提取

曾 理 李宗劍 劉長江

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

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

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

1 引言

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

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

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

多尺度幾何分析是近十多年來，圖像稀疏分解發(fā)展的一個分支，其非常適合檢測、表示、處理某些類型的高維數(shù)據(jù)。Wedgelet變換作為多尺度幾何分析方法家族中的一員[7]，由斯坦福大學的Donoho于1999年提出。Wedgelet以區(qū)域基的方式分析數(shù)據(jù)，對以不同灰度區(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的邊界點集為V(如圖2黑點所示，其中，1小立方體塊代表1體素)。任取不在同一邊上的3點(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 立方體塊中的邊界點

2.2 快速3D-Wedgelet分解

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

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

圖3 多邊形Planelet

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

圖4 坐標系

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

圖5 三邊形Planelet

設f(x, y, z)=0是圖5(a)的一個三邊形，則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)的一個三邊形。圖5(a)中，三邊形的3頂點坐標分別為(x1,0,0)，(0,y2,0)，(0,0,z3)。

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

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

對于由相鄰4個面形成的四邊形Planelet，設f(x, y, z)=0是圖6(d)的一個四邊形，則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)的一個四邊形。圖6(d)中，頂點A，B，C的坐標分別為(0,0,z1)，(x2,0,N)，(0,N, z3)。

(3)五邊形Planelet 五邊形Planelet與立方體塊的5個面相交，換句話說，只與1個面不相交，據(jù)此可粗分為6種情形。而對于每一種情形，又可細分為4種。設f(x, y, z)=0是圖7的一個五邊形，則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分別是立方體中的一個不重復五邊形。圖7中，頂點A，B，C的坐標分別為(0,N,z1)，(x2,0,N)，(x3,0,0)。且x2≠x3，因為，若x2=x3，則為四邊形Planelet。

圖6 四邊形Planelet

(4)六邊形Planelet 六邊形Planelet與立方體塊的6個面均相交，且3組對邊兩兩平行，可以將其分為4種情形。設f(x, y, z)=0是圖8的一個六邊形，則f(N?y, x, z)=0，f(N?x, N?y, z)=0，f(y, N?x, z)=0分別是立方體中的一個不重復六邊形。圖8中，頂點A，B，C的坐標分別為(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，包含的體素個數(shù)分別為N1, N2，N1+ N2=N3。

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

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

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

重復上述步驟，依次計算完立方體內的三邊形、四邊形、五邊形和六邊形Planelet。

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

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

3.1 單尺度下的特征提取

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

3.2 多尺度下的特征提取

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

由文獻[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é)單尺度下的特征提取思想，給出多尺度下特征提取的判斷條件。

對體數(shù)據(jù)而言，當

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

對切片圖像而言，當

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

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

3.3 線特征的融合

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

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

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

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

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

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

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

4 實驗結果與分析

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

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

圖10 實驗結果

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

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

從計算時間上看，基于3D-Wedgelet方法的計算時間大致是基于Wedgelet方法的4.667倍。

5 結束語

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

[1] 曾理, 安貝貝, 馬睿. 脊波在工業(yè)CT圖像裂紋邊緣檢測中的應用. 儀器儀表學報, 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] 焦李成, 譚山. 圖像的多尺度幾何分析: 回顧和展望. 電子學報, 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] 張丹丹. 基于邊緣算子的結構損傷診斷方法研究. [碩士論文],湖南大學, 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] 馬睿, 曾理, 盧艷平. 改進的基于Facet模型的亞像素邊緣檢測. 應用基礎與工程科學學報, 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] 胡正磊, 孫進平, 袁運能等. 基于小波邊緣提取和脊線跟蹤技術的SAR圖像河流檢測算法. 電子與信息學報, 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] 中國科學院自動化研究所醫(yī)學影像研究室. 3DMed. http://www.3dmed.net.Institute of Automation Chinese Academy of Science, Medical Image Processing Group. 3DMed. http://www. 3dmed.net.