梁愷彬,管一弘
昆明理工大學 理學院,昆明 650500
醫學圖像分割是臨床醫療診斷的重要前提,也是醫學圖像分析和應用的重要基礎環節[1],而其中人腦圖像分割更是醫學圖像分割的一個重要分支。人腦圖像在人腦疾病的診斷中發揮著重要作用。核磁共振圖像(Magnetic Resonance Images,MRI)分割技術能幫助醫生對腦瘤、腦梗塞、帕金森綜合癥等腦部疾病患者進行快速而準確的診斷。利用腦影像檢查技術定性和定量地分析腦功能對有效診斷腦疾病有重要幫助[2]。MRI成像發明于1970年,是一種流行的醫療成像技術[3]。與其他成像方式相比,MRI成像技術具有更大的潛在優勢。MRI不會產生偽影,并且利用MRI技術進行檢測時不會對人體產生電離輻射,所以MRI成像技術對人體更加安全[4]。但是MRI成像的原理使得圖像表現出一定的模糊性,即部分容積效應(Partial Volume Effects,PVE),使得各個軟組織之間邊界不明確和不連續[5],同時在MRI成像的過程中,由于圖像的采集,傳輸,壓縮解碼,以及熱、電噪音,射頻線圈和磁場的非均勻性等影響,使得人腦MRI存在加性噪聲,這些都會影響人腦MRI的分割,增加了分割困難。
直到目前,人腦MRI的分割方法已經提出了很多,但是依然沒有一種分割方法能夠對各種目標進行有效的分割,臨床上使用的分割方法主要還是以人工分割為主。其中,陳允杰等人[2]為解決傳統的活動輪廓模型對強噪音或弱邊界的人腦MRI很難得到真實解的缺陷,提出了基于高斯混合模型的活動輪廓模型腦MRI分割,引入高斯混合模型并構造新的約束項,在新的約束項作用下減少噪音的影響,防止弱邊界泄漏。宋艷濤等人[5]為了解決高斯混合模型對噪聲十分敏感和馬爾科夫隨機場模型分割結果容易出現過平滑現象的缺陷,提出了一種新的基于圖像片權重的馬爾科夫隨機場圖像分割模型,對鄰域內的不同圖像片根據相似度賦予不同的權重,使其在克服噪聲影響的同時能保持圖像細節信息。Szilágyi等人[6]提出了一種改進的模糊C均值聚類算法(Enhanced Fuzzy C-Means,EnFCM),用灰度級數代替圖像大小重新定義隸屬度,使得算法運算速度有效提升。Cai等人[7]提出了一種快速的可泛化模糊C均值算法(Fast Generalized Fuzzy C-Means,FGFCM),它通過引入一種結合了空間相似度與灰度相似度的局部空間相似度,從而提升分割算法的分割準確性。Adhikari等人[8]提出了一種條件空間的模糊C均值聚類算法(Conditional Spatial Fuzzy C-Means,CSFCM),利用鄰域平均隸屬度來求解聚類中心,使得算法在抗噪性能方面有一定的提高。Xiong Taisong等人[9]提出了一種結合空間局部信息的雙重高斯混合模型(Double Gaussian Mixture Model,DGMM),有效地在高斯混合模型中融入了馬爾科夫隨機場(Markov Random Field,MRF)的空間局部關系的優點,得到了良好的分割效果。
在眾多分割方法當中,基于統計模型的分割方法也經常被使用。其中高斯模型是被運用的最廣泛的模型之一。而在基于MRF的分割方法中,經常使用高斯模型來建立圖像的條件概率密度函數從而簡化模型的建模過程。而基于有限高斯混合模型的分割方法,由于它的形式簡單和易于實現的特點而越來越被人們所關注。由于人腦MRI的強度可以看作一系列隨機變量,而這些隨機變量符合高斯分布特點,所以基于高斯模型的分割方法成為了人腦MRI分割的主要方法之一。但是傳統的高斯混合模型由于其忽略了空間信息以及高斯模型的線性加權結構,造成高斯混合模型的抗噪性能和魯棒性并不理想,對強噪聲或邊緣模糊的人腦MRI的分割效果較差。同時,傳統的高斯混合模型往往只針對觀測值的分布建模,忽略了分割結果的分布情況,使得模型不完整。針對這些缺點,本文提出一種基于隱高斯混合模型的人腦MRI分割方法。本文以真實人腦MRI以及模擬圖像同時進行實驗,實驗結果表明本文方法在魯棒性、可靠性和抗噪性能方面擁有很大的優勢。
高斯統計模型是用來逼近圖像直方圖的概率模型,在對灰度圖像的統計特性描述中,高斯模型是描述區域內灰度緩變的理想模型之一[10]。但是傳統的有限高斯混合模型具有抗噪性能和魯棒性較差的缺點,使其對噪聲圖像的分割往往得不到良好的分割效果。針對傳統的有限高斯混合模型的缺點,本文提出了一種隱高斯混合模型來對人腦MRI進行分割。隱高斯混合模型把分割結果的概率密度函數作為隱含數據引入到非線性結構的高斯混合模型中,運用非線性加權的混合高斯結構對人腦MRI進行建模,能夠增強抗噪性能和魯棒性。同時引入分割結果的概率分布到高斯函數項中,能使模型更趨于合理和完整。
假設人腦MRI的像素強度集為X={x1,x2,…,xN},xi為第i個像素的觀測值,并且xi(i=1,2,…,N)滿足獨立同分布的特性。設圖像所含有的目標類為C個,每個目標類都滿足已知的概率分布,記為N(xk|θi)。為了使得模型更加完整,把假設分割結果的概率分布作為隱含數據引入高斯混合模型當中,記為zik。同時為了簡化求解的過程,運用相乘形式的高斯混合形式來對模型進行構建,得到的隱高斯混合分布式為:

