張海燕,李海云
首都醫科大學生物醫學工程學院,北京,100069
核磁共振圖像腦組織自動提取方法
張海燕,李海云
首都醫科大學生物醫學工程學院,北京,100069
核磁共振圖像的腦組織提取是神經圖像處理研究中的一個重要步驟。將傳統的幾何活動輪廓模型與二值水平集函數相結合,提出了一種新型的二值水平集活動輪廓模型,并基于該模型提出了一種能夠自動、準確實現MRI腦組織提取的方法。該方法在腦組織內部自動設定最優初始輪廓曲線,將該演化曲線隱含地表示成一個高維函數的零水平集,零水平集在基于區域的圖像力驅動下不斷演化并達到待分割腦部圖像的邊緣。將基于該方法的腦組織提取結果與作為金標準的專家手動分割結果和其他流行算法相比較,結果表明提出的腦組織提取方法能夠自動、準確和快速地提取MRI腦組織,是一種魯棒性較好的MRI腦組織提取方法。
腦組織提取;灰度直方圖;活動輪廓模型;水平集;質量評估
核磁共振圖像(M agnetic Resonance Imaging,MRI)腦組織提取是將核磁共振序列腦圖像中的腦部組織與非腦組織分離,去除腦外組織,也稱為頭骨剝離。腦組織的準確提取是神經圖像處理中一個重要的步驟,也是其他圖像處理算法例如配準、腦組織分類或者灰度不均勻性校正等的預處理步驟[1-3],例如,只有準確地去除非腦組織才能夠實現腦內組織如腦白質、腦灰質和腦脊液的準確分割[4]。腦組織的準確提取不僅是腦組織體積測量以及實現三維重建的關鍵技術,在神經圖像處理和分析中也是最耗時的預處理步驟之一,因此國際上提出了許多的腦組織提取算法(Brain Extraction A lgorithm s,BEAs),目前許多研究致力于開發準確、自動的腦組織提取算法,這些算法在不同程度上提高了腦組織提取的精度或速度[5-6],但是各種算法的分割質量也大不相同,從而會影響后續的圖像分析[7]。Fennema-Notestine等人[8]比較了最常用的幾種頭骨剝離算法的性能,認為每種算法都有自己的優點和缺點,但是還沒有一種算法能夠適用于大規模圖像的處理與分析。
目前MRI腦組織提取方法大體可分為三類:基于圖像灰度的方法、基于形態學的方法和基于活動輪廓模型的方法。基于活動輪廓模型的方法是目前圖像分割中應用最廣泛的一種方法。Smith等人[9]提出了一個用于頭骨剝離的Brain Extraction Tool(BET)算法,該算法將一組基于形態學和基于圖像特征的力應用于演化曲面的切向和法向,驅動演化曲線到達待分割圖像邊緣。Zhuang等人[10]提出了一個用于頭骨分離的基于模型的水平集方法(Model-based Level Set,M LS),該模型將演化曲線作為高維函數的零水平集,利用基于腦皮層灰度特征的圖像力和曲率力控制曲線的演化。形變模型的主要優點是能夠直接產生閉合的參數曲線或曲面,并對噪聲和偽邊界有較強的魯棒性。與其他兩種方法相比,基于形變模型的方法更容易實現自動的頭骨剝離,但對初始輪廓線的位置較敏感、計算量大造成分割時間較長。
隨著MR圖像數量的增長及具有臨床應用意義的實時圖像處理的需要,開發全自動、準確和快速的腦組織提取算法變得尤為重要。本文基于二值水平集活動輪廓模型提出了一種新的MRI腦組織提取算法,使用該算法對臨床收集的MR真實腦圖像和網絡上下載的的MR腦圖像數據自動分割腦組織與非腦組織邊界,并將分割結果與專家手動分割結果相比較,實驗結果表明該算法能自動、準確地進行腦組織提取而且有很好的魯棒性。
2.1 幾何活動輪廓模型
傳統幾何活動輪廓模型(Geometric Active Contour Model,GACM)是由Caselles等人[11]與Malladi等人[12]分別獨立地提出的基于曲線演化和水平集方法的活動輪廓模型,該模型利用水平集函數來隱含地表示模型輪廓線,從而成功地解決了輪廓線拓撲變化的問題。傳統的幾何活動輪廓模型的曲線演化方程可表示成以下形式:


