基于幾何主動輪廓模型的圖像分割方法

張榮華，吳 益，王 亨，肖 磊，唐勁天

（1.天津工業(yè)大學 電氣工程與自動化學院，天津 300387；2.清華大學工程物理系，北京 100084；3.清華大學 粒子技術(shù)與輻射成像教育部重點實驗室，北京 100084）

變分水平集是幾何主動輪廓模型的重要代表，自O(shè)sher等[1]在研究速度與曲率相關(guān)的曲線演化時提出以來，因其諸多優(yōu)點被廣大學者研究，在圖像處理領(lǐng)域尤其是醫(yī)學圖像分割方面被廣泛改進和運用.該方法基本思想是用高維曲面來映射低維曲線，曲線輪廓演化的方程式被隱含表示，所以水平集函數(shù)的演化方程能夠從高階閉合曲面或曲線的演化方程中推導出來，內(nèi)嵌的演化曲面或曲線就是其零水平集.因此，只需要設(shè)計零水平集，然后就可以計算出演化曲面或曲線的結(jié)果[2-4].該方法的最大優(yōu)點就是不依賴主動輪廓的參數(shù)化，而且能夠合理地處理拓撲結(jié)構(gòu)發(fā)生變化的曲線.

傳統(tǒng)水平集方法必須定期地初始化水平集函數(shù)作為符號距離函數(shù)，以此來縮小迭代過程帶來的計算錯誤，Li等[5]引入一個內(nèi)部能量項到能量泛函來糾正水平集函數(shù)與符號距離函數(shù)的偏差，提出一種無需初始化的演化方法.在此基礎(chǔ)上，Xu[6]又加入新的符號距離函數(shù)正則項來糾正偏差，提出距離正則水平集.近些年，許多其他有價值的分割模型也結(jié)合變分水平集相繼被提出.CV模型將目標邊緣終止函數(shù)消除掉，是一種無需邊界的主動輪廓模型[7].Gao等[8]用局部魯棒統(tǒng)計（LRSS）來對抗圖像噪聲，并利用稀疏場算法讓能量項的計算更加高效.Brox等[9]針對灰度不均勻分割圖像提出分段常值模型（PS）.Li等[10-11]提出局部二值方法（LBF），該模型將圖像的局部信息嵌入其中，能夠分割出灰度不均勻的血管圖像.

上述方法聚焦于水平集演化方向和速度的控制，忽略了終止函數(shù)對分割最終效果的影響，不能在尖角邊界終止演化，對初始演化區(qū)域敏感并且不能正確識別多個目標.本文構(gòu)造了新的終止函數(shù)和自適應性選擇系數(shù)，提出一種新的邊緣主動輪廓模型的自適應算法：一方面，圖像信息用來決定水平集函數(shù)的演化方向，因而演化結(jié)果不依賴于初始輪廓的位置；另一方面，外部能量項與終止速度函數(shù)相互獨立，收斂速度也得到提高.

1 水平集原理

考慮曲線在有界區(qū)域R2內(nèi)跟隨時間變化，設(shè)定C（x，y，t），x，y∈R2是一條閉合曲線；水平集函數(shù)設(shè)為U（x，y，t）∶R2× [0，t]→R，將曲線嵌入讓水平集函數(shù)變?yōu)殡[函數(shù)Uc（x，y，t）=d，d∈R為常數(shù).基于曲線演化理論有：

式中：F為速度函數(shù)，決定曲線C上每個點的演化速度；N為曲線的單位內(nèi)法向量，具體表示為N=-ΔU/|ΔU|，Δ為梯度算子.

結(jié)合（1）式對隱函數(shù)求其關(guān)于t的微分有如下結(jié)果：

因此 C（x，y，t）的演變轉(zhuǎn)變?yōu)樗郊瘮?shù)的偏微分方程可表示為：

圖1所示為水平集演化思想的示意圖.圖1中：曲面代表水平集函數(shù)U（x，y，t）；黑色加粗曲線表示演變的輪廓曲線C.

圖1 水平集思想示意圖Fig.1 Illustration of level set principle

