融合圖像局部能量和梯度的水平集分割方法

包立君，劉宛予，浦昭邦

(哈爾濱工業(yè)大學(xué)HIT-INSA中法生物醫(yī)學(xué)圖像聯(lián)合研究中心，150001哈爾濱，baolij@gmail.com)

水平集方法易于處理零水平集拓?fù)浣Y(jié)構(gòu)的變化，易于擴(kuò)展，對(duì)于結(jié)構(gòu)復(fù)雜的圖像和多目標(biāo)圖像的分割表現(xiàn)良好，常用于生物醫(yī)學(xué)圖像的分割［1-2］.水平集分割方法可分為基于邊緣信息的活動(dòng)輪廓模型和基于區(qū)域的水平集模型.基于邊緣的活動(dòng)輪廓模型［3-6］依賴于圖像的邊緣信息，分割存在裂縫和缺口的圖像不會(huì)產(chǎn)生邊緣斷裂.但當(dāng)圖像的目標(biāo)與背景對(duì)比度不足，邊界模糊或噪聲較強(qiáng)時(shí)，這類方法會(huì)導(dǎo)致水平集演化抖動(dòng)或從弱邊界泄露.基于區(qū)域的水平集模型［7］將邊緣信息與全局同質(zhì)區(qū)域信息相結(jié)合，適用于邊緣模糊或邊緣不連續(xù)的情況，并對(duì)噪聲具有一定的抑制作用;但其分割灰度不均勻圖像時(shí)，會(huì)將與背景相似的目標(biāo)區(qū)域錯(cuò)分為背景.研究者們相繼提出各種改進(jìn)模型［8-10］，但這些模型存在能量函數(shù)依靠全局特征的不足.Li提出了一種基于圖像局部能量最小化的變分水平集方法［11］，這里簡(jiǎn)稱RSF (region-scalable fitting)算法.該方法根據(jù)圖像局部灰度信息提出了局部能量項(xiàng)，將其作為水平集的外部驅(qū)動(dòng)能量，更適用于分割灰度不均勻圖像.但是，受零水平集初始輪廓的影響，局部能量項(xiàng)會(huì)出現(xiàn)計(jì)算偏差，使得算法的穩(wěn)定性受到限制，并且算法的收斂速度和水平集的初始化方法有待進(jìn)一步提高.

本文將圖像的梯度信息和局部能量相結(jié)合，提出了一種基于局部能量和梯度敏感性的自動(dòng)初始化水平集算法.梯度敏感項(xiàng)依據(jù)圖像特征，自適應(yīng)地確定對(duì)水平集的驅(qū)動(dòng)方向.該能量函數(shù)能夠矯正由局部能量項(xiàng)計(jì)算偏差導(dǎo)致的演化錯(cuò)誤，提高算法的穩(wěn)定性，加快了水平集的演化.并且，采用圖像的區(qū)域極值初始化灰度擬合值，不需要交互式的初始化，提高了算法的效率.

1 RSF算法

RSF算法是一種基于區(qū)域的水平集分割方法.設(shè)定圖像域?yàn)棣浮蔙，已知灰度圖像I，閉合輪廓線C將圖像域劃分為外部區(qū)域Ω1和內(nèi)部區(qū)域Ω2.符號(hào)距離函數(shù)u初始化為:曲線C(零水平集)處u= 0，曲線C外部區(qū)域?yàn)棣?:u=c0，內(nèi)部區(qū)域?yàn)棣?: u=-c0，c0=2.RSF算法的能量泛函定義為

式中:L(u)是弧長(zhǎng)約束項(xiàng)，v、μ為權(quán)值;P(u)是自適應(yīng)的距離保持函數(shù)，避免了水平集演化過程中需重新初始化帶來的繁雜計(jì)算［6，11］;ξ(u，f1，f2)是水平集演化的主要驅(qū)動(dòng)能量，是圖像域中各點(diǎn)的局部能量項(xiàng)ξx(u，f1(x)，f2(x))在圖像域的積分.局部能量項(xiàng)表征了圖像域的局部灰度值與相應(yīng)的灰度擬合值f1(x)和f2(x)的近似度，其水平集函數(shù)定義為