其中,zik為隱含數據的概率密度分布,wik為先驗概率,并滿足約束條件,參數集 Θ={θ1,θ2,…,θC}且 θi={μi,| Σi|},C 為類別總數,N(xk|θi)為單高斯項,μi與 | Σi|分別為類均值和類方差。
為了在分割模型中考慮到空間鄰域信息,定義了一個含有空間鄰域信息的高斯權重置信指數δi如下:

其中,β為平滑系數,?n表示以xn為中心的二階鄰域,為了增加分割準確性與抗噪性能,本文引入了含空間信息與平滑系數的高斯權重置信指數δi,提出了新的貝葉斯先驗概率形式wik如下:

由wik的定義可以看出,它是一個動態變化的值,通過參數μi和Σi以及平滑系數β的迭代更新實現其優化的過程。
隱含數據zik為高斯函數項的補充信息,根據貝葉斯公式,可以計算像素xi屬于每一個類的概率zik,則隱含數據zik概率密度函數如下:

這種新的表達形式(1),在迭代過程中,能夠將當前像素類的分布根據其鄰域的分布信息進行修改,并通過加權的表達形式,在分類過程中加入鄰域信息。
假設人腦MRI的像素強度集為X={x1,x2,…,xN},xi為第i個像素的觀測值,并且xi(i=1,2,…,N)滿足獨立同分布的特性,則有數據集X={x1,x2,…,xN}的似然估計表達式為:

由此可得其對數似然估計表達式為:

最大似然估計是常用的求解概率分布參數的方法,為了使得似然估計最大化,使用期望最大化(Expectationmaximization,EM)算法來對類均值 μi與類方差 | Σi|分別進行求解。EM算法是一種迭代算法,算法的每一次迭代都能保證似然函數值增加,并且收斂到一個局部極大值[11-12]。為了使似然函數能夠得到極大值,使夠得到參數μi和Σi的迭代公式如下:

為了對平滑系數β進行求解,使用牛頓迭代法來對β進行尋優。牛頓迭代法是一種使用導數的尋優算法,它每一步的迭代方向都是沿著當前點函數值下降的方向,具有收斂速度快的特性[13]。根據上述分析,運用牛頓迭代法得出本文提出方法的參數估計迭代更新公式如下:

因為梯度下降法按照固定的步長尋優,步長由人經驗設定很難取得良好的尋優效果;同時它容易陷入局部最優解。鑒于這些缺點,本文選用牛頓迭代法對平滑系數β進行尋優。牛頓迭代法的最大優點就是尋優速度為二階收斂,無需人為設定步長,尋優具有自適應性等特點,能夠很好地繞過局部最優解從而尋得一個較好解。
最后根據最大后驗概率準則,對圖像進行分割:

根據上述的分析,可以得到本文提出方法的參數估計迭代更新公式如下:

其中