曲線的演變過程其實是曲線上不同力的作用過程，也即能量的作用.分割目標區(qū)域趨于能量最小而存在，所以，建立能量模型然后按照一定的約束條件求能量最小，即可分割出目標輪廓[12].水平集方法用于圖像分割的幾何主動輪廓模型的一般過程可以概括為：

Step1：圖像輸入與預處理；

Step2：建立能量模型后利用變分法取得Euler-Lagrange方程；

Step3：利用水平集方法獲取曲線演化的偏微分方程；

Step4：設(shè)定初始水平集后選取PDE數(shù)值算法進行迭代運算；

Step5：利用終止函數(shù)判斷是否滿足停止演化條件，不滿足條件時繼續(xù)迭代；

Step6：演化終止，輸出水平集函數(shù)，得到最終的分割曲線.

2 本文模型

2.1 邊緣終止函數(shù)

在圖像分割時，主動輪廓模型可被分為基于邊界模型[3，10-12]和基于區(qū)域的模型[4，13-15]兩大類.基于邊界的水平集利用圖像的梯度來構(gòu)建邊緣檢測函數(shù)，以此來終止曲線演化，勾畫出目標邊界.通常，速度終止函數(shù)有如下兩類選取方法：

式中表示對圖像進行高斯濾波預處理.圖2反映了多種速度終止函數(shù)的特征.

圖2 多種速度終止函數(shù)特征曲線Fig.2 Variety function characteristic curves of velocity termination

圖2中，L=1時，g（s1）在A點以后對梯度變化很敏感，快速趨近于0，g（s2）則相對平緩.L=2時，g（s1）和g（s2）在s值較小區(qū)域下降緩慢，這將影響曲線收斂的時間，降低分割的速度.對于L的選取，方法不一，如將梯度取絕對平均值、將絕對偏差取中值、引入Canny噪聲估計[13]等，但所有方法都是以分離目標邊緣和噪聲作為出發(fā)點.當分割目標處于低信噪比圖像中時，很難區(qū)分梯度的變化是由噪聲還是邊緣信息產(chǎn)生，因此難以權(quán)衡噪聲的平滑和弱邊緣信息的保留.綜合以上分析，一般情況下：

由式（4）可知，函數(shù)g在演化過程中依賴圖像信息，梯度的模值起決定作用，在目標邊界時g接近0.弱邊界時，此邊緣檢測函數(shù)往往不連續(xù)并且不完全等于0，因此不能充分地終止水平集演化.再者，g始終為正值，曲線朝單一方向收斂，目標點一旦泄露便不能返回.

因此，本文提出L函數(shù)化的方案，將其定義為：

式中：I（x，y）為圖像灰度信息；m為分割前輸入?yún)?shù)；h1、h2為2個卷積函數(shù)的商，分別近似于圖像在高斯窗內(nèi)和在高斯窗外的平均灰度值.具體定義為：

式中：F 是修改的 Heaviside 函數(shù)[11，14]；Kσ為高斯核，核函數(shù)中的尺度參數(shù)σ＞0，兩者的函數(shù)表達式為：

因此得出如下新的邊緣終止函數(shù)：

當水平集演化到目標邊界時，式（5）中分子分母均接近于1，故L（I，h1，h2）≈1.此時，|（ΔGσ·I）|很大，由式（10）可知，g函數(shù)值近似零，能量模型陷入極小；當目標點被水平集包含但離邊界距離比較遠時，式（5）中分子項遠大于 1，而分母項近似 1，所以 L（I，h1，h2）＞＞1，g（I，L）≈1，模型能量未達極小，水平集在非目標邊界處繼續(xù)演化.

2.2 可變權(quán)重系數(shù)

可變權(quán)重系數(shù)ω（I）能讓算法根據(jù)圖像信息自動的改變演化方向，所以零水平集自適應演化，克服了對初始位置敏感的缺陷.可變權(quán)重系數(shù)ω（I）定義為：

式中：I（x，y）為圖像的灰度值；I（x，y）為反向灰度值；sgn（·）象征一個符號函數(shù)；α為大于零的常數(shù).

