李淑玲，李小林

重慶師范大學(xué)數(shù)學(xué)學(xué)院，重慶 401331

全局正定徑向基函數(shù)在圖像分割中的應(yīng)用

李淑玲，李小林

重慶師范大學(xué)數(shù)學(xué)學(xué)院，重慶 401331

將全局正定徑向基函數(shù)和圖像分割中基于偏微分方程水平集方法的發(fā)展方程相結(jié)合，提出了一種基于全局正定徑向基函數(shù)的圖像分割算法。用全局正定徑向基函數(shù)插值發(fā)展方程中的水平集函數(shù)，得到的插值函數(shù)具有較高的精度和光滑性，克服了傳統(tǒng)水平集方法中復(fù)雜費(fèi)時(shí)的重新初始化過程和水平集對(duì)初始輪廓位置敏感等缺點(diǎn)，非線性發(fā)展方程最終被轉(zhuǎn)化成常微分方程組并用Euler法求解。實(shí)驗(yàn)結(jié)果表明該算法不需要重新初始化過程，并且在沒有初始輪廓時(shí)也能夠快速正確地分割圖像。

圖像分割；徑向基函數(shù)；偏微分方程；發(fā)展方程；水平集

1 引言

圖像分割是數(shù)字圖像處理和圖形分析中最基礎(chǔ)的關(guān)鍵步驟，其任務(wù)是從圖像中分割出目標(biāo)并得到相應(yīng)的邊界曲線。水平集方法是一種重要的圖像分割算法，它把邊界曲線當(dāng)成零水平集嵌入到高一維水平集函數(shù)中，通過更新水平集函數(shù)來演化曲線，從而將曲線的演化轉(zhuǎn)化為偏微分方程的求解。與傳統(tǒng)的圖像分割算法相比，水平集方法拓?fù)溥m應(yīng)性強(qiáng)，可以分割結(jié)構(gòu)復(fù)雜的物體，因此已成為當(dāng)今圖像分割領(lǐng)域的研究熱點(diǎn)之一[1-3]。

由水平集方法得到的偏微分方程依賴于時(shí)間，是一類非線性的水平集發(fā)展方程，通常需用迎風(fēng)有限差分法數(shù)值求解，這種差分法計(jì)算量大，對(duì)水平集初始輪廓的位置比較敏感，同時(shí)為了保證水平集函數(shù)的規(guī)則性和穩(wěn)定性，在演化過程中需要周期性地重新初始化水平集函數(shù)。重新初始化算法復(fù)雜，計(jì)算耗時(shí)，而且數(shù)值誤差會(huì)使零水平集偏離目標(biāo)位置。近年來出現(xiàn)了一些不需要重新初始化的水平集方法[4-7]，這些方法需要添加一些額外項(xiàng)來避免重新初始化，增加了相應(yīng)算法的復(fù)雜度和計(jì)算量。另外，傳統(tǒng)水平集方法需要人工恰當(dāng)定義初始輪廓，而初始輪廓的位置、大小和形狀都可能導(dǎo)致不同甚至完全錯(cuò)誤的分割結(jié)果，如何有效地解決輪廓初始化問題仍具挑戰(zhàn)性。