傳統幾何活動輪廓模型是通過檢測圖像邊緣來提取目標邊界的一種方法,對有較好對比度的圖像分割效果較好。
2.2 二值水平集活動輪廓模型
傳統水平集活動輪廓模型雖然具有自動處理輪廓線拓撲變化的優點,但為了保證水平集方法在數值實現過程中的穩定性和精確性,水平集函數在演化過程中需要反復地被重新初始化為其零水平集的符號距離函數,而這個初始化過程是個非常耗時的運算過程。Zhang等人[13]提出了一種二值水平集活動輪廓模型。在該活動輪廓模型中,水平集函數僅取-1和1這兩個數值。在二值水平集函數更新過程中,將該二值水平集函數與一個二值圖像相互進行轉化,并使用形態學中的腐蝕與膨脹算子對該二值圖像進行操作,使其模擬曲線的演化過程。與傳統水平集活動輪廓模型相比,二值水平集活動輪廓模型在運行效率上得到了極大的提高,同時保持了自動處理輪廓線的拓撲變化的能力,但卻喪失了曲線演化的漸進性,從而使得分割效果類似于簡單的閾值分割法,不具有針對性。
2.3 一種新型的二值水平集活動輪廓模型
由于MR圖像上腦組織與非腦組織的邊緣較弱,腦組織之間邊界模糊以及圖像噪聲的影響,導致MR圖像的分割問題非常復雜和困難。基于傳統幾何輪廓模型的分割方法使用圖像的邊界梯度信息,邊緣定位較精確,性能上優于傳統的基于邊緣檢測的分割方法,但是也存在一定的局限性,即對噪聲和初始條件很敏感,檢測腦組織和非腦組織的邊界時輪廓曲線容易發生邊緣泄露或者過分割,而且在輪廓線的初始化問題上,初始輪廓線必須被完全地設置在目標邊界的內部或者外部,否則無法得到正確的分割結果。本文在傳統幾何活動輪廓模型的框架下,將圖像的邊界信息和區域統計信息相結合,提出了一種新型的二值水平集活動輪廓模型。首先構造具有圖像統計特征信息的SPF函數如下:

其中Gσ是標準偏差為σ的高斯核函數,*代表卷積算子。f1和f2為兩個平滑常數,定義如下:

其中是正則化的Heaviside函數,Gσ是方差為σ的高斯核函數。
方程(3)中的SPF函數也可以稱為區域函數[14],它的取值范圍為[-1,1]。區域函數可以利用區域統計信息調節壓力的符號,使得當演化曲線位于感興趣區域內部時演化曲線膨脹或位于外部時收縮,因此SPF函數可以解決所謂的由弱邊緣引起的邊緣泄露問題。用方程(3)中的SPF函數代替方程(2)的邊緣檢測函數g,得到本文提出的用于腦組織提取的水平集演化方程如下:

基于以上提到的新型二值水平集活動輪廓模型,本文提出的MR圖像腦組織提取方法主要包括三個步驟:首先根據圖像灰度直方圖估計圖像灰度參數以得到頭部的二值圖像;然后在腦區內部自動設定自適應的初始輪廓曲線,并將該輪廓曲線隱含地表示成高維水平集函數的零水平集,零水平集在圖像力的驅動下不斷演化并達到待分割腦部圖像的邊緣。最后介紹了該方法中的人腦MR圖像分割的簡化原則和兩頭分割方案。
3.1 圖像灰度參數估計及頭部的二值化圖像
采用Sm ith[9]提出的方法來估計灰度的有效范圍。通過分析人腦MR圖像的直方圖特征,選取和設定高低兩個閾值t1和t2構成灰度的有效范圍[t1,t2],其中較低的閾值t1用來去除代表背景,顱骨,空氣和其他較低灰度級噪點的立體像素;較高的t2閾值用來消去類似脂肪和血管系統的較亮的非腦組織。然后根據學習經驗,確定腦部區域和非腦區域的分割閾值tm:

