水平集理論在磁共振腦圖像分割中的模型研究

田 丹，吳靜飛，范立南

（沈陽(yáng)大學(xué) 信息工程學(xué)院，遼寧 沈陽(yáng) 110044）

在醫(yī)學(xué)圖像處理領(lǐng)域，磁共振腦圖像分割技術(shù)已經(jīng)成為腦疾病診斷的重要輔助手段.其主要優(yōu)勢(shì)在于，腦磁共振成像具有較高的軟組織對(duì)比度、非侵入性特性和較高的空間分辨率.人腦的軟組織大致可以分為幾個(gè)主要區(qū)域，如灰質(zhì)（GM）、白質(zhì)（WM）和腦脊液（CSF）等.魯棒的腦組織分割是一個(gè)極具挑戰(zhàn)性的任務(wù)，因?yàn)榛颊叩囊苿?dòng)、采集時(shí)間的有限性和軟組織的邊界灰度不均勻性，均會(huì)帶來(lái)大量的數(shù)據(jù)處理問(wèn)題.

為了實(shí)現(xiàn)魯棒分割，一些經(jīng)典方法中引入了形狀和灰度分布等先驗(yàn)知識(shí).這些方法主要包括水平集方法［1-3］、神經(jīng)網(wǎng)絡(luò)方法［4］和統(tǒng)計(jì)法［5］等.其中，水平集方法的基本思想是，在圖像區(qū)域中指定一個(gè)初始輪廓，然后根據(jù)速度函數(shù)迭代更新該輪廓的形狀和位置，使其最終定位到目標(biāo)邊界.由于該方法能提取連續(xù)的目標(biāo)拓?fù)湫螤?，目前在醫(yī)學(xué)圖像分割領(lǐng)域應(yīng)用很廣泛.例如，Chan等人［6］提出了一種著名的基于區(qū)域特性的水平集方法，即CV模型.該模型假設(shè)圖像中的目標(biāo)和背景區(qū)域均具有灰度均勻性.而在實(shí)際應(yīng)用中，數(shù)字圖像通常存在灰度不均勻性問(wèn)題，特別是醫(yī)學(xué)圖像.例如，由于射頻場(chǎng)的不均勻性等因素，導(dǎo)致腦磁共振圖像的灰度均勻性變差，表現(xiàn)為同一組織的像素灰度沿空間呈緩慢平滑的變化.為了解決該問(wèn)題，國(guó)內(nèi)外學(xué)者提出了一些改進(jìn)方法.Li等人［7］提出一種能同時(shí)分割和去除偏移場(chǎng)的活動(dòng)輪廓模型.Tsai等人［8］提出了一種分段光滑模型，但該模型計(jì)算量較大，在每次迭代水平集函數(shù)時(shí)要額外求解兩個(gè)偏微分方程.本文在CV模型的基礎(chǔ)上，引入一種可控的局域算子，使水平集曲線的進(jìn)化基于圖像的局部信息，從而克服了CV模型的圖像均勻特性限制.

1 問(wèn)題描述

1.1 能量準(zhǔn)則

假設(shè)Ω?R2是圖像域，I：Ω→R是一個(gè)指定的灰度圖像.在水平集描述中，由水平集曲線進(jìn)化形成的進(jìn)化平面?？捎闪闼郊瘮?shù)φ表示，滿足如下關(guān)系［9］：

式中，Ωin表示Ω 的內(nèi)部區(qū)域；Ωout可定義為Ω＼Ωin，表示Ω的外部區(qū)域.

圖像分割問(wèn)題可以通過(guò)水平集曲線C的進(jìn)化實(shí)現(xiàn)，這一進(jìn)化過(guò)程等效于指定能量函數(shù)的最小化過(guò)程，其進(jìn)化終了的穩(wěn)定狀態(tài)對(duì)應(yīng)目標(biāo)邊界.著名的CV模型的能量函數(shù)定義如下：

