，，， ，

(1.中原工學(xué)院，鄭州 450007；2.河南省人民醫(yī)院，鄭州 450003；3.劍橋大學(xué)，Cambridge CB2 1TN)

圖像分割在醫(yī)學(xué)圖像分析中占據(jù)至關(guān)重要的位置，準確的圖像分割可以為醫(yī)學(xué)圖像分析、模式識別以及機器視覺等高層應(yīng)用奠定良好的基礎(chǔ)。

在眾多的圖像分析方法中，主動輪廓模型(Active Contour Model，也被稱為"蛇"模型) 使用曲線或者曲面，在外力場和內(nèi)力場的作用下使之進行連續(xù)的形變，從而達到捕捉物體輪廓線的目的。由于該模型具有針對形狀變化自然處理的能力，因此被廣泛應(yīng)用于物體形狀的提取[1-3]。

主動輪廓模型可以分為顯式模型(explicit model)和隱式模型(implicit model)兩種。顯式模型使用參數(shù)方程對輪廓進行連續(xù)的形變。顯式模型利用輪廓上的點在時間上進行連續(xù)變化，計算量小，適用于時間敏感的實時輪廓提取任務(wù)。隱式模型通常是基于水平集(level set)方法[4]。這類方法將輪廓模型表示為一個高維標量函數(shù)的水平集，輪廓線的形變是基于曲率或者法向量等幾個量的變化，使輪廓線的演變獨立于參數(shù)方程。

本文主要研究顯式主動輪廓模型在改進的外力場模型下輪廓演變的計算。顯式主動輪廓模型由于其計算復(fù)雜度小而被廣泛應(yīng)用在視覺和圖像分析任務(wù)中。和其他特征提取技術(shù)相比較，在某種程度上其克服了傳統(tǒng)技術(shù)的局限性。圖像的初始輪廓估計、目標輪廓特征以及各種約束條件都集成在一個特征提取的過程中。經(jīng)過初始輪廓的估計后，顯式的輪廓線能夠自動收斂到能量的極小狀態(tài)，最終運動到特征邊緣附近，從而完成圖像的分割。輪廓線的形變受到內(nèi)力和外力的控制，其中內(nèi)力用于影響輪廓線的形狀，外力控制輪廓線的位移。近年來，研究者提出了各種新方法，對輪廓線內(nèi)力和外力的具體形式進行了不斷改進[1-2,5-11]。

本文針對腦部MRI圖像的目標區(qū)域進行輪廓提取。首先對腦部MRI圖像的二維直方圖進行處理以提取閾值，對圖像進行二值處理后，求得二值圖像的梯度邊緣。經(jīng)過上述處理后，基于圖像像素的外力場建模的復(fù)雜度大大降低。本文使用靜電力模型對預(yù)處理后的目標圖像計算初始輪廓線所受的力，輪廓線所受的力使目標圖像不斷收縮和位移，最終達到收斂的狀態(tài)，完成分割MRI 圖像。

1 主動輪廓模型

1.1 基本模型

顯式動態(tài)輪廓線的參數(shù)方程v(s)=[x(s),y(s)]以歸一化的弧長s為參數(shù)，其中s∈[0,1]。輪廓線的能量定義如下：

(1)

其中Eint(s)、Eext(s)、Econt(s)分別為輪廓線的內(nèi)能強度、外能強度以及人為交互力的強度。如果忽略外部干預(yù)的話，則可以令Econt(s)=0。

在Econt(s)=0的條件下，輪廓線的內(nèi)能強度Eint(s)由下面公式給出：

(2)

式(2)中，Eint(s)表達式的第一項為彈性勢能，第二項為剛性能量。α(s)和β(s)為相應(yīng)的控制系數(shù)，當(dāng)兩參數(shù)分別為零時，允許輪廓線上的點(也稱為蛇點)在此處的彈性能量和剛性能量取一個極大值，即允許輪廓線出現(xiàn)不連續(xù)或曲率較大[11]。

1.2 圖像的外力

圖像的外力一般可以?。?/p>

Eext(s)=-r(s)|(Gσ(x,y)*I(x,y))|2

(3)

其中I(x,y)為圖像在(x,y)位置處的灰度，Gσ(x,y)為二維高斯函數(shù)，r(s)為加權(quán)系數(shù)。

在二維高斯函數(shù)取極值的情況下簡化外力為：

(4)

取梯度的模值，可以將式(3)簡化為下面的形式：

Eext(s)=-r(s)|I(x,y)|

(5)

利用主動輪廓模型進行圖像分割的目標就是使式(1)最小化。設(shè)輪廓線可以由一系列離散點構(gòu)成的集合來表示，V={v1,…,vn}，其中vi=[xi,yi]，離散化的輪廓線能量可以寫成：

(6)

具有最小能量的輪廓線可以表示為：

(7)

離散化的內(nèi)能項可以寫成：

(8)