其中局部閾值Th為一常數,通常該常數取值由圖像數據的學習經驗確定,且與圖像的掃描方向有關,即橫軸位,冠狀位和矢狀位三個方向的取值不同。分割閾值tm確定后,以該閾值將腦部組織和非腦組織進行粗略分割,得到頭部二值化圖像。
3.2 自適應的水平集初始輪廓曲線的自動設定
水平集初始輪廓曲線的自動設定不僅是保證腦組織提取完全自動化的首要步驟,同時由于分割結果的準確性對初始輪廓曲線極為敏感,所以初始輪廓線應盡可能地靠近實際的輪廓才能得到最好的分割結果。另外由于一些非腦組織(例如皮膚和肌肉組織)具有和腦組織相似的灰度特性,所以如果初始輪廓曲線包含這些非腦組織分割時就會產生錯誤。現存的初始化方法通常存在準確性不高或者運算效率不高的缺點。Zhuang等人[10]將腦部區域的中心作為初始化水平集的圓心,取一個較小的半徑作為初始化水平集的半徑,以保證初始化水平集完全在腦區內部。這種水平集的初始化方法雖然能提高分割結果的正確性,但是由于初始水平集的半徑很小,降低了運算效率。本文提出了一種自適應的初始化水平集的自動設定方法,取二值化頭部圖像腦區的中心為初始化水平集的圓心,但是在保證初始化水平集在腦區內部的前提下半徑盡可能的大,從而提高了運算效率,如圖1所示。初始化方法如下:
(1)對原始圖像用閾值tm處理后得到腦部二值化圖像。
(2)對二值化后的圖像計算其中腦部區域的最左點和最右點,得到最左列和最右列的直線距離d。
(3)通過觀察頭顱MR圖像特點,發現頭顱大小在MR圖像的三個掃描方向上有所不同,因此根據學習經驗,對于橫軸位的MR圖像,可以將顱腦形狀近似為邊長等于d的正方形;而對于冠狀位和矢狀位的MR圖像,顱腦形狀都可近似為矩形,矩形的寬度都等于d,但是寬與高之比近似為5∶4。
(4)基于步驟3得到顱腦的近似形狀,初始化零水平集曲線為一個圓,圓心位于正方形或矩形的中心,半徑分別等于正方形邊長的三分之一(相對于橫軸位MR圖像)和矩形寬度的五分之一(相對于冠狀位和矢狀位MR圖像)。這種初始輪廓曲線的自動設定過程稱之為自適應的初始輪廓曲線自動設定方法,即根據輸入的顱腦MR圖像掃描方向或個體差異可自己調整初始輪廓線的位置及大小,從而保證初始水平集在腦區內部并盡可能的靠近圖像待分割邊緣。
3.3 二值水平集函數設定及分割過程中簡化原則和兩頭分割方案
為簡化水平集函數的初始化過程,使用二值水平集函數取代傳統水平集函數。初始化水平集函數如下:

通過對MR腦圖像進行研究,發現同一序列的前幾層和后幾層圖像上包含的腦組織區域很小甚至為零,更談不上分割腦組織的邊緣,所以前面和后面的幾層可以不予理會。根據學習經驗提出以下簡化原則:算法分割從每一圖像序列的第十分之一層開始,到第十分之九層結束。又因為同一序列顱腦MR圖像的相鄰層具有較強的相干性,所以采取如下的初始輪廓曲線設定方案:(1)從開始分割的層數到序列的中間層,由于腦區面積呈現增大的趨勢,所以當前層的最終演化輪廓線可以作為后面相鄰層的初始輪廓線,這樣既保證初始輪廓曲線在腦區內部又使初始輪廓曲線盡可能地靠近待分割邊緣;(2)從序列的中間層到分割結束的層數,由于腦區面積呈現逐漸減小的趨勢,所以為了保證初始輪廓曲線在腦區內部,采取反向分割順序,即先從最后一層開始分割,將其作為前面相鄰層的初始輪廓線,直到中間層結束。這種方案也可以稱為兩頭分割方案,即分割時分別從前、后兩頭向中間層過渡。實驗結果表明以上介紹的分割簡化原則和兩頭分割方案,不僅可以完全自動地在腦區內部設定與待分割邊緣靠近的初始輪廓線,而且大大提高了算法效率和分割結果的準確性。
本文采用M atlab 7.0作為軟件實驗平臺,算法測試的硬件配置為CPU:Intel雙核T6500 2.1 GHz,內存2 GB。測試數據為臨床收集的10組正常人的T1 MR腦圖像,每位受試者的體數據共176層,空間分辨率為256×256,層厚為1 mm,無間隔。圖2為應用本文算法對10例正常人腦MR圖像的部分腦組織提取結果。