式中，λ1和λ2為正常量；inside（C）和outside（C）分別表示輪廓線C的內(nèi)部和外部區(qū)域；c1和c2分別近似輪廓線外部和內(nèi)部區(qū)域的灰度平均值；｜C｜代表輪廓線C的長(zhǎng)度，是保證曲線平滑的調(diào)整項(xiàng).式（2）中，能量函數(shù)的前兩項(xiàng)為數(shù)據(jù)保真項(xiàng).CV模型的前提條件是假設(shè)圖像為分段常量，即在每個(gè)目標(biāo)區(qū)域內(nèi)灰度是均勻的.CV模型中能量函數(shù)的構(gòu)建是十分合理的.若曲線C在目標(biāo)外部，則能量函數(shù)的第一項(xiàng)近似于0，而第二項(xiàng)大于0.若曲線C在目標(biāo)內(nèi)部，則能量函數(shù)的第一項(xiàng)大于0，而第二項(xiàng)近似于0.若曲線C跨越目標(biāo)的內(nèi)部和外部區(qū)域，則能量函數(shù)的前兩項(xiàng)均大于0.不難看出，當(dāng)且僅當(dāng)曲線C在目標(biāo)邊界上時(shí)，能量函數(shù)才能達(dá)到最小值.

1.2 局域化處理

考慮到磁共振成像腦圖像存在灰度不均勻性問(wèn)題，這時(shí)依靠單一的全局灰度均值不能表征圖像特征，CV模型將區(qū)分不出目標(biāo)和背景.本文在CV模型基礎(chǔ)上，在局域范圍內(nèi)重新描述能量函數(shù)，從而使水平集輪廓進(jìn)化可以基于圖像的局域信息.基于局域信息的能量函數(shù)定義如下：

能量函數(shù)中引入了一個(gè)新的空間變量y，它獨(dú)立于變量x.換句話說(shuō)，x和y分別代表Ω內(nèi)獨(dú)立的空間點(diǎn).在函數(shù)L（x，y）的約束下，能量函數(shù)中的灰度I（y）可以被限制在以點(diǎn)x為中心的局部區(qū)域范圍內(nèi)，其局域尺寸可由半徑參數(shù)r控制.隨著點(diǎn)y逐漸遠(yuǎn)離中心點(diǎn)x，灰度I（y）對(duì)能量函數(shù)的作用將逐漸遞減并趨近于0.該局域區(qū)域可以是一個(gè)小的鄰域，也可以逐步擴(kuò)展到整個(gè)圖像區(qū)域.

為了獲取整個(gè)目標(biāo)邊界，能量函數(shù)應(yīng)該在整個(gè)圖像區(qū)域中作最小化處理.這可以通過(guò)最小化圖像區(qū)域中所有中心點(diǎn)x對(duì)應(yīng)的能量函數(shù)的積分來(lái)實(shí)現(xiàn)，定義如下：這里，能量函數(shù)的參量設(shè)定為水平集曲線C.為了處理拓?fù)浣Y(jié)構(gòu)的變化，將其轉(zhuǎn)變?yōu)樗郊Ｐ?

1.3 水平集模型

為了最優(yōu)化能量函數(shù)，進(jìn)化曲線C可以利用前面提到的零水平集函數(shù)φ表示.能量函數(shù)可以重新描述為

式中，M1（φ）＝H（φ），M2（φ）＝1-H（φ）；H（·）是Heaviside函數(shù)，定義為

這里，進(jìn)化曲線的長(zhǎng)度由Heaviside函數(shù)的積分求出.為了最小化該能量函數(shù)，需要滿足如下歐拉方程：

式中，δ（φ）是一維Dirac函數(shù)：

對(duì)于一個(gè)固定的水平集函數(shù)φ，c1和c2分別對(duì)應(yīng)進(jìn)化曲線外部區(qū)域和內(nèi)部區(qū)域的灰度平均值.