式中，λi是權(quán)系數(shù)(通常取λ1=λ2)，fi(x)是點(diǎn)x相對(duì)區(qū)域Ωi的灰度擬合函數(shù)，y是以x為中心的局部圖像域中任意點(diǎn)的坐標(biāo)，I(y)表示點(diǎn)y的灰度值.局部圖像域的尺度為σ，由高斯核函數(shù)Kσ定義［11］.Hi(u)為Heaviside函數(shù)，且H1(u)= H(u)，H2(u)=1-H(u).

從式(1)的歐拉-拉格朗日方程解得使ξx(u，f1(x)，f2(x))最小化的灰度擬合函數(shù)fi(x)為

給定中心點(diǎn)x，在輪廓線C演化到實(shí)際目標(biāo)邊界且f1(x)和f2(x)是C兩側(cè)區(qū)域Ω1和Ω2的最優(yōu)逼近時(shí)，ξx(u，f1(x)，f2(x))取得最小值.

灰度擬合函數(shù)fi(x)是u的函數(shù)，因此輪廓線位置直接關(guān)系到fi(x)的計(jì)算，進(jìn)而影響ξx(u，f1(x)，f2(x))的取值，從而水平集的演化對(duì)初始輪廓敏感.初始輪廓選取的不適當(dāng)往往會(huì)導(dǎo)致水平集收斂所需的迭代次數(shù)增加，甚至最終不能收斂到目標(biāo)邊界.特別當(dāng)初始輪廓線離目標(biāo)邊界較遠(yuǎn)，或者分割多目標(biāo)圖像時(shí)，局部能量項(xiàng)的計(jì)算易出現(xiàn)偏差，使得圖像域中本應(yīng)屬于同一個(gè)區(qū)域(Ω1和Ω2)的點(diǎn)被錯(cuò)誤分割到不同區(qū)域.

另一方面，局部能量項(xiàng)在零水平集附近且圖像灰度不均勻的區(qū)域具有很強(qiáng)的驅(qū)動(dòng)力.然而，在距離零水平集較遠(yuǎn)的目標(biāo)邊界處f1(x)≈f2(x)，局部能量項(xiàng)的驅(qū)動(dòng)力較弱.同樣，在圖像灰度均勻的區(qū)域f1(x)≈f2(x)，局部能量項(xiàng)的驅(qū)動(dòng)力也很弱，水平集演化緩慢甚至?xí)Ｖ沟椒沁吘壧?可見，RSF方法的穩(wěn)定性和水平集的演化速度有待進(jìn)一步提高.

2 本文算法

2.1 融合局部能量和梯度敏感性的水平集模型

圖像梯度是描述圖像特征的一個(gè)重要信息，在曲線演化過程中，具有重要的指導(dǎo)作用.為了提高水平集分割方法的穩(wěn)定性以及水平集的演化速度，本文采用梯度信息構(gòu)造梯度敏感項(xiàng)，并結(jié)合局部能量項(xiàng)，提出基于圖像局部能量和梯度敏感性的水平集分割方法，定義水平集能量泛函如下:

式中:Φ(u，▽I)為梯度敏感的外部能量函數(shù); G(u，▽I)為梯度敏感的內(nèi)部能量函數(shù);γ和ρ為Φ(u，▽I)和G(u，▽I)的權(quán)值，γ較大時(shí)，Φ(u，▽I)將零水平集拉向目標(biāo)邊界的作用增強(qiáng);ρ較大時(shí)，G(u，▽I)將加快零水平集離開灰度均勻區(qū)域的速度.

采用梯度下降法最小化能量泛函E(u)，得到水平集演化的速度函數(shù)

式中:ei由ξx(u，f1(x)，f2(x))導(dǎo)出，

Heaviside函數(shù)及其導(dǎo)數(shù)分別定義為

式中，參數(shù)ε控制了圖像域內(nèi)各點(diǎn)的驅(qū)動(dòng)能量在速度函數(shù)中的權(quán)重，實(shí)驗(yàn)中選取ε=1.

基于圖像局部能量和梯度敏感性的水平集模型中，水平集的運(yùn)動(dòng)受到2種外力的支配，1種是來自圖像局部灰度能量的外力，1種是來自圖像梯度特性的外力.在圖像的不同區(qū)域，都有相應(yīng)的力驅(qū)動(dòng)水平集運(yùn)動(dòng).下面對(duì)梯度敏感項(xiàng)做詳細(xì)介紹.

2.2 梯度敏感的外部能量函數(shù)