主動輪廓在目標區(qū)域外時，I（x，y）-I（x，y）＞ 0，可變權(quán)重系數(shù)引導水平集向目標區(qū)域內(nèi)靠近，相反，輪廓位于區(qū)域內(nèi)時，I（x，y）-I（x，y）＜ 0，ω（I）則引導水平集朝區(qū)域外演化.換句話說，零水平集能自適應的決定演化方向，有利于擺脫初始位置的干擾.與此同時，α的不同取值也能在復雜圖像分割時，對不同輪廓層起識別作用.

2.3 數(shù)值描述與算法實現(xiàn)

基于以上討論，本文模型的能量泛函為：

式中：β＞0為常數(shù)；Mq（U）用來修正水平集函數(shù)的內(nèi)部結(jié)構(gòu)，定量的描述水平集函數(shù)與符號距離函數(shù)的偏差程度，定義如式（13）.其中的能量密度函數(shù) q（s）按照（14）式分段選取，為保證演化過程中符號距離特性只在零水平集附近成立，q（s）在s=1和s=0時最小，避免水平集受困于局部最優(yōu)解，水平集演化的穩(wěn)定性和有效性得到保證.

水平集定義在區(qū)域Ω內(nèi)，初始水平集U0（x，y）設(shè)置為分段常值函數(shù)[17]，若點在區(qū)域內(nèi)等于-2，否則設(shè)為2.能量函數(shù)中Rg（U）的用來計算輪廓的曲線積分，Ag（U）表示區(qū)域Ω-U={x ∶U（x）＜ 0}的加權(quán)面積.

水平集方法實現(xiàn)過程中，平滑函數(shù)δ（x）和H（x）分別用（17）式和（18）式中的Dirac函數(shù)和Heaviside函數(shù)[15]來代替，式中H′（x）=δ（x）.

運用變分方法，E（U）的最小值可以由如下的梯度流方程式[16]算出：

當t=0時，U（x，y，0）表示初始水平集U0（x，y），?E（U）/?U代表求E（U）的導數(shù).結(jié)合上述推導，偏微分方程式轉(zhuǎn)化為：

上述偏微分方程綜合諾伊曼邊界條件和初始輪廓曲線得出最終結(jié)果，式（20）的離散方式采用有限中心差分方法.

3 實驗與應用

為驗證算法的準確性和有效性，本文選取各具特征的圖像進行試驗，試驗平臺為i5-2320CPU，64bit win7企業(yè)版操作系統(tǒng)的計算機，編程環(huán)境為R2013a版本的MATLAB.

對比Li等[5-6]的算法，本文描述的幾何主動輪廓水平集方法有其優(yōu)點.首先，演化方向可以根據(jù)圖像信息自行決定，不再是單向的演化；然后，帶尖角的目標區(qū)域、多目標對象（灰度值可不同）均能被準確識別；再者，也能分割帶內(nèi)腔的多目標并且計算時間短.

不同初始區(qū)域分割多目標如圖3所示.為證明算法的自適應演化能力，圖3中初始區(qū)域和目標輪廓分別內(nèi)部包含、相互交叉和互不相交.選用不同的初始位置對多目標圖像進行分割實驗，該算法不僅能準確、快速地識別出所有目標，而且自適應的決定演化方向.然而，用Li等[5-6]提出的方法即使迭代580次（迭代次數(shù)已超過本文算法的3倍），變分水平集仍不能正確的收斂，出現(xiàn)漏分割和邊界識別不準確等分割結(jié)果不理想的現(xiàn)象.另外，它的分割方法太依賴初始輪廓所在的位置來人工調(diào)整曲線的演化方向.本文提出的算法用更少的迭代時間得到了有效的分割結(jié)果，證明了算法的有效性.