對(duì)于固定的c1和c2，采用標(biāo)準(zhǔn)的梯度下降法最小化能量函數(shù).水平集進(jìn)化的偏微分方程為

式中，div（）表示散度.在進(jìn)化方程中，第一項(xiàng)負(fù)責(zé)驅(qū)動(dòng)曲線C逐漸進(jìn)化到目標(biāo)邊界，第二項(xiàng)負(fù)責(zé)縮短和光滑進(jìn)化曲線.為了獲取整個(gè)目標(biāo)邊界，應(yīng)該最小化圖像區(qū)域中所有中心點(diǎn)x的能量函數(shù)的積分.

2 數(shù)值實(shí)現(xiàn)

針對(duì)本文提出的局域化水平集模型，需要計(jì)算進(jìn)化曲線上每個(gè)點(diǎn)的局域統(tǒng)計(jì)量.為了加快水平集的進(jìn)化速度，僅在由半徑r限制的零水平集附近的窄帶范圍內(nèi)迭代計(jì)算水平集函數(shù).而局域半徑r的選取應(yīng)基于感興趣目標(biāo)的尺寸和周圍鄰近雜斑的情況.例如，當(dāng)分割具有鄰近雜斑的小目標(biāo)時(shí)，應(yīng)選取較小的局域半徑；而當(dāng)分割具有較少鄰近雜斑的大目標(biāo)時(shí)，應(yīng)選取較大的局域半徑.

水平集函數(shù)的偏微分項(xiàng)可以由如下有限差分［10］機(jī)制作離散化處理：

進(jìn)化方程中的散度項(xiàng)記為L(zhǎng)，可作如下離散化處理：

水平集進(jìn)化過(guò)程的計(jì)算可分為兩步：初始化和迭代更新.而最終的圖像分割效果對(duì)初始水平集曲線并不敏感，但要求迭代過(guò)程中重新初始化水平集，以保證其滿足符號(hào)距離函數(shù)特性.這里如下定義初始水平集函數(shù)：

即在局域范圍內(nèi)部將其賦值為正整數(shù)3，而在局域范圍外部，將其賦值為負(fù)整數(shù)-3.這里所涉及的整數(shù)運(yùn)算與傳統(tǒng)的符號(hào)距離初始化函數(shù)相比能加快計(jì)算速度.而重新初始化水平集的目的是為了重構(gòu)水平集.因?yàn)樵谶M(jìn)化過(guò)程中，水平集函數(shù)的梯度方向可能會(huì)太陡或者太平，這可能導(dǎo)致數(shù)值計(jì)算不精確.當(dāng)初始化像素與進(jìn)化曲線交叉時(shí)，均需重新作初始化處理.本文提出的局域法的缺點(diǎn)是它要處理所有的局部數(shù)據(jù)統(tǒng)計(jì)量，計(jì)算量大.

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

為了實(shí)現(xiàn)魯棒圖像分割，多種方法在模型中引入了形狀和灰度分布等先驗(yàn)知識(shí).例如，水平集方法和神經(jīng)網(wǎng)絡(luò)方法.其中，水平集方法利用一個(gè)閉合曲線的進(jìn)化實(shí)現(xiàn)目標(biāo)檢測(cè)，所以能夠獲得連續(xù)的目標(biāo)邊界.下面對(duì)這兩種方法進(jìn)行性能比較.

對(duì)于具有灰度均勻性的近似分段常值圖像，著名的CV模型能得到連續(xù)且精確的分割效果.圖1中給出了基于CV模型分割一個(gè)兩目標(biāo)合成圖像的進(jìn)化過(guò)程.圖像大小為128×128，水平集參數(shù)υ＝0.2×255×255.圖1中分別用水平集初始輪廓、中間輪廓和最終收斂輪廓描述其進(jìn)化過(guò)程.

圖1 兩目標(biāo)合成圖像的分割結(jié)果Fig.1 Segmentation of a synthetic two-object image