為了矯正局部能量項(xiàng)導(dǎo)致的水平集演化錯(cuò)誤，加快零水平集向目標(biāo)邊界的運(yùn)動(dòng)速度，提出如下的梯度敏感的外部能量函數(shù):

式中:sgn(·)是符號(hào)函數(shù)，ΔI是由拉普拉斯算子定義的圖像的二階導(dǎo)數(shù).|▽I|是圖像梯度的模，φ(▽I)是歸一化的梯度響應(yīng)函數(shù)，b是梯度敏感因子.圖像灰度值在［0，255］區(qū)間時(shí)，b取值為［5，15］.

利用圖像二階導(dǎo)數(shù)在邊界低灰度值的一側(cè)為正，高灰度值一側(cè)為負(fù)的過零點(diǎn)性質(zhì)，定義能量函數(shù)Φ(u，▽I)的符號(hào)為sgn(ΔI).從而，梯度敏感的外部能量函數(shù)可以依據(jù)圖像的梯度信息自適應(yīng)地判定對(duì)水平集的驅(qū)動(dòng)方向，避免了人工設(shè)定符號(hào)對(duì)初始位置的依賴，以及由此帶來的應(yīng)用局限性和邊界泄露問題，如GAC模型中的面積項(xiàng)［5］.

Φ(u，▽I)的強(qiáng)度由φ(▽I)確定，而敏感因子b決定了φ(▽I)對(duì)不同梯度值的響應(yīng)度.圖1給出了當(dāng)max(|▽I|)=50時(shí)，不同b值對(duì)應(yīng)的響應(yīng)曲線.由圖可知，φ(▽I)只對(duì)大于一定強(qiáng)度的梯度值有響應(yīng)，對(duì)小于這個(gè)強(qiáng)度的梯度不響應(yīng)或者響應(yīng)度趨于零.隨著b的增大，梯度響應(yīng)閾值增高，響應(yīng)范圍縮小.調(diào)節(jié)φ(▽I)的響應(yīng)特性，能夠?qū)崿F(xiàn)不同強(qiáng)弱程度邊緣的提取.例如，降低b的取值，可以提取出圖像中的弱邊緣.

圖1 不同b值時(shí)φ(▽I)的梯度響應(yīng)曲線

可見，梯度敏感的外部能量函數(shù)的取值主要與圖像的梯度信息相關(guān)，在灰度不均勻區(qū)域具有顯著的驅(qū)動(dòng)力.符號(hào)距離函數(shù)u對(duì)Φ(u，▽I)的影響遠(yuǎn)小于對(duì)ξx(u，f1(x)，f2(x))的影響.因此，Φ(u，▽I)的引入能夠提高水平集分割算法的穩(wěn)定性，加速零水平集的收斂.

2.3 梯度敏感的內(nèi)部能量函數(shù)

為了提高水平集在圖像灰度均勻區(qū)域的演化速度，設(shè)計(jì)了梯度敏感的內(nèi)部能量函數(shù)

式中:k=div(▽u/|▽u|)是演化曲線的曲率; g(▽I)是歸一化的梯度響應(yīng)函數(shù);a是梯度敏感因子，圖像灰度值在［0，255］區(qū)間時(shí)，a取值為［1，10］.

梯度敏感的內(nèi)部能量函數(shù)G(u，▽I)對(duì)水平集的驅(qū)動(dòng)方向由sgn(k)決定.曲率為正的水平集曲線向內(nèi)部收縮，曲率為負(fù)的水平集曲線向外部擴(kuò)張.系數(shù)a控制著g(▽I)收斂到零的速度，如圖2所示.G(u，▽I)主要作用于灰度均勻區(qū)域，具有推動(dòng)零水平集離開平坦區(qū)域的作用，與梯度敏感項(xiàng)的外部能量項(xiàng)互為補(bǔ)充.同時(shí)該項(xiàng)還能平滑輪廓曲線，約束新輪廓的產(chǎn)生.

圖2 不同a值時(shí)g(▽I)的梯度響應(yīng)曲線

2.4 水平集自動(dòng)初始化方法

