工業CT三維圖像的表面平滑方法

王佳熙

(成都大學 計算機學院，四川 成都 610106)

0 引 言

工業CT(computed tomography，CT)技術能夠在不破壞工業部件的前提下得到具有復雜結構工業部件的三維圖像.但受檢測環境或設備自身等因素的影響，其獲得的測量數據在采集和傳輸過程中會受到噪聲的干擾，進而影響到經圖像重建后得到的工業部件的三維圖像質量.與一般的檢測方式不同，工業CT技術獲得的并不是工業部件表面的點云數據，而是被掃描工業部件的三維圖像灰度體數據，其需要通過圖像分割或其他的方法才能得到工業部件的三維圖像表面.但是，這些方法都會帶來一些偽影，比如圖像表面出現毛刺或部分不連續的情況，這會使三維圖像表面的網格模型的三角面片數目過多，進而影響后續的逆向設計和改造等圖像后處理過程.針對以上問題，科研人員進行了相關研究并取得了一定進展[1-6].針對噪聲問題，劉玲慧等[1]提出了三維魯棒Chan-Vese(3-dimension robust chan-Vese,3D-RCV)算法，其通過在3D-CV算法的基礎上引入局部灰度信息來解決工業CT三維圖像中存在的噪聲問題.其研究表明，3D-RCV算法能較好地獲得用水平集函數表示的工業部件的三維圖像表面.而針對圖像的平滑方法，Buades等[5]提出了非局部均值算法，其把圖像中出現的很多周期性冗余信息加以利用，在全局中搜索灰度相似塊并充分利用像素之間的相關性，從而更好地保護圖像的邊緣特征.董斌等[6]提出了基于水平集的非局部表面恢復算法，其利用水平集函數來表示三維圖像表面，從而使得三維圖像的表面在邊界得到保護的同時減少了大量的噪聲，而且在表面有毛刺和部分不連續的地方進行了重造，圖像整體變得更為平滑.在此基礎上，本研究針對工業CT三維圖像所含噪聲過多而導致得到的工業部件三維圖像表面有毛刺和部分不連續的缺點，將3D-RCV算法和基于水平集的非局部恢復算法結合用于工業部件三維圖像的表面平滑.同時，通過網上公開的工業部件的灰度體數據測試了算法的可行性與有效性.

1 工業CT三維圖像表面平滑方法

本研究探討的工業CT三維圖像的表面平滑方法大致可分為兩個步驟：其一，利用3D-RCV算法對工業部件的三維圖像灰度體數據進行三維分割，獲得用水平集函數φ表示的工業部件的三維圖像表面；其二，采用基于水平集的非局部表面恢復算法對獲得的結果進行三維圖像表面平滑處理.

1.1 3D-RCV算法

3D-RCV算法的主要思路是：對于圖像中的每一個像素點，它與其鄰域內的所有像素點的灰度信息都具有一定的相關性.例如，圖像中某個像素點x，在它的鄰域內一定存在沒有被噪聲污染的像素點.即使有些點被噪聲污染，仍然可以利用它的鄰域內沒被噪聲污染的像素點的信息來幫助曲面進行相對正確的演化，這樣就不會被噪聲點完全影響而導致錯誤的演化.

1.1.1 模型建立

設I：Ω→R是圖像域Ω?R3上的一幅灰度圖像，Ω1和Ω2是Ω被一個二維曲面C劃分成的兩個互不相交的同質子區域，即，Ω=Ω1∪C∪Ω2，Ω1∩Ω2=φ.給定圖像中的一個點x(x,y,z)∈Ωi,i=1,2，定義該點的局部能量如下,

(1)

式中，常數ci是子區域Ωi的平均灰度，I(y)是點x的鄰域N(x)內的像素點y(x,y,z)的灰度值，鄰域N(x)的大小通過核函數Kσ來控制.這里，選取截斷的高斯核函數作為Kσ，定義如下，

