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

李淑玲，李小林

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

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

李淑玲，李小林

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

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

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

1 引言

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

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

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

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

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

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

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

當系數(shù)矩陣：

可逆時，方程組（2）有唯一解。因為徑向基函數(shù)滿足ri(xj)=rj(xi)，所以R對稱，故R正定時，即徑向基函數(shù)嚴格正定時，R可逆。然而，許多全局徑向基函數(shù)僅僅條件正定。

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

是全局的和嚴格正定的[8]，本文選用該函數(shù)，此時系數(shù)矩陣R對稱正定，因此可逆，從而可由式（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ū)域：目標區(qū)域Ω1和背景區(qū)域Ω2。在水平集方法中，C、Ω1和Ω2滿足：

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

其中V(x，t)是速度函數(shù)，決定曲線C上每一點的運動速度，φ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）中，右邊第一項是通過最小化曲線C的長度得到的，由于用徑向基函數(shù)插值得到的近似函數(shù)全局光滑，該項可略去，故令γ=0；右邊第二項和第三項分別是通過最小化曲線C內(nèi)部和外部的灰度值得到的，參數(shù)λ1和λ2是這兩項的權(quán)系數(shù)，故可令λ1=λ2=1。最終，速度函數(shù)可表示為：

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

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

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

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

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

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

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

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

在初始演化時刻t=0，令式（8）在N個節(jié)點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是時間步長。

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

最后式（12）等價于：

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

5 實驗結(jié)果

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

5.1 算法驗證

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

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

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

表2 表1中分割的計算時間s

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

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

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

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

表4 表1中分割的DSC值

表5 表1中分割的RSE值

5.2 算法應(yīng)用

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

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

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

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

表7 表6中分割的計算時間s

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

6 結(jié)束語

本文提出了一種基于全局正定徑向基函數(shù)和偏微分方程水平集格式的圖像分割算法。該算法有效地解決了傳統(tǒng)水平集方法中需要不斷重新初始化水平集函數(shù)和對水平集初始輪廓位置敏感等問題。數(shù)值實驗結(jié)果表明，本文算法允許常值初始化方案，消除了水平集演化對初始輪廓的需要。另外，與基于徑向基函數(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].計算機工程與應(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].軟件學報，2008，19（12）：3161-3169.

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

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

[8]張雄，劉巖.無網(wǎng)格方法[M].北京：清華大學出版社，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]李慶揚，王能超，易大義.數(shù)值分析[M].5版.北京：清華大學出版社，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.

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

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

2013-09-10

2013-11-07

1002-8331（2014）06-0139-05