依據本文提出的基于隱高斯混合模型,實現人腦MRI分割的算法步驟如下所示:
初始化
1.根據公式(2)計算每個像素的高斯分布;
2.根據公式(15)計算高斯權重置信度;
3.根據公式(14)計算高斯權重;
4.根據公式(13)計算隱含數據的雙高斯分布;
5.根據公式(16)更新高斯函數的類均值;
6.根據公式(17)更新高斯函數的類方值;
7.根據公式(18)更新平滑系數;
判斷:如果滿足條件迭代停止;否則繼續過程(1)到(7)。
輸出:根據最大后驗概率準則求得圖像分割后的標記場。
為了驗證算法的有效性,本文采用EnFCM、FGFCM、CSFCM、文獻[4]提出的基于模擬退火(Simulated Annealing,SA)與高斯混合模型的MRF算法(GMMMRF-SA)以及DGMM模型與本文方法進行比較。實驗分為兩步進行:第一步,采用EnFCM、FGFCM、CSFCM、GMM-MRF-SA、DGMM以及本文方法對模擬圖像進行分割,以錯分率(Misclassification Ratio,MCR)為指標比較它們的分割準確性;第二步,用EnFCM、FGFCM、CSFCM、GMM-MRF-SA、DGMM以及本文方法對真實的人腦MRI進行分割,比較它們的魯棒性、收斂性、可靠性以及區域內均勻性,并對比它們的分割效果。本文的所有實驗都是在64位Windows 8.1操作系統,i7處理器,2.4 GHz主頻,4 GB內存,Matlab2010a的平臺上完成。
為了對分割算法的分割準確性進行客觀的評價,本文采用三幅類別數分別為三類(k=3)、四類(k=4)和六類(k=6)的模擬圖像。由于模擬圖像為非圖像傳感器獲取的圖像,為了使得圖像環境更加接近真實的圖像環境,本文選擇高斯白噪聲與椒鹽噪聲的混合噪聲添加到模擬圖像當中,為每幅模擬圖像加入六種強度遞增的高斯噪聲與椒鹽噪聲的混合噪聲,在混合噪聲中的分割結果更能體現真實環境中的分割效果,也更能體現出分割模型的真實抗噪性能。分割結果與標準分割結果進行比較,采用MCR作為分割準確性的衡量指標。MCR被定義為錯分像素數與圖像總像素數的百分比,即

其中,NP代表分割錯誤的像素數,TP代表像素總數。實驗結果如圖1~6所示。
在上述實驗當中,圖1、圖3和圖5分別展示了EnFCM、FGFCM、CSFCM、GMM-MRF-SA、DGMM以及本文方法分別對k=3,k=4和k=6的噪聲圖像分割的效果。而圖2、圖4和圖6分別展示了三幅圖像在不同噪聲等級下分割的錯分率曲線。由錯分率曲線可以看出,與另外四種分割方法相比,本文分割方法對k=3,k=4和k=6的模擬圖像進行分割,在噪聲等級遞增的條件下依然保持穩定的低MCR值,這表明本文分割方法具有很好的分割準確性和抗噪聲能力。

圖1 模擬圖像的分割結果

圖2 k=3的模擬圖像的錯分率曲線

圖3 模擬圖像的分割結果
在對k=6的模擬圖像進行分割時,由于需要分割的類越多,各類的相似度就越高,從而導致各類的區分度下降,分割難度增加,再加上圖像受噪聲的影響,從而導致對分割難度的增加。從圖5可以看出EnFCM、FGFCM、CSFCM以及GMM-MRF-SA的分割結果依然帶有許多噪聲,而DGMM的分割結果雖然噪聲含量較低,但是圖像中間矩形部分被錯誤分割,而本文方法的分割依然能夠保持著良好的分割效果。從圖6可以看出,隨著噪聲強度的增加,DGMM展現出了不穩定性,而本文方法雖然也隨著噪聲強度增加錯分率有所上升,但是其錯分率一直保持較低,這證明了本文分割方法良好的分割能力和分割穩定性。

圖4 k=4的模擬圖像的錯分率曲線

圖5 模擬圖像的分割結果

