郭元卡

(西安電子科技大學(xué)電子工程學(xué)院，陜西西安 710071)

圖像分割是由圖像處理到圖像分析的關(guān)鍵步驟，分割的準(zhǔn)確性直接影響后續(xù)任務(wù)的有效性，因此具有重要意義。變分水平集方法［1］是直接通過(guò)極小化關(guān)于水平集函數(shù)自身的能量函數(shù)而得到演化的偏微分方程(PDE)，故具有參數(shù)化Snake模型的優(yōu)點(diǎn)，又具有水平集方法的優(yōu)點(diǎn)。相對(duì)于傳統(tǒng)的由純粹PDE驅(qū)動(dòng)的水平集圖像分割方法，基于變分水平集的圖像分割方法可以在能量函數(shù)中自然地融入附加約束信息，因而可以產(chǎn)生魯棒性更強(qiáng)的結(jié)果。

1 Level Set方法

水平集的基本思想是將曲線(xiàn)看成高一維空間中某一函數(shù)φ的零水平集，同時(shí)曲線(xiàn)的演化也擴(kuò)充到高一維的空間中。將水平集函數(shù)按照它所滿(mǎn)足的演化方程進(jìn)行演化或迭代，由于水平集函數(shù)不斷進(jìn)行演化，所以對(duì)應(yīng)的零水平集也在不斷變化，當(dāng)水平集演化趨于平穩(wěn)時(shí)，演化停止，得到曲線(xiàn)形狀。

考慮零水平集x(t)所對(duì)應(yīng)的水平集函數(shù)φ，則有

對(duì)方程兩邊求關(guān)于時(shí)間的偏導(dǎo)數(shù)，有

假設(shè)F為外法向方向的速度，那么

基本方程式是水平集函數(shù)及相應(yīng)的水平集在法向力F推動(dòng)下的演化方程。

2 算法原理

傳統(tǒng)的水平集方法在每次迭代后，都需要對(duì)水平集函數(shù)進(jìn)行重新初始化，以此修正迭代造成的水平集位置偏離。雖然這樣保證了演化計(jì)算的穩(wěn)定性，但花費(fèi)時(shí)間長(zhǎng)，而且重新初始化較復(fù)雜，嚴(yán)重阻礙了水平集方法的廣泛應(yīng)用。

Gomes J與 Faugeras O［2］指出，重新初始化水平集函數(shù)是水平集理論與實(shí)踐的矛盾之處。理論上，重新初始化沒(méi)有必要;實(shí)踐中，為保證水平集在演化過(guò)程中不會(huì)偏離符號(hào)距離函數(shù)(SDF)太遠(yuǎn)，時(shí)間步長(zhǎng)必須選得足夠小，這樣做又增加了迭代次數(shù)，增加了演化的時(shí)間。下面介紹一種無(wú)需重新初始化的變分水平集算法［3］。

(1)內(nèi)部能量。設(shè)φ為水平集函數(shù)，SDF滿(mǎn)足一個(gè)非常好的性質(zhì)，即反之，滿(mǎn)足1的任意函數(shù)φ都可表示成一個(gè)SDF加一個(gè)常數(shù)［4］。所以引出如下積分

作為函數(shù)φ與Ω中SDF接近程度的一個(gè)衡量標(biāo)準(zhǔn)。

根據(jù)上述定義的函數(shù)P(φ)，提出以下變分公式

式中，μ＞0為常數(shù)，它控制著“懲罰”水平集φ偏離SDF的力度，P(φ)稱(chēng)為函數(shù)φ的內(nèi)部能量;基于圖像數(shù)據(jù)的能量εm(φ)用于驅(qū)動(dòng)φ的零水平集曲線(xiàn)向圖像的邊緣演化，稱(chēng)為外部能量。

他是最小化函數(shù)ε的梯度流［5］。

函數(shù)φ根據(jù)式(7)表示的梯度流進(jìn)行演化，在φ的演化過(guò)程中，零水平集曲線(xiàn)通過(guò)外部能量εm(φ)進(jìn)行演化;與此同時(shí)，由于內(nèi)部能量“懲罰”的作用，演化函數(shù)將自動(dòng)保持為一個(gè)近似的SDF。因此，在提出的公式中重新初始化步驟被完全消除了。

(2)外部能量。

設(shè)I為原始圖像，g為邊緣檢測(cè)函數(shù)，定義如下

式中，Gσ是標(biāo)準(zhǔn)差為σ的高斯核。當(dāng)輪廓線(xiàn)位于圖像的邊緣處時(shí)，圖像的高斯梯度較大，邊緣檢測(cè)函數(shù)g較小，曲線(xiàn)的演化速度幾乎為零，輪廓線(xiàn)就會(huì)停止在圖像的邊緣處。

定義函數(shù)φ(x，y)的外部能量如下

式中，λ＞0和v均為常數(shù)。Lg(φ)和Ag(φ)定義如下