選用一幅人造圖像作為實驗對象，該圖像由3個灰度值不同的圖形組成，其中有一個凹形和帶封閉內(nèi)腔的圖形.在初始輪廓位置相同的情況下試驗，結(jié)果如圖4所示.從左至右依次為：初始區(qū)域、Li[5]的方法分割結(jié)果、本文分割算法結(jié)果.提出的算法不僅快速識別出外部邊界和凹形目標區(qū)域，還分割出圖形的內(nèi)部腔體邊界，但是Li[5]的方法分割效果不佳.

圖3 不同初始區(qū)域分割多目標Fig.3 Multiple targets segmentation with different initial regions

圖4 分割復雜的人工圖像Fig.4 Segmentation of complex artificial images

尖角邊界識別如圖5所示.為了證明算法識別尖角邊界的能力，圖5中用來進行對比試驗的樹葉輪廓有多個尖角，仍然是相同的初始輪廓下，本文算法準確、快速地分割出目標區(qū)域的所有尖角.然而，用Li等[5]提出的方法進行試驗，即使580次迭代后，尖角依然很難被識別，而且改變初始位置或者增加迭代次數(shù)繼續(xù)分割也沒有好轉(zhuǎn)的趨勢.

圖5 尖角邊界識別Fig.5 Recognition of sharp angle boundary

CT肝臟腫瘤識別如圖6所示.為進一步驗證算法的可行性，圖6為將算法應用到CT圖像識別人體肝臟腫瘤位置的實驗結(jié)果圖.由圖6可知，提出的算法能夠有效地識別出腫瘤位置，可用于醫(yī)學圖像的組織器官分割[17-18].

圖6 CT肝臟腫瘤識別Fig.6 Recognition of CT liver tumor

算法開始前的輸入?yún)?shù)β=1.5，λ=1.5跟Li等[5]的方法保持一致，另外2個關(guān)鍵輸入?yún)?shù)m和α由實驗得出，實驗發(fā)現(xiàn)m在范圍[1，2.5]取值效果較好.表1給出了以上實驗圖像的信息（按出現(xiàn)順序標號從小到大）、分割參數(shù)的選擇以及各實驗的迭代次數(shù)和用時，不難看出本文提出算法的分割效率較高.

表1 本文算法分割參數(shù)及效率Tab.1 Segmentation parameters and efficiency of this algorithm

分割評價用 Jaccard（J）系數(shù)、Dice（D）系數(shù)進行綜合評定.J系數(shù)和D系數(shù)是分割中基于相似度評測的常用量化標準，其值越接近1表明結(jié)果與標準相似度越高，兩大系數(shù)分別定義如下：

式中：S1為標準圖像輪廓的像素集合；S2為分割算法結(jié)果圖的像素集合.

本文算法分割結(jié)果的量化評估如表2所示.其中圖像No.1（即圖3實驗圖像）對應的3次分割結(jié)果的評估參數(shù)均近似0.95.同時，用Li[5]的算法進行評估參數(shù)對比，圖3中不同初始區(qū)域下，Li[5]等的方法J系數(shù)分 別 為 ：0.89、0.73、0.56；D系數(shù)分別為：0.90、0.76、0.61.圖4和圖5的對比試驗中，其方法的評估系數(shù)同樣低于本文方法，具體為：圖4中J系數(shù)為0.72，D系數(shù)為0.75；圖5中J系數(shù)為0.83，D系數(shù)為0.84.

表2 本文算法分割系數(shù)Tab.2 Segmentation coefficient of this algorithm

4 結(jié)論

為克服傳統(tǒng)水平集方法的缺陷，促進水平集分割方法在醫(yī)學圖像上的應用，本文結(jié)合幾何主動輪廓模型提出了一種新的改進方法.對比傳統(tǒng)和現(xiàn)行的算法，該算法具有一定的優(yōu)勢.

（1）算法可以自我調(diào)整演化方向，讓分割過程不依賴曲線所在的初始位置.

（2）算法能準確分割出多目標、尖角弱邊界和凹形、帶內(nèi)腔圖案且分割結(jié)果的J系數(shù)在0.91以上，D系數(shù)不低于0.95.最后，該算法在更為復雜的人體器官CT圖像分割中也表現(xiàn)良好，能較好的分割出目標腫瘤輪廓.算法的準確性、快速性和可靠性得到驗證.