(2)

式中，σ是高斯函數的標準差，a是正則化常數，鄰域N(x)的半徑記為r.

x,y∈Ω

(3)

式中，λ1>0,λ2>0,為固定的參數.

下面引入面積正則項，則該算法的總體能量定義為，

F(C,c1,c2)

c2|2dy)dx,x,y∈Ω

(4)

在水平集方法中，演化曲面C?Ω可以用Lipschitz函數φ：Ω→R(也稱為水平集函數)的零水平集來表示，具體形式如下,

(5)

按文獻[7]的變分水平集表示方法，分別定義正則化Heaviside函數Hε(z)和Dirac函數δε(z)如下，

(6)

(7)

為了保持水平集函數的光滑性，使用文獻[8]定義的水平集函數正則項P，

(8)

1.1.2 模型求解

綜上，總體能量F(C,c1,c2)可以被改寫成如下形式，

Fε(φ,c1,c2)

|I(y)-c2|2dy)(1-Hε(φ(x)))dx,

x,y∈Ω

(9)

式中，μ≥0,ν≥0,分別是面積項和水平集正則項的參數.

該方程可以采用變量交替迭代和梯度下降的方法進行求解.由此，使得總體能量最小化的各個變量的解析表達式如下：

首先，固定ci,i=1,2，Fε(φ,c1,c2)對水平集函數φ的導數可以轉化為如下梯度下降流方程，

ν·div(dp(|▽φ|)▽φ)

(10)

然后，固定水平集函數φ，令泛函F(φ,c)對ci求導，可得到使函數F(φ,c)最小時的ci解析式為,

(11)

式中，Mi(φ(x))定義為，

(12)

1.2 基于水平集的非局部表面恢復算法

1.2.1 模型建立

近年來，圖像平滑算法已延拓到了三維圖像的表面平滑，其主要采取兩種方式來表示三維圖像的表面：其一是用三角形網格；其二是用比較隱式的方式，比如水平集函數.后者的優點是數值計算簡單與拓撲結構靈活多變.而拓撲結構的靈活多變是很重要的，它不僅能去噪，還能進行拓撲修正[8].

非局部均值方法是Buades等[5]提出來的，其主要用來對圖像進行去噪，具體過程為，

NL(u)(x)

(13)

式中，c(x)是一個標準化因子.

da(u(x),u(y))

(14)

式中，Ga是一個標準偏差為a的高斯函數，用非局部均值來對圖像進行去噪可以得到較好的效果，其能量函數[6]可以寫成，

(15)

該函數對應的梯度下降流是一個非局部熱方程，

ut=Δωu：

(16)

式中，x∈Ω.

式(16)經過有限差分法可得到其離散版本，

(17)

式中,uj表示u在網格點j處的值，j取遍計算域的所有的網格點.Nj是網格點j的鄰域使得當l∈Nj時候，ω(l,j)>0.CFL限制條件使得dt滿足,

(18)

1.2.2 模型求解

基于水平集的非局部表面的恢復算法具體步驟如下：

(1)用水平集函數φ來表示圖像的表面，

(19)

式中，∑表示一個區域，S表示該區域的邊界.

(2)把φ的0水平集附近的一個狹窄的區域記作∑η，η是該區域的寬度.

(3)計算權重ω(x,y)和相似函數D(x,y)，

(20)

x∈∑η,y∈Nx

(21)

式中,Nx是∑η中x的鄰域，φ[x]是φ以x為中心處的3D塊.

(4)對能量函數J(u)的梯度下降流使用有限差分法得到其離散的迭代格式如下，

(22)

式中，ωjl可以由上式計算得出，令dt為,

(23)

式中，j遍歷完∑η中所有的網格點.

(5)設置算法的迭代次數k.

2 結果與分析