但在實(shí)際應(yīng)用中，特別是對(duì)于醫(yī)學(xué)圖像，通常存在灰度不均勻性.圖2中分別給出了基于CV模型和局域化CV模型分割一人腦MR切片圖像的實(shí)驗(yàn)結(jié)果.圖像大小為128×128，水平集參數(shù)υ＝0.35×255×255.

為了分析和檢驗(yàn)當(dāng)圖像存在灰度不均勻性和隨機(jī)噪聲時(shí)分割算法的穩(wěn)定性，從brainweb數(shù)據(jù)庫(kù)中選取了多幅仿真腦磁共振圖像做仿真實(shí)驗(yàn).這里以其中兩幅正常腦解剖T1加權(quán)像為例.其噪聲分別是0和5%，灰度不均勻性分別為0和20%.圖像切片厚度為1.0mm，大小為258×258.

針對(duì)神經(jīng)網(wǎng)絡(luò)方法，采用自組織映射（SOM）網(wǎng)絡(luò)模型.為了得到一個(gè)滿意的映射關(guān)系，需要設(shè)定好網(wǎng)絡(luò)的拓?fù)浣Y(jié)構(gòu).實(shí)驗(yàn)結(jié)果表明，六角網(wǎng)格的拓?fù)浣Y(jié)構(gòu)效果最好.另一個(gè)影像分割效果的參數(shù)就是網(wǎng)絡(luò)尺寸.分別采用10×10，11×11，12×12，13×13，14×14，15×15的映射尺寸作測(cè)試.隨著尺寸的增加，分割效果得到改善.但當(dāng)網(wǎng)絡(luò)尺寸大于12×12時(shí)，分割效果不再有明顯的改善.因此，為了得到更精確的分割結(jié)果，同時(shí)減少網(wǎng)絡(luò)運(yùn)算時(shí)間，采用大小為12×12的神經(jīng)網(wǎng)絡(luò)［11］.

圖2 人腦MR切片分割結(jié)果Fig.2 Segmentation of a brain MR image

SOM神經(jīng)網(wǎng)絡(luò)方法包括兩個(gè)主要步驟：特征提取和像素分類.在特征提取過(guò)程中，提取一個(gè)二維歸一化的統(tǒng)計(jì)向量X＝（x1，x2）作為特征向量.該特征向量由像素及其八鄰域的歸一化均值和方差組成.該向量作為SOM神經(jīng)網(wǎng)絡(luò)的輸入.對(duì)于正常的腦磁共振圖像，網(wǎng)絡(luò)輸出節(jié)點(diǎn)數(shù)為4，分別對(duì)應(yīng)白質(zhì)、灰質(zhì)、腦脊液和背景的分類.為了提高網(wǎng)絡(luò)對(duì)噪聲的魯棒性，訓(xùn)練數(shù)據(jù)由1 000個(gè)隨機(jī)采樣的像素點(diǎn)構(gòu)成.網(wǎng)絡(luò)迭代數(shù)設(shè)定為1 000.

針對(duì)局域水平集方法，研究了半徑參數(shù)r對(duì)分割結(jié)果的影響.實(shí)驗(yàn)結(jié)果表明，引入局部半徑r能增強(qiáng)輪廓內(nèi)部和外部統(tǒng)計(jì)特征的平滑性.半徑越大，局部統(tǒng)計(jì)特性越平滑.但當(dāng)提取小目標(biāo)時(shí)，應(yīng)該選擇一個(gè)較小的半徑.這里令r＝7.

水平集模型中的其他參數(shù)設(shè)置如下：λ1＝λ2＝1.0，υ＝0.001×255×255 ，時(shí)間步長(zhǎng) Δt＝0.1.圖3和圖4分別給出了原始腦圖像和兩種算法的分割結(jié)果.

圖3 噪聲和灰度不均勻性均為0%的仿真MR腦圖像的分割結(jié)果Fig.3 Segmentation results of a simulated MR image with 0noise level and 0inhomogeneity effect