徑向基函數(shù)是一類以源點(diǎn)和場點(diǎn)之間的距離為自變量的函數(shù)，已成為數(shù)值求解科學(xué)工程問題的有力工具[8-10]。近年來，Wang等[11]用全局徑向基函數(shù)求解了結(jié)構(gòu)拓?fù)鋬?yōu)化中的水平集發(fā)展方程。在文獻(xiàn)[11]的基礎(chǔ)上，Xie等[12]用全局徑向基函數(shù)求解了可變形模型中的水平集發(fā)展方程，Gelas等[13]用緊支徑向基函數(shù)求解了圖像分割中的水平集發(fā)展方程。全局徑向基函數(shù)的插值精度高，但文獻(xiàn)[11]和文獻(xiàn)[12]使用的全局徑向基函數(shù)是條件正定的[8]，因此需要添加額外的方程來保證唯一性，增加了計(jì)算難度。文獻(xiàn)[13]使用的緊支徑向基函數(shù)雖然嚴(yán)格正定，但比全局徑向基函數(shù)的插值精度要低一些[8]，致使文獻(xiàn)[13]的圖像分割算法收斂速度較慢。另外，文獻(xiàn)[13]中的算法求解的是基于變分水平集方法的發(fā)展方程，涉及正則化Dirac函數(shù)δε(φ)。因?yàn)棣摩?φ)具有緊支集，所以發(fā)展方程的控制作用是局部的，從而限制了水平集函數(shù)演化邊界曲線的效率[1]。為此，可將δε(φ)用‖▽?duì)铡?，得到更高效的基于偏微分方程水平集方法的發(fā)展方程[1]。

針對(duì)圖像分割中基于偏微分方程水平集方法的發(fā)展方程，本文提出了一種基于全局正定徑向基函數(shù)的圖像分割算法。在該算法中，通過用徑向基函數(shù)插值水平集函數(shù)，原來的發(fā)展方程被轉(zhuǎn)化為常微分方程組。由于使用的徑向基函數(shù)是全局的和嚴(yán)格正定的，從而避免了文獻(xiàn)[11]中的額外方程，且比文獻(xiàn)[13]的算法收斂快。另外，由于徑向基函數(shù)插值的精度和光滑性都較高，所以與傳統(tǒng)水平集方法相比，本文算法不需要復(fù)雜費(fèi)時(shí)的重新初始化過程，并且對(duì)水平集初始輪廓的位置也不敏感。實(shí)驗(yàn)結(jié)果表明本文算法比基于有限差分格式的水平集算法[2]和文獻(xiàn)[13]中的算法具有更快的收斂速度。

2 徑向基函數(shù)插值

設(shè)Ω是平面上的有界區(qū)域，平面上的任一點(diǎn)記為x=(x，y)。將Ω用N個(gè)節(jié)點(diǎn)xj(j=1，2，…，N)離散，則函數(shù)f(x)在Ω內(nèi)的徑向基函數(shù)插值為：

其中ri(x)是以節(jié)點(diǎn)xi為中心的徑向基函數(shù)，rT(x)= [r1(x),r2(x),…,rN(x)]，ai是待定系數(shù)，a=[a1,a2,…，aN]T。令式（1）在xj滿足，可得：

當(dāng)系數(shù)矩陣：

可逆時(shí)，方程組（2）有唯一解。因?yàn)閺较蚧瘮?shù)滿足ri(xj)=rj(xi)，所以R對(duì)稱，故R正定時(shí)，即徑向基函數(shù)嚴(yán)格正定時(shí)，R可逆。然而，許多全局徑向基函數(shù)僅僅條件正定。

全局徑向基函數(shù)插值具有很高的精度，但文獻(xiàn)[11-12]使用的全局徑向基函數(shù)是條件正定的[8]，所以需要添加方程來保證唯一性，這增加了計(jì)算的難度。文獻(xiàn)[13]使用的緊支徑向基函數(shù)是嚴(yán)格正定的，但與全局徑向基函數(shù)相比，緊支徑向基函數(shù)的插值精度要低一些，導(dǎo)致相應(yīng)圖像分割算法的收斂速度較慢。由于逆MQ徑向基函數(shù)：

是全局的和嚴(yán)格正定的[8]，本文選用該函數(shù)，此時(shí)系數(shù)矩陣R對(duì)稱正定，因此可逆，從而可由式（2）解出系數(shù)向量a，最后代入式（1）可得：

其中f=[f(x1)，f(x2)，…，f(xN)]T，rT(x)R-1是徑向基函數(shù)插值得到的形函數(shù)向量。

3 基于偏微分方程水平集方法的分割模型

3.1 水平集發(fā)展方程

