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

包立君，劉宛予，浦昭邦

(哈爾濱工業大學HIT-INSA中法生物醫學圖像聯合研究中心，150001哈爾濱，baolij@gmail.com)

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

本文將圖像的梯度信息和局部能量相結合，提出了一種基于局部能量和梯度敏感性的自動初始化水平集算法.梯度敏感項依據圖像特征，自適應地確定對水平集的驅動方向.該能量函數能夠矯正由局部能量項計算偏差導致的演化錯誤，提高算法的穩定性，加快了水平集的演化.并且，采用圖像的區域極值初始化灰度擬合值，不需要交互式的初始化，提高了算法的效率.

1 RSF算法

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

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

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

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

給定中心點x，在輪廓線C演化到實際目標邊界且f1(x)和f2(x)是C兩側區域Ω1和Ω2的最優逼近時，ξx(u，f1(x)，f2(x))取得最小值.

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

另一方面，局部能量項在零水平集附近且圖像灰度不均勻的區域具有很強的驅動力.然而，在距離零水平集較遠的目標邊界處f1(x)≈f2(x)，局部能量項的驅動力較弱.同樣，在圖像灰度均勻的區域f1(x)≈f2(x)，局部能量項的驅動力也很弱，水平集演化緩慢甚至會停止到非邊緣處.可見，RSF方法的穩定性和水平集的演化速度有待進一步提高.

2 本文算法

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

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

式中:Φ(u，▽I)為梯度敏感的外部能量函數; G(u，▽I)為梯度敏感的內部能量函數;γ和ρ為Φ(u，▽I)和G(u，▽I)的權值，γ較大時，Φ(u，▽I)將零水平集拉向目標邊界的作用增強;ρ較大時，G(u，▽I)將加快零水平集離開灰度均勻區域的速度.

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

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

Heaviside函數及其導數分別定義為

式中，參數ε控制了圖像域內各點的驅動能量在速度函數中的權重，實驗中選取ε=1.

基于圖像局部能量和梯度敏感性的水平集模型中，水平集的運動受到2種外力的支配，1種是來自圖像局部灰度能量的外力，1種是來自圖像梯度特性的外力.在圖像的不同區域，都有相應的力驅動水平集運動.下面對梯度敏感項做詳細介紹.

2.2 梯度敏感的外部能量函數

為了矯正局部能量項導致的水平集演化錯誤，加快零水平集向目標邊界的運動速度，提出如下的梯度敏感的外部能量函數:

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

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

Φ(u，▽I)的強度由φ(▽I)確定，而敏感因子b決定了φ(▽I)對不同梯度值的響應度.圖1給出了當max(|▽I|)=50時，不同b值對應的響應曲線.由圖可知，φ(▽I)只對大于一定強度的梯度值有響應，對小于這個強度的梯度不響應或者響應度趨于零.隨著b的增大，梯度響應閾值增高，響應范圍縮小.調節φ(▽I)的響應特性，能夠實現不同強弱程度邊緣的提取.例如，降低b的取值，可以提取出圖像中的弱邊緣.

圖1 不同b值時φ(▽I)的梯度響應曲線

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

2.3 梯度敏感的內部能量函數

為了提高水平集在圖像灰度均勻區域的演化速度，設計了梯度敏感的內部能量函數

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

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

圖2 不同a值時g(▽I)的梯度響應曲線

2.4 水平集自動初始化方法

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

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

2.5 算法實現

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

圖3 算法流程

3 實驗結果

為了驗證算法的有效性，對本文算法和RSF算法進行實驗對比.除特殊說明，設定實驗參數為:σ=3，v=0.003×255×255，a=10，b=5，γ=0.03×255×255，ρ=0.006×255×255.

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

圖4 米粒圖像分割結果對比

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

圖5 大腦MRI圖像分割結果

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

圖6 米粒圖像水平集的自動初始化

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

表1 本文算法與RSF算法的收斂時間對比 s

4 結論

1)本文融合了圖像的局部灰度信息和梯度信息，提出了一種基于圖像局部能量和梯度敏感性的自動初始化水平集分割算法，可以實現灰度不均勻圖像、多目標圖像、結構復雜的多層輪廓圖像的分割.

2)局部能量項強調圖像的局部信息，是對傳統C-V模型全局優化思想的改進.梯度敏感的外部能量函數的提出有助于提高水平集演化的穩定性.

3)通過調節梯度敏感因子，可以控制算法對不同強弱邊緣的響應度，準確提取出弱邊緣.通過初始化灰度擬合值，實現了水平集的自動初始化，進一步提高了算法的效率.

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

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

［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］何傳江，李夢，詹毅.用于圖像分割的自適應距離保持水平集演化［J］.軟件學報，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］原達，張彩明，李晉江，等.基于Mumford-Shah模型的高精度MR圖像輪廓提取算法［J］.計算機學報，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.