對于目標函數(shù)(7)的最小化，可以使用Amini等提出的動態(tài)規(guī)劃法[1]，或者Williams等提出的貪婪法[5]，其中貪婪法收斂速度較快。如果設(shè)當(dāng)前蛇點為vi，其鄰域點集為Pi，使用貪婪法移動一次Snake輪廓線的算法如下：

Algorithm 1:Moving Snake1:Procedure Moving Snake(v1,…,vn)2: for all vi∈{v1,…,vn} do3: Emin←Esnake(vi);4: for all p∈Pi do5: if Emin>Esnake(p) then6: Emin←Esnake(p);7: vi←p8: end if9: end for10: end for11:end procedure

2 Snake算法的改進

2.1 蛇點移動鄰域

(9)

如圖1所示，對目標蛇點鄰域的搜索算法可以做如下改進(以蛇點的5×5鄰域為例)——當(dāng)前蛇點為vi，相鄰蛇點為vi-1和vi+1，兩個相鄰蛇點的中垂線為l，如果中垂線l和vi的鄰域相交的話，則將相交的鄰域點作為搜索點(圖1中陰影部分)，否則，將當(dāng)前點向中垂線附近的鄰域進行移動。

圖1 目標蛇點的搜索域

判斷鄰域點和l相交，可以采用下面的準則：

設(shè)ε為定義好的一個比較小的值，設(shè)vi的鄰域點集為Pi，對于v(v∈Pi)，計算

(10)

若

f(v,vi+1,vi-1)≤ε

(11)

則判定v在中垂線附近，并把v歸為搜索域P(vi)(其中P(vi)?Pi)中的一個點。

如果鄰域點均不滿足式(11)，則P(vi)中僅含一個點v，且滿足下式

(12)

綜上所述，給定蛇點vi，圖像鄰域Pi和常數(shù)ε<1，則蛇點vi的移動鄰域P(vi)可以由下面的算法(Algorithm 2)得到。

Algorithm 2:Get Search Neighborhood1:procedure2: GetSearch Neighborhood(vi,Pi,ε)3: P(vi)←?;4: for all v∈Pi do5: Calculate f(v,vi+1,vi-1) as in (10);6: if f(v,vi+1,vi-1)≤ε then7: P(vi)←P(vi)∪{v};8: end if9: end for10: if P(vi)=? then11: P(vi)←{argminv∈Pif(v)};12: end if13: return P(vi);14:end procedure

2.2 外力模型

對于主動輪廓的外力場可以用不同的模型進行描述，例如GVF (Gradient Vector Flow)模型[2]和基于靜電場的外力模型[6]。本文使用靜電場模型計算主動輪廓線的外力，首先對目標圖像進行預(yù)處理，達到簡化圖像的目的，然后在此基礎(chǔ)上計算圖像的外力。

設(shè)圖像中非蛇點的像素集為M={m1,…,mM}，該點集每個象素的位置放置負電荷ei=-I(mi)(i=1,2,…，M)。其中I(mi)是象素mi的灰度值。蛇點集合為V={v1,…,vn}，每個蛇點的正電荷為1，由于圖像邊緣的電場強度較大，蛇點將在電場的作用下向物體的邊緣進行移動。同時，由于蛇點之間的斥力作用，蛇點將趨于均勻分布。設(shè)蛇點的數(shù)量為n，則輪廓線在靜電場中的勢能為：

(13)

其中k1、k2為事先選定的常數(shù)。

這里整個輪廓線V={v1,…,vn}的能量snake(V)定義為：

(14)

綜上所述，移動整個輪廓線V={v1,…,vn}則需分別對各個蛇點vi(i=1,…,n)進行操作，分別計算各蛇點的移動鄰域P(vi)、snake(vi)。具體過程可以由下面算法(Algorithm 3)來實現(xiàn)。

Algorithm 3:Modified Moving Snake1:procedure2:Modified Moving Snake(v1,…,vn,ε)3: for all vi∈{v1,…,vn} do4: Emin←êsnake(vi)5: Caculate P(vi) using Algorithm 2;6: for all p∈P(vi) do7: Calculate êsnake(P) using (14);8: if Emin>êsnake(P) then9: Emin←êsnake(P);10: vi←p;11: end if12: end for13:end for14:end procedure

2.3 圖像的二值化處理

為了減少構(gòu)造外力的計算量，在使用主動輪廓模型計算腦部MRI圖像外力前，首先利用圖像的直方圖對原始圖像進行二值化處理。在圖2中，圖2(a)為原始圖像，圖2(b)是原始圖像的灰度直方圖，由于直方圖的雙峰間的峰谷過于平坦，使用一維的閾值分割方法往往難以得到理想的灰度閾值。因此，為了得到較優(yōu)的灰度閾值，可以使用二維直方圖，本文使用原始圖和其八鄰域平滑圖得到二維直方圖。

(a)原始圖像 (b)灰度直方圖

圖3 二維灰度直方圖

圖3中高亮度的區(qū)域?qū)?yīng)二維直方圖中取值大的部分，由圖3可見，灰度二維直方圖取值較大的部分集中在平面的對角線附近。取橫軸閾值T1，縱軸閾值T2，可以把二維直方圖分為如下4個部分：