圖1 自適應的初始輪廓曲線自動設定方法

圖2 MR TIW I真實圖像橫軸位、冠狀位和矢狀位的腦組織提取結果
為了對本文算法的分割效果進行質量評估,從網絡腦圖像標準分割庫(Internet Brain Segmentation Repository,IBSR)[15]下載10組正常人體大腦MR圖像(包括待分割圖像數據和對應的專家手動分割結果)。用提出的方法對這10組正常人體冠狀位MR腦圖像進行腦組織提取,分割結果如圖3所示。將算法自動分割結果與專家手動分割結果進行比較,采用靈敏度(sensitivity)、特異性(specificity)和FP_Rate系數來衡量算法自動分割結果與IBSR提供的專家手動分割結果的一致性:
靈敏度和特異性分別是算法識別出的腦組織和非腦組織的正確分類的像素的百分比:

FP_Rate系數是算法自動分割結果中把實際上是非腦組織的體素誤判為腦組織類的體素個數與專家手動分割結果中分類為腦組織體素個數的比值。因此FP_Rate系數的值越低,表明算法的分割結果越準確。

圖3 本文算法的分割結果與專家手動分割結果

其中TP表示判斷正確的腦部組織的體素;TN表示判斷正確的非腦組織區域;FN表示分割結果中把實際上是腦組織的體素誤判為非腦組織類的體素個數;FP表示分割結果中把實際上是非腦組織的體素誤判為腦組織類的體素個數。靈敏度和特異性的取值在0~1之間,通常來講,靈敏度的值越大,分割結果越準確。但是有一種特例,那就是如果分割算法比較保守(臨床應用時通常要采取保守措施),即寧肯將許多非腦組織誤判為腦組織而盡量避免將任何腦組織分類為非腦組織,那么此時FN就等于0,所以這種情況下靈敏度總是等于1。所以一個分割算法不能只是單一的使用靈敏度系數,必須綜合其他系數進行評估。
將本文提出的算法分割結果與專家手動分割結果比較,以專家手動分割結果為模板計算算法的靈敏度、特異性和FP_Rate系數,結果如表1所示。低靈敏度表示分割結果遺漏一些腦部區域,低特異性表示頭骨剝離后結果中腦部區域錯誤地包含非腦區域,而FP_Rate系數的值越低表明算法的誤判率越低。從表1可見,綜合靈敏度、特異性和FP_Rate系數三個參數,本文提出的方法與傳統的腦組織提取方法BET,BSE相比具有明顯的優越性,分割精度明顯提高,因為BET算法在分割過程中總是采取保守措施,所以它的靈敏度最高但也導致特異性和FP_Rate系數性能最差。BET算法如此設計的初衷可能是因為從臨床應用的角度出發,避免誤判任何的腦組織顯然要比去掉所有的非腦組織更重要。與M LS相比,本文的算法在分割精度上性能稍差,但FP_Rate和Specificity兩項系數性能略優于M LS;分割效率方面,由于本文提出的初始化水平集自動設定方法在很大程度上提高了運算效率,因此與其他三種算法相比,分割速度具有明顯的優越性,因此提出的腦組織提取方法可有效地實現MR圖像腦組織的自動和準確提取。