圖6 k=6的模擬圖像的錯分率曲線
為了驗證本文方法對真實人腦MRI分割效果,本次實驗采用高斯噪聲分別為0%,3%,6%,9%,12%和15%真實的人腦T2加權圖像進行分割。把人腦MRI分為灰質(Gray Matter,GM)、白質(White Matter,WM)、腦脊液(Cerebrospinal Fluid,CSF)和背景(Background,B G)四類,分別采用EnFCM、FGFCM、CSFCM、GMM-MRF-SA、DGMM以及本文方法進行分割。為了全面客觀評價各分割方法,并針對真實人腦MRI缺少分割參考標準這一缺陷,本文提出利用系統方差均值RB和似然函數收斂性CV分別作為魯棒性和收斂性的衡量指標,并采用分割可靠性P[14]和區域內均勻性UR[15]作為分割效果的衡量指標。系統方差均值RB和函數收斂性CV的定義如下所示:

其中,i為實驗數,M為實驗總數,φi為分割結果不同組織的百分比,E為組織百分比均值,M為實驗總次數,t為迭代次數。RB代表在不同噪聲等級下對各組織分割的穩定程度,RB越接近零,表示分割方法越穩定,受噪聲干擾程度越低,能反映分割方法的魯棒性。分割可靠性P和區域內均勻性UR則是反映分割效果的指標,分割可靠性P與RB相類似,反映的是不同噪聲下分割結果的可靠程度,所不同的是P是全局分割效果的測度;而區域內均勻性UR反映的則是類像素間的相似程度,由于人腦MRI相同組織的像素相似程度高,所以利用區域內均勻性UR可以作為分割準確性的側面評價標準。P與UR的定義如下:

其中,Sn為第n次實驗的分割結果,| S1?S2?…?Sn|表示n次實驗的分割結果交集的總像素數,|S1?S2?…?Sn|表示n次實驗的分割結果并集的總像素數。N為圖像總像素,C為分割區域數,| Ri|表示區域Ri中的像素數,x(s)表示位于s像素的強度。實驗結果如圖7、8所示。

圖7 函數收斂性曲線圖

圖8 人腦MRI分割結果圖像(k=4(256×256))
圖7為牛頓迭代法與梯度下降法的函數收斂性曲線。由曲線可以看出,梯度下降法在第四次迭代時趨于發散并在第五次迭代迅速收斂,這表明梯度下降法尋優陷入了局部最優當中,收斂得出的結果為局部最優解;而牛頓迭代法的收斂趨勢較為平緩,雖然在第六次迭代時也有趨于發散,但是并沒有在第七次迭代中迅速收斂,而是平穩的到第十次迭代才收斂,這說明牛頓迭代法跨越了局部最優解尋得收斂結果。
圖8為EnFCM、FGFCM、CSFCM、GMM-MRF-SA、DGMM以及本文方法對含3%、9%和15%高斯噪聲的人腦MRI的分割結果。而表1則是它們的迭代次數、運算時間以及分割結果統計特征的比較。從分割效果上看,EnFCM、FGFCM、CSFCM、GMM-MRF-SA、DGMM以及本文方法對含3%、9%和15%加性噪聲的人腦MRI的分割效果各有不同,可以明顯看出,本文方法分割的WM連通性明顯好于其他分割方法,并且對GM的細節明顯保持更好;由表1可以看出,本文方法迭代次數較少,說明它具有較好的收斂性,但是運算時間較長,時效性方面并不占優勢。從統計特征上的對比,結合不同噪聲等級下各組織百分比曲線圖9~圖11與表2,可以看出,本文方法在不同噪聲等級下分割的WM、GM以及CSF的像素數和百分比變化較小,WM、GM以及CSF的百分比曲線起伏較小,且通過各組織的RB值對比可以看出,本文分割方法對WM、GM的RB值最小,對CSF的RB值為次最小,說明本文分割方法對WM、GM以及CSF的分割結果均較為穩定;同時,本文分割方法擁有最高的P值,這說明本文方法分割結果的可靠性最高,這都表明了本文分割方法具有優良的分割效果和魯棒性。

圖9 不同噪聲等級下的WM百分比

圖10 不同噪聲等級下的GM百分比

圖11 不同噪聲等級下的CSF百分比

表1 含3%,9%和15%加性噪聲的人腦MRI分割結果比較