令Ω為圖像區(qū)域，圖像分割旨在尋找曲線C，把Ω分割成不相重疊的區(qū)域：目標(biāo)區(qū)域Ω1和背景區(qū)域Ω2。在水平集方法中，C、Ω1和Ω2滿足：

其中t是時(shí)間變量。水平集函數(shù)φ(x，t)是由曲線C生成的符號(hào)距離函數(shù)，滿足如下Cauchy問題[1]：

其中V(x，t)是速度函數(shù)，決定曲線C上每一點(diǎn)的運(yùn)動(dòng)速度，φ0(x)表示曲線C的初始位置。

3.2 C-V模型

構(gòu)建水平集發(fā)展方程（5）的模型很多。C-V模型[2]是最有代表性的模型之一，該模型的速度函數(shù)為：

其中γ≥0、λ1＞0和λ2＞0都是固定的參數(shù)，I(x)是圖像的灰度。

分別是曲線C內(nèi)部和外部的圖像灰度平均值。

是正則化的Heaviside函數(shù)，ε是參數(shù)。

在式（7）中，右邊第一項(xiàng)是通過最小化曲線C的長度得到的，由于用徑向基函數(shù)插值得到的近似函數(shù)全局光滑，該項(xiàng)可略去，故令γ=0；右邊第二項(xiàng)和第三項(xiàng)分別是通過最小化曲線C內(nèi)部和外部的灰度值得到的，參數(shù)λ1和λ2是這兩項(xiàng)的權(quán)系數(shù)，故可令λ1=λ2=1。最終，速度函數(shù)可表示為：

4 圖像分割模型的數(shù)值解法

水平集發(fā)展方程（5）是隨時(shí)間變化的偏微分方程。在傳統(tǒng)水平集方法中，數(shù)值求解Cauchy問題式（5）、（6）時(shí)需要適當(dāng)?shù)挠L(fēng)差分格式和周期性重新初始化算法，這增加了算法的復(fù)雜度和計(jì)算量，從而限制了水平集方法在圖像分割中的應(yīng)用。

本文通過用徑向基函數(shù)插值水平集函數(shù)φ(x，t)，偏微分方程Cauchy問題式（5）、（6）被轉(zhuǎn)換成極易求解的常微分方程Cauchy問題。另外，用徑向基函數(shù)得到的近似水平集函數(shù)全局光滑，能有效地避免復(fù)雜費(fèi)時(shí)的重新初始化過程和輪廓初始化問題，從而降低了算法的復(fù)雜度和計(jì)算量。

因?yàn)樗郊瘮?shù)φ(x，t)中的時(shí)間變量是人為添加的，所以根據(jù)文獻(xiàn)[11-13]，可以假定φ(x，t)在時(shí)間和空間上是可分離的，從而由式（1），φ(x，t)的徑向基函數(shù)插值可表示為：

其中a(t)=[a1(t)，a2(t)，…，aN(t)]T是僅依賴于時(shí)間變量t的未知向量。把式（8）代入式（5）可得：

令上式在節(jié)點(diǎn)xj滿足，得到常微分方程組：

式（9）在形狀和拓?fù)浣Y(jié)構(gòu)優(yōu)化時(shí)很有效[11]，但由于水平集函數(shù)很可能存在駐點(diǎn)，致使水平集函數(shù)在圖像區(qū)域內(nèi)的演化變慢，所以式（9）在圖像分割時(shí)可能失效。為此，可將V(xj，t)替換為規(guī)范化形式[12]，V(xj，t)/ ||▽(rT(xj)a(t))||，則式（9）簡化為：

其中R是式（3）給出的插值矩陣，a(t)是未知向量，

在初始演化時(shí)刻t=0，令式（8）在N個(gè)節(jié)點(diǎn)xj上得到滿足，再利用式（6）可得：

其中Φ0=[φ0(x1)，φ0(x2)，…，φ0(xN)]T已知，據(jù)此可得方程組（10）的初始條件為：