為進(jìn)一步提高算法的效率，本文根據(jù)局部能量模型的特點(diǎn)，提出利用區(qū)域極值初始化灰度擬合函數(shù)，實(shí)現(xiàn)水平集自動(dòng)初始化的方法.具體步驟為:首先對(duì)圖像的灰度分布做直方圖統(tǒng)計(jì)，并采用二次差分法繪出直方圖的包絡(luò)線.然后，對(duì)包絡(luò)線上的點(diǎn)再次求取二次差分得到一組極大值點(diǎn).將極值點(diǎn)排序后找出灰度值概率密度最大的3個(gè)點(diǎn)，用這3個(gè)點(diǎn)中灰度值最低的點(diǎn)初始化f1，灰度值最高的點(diǎn)初始化f2.最后，代入u=c0*sgn(-(e1-e2))，完成自動(dòng)初始化.

圖像的直方圖表征了圖像灰度的全局概率分布情況，因此包絡(luò)線的極大值點(diǎn)所對(duì)應(yīng)的灰度值代表著圖像的主要灰度分布區(qū)域的灰度中心.本文對(duì)f1和f2的初始化是圖像區(qū)域灰度中心估計(jì)與灰度擬合函數(shù)初始化的有機(jī)結(jié)合.可視為將局部域擴(kuò)展到整幅圖像的一種特例.自動(dòng)初始化完畢后，零水平集基本位于目標(biāo)邊界附近，這不但免去了交互式的手動(dòng)初始化，而且提高了算法的效率.

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

采用中心差分法迭代求解速度函數(shù)，時(shí)間步長(zhǎng)Δt=0.1.φ(▽I)和g(▽I)在迭代中可作為常量使用.求解速度函數(shù)前，需要先對(duì)灰度擬合函數(shù)f1和f2進(jìn)行更新.算法的計(jì)算量主要在于每次迭代時(shí)求解-δ(u)(e1-e2)、f1和f2的4次卷積運(yùn)算.本文算法的流程圖如圖3所示.

圖3 算法流程

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

為了驗(yàn)證算法的有效性，對(duì)本文算法和RSF算法進(jìn)行實(shí)驗(yàn)對(duì)比.除特殊說明，設(shè)定實(shí)驗(yàn)參數(shù)為:σ=3，v=0.003×255×255，a=10，b=5，γ=0.03×255×255，ρ=0.006×255×255.

實(shí)驗(yàn)1 采用多目標(biāo)的米粒圖像驗(yàn)證算法的穩(wěn)定性.手動(dòng)初始水平集如圖4(a)，本文算法迭代70次后準(zhǔn)確分割出每個(gè)米粒，結(jié)果如圖4(b)所示.RSF算法在迭代290次后，水平集收斂如圖4(c)，無法分割出初始輪廓附近的米粒.改變水平集初始輪廓如圖4(d)所示，本文算法和RSF算法的分割結(jié)果見圖4(e)和4(f).可以看出，本文方法對(duì)水平集初始輪廓的敏感性降低，算法的穩(wěn)定性得到提高.

圖4 米粒圖像分割結(jié)果對(duì)比

實(shí)驗(yàn)2 采用大腦MRI圖像驗(yàn)證算法的準(zhǔn)確性.大腦MRI圖像結(jié)構(gòu)復(fù)雜，含有內(nèi)、外多層輪廓.采用圖5(a)所示的自動(dòng)初始化零水平集，本文算法和RSF算法的分割結(jié)果分別見圖5(b)和圖5(c).由圖可知，本文算法和RSF算法都能分割出大腦的內(nèi)外層輪廓，而本文算法對(duì)弱邊緣和細(xì)節(jié)的提取更準(zhǔn)確，具體可參見2個(gè)ROI區(qū)域的局部放大圖5(d)～(g).實(shí)驗(yàn)選取a=2、b=5，此時(shí)梯度敏感的外部能量函數(shù)對(duì)弱邊緣有響應(yīng)，因此算法能夠提取出弱邊緣.可見，本文算法可以根據(jù)分割需要調(diào)節(jié)梯度敏感因子，實(shí)現(xiàn)對(duì)圖像不同強(qiáng)弱邊緣的提取.

圖5 大腦MRI圖像分割結(jié)果

實(shí)驗(yàn)3 對(duì)算法采用不同初始化方法所需的收斂時(shí)間與RSF算法進(jìn)行了對(duì)比，并給出一個(gè)自動(dòng)初始化實(shí)例.圖6(a)為米粒圖像的直方圖，從圖中包絡(luò)線的極值點(diǎn)(由實(shí)心圓點(diǎn)標(biāo)注)中選出3個(gè)極大值點(diǎn)，即為左起的第2、3、5點(diǎn).其中，灰度值最低點(diǎn)為(19，616)，最高點(diǎn)為(232，495).初始化f1=19，f2=232，相應(yīng)的水平集自動(dòng)初始化如圖6(b)所示.自動(dòng)初始化后，零水平集曲線基本位于目標(biāo)邊界的附近.在此初始化基礎(chǔ)上，水平集能夠快速收斂.