(15)

其中C和D區(qū)對應(yīng)圖像中的噪聲和圖像邊緣部分，概率較小，從而閾值劃分問題轉(zhuǎn)化為求閾值矢量(T1,T2)，使得熵函數(shù)φ(T1,T2)最小。

(16)

其中

(17)

即優(yōu)化目標為

(18)

對于式(15)中最優(yōu)參數(shù)，可以采用各種優(yōu)化算法進行計算[12-13]，如可以使用遺傳算法計算得到該類圖像灰度閾值，由腦部MRI圖像的先驗知識可知這類圖像的閾值分布情況，從而可以進一步縮小優(yōu)化算法的搜索空間，加快算法的收斂速度。

3 實驗結(jié)果及討論

為了突出對二值圖像邊緣的影響，本文還對目標圖像進行梯度運算，并以梯度模值表征圖像邊緣。圖4(a)是梯度運算的結(jié)果，圖4(b)是使用靜電場模型對圖4(a)的目標區(qū)域(矩形框中)進行計算后得到的圖像外力結(jié)果。

(a) (b)

由圖4(b)可以看出，在水腫區(qū)域的邊緣處形成較強的外力，即使在目標區(qū)域的凹陷區(qū)也有外力的存在。經(jīng)過上述算法進行圖像處理得到最終的實驗結(jié)果(見圖5)。

(a)初始分割結(jié)果

(b)使用原始圖像的主動輪廓線的分割結(jié)果

(c)使用簡化圖像并使用梯度力的分割結(jié)果

(d)使用靜電場模型得到的分割結(jié)果

4 結(jié) 語

本文討論了使用主動輪廓模型進行醫(yī)學(xué)MRI圖像分割的方法。針對預(yù)處理后的圖像采用了改進的輪廓點鄰域的搜索方法。利用圖像的二維直方圖對圖像二值化的閾值進行了優(yōu)化求解，目標圖像的二值化大大降低了后續(xù)處理的計算復(fù)雜度。同時，使用靜電場模型對預(yù)處理后的圖像的梯度場進行了外力計算。在該外力場的影響下，主動輪廓線可以較好地提取目標圖像的形狀。對腦部MRI圖像進行的分割實驗表明，該方法有效。

參考文獻:

[1] Amini A A，Weymouth T E，Jain R C.Using Dynamic Programming for Solving Variational Problem in Vision[J].IEEE Transactions on Pattern Analysis and Machine Intelligence，1990，12 (9) : 855-867.

[2] Xu C, Prince J L.Snakes，Shapes and Gradient Vector Flow[J].IEEE Transactions on Image Processing，1998，7(3): 359-369.

[3] Whitaker R.Modeling Deformable Surfaces with Level Sets[J].IEEE Computer Graphics and Applications, 2004,24(5):6-9.

[4] Malladi R, Sethian J A, Vemuri B C.Shape Modelling with Front Propagation: A Level Set Approach[J].IEEE Transactions on Pattern Analysis and Machine Intelligence, 1995,17(2):158-175.

[5] Williams D J，Shah M.A Fast Algorithm for Active Contours and Curvature Estimation[J].CVGIP: Image Understanding，1992，55(1): 14-26.

[6] Jalba A C，Wilkinson M H，Roerdink J B.CPM: A Deformable Model for Shape Recovery and Segmentation Based on Charged Particles[J].IEEE Trans.Pattern Analysis and Machine Intelligence，2004，26(10): 1320-1335.

[7] Xie X , Mirmehdi M.MAC: Magnetostatic Active Contour Model [J].IEEE Transactions on Pattern Analysis and Machine Intelligence，2008，30(4) :632-646.

[8] Wang T，Cheng I，Basu A.Fluid Vector Flow and Applications in Brain Tumor Segmentation [J].IEEE Transactions on Pattern Analysis and Machine Intelligence，2009，56(3): 781-789.

[9] Fatakdawala H，Xu J，Basavanhally A，et al.Expectation-maximization-driven Geodesic Active Contour with Overlap Resolution (EMaGACOR): Application to Lymphocyte Segmentation on Breast Cancer Histopathology [J].IEEE Transactions on Biomedical Engineering，2011，57(7): 1676-1689.

[10] Xie X.Active Contouring Based on Gradient Vector Interaction and Constrained Level Set Diffusion[J].IEEE Transactions on Image Processing，2010，19(1): 154-164.

[11] 廖亮,林土勝,張衛(wèi)東.基于靜電力方法的主動輪廓模型的腦部MRI圖像分割[J].生物醫(yī)學(xué)工程學(xué)雜志, 2008(4):26-29,45.

[12] 陳果，左洪福.圖像分割的二維最大熵遺傳算法[J].計算機輔助設(shè)計與圖形學(xué)學(xué)報，2002，14(6): 530-534.

[13] 周露芳，古樂野.基于量子遺傳算法的二維最大熵圖像分割[J].計算應(yīng)用，2005，25(8):1805-1807.