從而，偏微分方程Cauchy問題式（5）、（6）轉(zhuǎn)化為常微分方程Cauchy問題式（10）、（11）。

由于插值矩陣R可逆，式（10）可寫成：

利用向前Euler方法[14]，上式可離散為：

其中τ＞0是時(shí)間步長。

在式（12）中，直接求R-1的計(jì)算量較大。由于R對(duì)稱正定，故利用Cholesky分解法[14]可避免直接計(jì)算R-1，為此將R分解為下三角矩陣L及其轉(zhuǎn)置的乘積，即R=LLT，則在每次迭代時(shí)先求解三角形方程組：

最后式（12）等價(jià)于：

類似文獻(xiàn)[13]中性質(zhì)1，可證||a(tn)||會(huì)隨著迭代次數(shù)n緩慢增長。由于水平集函數(shù)φ(x，t)與任一非零實(shí)數(shù)的乘積不影響零水平集（即曲線C）的位置，因此在每次計(jì)算式（13）之后，可將a(tn)規(guī)范化，即令a(tn+1)=a(tn+1)/ ||a(tn+1)||，此時(shí)迭代終止條件[13]可取為||a(tn)-a(tn-1)||≤||a(tn-1)-a(tn-2)||。

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

分別采用本文的基于全局正定徑向基函數(shù)的偏微分方程水平集算法、文獻(xiàn)[13]的基于緊支徑向基函數(shù)的變分水平集算法以及文獻(xiàn)[2]的基于有限差分格式的C-V模型水平集算法對(duì)一些圖像進(jìn)行了分割實(shí)驗(yàn)。參數(shù)選取如下：對(duì)本文算法，Heaviside函數(shù)中的參數(shù)ε=1，時(shí)間步長τ=1；對(duì)文獻(xiàn)[13]的算法，Dirac函數(shù)和Heaviside函數(shù)中的參數(shù)ε=1，時(shí)間步長τ=1；對(duì)文獻(xiàn)[2]的算法，C-V模型參數(shù)ν=0，μ=0.5×2552，λ1=λ2=1，時(shí)間步長τ=0.1，并且水平集每更新10次重新初始化水平集函數(shù)1次。實(shí)驗(yàn)平臺(tái)是操作系統(tǒng)為Windows XP的PC（Intel?CoreTM2 Duo CPU，1.58 GHz，2.00 GB內(nèi)存），程序編寫使用Matlab7.0.1。

5.1 算法驗(yàn)證

表1給出了使用三種算法分割一幅人造圖像（表1（a3））的結(jié)果。結(jié)果表明，采用本文算法（表1（b1～b3））和文獻(xiàn)[13]的算法（表1（c1～c3）），水平集函數(shù)的初始輪廓曲線無論置于何處（表1（a1，a2）），甚至沒有給定初始輪廓時(shí)（表1（a3）），均能正確分割出目標(biāo)物體，而采用文獻(xiàn)[2]的算法（表1（d1～d3）），只有當(dāng)初始輪廓完全包圍目標(biāo)物體時(shí)，才能分割（表1（d1）），否則不能（表1（d2，d3））。

表1 三種算法對(duì)應(yīng)于不同初始輪廓的分割結(jié)果

表2和表3分別給出了三種算法的計(jì)算時(shí)間和迭代次數(shù)。對(duì)每種初始輪廓，本文算法所需迭代次數(shù)和迭代時(shí)間都最少，并且對(duì)初始輪廓最不敏感，因此更高效。

表2 表1中分割的計(jì)算時(shí)間s

為了定量地評(píng)估三種算法分割結(jié)果的精確性，計(jì)算了兩種區(qū)域交疊性度量：Dice相似性系數(shù)（DSC）[15]和錯(cuò)誤分割率（RSE）[16]，計(jì)算公式為：

表3 表1中分割的迭代次數(shù)