為了測試本研究所提出方法的可行性，利用網上公開的工業部件的灰度體數據進行相關測試實驗.軟件平臺為,MATLAB_R2016b分析工具.硬件平臺為,組裝臺式機，其主要配置為,Inter(R)Core(TM)i5-6500 CPU@3.20 GHz，8 GB內存，64位Windows操作系統.

2.1 參數選取

本研究的方法中存在一些參數，而參數的選取會影響實驗的結果，故需要先對這些參數進行分析和討論.

首先，在3D-RCV算法中，參數λ1、λ2和水平集正則項參數ν在實驗中可以設置為λ1=λ2=1,ν=1，而面積項參數μ則需要根據待分割圖像的不同類型進行調整.面積項參數調整的規則是：當μ取得比較小的時候可以分割出圖像中的細小物體；如果需要分割出較大的目標，μ就應取得相對大一點.另外，算法中高斯窗口的半徑大小的選取對分割結果也有一定的影響，它的選取主要和噪聲的強弱有關.根據經驗來說，如果圖像中的噪聲較大，高斯窗口的半徑就應取得大一些，反之亦然.在實驗中，高斯窗口選取的半徑為7×7×7.

其次，基于水平集的非局部表面恢復算法中，參數b1是兩個網格點之間距離的懲罰權重，b2是兩個3D塊之間的懲罰相似性權重.b1比較大的時候，更能利用遠處的信息；而b2比較大時候，更能保護尖銳的特征.本研究根據經驗選取，b1為1 000，b2為600，區域的寬度η選擇為2.因為圖像所含的噪聲種類不同，難以找到一個統一的停止準則，所以算法的停止迭代步數通常根據經驗選取，本次實驗選取停止迭代步數為300步.

2.2 結果與分析

實驗中，本研究通過網上公開的發動機的灰度體數據來測試方法的可行性.

首先，利用3D-RCV算法對發動機的灰度體數據進行三維分割得到用水平集函數表示的工業部件的三維圖像表面，結果如圖1(a)所示.然后，采用基于水平集的非局部表面恢復算法對其進行平滑處理，結果如圖1(b)所示.

從圖1(a)可以看出，經過3D-RCV算法分割后得到的發動機的三維圖像表面有很多毛刺和部分不連續的地方.由圖1(b)可以看出，而經過基于水平集的非局部表面恢復算法平滑處理后，發動機的三維圖像表面在邊界得到保護的同時減少了大量的噪聲，在表面有毛刺和部分不連續的地方還進行了重造，圖像整體變得更加平滑.

(a)3D-RCV分割后的圖像

下面再從圖1中選取2個局部區域進行放大來具體說明上述結論.

圖2是圖1中用箭頭①所示矩形框標出來的左下部分的局部放大圖.圖2(a)中的圖像表面凹凸不平，整體非常的不光滑，而經過平滑處理后，圖像的表面整體變得非常平滑(見圖2(b)).

(a)3D-RCV分割后的圖像

圖3是圖1中用箭頭②所示矩形框標出來的右下部分的局部放大圖.圖3(a)中圖像存在有毛刺和部分不連續的地方，經過平滑處理后，其在邊界的大體形狀得到保護的同時，還進行了重造，圖像整體變得更為平滑(見圖3(b)).

(a)3D-RCV分割后的圖像

3 結 論

針對工業CT圖像所含噪聲過多而導致所得到的工業部件三維圖像表面有毛刺和部分不連續的情況，本研究將3D-RCV算法和基于水平集的非局部表面恢復算法相結合用于工業部件三維圖像的表面平滑處理.實驗測試結果表明，經過本方法平滑處理后的工業部件的三維圖像表面在邊界得到保護的同時減少了大量的噪聲，而且在表面有毛刺和部分不連續的地方進行了重造，圖像整體變得更為平滑.本研究為工業CT三維圖像的可視化軟件的逆向設計和改造提供了高質量的三維圖像，減少了后期處理的計算量和存儲量.