，，，,2

(1.中國計量大學 信息工程學院，杭州 310018; 2. 北卡羅來納大學教堂山分校 醫(yī)學院，美國)

基于灰度直方圖的MR腦組織的提取方法

朱冠菲1，徐永秋1，米紅妹1，朱建明1,2

(1.中國計量大學信息工程學院，杭州310018; 2.北卡羅來納大學教堂山分校醫(yī)學院，美國)

根據腦磁共振圖像(magnetic resonance imaging，MRI)的灰度直方圖中不同峰值對應不同的腦組織灰度分布的特點，提出了一種基于灰度直方圖提取MR圖像中腦組織的方法；首先，為了克服傳統(tǒng)方法主觀選擇門限閾值等方面的不足，利用多項式曲線擬合灰度直方圖獲取區(qū)域生長門限閾值確定最優(yōu)種子點，并結合形態(tài)學重構方法進行顱骨分離，對腦MR圖像進行了預處理；然后，結合K均值聚類算法通過對灰度直方圖多峰值的選取確定初始聚類中心，將顱骨剝離后的腦組織圖像高效、精確地細化分割出腦脊液、灰質、白質；文中分別使用了模擬腦MRI數(shù)據以及真實腦部MRI數(shù)據進行測試，對分類結果進行定性、定量的分析，并與模糊C均值算法進行比較；結果表明，該方法提高了提取腦組織的工作效率和準確度。

灰度直方圖；曲線擬合；區(qū)域生長；K均值；腦組織提取

0 引言

核磁共振圖像(magnetic resonance imaging，MRI)腦組織提取在神經圖像分析中是配準、腦組織分類等的預處理步驟，起著重要作用[1]。在實際應用中，受圖像噪聲、灰度不均勻性等影響，給腦MR圖像的精準分割造成困難和挑戰(zhàn)[2-3]。