其中N(·)表示某閉合區(qū)域中像素個(gè)數(shù)，Re表示精確的目標(biāo)區(qū)域，Rn表示數(shù)值算法獲得的目標(biāo)區(qū)域。顯然，DSC值越接近1，同時(shí)RSE值越接近0，分割結(jié)果越精確。

表4和表5給出了三種算法的量化評(píng)估結(jié)果。對(duì)三種初始輪廓曲線，本文算法的DSC值接近1，而RSE值接近0，說明本文算法的分割結(jié)果很精確，文獻(xiàn)[13]的算法具有同樣的精度，但文獻(xiàn)[2]算法的精度要差一些。

表4 表1中分割的DSC值

表5 表1中分割的RSE值

5.2 算法應(yīng)用

表6給出了分割五幅測試圖像的結(jié)果。原始圖像放在第一行，它們分別是：含有孔洞區(qū)域的圖像（表6（a1）），部分目標(biāo)跨越邊界的多目標(biāo)米粒圖像（表6（a2）），目標(biāo)邊緣較弱的多目標(biāo)細(xì)胞圖像（表6（a3）），含深度凹陷的圖像（表6（a4））和鞍馬形深度圖像（表6（a5））。對(duì)表6給出的所有分割結(jié)果，水平集演化都開始于常值函數(shù)，因此沒有初始輪廓。從視覺上觀察，本文算法和文獻(xiàn)[13]的算法都得到了很好的分割結(jié)果。在實(shí)驗(yàn)中發(fā)現(xiàn)文獻(xiàn)[2]的C-V模型算法在迭代2 500次后仍然無法分割出目標(biāo)物體。

表6 比較兩種算法的分割結(jié)果

值得注意的是，對(duì)表6（a2，a3）給出的部分目標(biāo)跨越邊界的多目標(biāo)圖像，無法選取一條初始輪廓曲線將所有目標(biāo)都包含在內(nèi)，所以那些需要事先定義初始輪廓的算法，只使用一個(gè)初始輪廓很難分割這類圖像，而本文算法在沒有給定初始輪廓時(shí)，經(jīng)過幾次迭代便能準(zhǔn)確分割所有目標(biāo)。另外，表6（a5）給出的深度圖像的灰度分布較為復(fù)雜，鞍馬的底邊緣呈階梯狀，而側(cè)邊緣呈屋頂狀，基于常值逼近的圖像分割算法很難分割這類深度圖像，但本文算法很好地提取出了包括階梯狀和屋頂狀邊緣在內(nèi)的所有目標(biāo)邊緣。

表7和表8分別給出了兩種算法的計(jì)算時(shí)間和迭代次數(shù)。對(duì)所有圖像，本文算法在迭代次數(shù)和計(jì)算時(shí)間方面都優(yōu)于文獻(xiàn)[13]的算法，因此本文算法具有更快的收斂速度和更高的計(jì)算效率。

表7 表6中分割的計(jì)算時(shí)間s

表8 表6中分割的迭代次數(shù)

6 結(jié)束語

本文提出了一種基于全局正定徑向基函數(shù)和偏微分方程水平集格式的圖像分割算法。該算法有效地解決了傳統(tǒng)水平集方法中需要不斷重新初始化水平集函數(shù)和對(duì)水平集初始輪廓位置敏感等問題。數(shù)值實(shí)驗(yàn)結(jié)果表明，本文算法允許常值初始化方案，消除了水平集演化對(duì)初始輪廓的需要。另外，與基于徑向基函數(shù)和變分水平集格式的圖像分割算法以及基于有限差分格式的水平集方法相比，本文算法的迭代次數(shù)少、計(jì)算時(shí)間短，并且在沒有水平集初始輪廓時(shí)也能正確地分割圖像。

[1]Tsai R，Osher S.Level set methods and their applications in image science[J].Comm Math Sci，2003，1：623-656.

[2]Chan T，Vese L.Active contours without edges[J].IEEE T Image Process，2001，10：266-277.