圖4 噪聲為5%、灰度不均勻性為20%的仿真MR腦圖像的分割結(jié)果Fig.4 Segmentation results of a simulated MR image with 5%noise level and 20%inhomogeneity effect

仿真結(jié)果表明，水平集方法可以提取連續(xù)的目標(biāo)邊界，并且在隨機(jī)噪聲和灰度不均勻性存在的情況下魯棒性更強(qiáng).

4 結(jié) 語(yǔ)

本文提出了一種局域水平集模型用于磁共振腦圖像的魯棒分割.因?yàn)樗郊椒ɡ瞄]合曲線的進(jìn)化獲取目標(biāo)邊界，所以其固有的特性就是能提取連續(xù)的目標(biāo)邊界.通過(guò)局域化處理，可以提高水平集方法對(duì)噪聲的魯棒性，并適用于灰度不均勻性圖像的處理.仿真結(jié)果表明，與傳統(tǒng)的引入先驗(yàn)知識(shí)的方法（例如神經(jīng)網(wǎng)絡(luò)方法）相比，局域水平集方法在魯棒去噪和目標(biāo)邊界提取的連續(xù)性方面效果更佳.

［1］ Li C，Xu C，Gui C，et al.Distance Regularized Level Set Evolution and its Application to Image Segmentation［J］.IEEE Trans.Image Processing， 2010，19 （12）：3243-3254.

［2］ Sanjay Kumar，Santosh Kumar Ray，Peeyush Tewari.A Combined Approach Using Fuzzy Clustering and Local Image Fitting Level Set Method for Global Image Segmentation［J］.Canadian Journal on Image Processing and Computer Vision，2012，3（1）：1-5.

［3］ Antunes S，Silva J，Santos J，et al.Phase Symmetry Approach Applied to Children Heart Chambers Segmentation：A Comparative Study［J］.IEEE Trans.Biomedical Engineering，2011，58（8）：2264-2271.

［4］ Awad M.An Unsupervised Artificial Neural NetworkMethod for Satellite Image Segmentation ［J］. The International Arab Journal of Information Technology，2010（7）：199-205.

［5］ Nadjib Nasri，Karim Mokrani，Abdenour Mekhmoukh.MRI Images Segmentation by FCM and Neighbor’s Statistical Characteristics［J］.International Journal of Research and Reviews in Soft and Intelligent Computing，2012，2（1）：127-130.

［6］ Chan T，Vese L.Active Contours without Edges［J］.IEEE Trans.Image Processing，2001（10）：266-277.

［7］ Li C，Huang R，Ding Z，et al.A Level Set Method for Image Segmentation in the Presence of Intensity Inhomogeneities with Application to MRI［J］.IEEE Trans.Image Processing，2011，20（7）：2007-2016.

［8］ Tsai A，Yezzi A， Willsky A S.Curve Evolution Implementation of the Mumford-Shah Functional for Image Segmentation， Denoising，Interpolation，and Magnification［J］.IEEE Trans.Image Processing，2001，10（8）：1169-1186.

［9］ Bernard O，F(xiàn)riboulet D，Thevenaz P，et al.Variational B-spline Level-set：A Linear Filtering Approach for Fast Deformable Model Evolution［J］.IEEE Trans.Image Processing，2009（18）：1179-1191.

［10］ 郝哲，朱一飛，王鐵男.基于差分法的土石壩數(shù)值模擬研究［J］.沈陽(yáng)大學(xué)學(xué)報(bào)，2010，22（2）：23-27.（Hao Z，Zhu Y F，Wang T N.Numerical Simulation Study on Earth-Rock Dam Based on Calculus of Difference［J］.Journal of Shenyang University，2010，22（2）：23-27.）

［11］ Dan Tian，Linan Fan.MR Brain Image Segmentation Based on Wavelet Transform and SOM Neural Network［C］∥Chinese Control and Decision Conference.Suzhou，2010：4243-4246.