目前，腦組織分類和提取有大體可以分為三類：基于活動輪廓模型，基于形態(tài)學和基于圖像灰度?；诨顒虞喞挠校篠udipta R[4]等使用水平集產生遞階結構進行腦組織的三類分割和提取?；谛螒B(tài)學的有：陳筱[7]等和王大溪[8[9]等開發(fā)了非局部均值擴散(nonlocal means diffusion,NLD)算法, 進行組織灰度相似性統(tǒng)計，田換[10]等提出圖割和K均值相結合的算法對腦MRI圖灰度不均勻性進行了處理并細化達到分割腦白質和腦灰質的目的。其中基于灰度分割的主流方法是上述的聚類方法。傳統(tǒng)的模糊C均值(Fuzzy C Means, FCM)算法對參數(shù)初始值十分敏感，隨機產生初始聚類中心，用戶選定的種子點主觀性強，導致腦MR圖像分割結果不穩(wěn)定。為了克服FCM對參數(shù)初始值敏感的問題，本文提出了一種基于灰度直方圖的提取磁共振圖像腦組織的方法。首先，對磁共振腦圖進行交互式操作，利用多項式曲線擬合灰度直方圖得到選取區(qū)域生長全局閾值以確定最優(yōu)種子點進行顱骨剝離，并記錄各峰谷值，再通過對灰度直方圖多峰值的選取確定初始聚類中心，最后采用K均值聚類方法達到分割腦白質、腦灰質和腦脊液的目的。實驗表明，該算法能有效準確的從腦MRI圖中分割出腦白質、腦灰質和腦脊液。

1 基于灰度直方圖的圖像預處理

直方圖在醫(yī)學處理中有廣泛的應用，MR腦圖中多個峰值對應不同的腦組織灰度分布[11]。圖像的灰度直方圖作為對圖像各灰度值進行統(tǒng)計所得到的關于灰度值的函數(shù)，反映了圖像中每種灰度值出現(xiàn)的頻數(shù)。最簡單直接的閾值分割是根據灰度直方圖不同峰值進行初步預選，主觀性強。腦組織主要包括腦白質、腦灰質以及腦脊液，3種組織在圖中的灰度差異比較明顯。一幅典型無噪聲T1加權腦MRI圖(見圖1)的灰度直方圖如圖2所示?？梢姵ケ尘盀榱愕姆?，有3種峰值灰度差異較大，其相對應腦白質、腦灰質以及腦脊液3種腦組織。

圖1 無噪聲T1加權模擬腦MRI圖

圖2 圖1腦MRI圖對應的灰度直方圖

1.1 顱骨分離

由于一幅腦MRI圖像除了要分割的組織成分外，還包含顱骨和頭皮等成分，這些非腦組織的成分會直接影響腦組織的分割效果，所以需要進行顱骨剝離。采用區(qū)域生長法可分割出顱骨外皮脂和顱內腦組織，該算法判定準則中需要考慮的像素與初始種子點的像素灰度值差的絕對值小于某個門限值，即全局閾值。所以全局閾值的選取非常重要。本文在得到最佳全局閾值后，結合區(qū)域生長函數(shù)及形態(tài)學處理，將顱骨去除后的腦組織圖像作為先驗信息，為后續(xù)提取分類腦組織做準備。

1.1.1 最佳全局閾值獲取

灰度直方圖存在多個峰谷，主觀選取單一閾值難以達到分割的最好效果，熊平[12]等提出了一種根據CT圖像的灰度直方圖進行曲線擬合得到局部極小值對應的最大閾值對CT圖像分割。據此，本文提出多項式曲線擬合獲取局部極大值作為最佳全局閾值。最小二乘法多項式曲線擬合，是根據給定的m個像素點,并不要求這條曲線精確地經過這些點。通過給定數(shù)據點pi(xi,yi)，其中xi為像素灰度，yi為像素統(tǒng)計個數(shù)，i為下標，i=1,2,...,m，m為像素點個數(shù)，pi(xi,yi)即為灰度直方圖上橫軸灰度與縱軸統(tǒng)計個數(shù)所對應的點，求近似曲線y=f(x)。設擬合多項式為：

y=a0+a1x+...+akxk

(1)

式中，ai是系數(shù)，k為冪次。按偏差平方和最小的原則選取擬合曲線，并且采取二項式方程為擬合曲線的方法，得到多項式系數(shù)ai，同時我們就得到了擬合曲線。

擬合函數(shù)f(x)為一連續(xù)函數(shù)且可導。令f(x)′=0一階導便可求得像素灰度分布的局部極值。再根據二階導f(x)′原理及其凹凸性，判斷局部極大值。由系數(shù)矩陣ai求得f(x)′=0時對應的灰度x值，作為區(qū)域生長函數(shù)的最佳全局閾值。

1.1.2 基于區(qū)域生長和形態(tài)學的顱骨剝離

為了得到顱骨剝離得到腦組織的先驗信息，提出結合區(qū)域生長和形態(tài)學方法對整幅腦MRI圖像進行顱骨剝離。

區(qū)域生長是根據預先定義的生長準則來把像素或子區(qū)域集合成較大區(qū)域的處理方法。我們先將獲取的最佳全局閾值結合區(qū)域生長方法應用到腦MRI圖上。圖像區(qū)域內像素的相似性度量以像素灰度為準，以|Zi-S|≤T為生長準則[13]，進行不同成分之間的區(qū)域分割。Zi是被處理腦MR圖的像素灰度值，參數(shù)S這里定義為顱骨外某皮脂像素的像素值，T為1.1.1章節(jié)得到的全局閾值，此原則用來測試圖像中的像素是否與皮脂像素灰度(種子點)足夠相似。由于骨質中無氫原子，顱骨外的皮脂等對病理研究沒有意義，所以在本文算法中，顱骨與腦組織之間最大類別差值即為圖像中有著該像素值的點的種子。生長分割完畢，轉化為二值圖像，再利用形態(tài)學重構進行腦組織提取。

數(shù)學形態(tài)學中，二值圖像F=f(x,y)與結構元素B都是定義在二維笛卡爾網格上的集合。形態(tài)學算子有腐蝕、膨脹、開運算和閉運算4種[14]。利用各算子的形態(tài)學處理和閾值分割，找到腦部MRI圖像除背景之外的最大的連通區(qū)域，結合區(qū)域生長方法，進行組織分割及形態(tài)學重構。重構操作需要一個標記(marker)圖像h、一個掩模(mask)圖像g和一個結構元素B，其中標記圖像h必須是掩模圖像g的子集。使用8鄰域連通分量，得到模板二值圖像圖像J，對原圖像f(x,y)進行點乘運算，便可得到顱骨分離后的腦組織圖像temp。

1.2 多峰值選取

K均值的算法將圖像中的像素分割成k個簇，根據某種距離原則，再分配到相似的有限元中。K均值的k個聚類單元影響整個算法的計算結果，所以找到合適的k個初始聚類中心十分重要。3種不同的腦組織的“波峰”不同，選取灰度直方圖的“波峰”作為K均值的初始聚類中心方法[15]，據此，我們提出使用雙邊濾波及高斯分布方法[16]判斷并找出圖像灰度直方圖的峰值與谷值。采用相似性濾波核函數(shù)Sa[16]：

(2)

其中，p為隨機給定的某像素點位置，q點為核函數(shù)對應的像素點位置，Ip為q點對應的像素值大小，Np為高斯分布中最近的峰值，σp為峰值和谷值之間的一半距離(高斯分布95%的置信區(qū)間)。循環(huán)比較各像素點灰度值與峰值對應灰度值大小，控制高斯分布參數(shù)中心和標準偏差，實現(xiàn)找到圖像局部極值目的。除去背景的峰值，所得到的3個峰值對應3種MR腦組織。因此，我們將得到的3種峰值作為K均值的3個初始聚類值進行聚類計算。

2 基于K均值聚類算法的腦組織分類提取

K均值聚類算法由J.B.MacQueen于1967年提出的一種基于劃分的經典聚類算法，算法簡單高效、收斂速度快[17]。該方法可以應用在不同的領域，通常用于數(shù)據挖掘，包括圖像處理，它可以用來作為基于聚類的圖像分割方法。

各部分峰值為K均值算法的輸入值。將1.2章節(jié)得到的局部最大峰值設定為“波峰”，作為K均值算法的初始聚類中心，并將經過顱骨分離的預處理后的腦組織圖temp作為輸入對象，進行腦組織的分類和提取?；贙均值聚類算法的腦組織分類提取的基本步驟描述如下：

第1步：自動獲取閾值進行顱骨剝離的腦組織圖temp，并得到各“波峰”對應的灰度值R(k)；

第2步：將峰值從小到大排列，作為k個數(shù)據對象依次對應的初始聚類中心；

第3步：分別計算各個數(shù)據對象即二維腦MRI灰度圖到各個聚類中心的曼哈頓距離(Manhattan Distance)，根據最近鄰原則min賦值逐個劃分到最近的聚類中心所代表的簇中，計算誤差平方和準則函數(shù)E值，即均方差總和；曼哈頓距離又稱絕對距離，數(shù)學定義為：

(3)

式中，d(xi,yj)表示數(shù)據對象xi,yj之間的距離，p表示數(shù)據空間中相應的一個數(shù)據點，每個數(shù)據對象有p個屬性，k=1,2,…,p。

誤差平方和準則函數(shù)E定義為：

(4)

其中：E表示數(shù)據集中所有數(shù)據對象的均方差總和，mi為聚類子集xi的聚類中心點。

第4步：分別計算各個簇中所有數(shù)據對象的平均值作為新中心，重新分配數(shù)據對象，并重復2、3步驟，直到計算得到的E值和前一次得到的E值小于等于預先設定的閾值，則收斂。

第5步：輸出滿足終止條件的結果，即腦灰質、腦白質以及腦脊液的圖。

3 實驗結果與分析

本文分別使用模擬數(shù)據以及真實腦部MRI數(shù)據進行測試，軟件實驗平臺采用Matlab2013b，算法測試硬件配置為酷睿i3處理器，2G內存。

3.1 模擬數(shù)據實驗

本實驗對來自BRAINWEB[18]數(shù)據庫的仿真T1加權模擬腦圖進行分割和組織提取。體素大小為181×217×181，無噪聲和偏移場作用，層厚1 mm。選取其中一片作為示例(即原模擬腦圖圖1)進行灰度直方圖的多項式曲線擬合并進行顱骨分割以提取腦組織，其灰度直方圖多項式擬合曲線、峰值分布如圖3和圖4。

圖3 模擬腦MRI(圖1)多項式擬合曲線圖

圖4 模擬腦MRI圖(圖1)峰值分布

多項式擬合全局閾值結果以及各峰值見表1：

表1 模擬腦MRI圖(圖1)的灰度值參數(shù)

全局閾值代入生長函數(shù)后并進行形態(tài)學重構(詳見1.1.2章節(jié))得到剝離顱骨的腦組織圖像temp，如圖5(a)。將各峰值作為輸入數(shù)據，灰度直方圖擬合曲線得到的系列極值進行順序排列并作為k個輸入初始聚類中心，得到3種腦組織如圖5(b)(c)(d)。

圖5 MRI腦原圖(圖1)顱骨剝離后的腦組織部分

3.2 真實數(shù)據實驗

測試真實腦MRI圖像數(shù)據為醫(yī)院收集的15組病人腦部磁共振T1加權圖，掃描序列為快速自旋回波，1.5T的通用GE系統(tǒng)采集，空間分辨率為512×512。

真實腦MRI圖一定存在噪聲。強場越高，噪聲越少。磁共振圖像大多符合瑞利分布的加性噪聲[19]。所以，首先對這一些列測試數(shù)據進行平滑預處理，本文采用常用的中值濾波，3×3模板進行預處理。圖6展示了其中一組序列中某一病人經室間孔的橫斷層面，對去噪前后進行了定性對比。

圖6 某一病人經室間孔的橫斷層面

本文同時對去噪預處理的效果進行了定量評價。均方誤差(Mean Squared Error, MSE)是衡量“平均誤差”的一種較方便的方法，可以評價圖像數(shù)據的變化程度。峰值信噪比(Peak Signal to Noise Ratio, PSNR) 是最普遍，最廣泛使用的評鑒畫質的客觀量測法。均方誤差越小，歸一化均方誤差(Normalized Mean Square Error，NMSE)[20]越小，峰值信噪比越大，說明圖像質量越好，實驗數(shù)據具有更好的精度。從15組病人選取出共30張經室間孔的橫斷層面的腦MR圖進行中值濾波操作，得到的實驗數(shù)據去噪后具體評價指數(shù)見表2。

表2 去噪后評價指數(shù)

由表2可知，去噪后MRI腦圖的均方誤差約4.2，歸一化均方誤差極小，而峰值信噪比達到71，由均方誤差越小，峰值信噪比越大，說明圖像質量越好可知，中值濾波去噪效果良好。

去噪后得到的灰度直方圖、多項式擬合曲線圖如圖7所示，各峰值分布圖如圖8所示。

圖7 真實腦MR圖(圖6(a))灰度直方圖及擬合曲線 圖8 真實腦MR圖(圖6(a))峰值分布

圖6(a)真實腦MR數(shù)據的全局閾值及各峰值見表3。

表3 真實腦MR圖(圖6(a))灰度值參數(shù)

同模擬實驗步驟，將全局閾值作為輸入數(shù)據，得到二值模板J和顱骨剝離后腦組織分別如圖9(a)、(b)。

本文采用商業(yè)軟件SPM分割的圖作為金標準(專家醫(yī)生手動分割的結果 )并增加了結合傳統(tǒng)的FCM經典算法作為對比，3種方法對圖6(a)提取并分類腦組織的結果分別如圖10～12。

在進行顱骨剝離等預處理操作后，將圖9中提取的腦組織圖(b)作為輸入數(shù)據，得到的各峰值作為k個輸入初始聚類中心，結合K均值算法進行腦組織的分類提取。

圖9 真實腦MR圖(圖6(a))的二值模板J(a)，提取出的腦組織部分(b)

為了對本文算法的分割效果進行質量評估，使用相似性測度KI ( Kappa index ) 指數(shù)[21-22]來定量地評價算法的分割性能，KI 指數(shù)的定義如下：

(5)

圖10 SPM軟件分割的腦MR圖分割提取后的 圖11 結合FCM算法腦MR圖分割提取后的 圖12 本文方法腦MR圖分割提取后的

表4 算法性能對比

由表4定量分析，可以看出結合K均值的方法進行提取分類腦白質，腦灰質和腦脊液的KI指數(shù)均比結合傳統(tǒng)FCM算法高，尤其是腦白質的分類效果上尤為顯著。另外，不考慮顱骨剝離等圖像處理操作情況下，在后處理運行時間上，K均值方法較于FCM算法平均時間快了近百分之三十八，本文方法在精度和效率上都提高了。

4 結束語

本文提出了一種基于灰度直方圖的提取磁共振圖像腦組織的方法。利用多項式曲線擬合灰度直方圖得到自適應選取區(qū)域生長最優(yōu)種子點并結合了形態(tài)學二值化模板，對腦組織進行提取，并優(yōu)化k均值聚類方法細化分類達到分割腦白質、腦灰質和腦脊液的目的。與傳統(tǒng)模糊C均值算法對比，具有計算時間快，精確度高的優(yōu)點。本文算法也存在一些缺陷，由于受噪聲等影響，顱骨剝離提取腦組織的圖像邊緣還是存在一些毛刺，穩(wěn)定性不高，分割精度會降低；此算法針對非眼球部位的腦組織有廣泛適用性，但是對于包含眼球的腦MRI圖具有局限性，需要進一步深入研究。

[1] Zhang Haiyan, Li Haiyun. Automated MRI brain tissue extraction method[J].Computer Engineering and Applications, 2014,50(16): 168-172.

[2]Parogios N, Deriche. Geodesic active regions and level set methods for supervised texture segmentation[J]. International Journal of Computer Vision, 2002,46(3): 223-247.

[3]Ginneken B, Alejandro F. F, Joes J. Active shape model segmentation with optimal features[J]. IEEE Transactions on Medical Imaging,2002,21(8): 924-933.

[4]Sudipta Roy, Samir Kumar B. A new method of brain tissues segmentation from MRI with accuracy estimation[J]. Procedia Computer Science,2016,85: 362-369.

[5]陳志彬, 邱天爽. 一種基于FCM和Level Set 的MRI醫(yī)學圖像分割方法[J].電子學報.2008,36(9): 1733-1736.

[6]王 克. 基于支持向量機和水平集方法的核磁共振腦圖像分割[D].西安:西安理工大學,2011.

[7]陳 筱. 基于區(qū)域生長和數(shù)學形態(tài)學的MRI圖像處理研究[D]. 武漢:中南民族大學,2012.

[8]王大溪. 基于形態(tài)學的腦部MRI圖像顱骨剝離算法[J]. 計算機技術與發(fā)展,2015,25(12): 206-209.

[9]Xing Xiuxi, Zhou Youlong, Jonathan S, et al. PDE based spatial smooting : a practical demonstrations of impacts on MRI brain extraction, tissue segmentation and registration[J]. Magnetic Resonance Imaging, 2011, 29: 731-738.

[10]田 換,覃 曉,元昌安,等. 基于K-means和圖割的腦部MRI分割算法[J].數(shù)據采集與處理,2016,31(5): 974-982.

[11]羅述謙, 周果宏. 醫(yī)學圖像處理與分析[M]. 北京:科學出版社,2003.

[12]熊 平, 黎 妲, 徐 平.一種基于MATLAB的CT腦組織圖像提取方法[J].生物醫(yī)學工程學進展,2009,30(1): 17-19

[13]岡薩雷斯. 數(shù)字圖像處理( matlab版)[M]. 北京:電子工業(yè)出版社,2005.

[14]楊 杰. 數(shù)字圖像處理及MATLAB實現(xiàn)(第2版) [M]. 北京:電子工業(yè)出版社,2013.

[15]兆 學, 喻海中, 陳 浩. 基于灰度直方圖多峰值選取的腦組織MRI圖像K-means聚類分割方法研究[J]. 生物醫(yī)學工程學雜志,2013,30(6): 1164-1170.

[16]De Silva, Fernando, Kodikaraarachchi, et al.Adaptive sharpening of depth maps for 3D-TV[J]. Electronics Letters, 2010,46(23): 1546-1548.

[17]Milan J, Petr S, Tomá? K, et al. Implementation of K-means segmentation algorithm on Intel Xeon Phi and GPU: Application in medical imaging[J]. Advances in Engineering Software, 2017,10(3): 21-28.

[18]Brain Web:Simulated Brain Database[DB/OL].http://brainweb.bic.mni.mcgill.ca/cgi/brainweb1.

[19]Crum WR, Camara O, Hill DL. Ceneralized overlap measures for evaluation and validation in medical image analysis[J]. IEEE Transactions on Medical Imaging,2006, 25(11): 1451-1461.

[20]陳藍鈺,常 嚴,王 雷,等.基于正則化迭代的并行磁共振圖像重建算法[J]. 計算機測量與控制,2015,23(12): 4177-4179.

[21]陳志彬. 非參數(shù)變形模型結合模糊技術的MRI圖像分割[D]. 大連:大連理工大學,2010.

[22]王新寧, 林相波, 袁 珍. 基于FCM聚類算法的MRI腦組織分割方法比較研究[J].北京生物醫(yī)學工程,2015,34(3): 221-228.

GrayHistogram-basedMethodofMRBrainTissueExtraction

Zhu Guanfei1, Xu Yongqiu1,Mi Hongmei1,Zhu Jianming1,2

(1.Department of Information Engineering, China JiliangUniversity, Hangzhou 310018, China; 2. College of Medicine, University of North Carolina(Chapel Hill), North Carolina,USA)

This paper presents a tissue segmentation method based on the fact that the peak distribution of gray scale histogram in MR brain images corresponds to different brain tissue distribution. Firstly, in order to separate the skull, we propose a method which combines optimal threshold selection and adaptive regional growing algorithm. The optimal threshold is obtained by using polynomial curve fitting to histogram data, overcoming the drawbacks of subjective threshold selection associated with traditional methods. We then apply adaptive regional growing algorithm to complete the separation of the skull. Secondly, the initial cluster centers are determined by selecting the peak values of gray histogram, and then the K-means clustering algorithm is used to refine the segmentation of cerebrospinal fluid, gray matter, white matter. Finally, this method is tested on both simulated MRI data and human brain MRI images. We perform both qualitative and quantitative analyses in comparison with other image segmentation algorithms. Results show that the proposed algorithm can improve the efficiency and accuracy brain tissue segmentation.

gray histogram; curve fitting; region growing; K- means; brain tissue extraction

2017-04-24；

2017-05-15。

朱冠菲(1992-)，女，安徽安慶人，碩士研究生，主要從事醫(yī)學圖像處理方向的研究。

朱建明(1963-)，男，美籍華人，博士，教授，主要從事醫(yī)學圖像處理、信號處理方向的研究。

1671-4598(2017)11-0170-04

10.16526/j.cnki.11-4762/tp.2017.11.043

TP391;R445.2

A