基于靜電場模型和二維灰度直方圖的腦部磁共振圖像的分割

，，， ，

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

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

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

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

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

本文針對腦部MRI圖像的目標(biāo)區(qū)域進(jìn)行輪廓提取。首先對腦部MRI圖像的二維直方圖進(jìn)行處理以提取閾值，對圖像進(jìn)行二值處理后，求得二值圖像的梯度邊緣。經(jīng)過上述處理后，基于圖像像素的外力場建模的復(fù)雜度大大降低。本文使用靜電力模型對預(yù)處理后的目標(biāo)圖像計(jì)算初始輪廓線所受的力，輪廓線所受的力使目標(biāo)圖像不斷收縮和位移，最終達(dá)到收斂的狀態(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)能強(qiáng)度、外能強(qiáng)度以及人為交互力的強(qiáng)度。如果忽略外部干預(yù)的話，則可以令Econt(s)=0。

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

(2)

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

1.2 圖像的外力

圖像的外力一般可以取：

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)

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

(6)

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

(7)

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

(8)

對于目標(biāo)函數(shù)(7)的最小化，可以使用Amini等提出的動態(tài)規(guī)劃法[1]，或者Williams等提出的貪婪法[5]，其中貪婪法收斂速度較快。如果設(shè)當(dāng)前蛇點(diǎn)為vi，其鄰域點(diǎn)集為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算法的改進(jìn)

2.1 蛇點(diǎn)移動鄰域

(9)

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

圖1 目標(biāo)蛇點(diǎn)的搜索域

判斷鄰域點(diǎn)和l相交，可以采用下面的準(zhǔn)則：

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

(10)

若

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

(11)

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

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

(12)

綜上所述，給定蛇點(diǎn)vi，圖像鄰域Pi和常數(shù)ε<1，則蛇點(diǎn)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 外力模型

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

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

(13)

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

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

(14)

綜上所述，移動整個輪廓線V={v1,…,vn}則需分別對各個蛇點(diǎn)vi(i=1,…,n)進(jìn)行操作，分別計(jì)算各蛇點(diǎn)的移動鄰域P(vi)、snake(vi)。具體過程可以由下面算法(Algorithm 3)來實(shí)現(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)造外力的計(jì)算量，在使用主動輪廓模型計(jì)算腦部MRI圖像外力前，首先利用圖像的直方圖對原始圖像進(jìn)行二值化處理。在圖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)化目標(biāo)為

(18)

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

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

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

(a) (b)

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

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

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

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

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

4 結(jié) 語

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

參考文獻(xiàn):

[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].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2002，14(6): 530-534.

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