圖6 米粒圖像水平集的自動(dòng)初始化

表1給出了針對(duì)4幅不同的圖像，本文算法與RSF算法在分別采用水平集手動(dòng)初始化和自動(dòng)初始化的收斂時(shí)間.2種算法使用同樣的初始輪廓，算法收斂所需的CPU時(shí)間和圖像尺寸如表1所示.程序運(yùn)行環(huán)境為Intel Pentium Dual E2140 1.60 GHz，內(nèi)存2.00 GB，操作系統(tǒng)為Windows XP，應(yīng)用Matlab2006a編程.由表可知，無論采用哪種初始化方法，本文算法的收斂時(shí)間都小于RSF算法的1/3倍.而采用自動(dòng)初始化方法，算法的效率較手動(dòng)初始化也得到顯著的提高.

表1 本文算法與RSF算法的收斂時(shí)間對(duì)比 s

4 結(jié)論

1)本文融合了圖像的局部灰度信息和梯度信息，提出了一種基于圖像局部能量和梯度敏感性的自動(dòng)初始化水平集分割算法，可以實(shí)現(xiàn)灰度不均勻圖像、多目標(biāo)圖像、結(jié)構(gòu)復(fù)雜的多層輪廓圖像的分割.

2)局部能量項(xiàng)強(qiáng)調(diào)圖像的局部信息，是對(duì)傳統(tǒng)C-V模型全局優(yōu)化思想的改進(jìn).梯度敏感的外部能量函數(shù)的提出有助于提高水平集演化的穩(wěn)定性.

3)通過調(diào)節(jié)梯度敏感因子，可以控制算法對(duì)不同強(qiáng)弱邊緣的響應(yīng)度，準(zhǔn)確提取出弱邊緣.通過初始化灰度擬合值，實(shí)現(xiàn)了水平集的自動(dòng)初始化，進(jìn)一步提高了算法的效率.

4)與RSF算法相比，本文算法穩(wěn)定性好、收斂速度快，分割結(jié)果更準(zhǔn)確，且不需要交互操作.

5)理論上，本文算法可以推廣到多相水平集，相應(yīng)的多相水平集模型及自動(dòng)初始化方法有待進(jìn)一步的研究.

［1］CHENOUNE Y，DELECHELLE E，PETIT E，et al. Segmentation of cardiac cine-MR images and myocardial deformation assessment using level set methods［J］. Computerized Medical Imaging and Graphics，2005，29 (8):607-616.

［2］DRAPACA CS，CARDENAS V，STUDHOLME C.Segmentation of tissue boundary evolution from brain MR image sequences using multi-phase level sets［J］. Computer Vision and Image Understanding，2005，100 (3):312-329.

［3］XU C，PRINCE J L.Snakes，shapes，and gradient vector flow［J］.IEEE Trans on Image Processing，1998，7(3):359-369.

［4］CASELLES V，CATTE F，COLL T，et al.A geometric model for active contours in image processing［J］.Numerische Mathematik，1993，66(1):1-31.

［5］CASELLES V，KIMMEL R，SAPIRO G.Geodesic active contours［J］.International Journal of Computer Vision，1997，22(1):61-79.

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

［7］CHAN T，VESE L.Active contours without edges［J］. IEEE Trans on Image Processing，2001，10(2):266-277.

［8］原達(dá)，張彩明，李晉江，等.基于Mumford-Shah模型的高精度MR圖像輪廓提取算法［J］.計(jì)算機(jī)學(xué)報(bào)，2009，32(2):268-274.

［9］VESE L，CHAN T.A multiphase level set framework for image segmentation using the Mumford and Shah model［J］.International Journal of Computer Vision，2002，50(3):271-293.

［10］CHEN S，SOCHEN N A，ZEEVI Y.Integrated active contours for texture segmentation［J］.IEEE Trans on Image Processing，2006，15(6):1633-1646.

［11］LI C，KAO C，GORE J，et al.Minimization of regionscalable fitting energy for image segmentation［J］. IEEE Trans on Image Processing，2008，17(10): 1940-1949.