王佳熙

(成都大學(xué) 計(jì)算機(jī)學(xué)院，四川 成都 610106)

0 引 言

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

1 工業(yè)CT三維圖像表面平滑方法

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

1.1 3D-RCV算法

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

1.1.1 模型建立

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

(1)

式中，常數(shù)ci是子區(qū)域Ωi的平均灰度，I(y)是點(diǎn)x的鄰域N(x)內(nèi)的像素點(diǎn)y(x,y,z)的灰度值，鄰域N(x)的大小通過核函數(shù)Kσ來控制.這里，選取截?cái)嗟母咚购撕瘮?shù)作為Kσ，定義如下，

(2)

式中，σ是高斯函數(shù)的標(biāo)準(zhǔn)差，a是正則化常數(shù)，鄰域N(x)的半徑記為r.

x,y∈Ω

(3)

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

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

F(C,c1,c2)

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

(4)

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

(5)

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

(6)

(7)

為了保持水平集函數(shù)的光滑性，使用文獻(xiàn)[8]定義的水平集函數(shù)正則項(xiàng)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,分別是面積項(xiàng)和水平集正則項(xiàng)的參數(shù).

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

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

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

(10)

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

(11)

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

(12)

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

1.2.1 模型建立

近年來，圖像平滑算法已延拓到了三維圖像的表面平滑，其主要采取兩種方式來表示三維圖像的表面：其一是用三角形網(wǎng)格；其二是用比較隱式的方式，比如水平集函數(shù).后者的優(yōu)點(diǎn)是數(shù)值計(jì)算簡單與拓?fù)浣Y(jié)構(gòu)靈活多變.而拓?fù)浣Y(jié)構(gòu)的靈活多變是很重要的，它不僅能去噪，還能進(jìn)行拓?fù)湫拚齕8].

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

NL(u)(x)

(13)

式中，c(x)是一個標(biāo)準(zhǔn)化因子.

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

(14)

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

(15)

該函數(shù)對應(yīng)的梯度下降流是一個非局部熱方程，

ut=Δωu：

(16)

式中，x∈Ω.

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

(17)

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

(18)

1.2.2 模型求解

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

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

(19)

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

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

(3)計(jì)算權(quán)重ω(x,y)和相似函數(shù)D(x,y)，

(20)

x∈∑η,y∈Nx

(21)

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

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

(22)

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

(23)

式中，j遍歷完∑η中所有的網(wǎng)格點(diǎn).

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

2 結(jié)果與分析

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

2.1 參數(shù)選取

本研究的方法中存在一些參數(shù)，而參數(shù)的選取會影響實(shí)驗(yàn)的結(jié)果，故需要先對這些參數(shù)進(jìn)行分析和討論.

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

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

2.2 結(jié)果與分析

實(shí)驗(yàn)中，本研究通過網(wǎng)上公開的發(fā)動機(jī)的灰度體數(shù)據(jù)來測試方法的可行性.

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

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

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

下面再從圖1中選取2個局部區(qū)域進(jìn)行放大來具體說明上述結(jié)論.

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

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

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

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

3 結(jié) 論

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