式中，δ為單變量Dirac函數(shù)，H為Heaviside函數(shù)。

假設(shè)φ的零水平集用可微的參數(shù)化曲線(xiàn)C(p)，p∈［0，1］表示，那么Lg(φ)表示的是在正形投影公式中φ的零水平集曲線(xiàn)的長(zhǎng)度［6］。

能量函數(shù)Ag(φ)用來(lái)加速曲線(xiàn)的演化，可以看出，當(dāng)邊緣檢測(cè)函數(shù)g為常數(shù)1時(shí)，Ag(φ)即為Ωφ=區(qū)域的加權(quán)面積。一般來(lái)說(shuō)，當(dāng)初始輪廓線(xiàn)位于目標(biāo)外部時(shí)，v為正值，這樣可以促使輪廓線(xiàn)快速收縮;反之，當(dāng)初始輪廓線(xiàn)位于目標(biāo)內(nèi)部時(shí)，v為負(fù)值，可以加快輪廓線(xiàn)的膨脹速度。

于是，得到完整的能量函數(shù)

式中，Δ為L(zhǎng)aplacian算子;▽?duì)帐铅盏奶荻?div為散度。這一梯度流即是零水平集函數(shù)的演化方程。第2、3項(xiàng)由外部能量而來(lái)，用于驅(qū)動(dòng)φ的零水平集曲線(xiàn)向目標(biāo)邊緣運(yùn)動(dòng)。第1項(xiàng)中的

3 算法實(shí)現(xiàn)

(1)偏微分方程的求解。該算法采用有限差分的方法來(lái)求解式(13)。在實(shí)際應(yīng)用中，Dirac函數(shù)采用正則化的 δε(x)進(jìn)行計(jì)算，δε(x)定義如下

文中所有實(shí)驗(yàn)選取ε=1.5。

(2)時(shí)間步長(zhǎng)的選擇。在實(shí)現(xiàn)水平集方法時(shí)，可以選擇的時(shí)間步長(zhǎng)τ的范圍明顯大于在傳統(tǒng)水平集方法中使用的時(shí)間步長(zhǎng)范圍，在這個(gè)算法中，τ可以設(shè)在0.1～100之間。實(shí)驗(yàn)中，為保持水平集演化的穩(wěn)定性，時(shí)間步長(zhǎng)τ和系數(shù)μ必須滿(mǎn)足

采用較大的時(shí)間步長(zhǎng)可以加速演化，但是，如果時(shí)間步長(zhǎng)的選擇過(guò)大，在邊緣定位時(shí)可能出錯(cuò)。一般地，對(duì)于大多數(shù)圖像選擇τ≤10.0。

(3)水平集函數(shù)的初始化。在傳統(tǒng)的水平集方法中，必須將水平集函數(shù)φ初始化為符號(hào)距離函數(shù)φ0，如果初始水平集函數(shù)明顯不同于SDF，那么重新初始化步驟也不能將這個(gè)函數(shù)重新初始化為一個(gè)SDF。在算法中，不僅重新初始化步驟得以完全消除，而且水平集函數(shù)也不再需要初始化為SDF。

這里提出以下函數(shù)作為初始函數(shù)φ0。令Ω0是圖像區(qū)域Ω的一個(gè)子集，?Ω0是Ω0的所有邊界點(diǎn)的集合，它們能夠被一些簡(jiǎn)單的形態(tài)學(xué)操作有效識(shí)別。于是，初始函數(shù)φ0定義如下

式中，ρ＞0是一個(gè)常數(shù)。ρ一般選取ρ≥2ε，ε為正則化Dirac函數(shù)δε(x)定義中的寬度。

與傳統(tǒng)由輪廓計(jì)算得來(lái)的SDF不同，上述初始水平集函數(shù)是從圖像區(qū)域Ω中任意一個(gè)子區(qū)域Ω0中計(jì)算而來(lái)，水平集函數(shù)的這種基于區(qū)域的初始化方法不僅計(jì)算高效，而且在一些情況下應(yīng)用較為靈活。例如，如果感興趣區(qū)域可以通過(guò)一些途徑粗略和自動(dòng)獲取，那么可以使用這些粗略獲得的區(qū)域作為區(qū)域Ω0來(lái)構(gòu)造初始水平集函數(shù)φ0。接下來(lái)，初始水平集函數(shù)將根據(jù)演化方程穩(wěn)定演化，同時(shí)它的零水平曲線(xiàn)可以收斂于感興趣區(qū)域的精確邊緣。

雖然這種初始函數(shù)φ0明顯偏離了SDF，但由于能量函數(shù)中的懲罰函數(shù)P(φ)的調(diào)節(jié)作用，使水平集函數(shù)φ在其零水平集附近能夠近似等于SDF，這對(duì)于邊緣的定位沒(méi)有太大影響。

(4)水平集函數(shù)收斂的判斷。水平集函數(shù)收斂的判斷條件為

算法的實(shí)現(xiàn)步驟為:

(1)初始化:設(shè)置初始輪廓，可以設(shè)置為圓形、矩形或是一條非閉合曲線(xiàn)，根據(jù)式(18)計(jì)算水平集函數(shù)的初始值φ0。

(2)利用迭代公式(17)，計(jì)算φn+1。

(3)檢查本次迭代是否收斂，如果收斂，則停止計(jì)算;否則，轉(zhuǎn)至步驟(2)，繼續(xù)計(jì)算。

(4)更新輪廓線(xiàn)，獲得分割圖像。

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

用多幅圖像對(duì)上述算法進(jìn)行了實(shí)驗(yàn)，獲得了較好的分割效果。所有實(shí)驗(yàn)都采用式(18)所定義的初始化方法。

圖1中，參數(shù)λ=5.0，μ=0.04，v=1.5，時(shí)間步長(zhǎng)τ=5給出了一幅像素大小為83×65的顯微鏡下觀測(cè)到的兩個(gè)細(xì)胞圖像的輪廓演化過(guò)程。可以看到，兩個(gè)細(xì)胞的某些邊緣非常模糊，演化曲線(xiàn)自動(dòng)分裂，將兩個(gè)細(xì)胞分割出來(lái)，實(shí)驗(yàn)結(jié)果表明這種方法在處理弱目標(biāo)邊緣時(shí)有較強(qiáng)的魯棒性。實(shí)驗(yàn)中，選擇了一個(gè)矩形區(qū)域作為區(qū)域Ω0來(lái)計(jì)算初始水平集函數(shù)φ0，如圖1(a)中所示。

圖1 細(xì)胞圖像(83×65)分割結(jié)果

圖2中，參數(shù)λ=5.0，μ=0.04，v=3.0，時(shí)間步長(zhǎng)τ=5給出了一幅像素大小為84×84的兩個(gè)目標(biāo)物圖像的輪廓演化過(guò)程。這個(gè)時(shí)間步長(zhǎng)明顯大于在傳統(tǒng)水平集方法中使用的時(shí)間步長(zhǎng)。在傳統(tǒng)的水平集方法中，初始輪廓線(xiàn)一般選取一條簡(jiǎn)單閉合曲線(xiàn)。而圖2(a)中使用了一條非閉合的直線(xiàn)作為初始輪廓線(xiàn)，演化曲線(xiàn)將兩個(gè)目標(biāo)物都提取了出來(lái)，這是傳統(tǒng)水平集方法做不到的。

圖2 兩個(gè)目標(biāo)物圖像(84×84)分割結(jié)果

圖3與圖2參數(shù)相同給出了一幅像素大小為128×128的腦部CT圖像的去骨分割過(guò)程。

圖3 腦部CT圖像(128×128)的去骨分割結(jié)果

這種方法具有幾個(gè)顯著優(yōu)點(diǎn):它使得真正影響水平集計(jì)算量的重新初始化步驟得以完全消除;靈活的水平集函數(shù)初始化方法，使得初始輪廓線(xiàn)的選擇更加自由，并且計(jì)算也簡(jiǎn)便了很多。但該模型中的初始輪廓線(xiàn)的選取依然受到很大限制，初始輪廓線(xiàn)需包圍所有待分割的物體，而且不能跨越同質(zhì)區(qū)域。

5 結(jié)束語(yǔ)

本文提出的無(wú)需重新初始化的變分水平集模型，獲得了滿(mǎn)意的分割效果。但由于這個(gè)模型僅利用了圖像的邊緣梯度信息，而沒(méi)有考慮圖像的全局區(qū)域信息，使得模型對(duì)邊緣模糊圖像、噪聲圖像、紋理圖像的分割效果不理想。此外，該模型的迭代次數(shù)過(guò)多，初始輪廓線(xiàn)的選取依然受到很大限制。

［1］ ZHAO H K，CHAN T F，MERRIMAN B，et al.A variational level set approach to multiphase motion［J］.Journal Of Computational Physics，1996，167:179 －195.

［2］ GOMES J，F(xiàn)AUGERAS O.Reconciling distance functions and level sets［J］.J.Visiual Communic and Image Representation，2000，11(2):209－223.

［3］ LI Chunming，XU Chengyang，GUI Changfeng，et al.Level set evolution without reinitialization:a new variational formulation［C］.San Diego:IEEE International Conference on Computer Vision and Pattern Recognition(CVPR)，2005:430－436.

［4］ ARNOLD V I.Geometrical methods in the theory of ordinary differential equations［M］.New York:Springer－Verlag，1983.

［5］ EVANS L C.Partial differential equations［M］.Providence R I:American Mathematical Society，1998.

［6］ VEMURI B，CHEN Y.Joint image registration and segmentation［M］.New York:Geometric Level Set Methods in Imaging，2003.

［7］ OSHER S，F(xiàn)EDKIW R.Level set methods and dynamic implicit surfaces［M］.New York:Spring－Verlag，2002.