表1 各種算法性能比較
本文提出了一種新型的二值水平集活動輪廓模型用于腦組織提取。與傳統的腦組織提取方法相比,該方法具有以下優點:初始化方法使用二值水平集函數代替傳統水平集函數,避免了傳統算法中費時的重新初始化步驟,而且自適應的水平集輪廓曲線能夠盡可能地靠近待分割的邊界,從而大大提高了算法效率;第二,提出的活動模型利用圖像的統計信息作為約束條件,對具有弱邊界的腦組織能夠獲得較好的分割效果;第三,可以在腦區內部自動設定初始水平集,實現了腦組織提取的全自動化;最重要的是,該方法非常簡單,對原始MR圖像直接分割,無需任何預處理步驟,方便使用。因此本文提出的腦組織提取方法可有效地實現腦組織和非腦組織的自動和準確分割。然而,本文提出的模型對一些含有病變組織的MR腦圖像,分割精度降低,因此可繼續深入研究,使算法能夠適用于大規模MR腦圖像的處理與分析。
[1]Woods R P,Dapretto M,Sicotte N L,et al.Creation and use of a talairach-compatible atlas for accurate,automated,nonlinear intersubject registration,and analysis of functional imaging data[J].Human Brain Mapping,1999,8:73-79.
[2]Shattuck D W,Sandor-leahy S R,Schaper K A,et al.Magnetic resonance image tissue classification using a partial volume model[J].Neuro Image,2001,13:856-876.
[3]Strother S,La Conte S,Kai Hansen L,et al.Optimizing the fMRI dataprocessing pipeline using prediction and reproducibility performance metrics:I.A preliminary group analysis[J].Neuro Image,2004,23:196-207.
[4]Zhang Y Y,Brady M,Smith S.Segmentation of brain MR images through a hidden Markov random field model and the expectation-maximization algorithm[J].IEEE Trans on Med Imaging,2001,20:45-57.
[5]M arroquin J,Vemuri B.An accurate and efficient Bayesian method for automatic segmentation of brain MRI[J].IEEE Trans on M ed Imaging,2002,21:934-945.
[6]Liu J X,Chen Y S,Chen L F.Accurate and robust extraction of brain regions using a deformable model based on radial basis functions[J].Journal of Neuroscience Methods,2009,183:255-266.
[7]Boesen K,Rehm K,Schaper K,et al.Quantitative comparison of four brain extraction algorithms[J].Neuro Image,2004,22(3):1255-1261.
[8]Fennema-notestine C,Ozyurt I B,Clark C P,et al.Quantitative evaluation of automated skull-stripping method applied to contemporary and legacy images,effects of diagnosis,bias correction,and slice location[J].Hum Brain Mapp,2006,27(2):99-113.
[9]Sm ith S M.Fast robust automated brain extraction[J]. Hum Brain Mapp,2002,17:143-155.
[10]Zhuang A H,Valentino D J,Toga A W.Skull-stripping magnetic resonance brain images using a model-based level set[J].Neuro Image,2006,32:79-92.
[11]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.
[12]Malladi R,Sethian J,Vemuri B.Shape modeling with front propagation:a level set approach[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,1995,17(2):158-175.
[13]Zhang S.A fast level set implementation method for image segmentation and object tracking[C]//Proc SPIE,Applications of Digital Image Processing XXIX,63121R, doi:10.1117/12.682228,2006-08-24.
[14]Zhang K H,Zhang L,Song H H.Active contours with selective local or global segmentation:a new formulation and level set method[J].Image and Vision Computing,2010,28:668-676.
[15]IBSR was developed by the center for Morphometric Analysis at Massachusetts General Hospital.Internet Brain Segmentation Repository[DB/OL].[2012-05-12].http://www. cma.mgh.harvard.edu/ibsr.
ZHANG Haiyan,LI Haiyun
College of Biomedical Engineering,Capital Medical University,Beijing 100069,China
The segmentation of brain region from non-brain region in magnetic resonance images,also referred to brain tissue extraction,is an important and difficult task due to the convoluted brain surface and weak boundaries between brain and non-brain tissue.This paper proposes an accurate and automatic method of brain tissue extraction based on a novel binary level set active contour model.An initial contour,which is automatically set inside the brain,is driven by a region-based image force until it well fit the brain and non-brain boundaries.Quantitative evaluation of the new method is performed by comparing the result of the new method to those obtained using manual segmentation in 10 sets normal adult MR brain images.Experimental results demonstrate that the proposed method can provide accurate and automatic brain tissue extraction results.
brain tissue extraction;intensity histogram;active contour model;level set;quality evaluation
A
TP391
10.3778/j.issn.1002-8331.1209-0273
ZHANG Haiyan,LI Haiyun.Automated MRI brain tissue extraction method.Computer Engineering and Applications,2014,50(16):168-172.
北京市自然科學基金(No.4122018)。
張海燕(1977—),女,博士,副教授,研究方向為醫學圖像處理,生物醫學工程。E-mail:zhangxyz@ccmu.edu.cn
2012-09-24
2013-01-16
1002-8331(2014)16-168-05
CNKI網絡優先出版:2013-01-25,http://www.cnki.net/kcms/detail/11.2127.TP.20130125.1658.002.htm l