表2 各組織的系統方差均值RB與分割可靠性P
表3為不同噪聲等級下的區域內均勻性UR。從表中縱向對比可以看出,與其他分割方法相比,本文分割方法,對含高斯噪聲分別為0%,3%,6%,9%,12%和15%的人腦MRI分割,區域內均勻性UR均為最高;而從橫向對比可以看出,本文分割方法的區域內均勻性UR穩定保持在99.5%左右,這些都側面反映了本文分割方法良好的分割準確性以及分割穩定性。但是從UR的定義以及EnFCM、FGFCM、CSFCM、GMM-MRFSA和DGMM的橫向比較可以看出,UR測度也存在一定的不足。由像素強度最大值與最小值的差值決定UR測度的分母,兩點間像素的強度差使得測度忽略了錯分像素的數量,容易導致測量結果過于片面而忽略整體效果,這也是EnFCM、FGFCM、CSFCM、GMM-MRFSA和DGMM的UR值隨噪聲等級的增加而有所增加的原因,同時分割方法的穩定性與魯棒性較差也是造成UR值隨噪聲等級的增加而有所增加的原因之一,這與分割方法的可靠性P值的高低保持了一致的反映,總的來說UR測度仍然有一定的參考意義。

表3 不同噪聲等級下的區域內均勻性UR %
傳統的基于高斯混合模型的分割方法由于忽略了空間信息和未考慮分割結果的分布情況導致模型不完整,從而導致對人腦MRI的分割效果并不理想,魯棒性與抗噪性能不強。本文運用加入了空間鄰域信息與分割結果的分布情況的非線性加權隱高斯混合模型作為人腦MRI主要分割模型,有效提高了分割的魯棒性和抗噪能力;同時,本文選用牛頓迭代法代替梯度下降法來進行尋優,使尋優能力大大增強,提升了尋優效率與尋優效果;雖然本文方法由于計算的信息較多且復雜程度較高使得計算耗時較長,時效性方面并不是最優的,但是所消耗的時間能夠使得分割的抗噪性能和魯棒性大大提高,這仍然是可以接受的。
:
[1]Duncan J S,Ayache N.Medical image analysis:progress overtwo decadesand the challengesahead[J].IEEE Trans on Pattern Anal Mach Intell,2000,22:85-106.
[2]陳允杰,張建偉,韋志輝,等.基于高斯混合模型的活動輪廓模型腦MRI分割[J].計算機研究與發展,2007(9):1595-1603.
[3]Balafar M A.Gaussian mixture model based segmentation methods for brain MRI images[J].Artificial Intelligence Review,2014,41(3):429-439.
[4]韓宇.基于馬爾科夫隨機場和貝葉斯理論的人腦MR圖像的分割研究[D].昆明:昆明理工大學,2014.
[5]宋艷濤,紀則軒,孫權森.基于圖像片馬爾科夫隨機場的腦MR圖像分割算法[J].自動化學報,2014(8):1754-1763.
[6]Szilágyi L,Benyó Z,Sizlágyi S M,et al.MR brain image segmentation using an enhanced fuzzy C-means algorithm[C]//Proceedings of the 25th Annual International Conference of the IEEE EMBS,Cancun,Mexico,2003,1(1):724-726.
[7]Cai Weiling,Chen Songcan,Zhang Daoqiang.Fast and robust fuzzy c-means clustering algorithms incorporating local information for image segmentation[J].Pattern Recognition,2007,40(3):835-838.
[8]Adhikari S K,Sing J K,Basu D K,et al.Conditional spatial fuzzy C-means clustering algorithm for segmentation of MRI images[J].Applied Soft Computing,2015(34):758-769.
[9]Xiong Taisong,Zhang Lei,Zhang Yi.Double Gaussian mixture model for image segmentation with spatial relationships[J]Journal of Visual Communication and Image Representation,2016,34:135-145.
[10]Udupa J K,Leblanc V R,Ying Zhuge,et al.A framework for evaluating image segmentation algorithms[J].Computerized Medical Imaging&Graphics,2006,30(2):75-87.
[11]邱藤.基于高斯混合模型的EM算法及其應用研究[D].成都:電子科技大學,2015.
[12]周鑫.基于EM算法的G0分布參數最大似然估計[J].電子學報,2013,41(1):178-184.
[13]姚金杰,韓焱.基于粒子群和牛頓迭代法的目標定位方法研究[J].計算機應用研究,2010(5):1700-1701.
[14]Chabrier S,Emile B,Rosenberger C,et al.Unsupervised performance evaluation of image segmentation[J].EURASIP Journal on Advances in Signal Processing,2006(1):1-12.
[15]Chabrier S,Rosenberger C,Laurent H,et al.Evaluating the segmentation result of a gray-level image[C]//Signal Processing Conference,2004.