[3]楊智鵬，楊玲.自適應(yīng)水平集模型在云圖弱邊界分割的應(yīng)用[J].計(jì)算機(jī)工程與應(yīng)用，2013，49（12）：116-120.

[4]Zhang K H，Zhang L，Song H H，et al.Re-initialization free level set evolution via reaction diffusion[J].IEEE T Image Process，2013，22：258-271.

[5]何傳江，李夢，詹毅.用于圖像分割的自適應(yīng)距離保持水平集演化[J].軟件學(xué)報(bào)，2008，19（12）：3161-3169.

[6]張世征，何傳江，原野.結(jié)合局部熵的無需重新初始化水平集演化[J].計(jì)算機(jī)工程與應(yīng)用，2010，46（18）：174-176.

[7]孫文杰，陳允杰，湯楊，等.一種改進(jìn)的活動(dòng)區(qū)域輪廓模型——無需水平集重新初始化[J].計(jì)算機(jī)工程與應(yīng)用，2008，44（2）：8-11.

[8]張雄，劉巖.無網(wǎng)格方法[M].北京：清華大學(xué)出版社，2004.

[9]Li X L，Zhu J L，Zhang S G.A hybrid radial boundary node method based on radial basis point interpolation[J]. Eng Anal Bound Elem，2009，33：1273-1283.

[10]Li X L，Zhu J L.The method of fundamental solutions for nonlinear elliptic problems[J].Eng Anal Bound Elem，2009，33：322-329.

[11]Wang S Y，Wang M Y.Radial basis functions and level set method for structural topology optimization[J].Int J Numer Meth Eng，2006，65：2060-2090.

[12]Xie X H，Mirmehdi M.Radial basis function based level set interpolation and evolution for deformable modelling[J].Image Vision Comput，2011，29：167-177.

[13]Gelas A，Bernard O，F(xiàn)riboulet D，et al.Compactly supported radial basis functions based collocation method for level-set evolution[J].IEEE T Image Process，2007，16：1873-1887.

[14]李慶揚(yáng)，王能超，易大義.數(shù)值分析[M].5版.北京：清華大學(xué)出版社，2008.

[15]Shattuck D W，Sandor-Leahy S R.Magnetic resonance image tissue classification using a partial volume model[J]. Neuroimage，2001，13：856-876.

[16]Liu B，Cheng H D，Huang J，et al.Probability density difference-based active contour for ultrasound image segmentation[J].Pattern Recognition，2010，43：2028-2042.

LI Shuling,LI Xiaolin

College of Mathematics Science,Chongqing Normal University,Chongqing 401331,China

A numerical algorithm based on globally supported and positive definite radial basis functions is developed in this paper to solve evolution equations which arise in image segmentation using partial differential equation based level set methods.In this algorithm,radial basis functions are used to interpolate level set functions with a high level of precision and smoothness,and then the nonlinear evolution equation is cast into ordinary differential equations and Euler’s scheme is used.Compared with conventional level set methods,this algorithm is free from the initial contour,and completely eliminates the need of the complex and costly re-initialization procedure.Experimental results indicate that the presented algorithm is free of re-initialization,and can segment images quickly even without any initial contour.

image segmentation;radial basis function;partial differential equation;evolution equation;level set

A

TP391.4

10.3778/j.issn.1002-8331.1309-0088

LI Shuling,LI Xiaolin.Application of image segmentation algorithm based on globally supported and positive definite radial basis functions.Computer Engineering and Applications,2014,50（6）：139-143.

國家自然科學(xué)基金（No.11101454）；重慶市教委科學(xué)技術(shù)研究項(xiàng)目（No.KJ130626）；重慶高校創(chuàng)新團(tuán)隊(duì)建設(shè)計(jì)劃資助項(xiàng)目（No.KJTD201308）。

李淑玲（1983—），女，碩士研究生，研究領(lǐng)域?yàn)閳D像處理；李小林（1983—），男，博士，副教授。E-mail：shuling1124@163.com

2013-09-10

2013-11-07

1002-8331（2014）06-0139-05