[1]OSHER Stanley，SETHIAN James A.Fronts propagating with curvature-dependent speed：Algorithms based on Hamilton-Jacobi formulations[J].Journal of Computational Physics，1988，79（1）：12-49.

[2]KIM W，KIM C.Active contours driven by the salient edge energy model[J].IEEE Transactions on Image Processing A Publication of the IEEE Signal Processing Society，2013，22（4）：1667-1673.

[3]呂曉琪，石靜，任曉穎，等．一種基于水平集的三維肝臟磁共振圖像混合分割方法[J]．數(shù)據(jù)采集與處理，2015（2）：350-358.LYU X Q，SHI J，REN X Y，et al.A 3D liver hybrid magnetic resonance image segmentation method based on level set[J].Data Acquisition and Processing，2015（2）：350-358（in Chinese）.

[4]鄭倩.醫(yī)學圖像分割方法研究及其應用[D].廣州：南方醫(yī)科大學，2014.ZHENG Q.Research and application of medical image segmentation method[D].Guangzhou：Southern Medical University，2014（in Chinese）.

[5]LI C，XU C，GUI C，et al.Level set evolution without re-ini－tialization：A new variational formulation[C]//Computer Vision and Pattern Recognition，2005.[s.l.]：IEEE Xplore，2005：430-436.

[6]XU C.Distance regularized level set evolution and its application to image segmentation[J].IEEE Transactions on Image Processing，2010，19（12）：3243.

[7]CHAN T F，VESE L A.Active contours without edges[J].IEEE Press，2001，10（2）：266-277.

[8]GAO Y，TANNENBAUM A，KIKINIS R.Simultaneous multiobject segmentation using local robust statistics and contour interaction[C]//Recognition Techni-ques and Applicationsin Medical Imaging.Berlin：Springer Berlin Heidelberg，2010：195-203.

[9]BROX T，CREMERS D.On the statistical interpretation of the piecewise smooth mumford-shah functional[C]//International Conference on Scale Space and Variational Methods in Computer Vision.Berlin：Springer-Verlag，2007：203-213.

[10]LI C，KAO C Y，GORE J C，et al.Implicit active contours driven by local binary fitting energy[C]//Computer Vision and Pattern Recognition，2007.[s.l.]：IEEE，2007：1-7.

[11]王茜，彭中，劉莉.一種基于自適應閾值的圖像分割算法[J].北京理工大學學報，2003，23（4）：521-524.WANG Q，PENG Z，LIU L.An image segmentation algorithm based on adaptive threshold[J].Journal of Beijing Institute of Technology，2003，23（4）：521-524（in Chinese）.

[12]ZHANG K，SONG H，ZHANG L.Active contours driven by local image fitting energy[J].Pattern Recognition，2010，43（4）：1199-1206.

[13]曹武.心臟CT序列圖像分割算法研究[D].北京：北京交通大學，2016.CAO W.Study of segmentation algorithms on cardiac CT sequence images[D].Beijing：Beijing Jiaotong University，2016（in Chinese）.

[14]ZHANG K，ZHANG L，ZHANG S.A variational multiphase level set approach to simultaneous segmentation and bias correction[J].IEEE International Conference on Image Processing，2010，119（5）：4105-4108.

[15]何小虎.圖像分割算法研究[J].福建電腦，2014，30（8）：3-4.HE X H.Research on image segmentation algorithm[J].Fujian Computer，2014，30（8）：3-4（in Chinese）.

[16]LU S，HUANG H，LIANG P，et al.Hepatic vessel segmentation using variational level set combined with non-local robust statistics[J].MagneticResonanceImaging，2017，36：180-186.

[17]ZHENG Q，DONG E，CAO Z，et al.Active contour model driven by linear speed function for local segmentation with robust initialization and applications in MR brain images[J].Signal Processing，2014，97：117-133.

[18]HUANG J，JIAN F，WU H，et al.An improved level set method for vertebra CT image segmentation[J].Biomedical Engineering Online，2013，12（1）：48.