王志剛, 馮云超
(湖南師范大學(xué) 信息科學(xué)與工程學(xué)院 ,湖南 長(zhǎng)沙 410081)
根據(jù)世界衛(wèi)生組織的最新統(tǒng)計(jì),腦腫瘤是世界上最常見的癌癥死亡類型之一[1],且腦腫瘤的發(fā)病率和死亡率近年呈上升趨勢(shì),對(duì)人類的生命健康造成了嚴(yán)重的威脅.在疾病診斷過程中,準(zhǔn)確獲得腦腫瘤的尺寸、形狀等信息能夠?yàn)槟X疾病的治療提供很大的幫助.因此,將腦部腫瘤從核磁共振成像(Magnetic Resonance Images,MRI)圖像中準(zhǔn)確完整地分割出來對(duì)于腦疾病的診斷具有重要的意義.
MRI對(duì)人體無電離輻射傷害,有較高的安全性[2],并且對(duì)人體軟組織顯示清晰,因此被廣泛地應(yīng)用于腦部疾病的診斷和治療中[3].然而,由于人類的腦部結(jié)構(gòu)非常復(fù)雜,且MRI成像存在邊緣模糊、灰度不均勻、易受噪聲干擾等問題[4],因此腦腫瘤MRI圖像分割一直是研究的熱點(diǎn)和難點(diǎn)問題.
核模糊C-均值聚類算法(KFCM)是一種無監(jiān)督聚類算法,具有效果穩(wěn)定、運(yùn)算速度快等特點(diǎn),因此被廣泛地應(yīng)用于醫(yī)學(xué)圖像處理領(lǐng)域[5].然而,KFCM算法對(duì)于噪聲較為敏感,且聚類精度不高,降低了算法的實(shí)用性.為此,學(xué)者們對(duì)該算法進(jìn)行了大量的研究和改進(jìn),王燕等[6]將核模糊聚類算法與馬氏距離相結(jié)合,提出一種FCM-KM算法,抗噪性得到明顯提升;汪敏等[7]將核函數(shù)與模糊局部信息聚類相結(jié)合,提出了一種基于方差系數(shù)加權(quán)的模糊聚類分割方法,具有較好的抗噪性和分割精度;Hua L等[8]將PSO算法與KFCM算法相結(jié)合,進(jìn)一步提高了算法的分割準(zhǔn)確率.但是,這些優(yōu)化算法大多數(shù)只考慮了像素的灰度信息,忽略了像素之間的鄰域信息,并且采用單個(gè)高斯核函數(shù)進(jìn)行數(shù)據(jù)映射,限制了算法的泛化和學(xué)習(xí)能力,因此仍然存在一些過分割和欠分割現(xiàn)象,并且對(duì)于噪聲的魯棒性也不高.
針對(duì)現(xiàn)有算法的不足,本文提出一種基于邊緣保持濾波和改進(jìn)核模糊聚類的腦腫瘤圖像分割方法.實(shí)驗(yàn)表明,該方法在去除噪聲的同時(shí)保留了更多的圖像邊緣信息,具有較好的噪聲魯棒性,并且分割精度也更高.
傳統(tǒng)的去噪算法雖然能有效去除圖像中的噪聲,但是卻模糊了圖像的邊緣,不利于醫(yī)學(xué)圖像的后處理,本文采用一種新型的加權(quán)引導(dǎo)濾波算法進(jìn)行圖像預(yù)處理,在去除噪聲的同時(shí)可以把圖像的邊緣較好地保留下來.
引導(dǎo)濾波是何凱明等[9]提出的一種新型濾波器,它假設(shè)輸入與輸出之間存在如下的線性關(guān)系:
qi=akIi+bk,?i∈wk,
(1)
其中:q為輸出圖像;I為引導(dǎo)圖像;i和k為像素索引;a、b是當(dāng)窗口w中心位于k時(shí)該線性函數(shù)的不變系數(shù).為了使輸入圖像p與輸出圖像q的差別最小,使其轉(zhuǎn)化為最優(yōu)化問題,其代價(jià)函數(shù)為:
(2)
其中:ε為調(diào)整系數(shù),是為了防止ak過大;p為輸入圖像;I可為任意圖像,也可以是輸入圖像本身.利用線性回歸,可解得系數(shù)a、b:
(3)
(4)
由于原始的引導(dǎo)濾波算法,在不同的窗口內(nèi)采用了相同的ε調(diào)整系數(shù),沒有考慮到不同窗口內(nèi)像素之間的差異性,使得邊緣保持效果不佳.Li等[10]提出一種基于方差的加權(quán)引導(dǎo)濾波算法,它使用3×3窗口內(nèi)的局部方差計(jì)算邊緣權(quán)重,利用邊緣權(quán)重對(duì)調(diào)整系數(shù)ε進(jìn)行懲罰,來增強(qiáng)算法的魯棒性,其定義如下:
(5)
由于并不是所有方差較大的區(qū)域都對(duì)應(yīng)圖像的邊緣區(qū)域,所以使用局部方差計(jì)算的邊緣權(quán)重并不是很好的邊緣懲罰因子.本文采用基于Canny算子的邊緣權(quán)重代替文獻(xiàn)[10]中的局部方差,Canny算子在邊緣監(jiān)測(cè)中精確度較高,但是傳統(tǒng)的Canny算子為了避免噪聲的影響會(huì)先進(jìn)行高斯濾波去噪操作,由于高斯濾波不能很好地保護(hù)圖像的邊緣,所以本文采用基于雙邊濾波的Canny算子,在降噪的同時(shí)還能很好地保留邊緣信息.改進(jìn)后的權(quán)重因子定義如下:
(6)
其中:CB(i)為雙邊濾波Canny算子;γ取(0.001×L)2,L為圖像的灰度值范圍;N為像素總數(shù).這樣一來,圖像的邊緣部分就會(huì)比平坦部分分配到更大的權(quán)重,將邊緣權(quán)重因子引入到代價(jià)函數(shù)中可得:
(7)
利用線性回歸,求解系數(shù)a、b得:
(8)
(9)
核模糊聚類算法的基本思想是將原始數(shù)據(jù)通過核函數(shù)映射到高維特征空間進(jìn)行處理[11],這樣可以突出不同類別樣本之間的特征差異,使得其在高維特征空間中變得線性可分,然后再進(jìn)行模糊C均值聚類得到分割結(jié)果.
設(shè)X={x1,x2,x3,…,xn}是原始空間Rs通過非線性變換φ(·)映射至特征空間Rp上的一個(gè)數(shù)據(jù)集,則KFCM算法在高維特征空間下的目標(biāo)函數(shù)為:
(10)
傳統(tǒng)的核模糊聚類算法(KFCM)采用單一的高斯核函數(shù)進(jìn)行數(shù)據(jù)映射,其定義為:
(11)
由于單一高斯核函數(shù)僅有一個(gè)可以調(diào)節(jié)的參數(shù)σ,這對(duì)聚類算法的學(xué)習(xí)和泛化能力起到了一定的限制作用[12],所以本文采用混合高斯核函數(shù)來替代原先的單一高斯核函數(shù),這樣可以通過調(diào)整多個(gè)σ參數(shù)來優(yōu)化數(shù)據(jù)在高維特征空間中的分布,可以有效提高算法的性能.
混合高斯核函數(shù)的定義如下:
(12)
為防止濾波后的殘余噪點(diǎn)對(duì)圖像分割造成干擾,本文將馬爾科夫隨機(jī)場(chǎng)的局部先驗(yàn)概率引入到算法中,由于其考慮了像素的鄰域信息,因此具有較強(qiáng)的抗噪性.
假設(shè)yij是(i,j)像素點(diǎn)的灰度值,yij=k表示該像素點(diǎn)屬于第k類,設(shè)該像素點(diǎn)的馬爾科夫模型為y={yij|yij∈K},其中K={1,2,3,…,k},k為類別數(shù),根據(jù)Hammersley-Clifford定理,吉布斯隨機(jī)場(chǎng)與馬爾科夫隨機(jī)場(chǎng)是等效的[13],可以得到馬爾科夫隨機(jī)場(chǎng)的先驗(yàn)概率計(jì)算公式:
(13)
其中:Ω表示標(biāo)號(hào)場(chǎng);y為標(biāo)號(hào)場(chǎng)中的元素;U(y)代表能量函數(shù),其定義如下式:
(14)
其中:c表示基團(tuán);C表示基團(tuán)集合;Vc(y)表示勢(shì)團(tuán)函數(shù),其定義如下式所示:
(15)
其中:yt為鄰域像素的標(biāo)號(hào);ys為中心像素標(biāo)號(hào);β為勢(shì)團(tuán)參數(shù).
本文在混合核函數(shù)聚類算法的目標(biāo)函數(shù)中添加一正則項(xiàng),將馬爾科夫隨機(jī)場(chǎng)的先驗(yàn)概率引入,增強(qiáng)算法的抗噪性,改進(jìn)的目標(biāo)函數(shù)如下式所示:
(16)
其中:uij為第j個(gè)數(shù)據(jù)對(duì)第i類的隸屬度;dij為第j個(gè)像素到第i個(gè)聚類中心在核空間中的歐式距離;Pij表示馬爾科夫隨機(jī)場(chǎng)的先驗(yàn)概率;λ為約束系數(shù).
根據(jù)拉格朗日乘數(shù)法,可以求得隸屬度迭代公式如式(17) 所示,聚類中心迭代公式如式(18)所示:
(17)
(18)
由于核模糊聚類算法的初始聚類中心為隨機(jī)確定,容易陷入局部極值造成算法的不穩(wěn)定性,所以對(duì)上述算法做出進(jìn)一步改進(jìn),利用粒子群算法初始化聚類中心,增強(qiáng)算法的穩(wěn)定性.
粒子群優(yōu)化 (Particle Swarm Optimization,PSO) 算法是由Kennedy和Eber-hart提出的優(yōu)化算法[14].該算法是受到鳥群搜尋食物飛行行為啟發(fā)提出的,具有易實(shí)現(xiàn)、需調(diào)節(jié)參數(shù)少、收斂速度快等優(yōu)點(diǎn).
在D維搜索空間中由N個(gè)粒子匯聚成一個(gè)群落,群落中的每個(gè)粒子都在尋找自己的最佳位置,第i個(gè)粒子由Xi=(xi1,xi2,…,xiD)表示位置,由Vi=(vi1,vi2,…,viD)表示其速度,搜尋過程中,個(gè)體的最優(yōu)位置定義為pi=(pi1,pi2,…,piD),記為pbest,群體中所有粒子在之前搜尋過程中的最佳位置用g表示,記為gbest.粒子的速度和位置更新公式如下所示:
(19)
(20)
其中:c1和c2為學(xué)習(xí)因子;r1和r2是兩個(gè)隨機(jī)值,取值范圍為[0,1];w是慣性權(quán)重.
算法的具體實(shí)現(xiàn)過程如下:
步驟1:設(shè)定加權(quán)引導(dǎo)濾波算法的相關(guān)參數(shù),規(guī)整化因子為0.01,半徑為3,對(duì)圖像進(jìn)行去噪預(yù)處理;
步驟2:設(shè)定粒子群算法的相關(guān)參數(shù),粒子群規(guī)模N=30,最大迭代次數(shù)tmax=30,學(xué)習(xí)因子c1=c2=2,然后通過粒子群算法求得算法的初始聚類中心V(0);
步驟3:設(shè)定聚類數(shù)目c=3,迭代終止條件ε=10-5,最大迭代次數(shù)Tm=100,高斯核數(shù)目g=2,高斯核參數(shù)σ1=150,σ2=20,核函數(shù)權(quán)重因子α1=0.7,α2=0.3,正則項(xiàng)約束系數(shù)λ=0.6;
步驟4:輸入由步驟2獲得的初始聚類中心V(0),根據(jù)式(17)和式(18)不斷更新聚類中心Vi和隸屬度矩陣U,如果‖J(i)-J(i-1)‖≤ε或者達(dá)到最大迭代次數(shù),則停止迭代,否則繼續(xù)進(jìn)行本步驟;
步驟5:根據(jù)最終得到的隸屬度矩陣把目標(biāo)像素指定為相應(yīng)的最佳聚類中心的值,輸出圖像;
步驟6:結(jié)合形態(tài)學(xué)對(duì)圖像進(jìn)行后處理,平滑腫瘤圖像的邊緣.
本次實(shí)驗(yàn)的腦腫瘤MRI圖像來自BraTS2018數(shù)據(jù)庫,該數(shù)據(jù)庫的MRI圖像共有四種形態(tài),分別是 T1、T1c、T2 和 FLAIR,由于FLAIR圖像常用于發(fā)現(xiàn)腫瘤[15],所以本次實(shí)驗(yàn)選擇FLAIR圖像.實(shí)驗(yàn)環(huán)境為Windows10、AMD Ryzen5 3500U處理器和16 GB內(nèi)存,實(shí)驗(yàn)平臺(tái)為Matlab 2018a.
為了證明加權(quán)引導(dǎo)濾波算法的有效性,從數(shù)據(jù)庫中隨機(jī)挑選某位患者的腦腫瘤MRI圖像,抽取其第65張切片和第75張切片為實(shí)驗(yàn)圖像,在圖像中添加方差為0.01的高斯噪聲,使用該算法進(jìn)行去噪,并且同原始引導(dǎo)濾波[9]、雙邊濾波[16]和基于方差的加權(quán)引導(dǎo)濾波[10]進(jìn)行對(duì)比,算法的相關(guān)參數(shù)均采用原文獻(xiàn)中的推薦參數(shù),其去噪效果如圖1~圖3所示.
圖1 第65張切片的去噪效果Fig.1 Denoising effect of the 65th slice
從圖1和圖2可以看出,本文算法在視覺效果上強(qiáng)于其他保邊濾波算法,去噪后的圖像無明顯噪點(diǎn);從圖3和圖4可以看出,經(jīng)本文算法去噪后,圖像邊緣清晰,跟原圖像最接近,其他算法的涂抹感較為嚴(yán)重.為了防止主觀評(píng)價(jià)因素的影響,本文選擇峰值信噪比(PSNR)和結(jié)構(gòu)相似性(SSIM)兩個(gè)評(píng)價(jià)指標(biāo)對(duì)算法進(jìn)行評(píng)估,兩者的值越大說明算法降噪性能越好.
圖2 第75張切片去噪效果Fig.2 Denoising effect of the 75th slice
圖3 切片65細(xì)節(jié)放大圖Fig.3 Enlargement of section 65 details
圖4 切片75細(xì)節(jié)放大圖Fig.4 Enlargement of section 75 details
峰值信噪比定義如下:
(21)
其中:MAX表示圖像的灰度級(jí);MSE表示原圖像與處理圖像之間均方誤差.
結(jié)構(gòu)相似性定義如下:
(22)
以上面兩幅圖像為例,對(duì)四種算法的PSNR值和SSIM值進(jìn)行比較,比較結(jié)果如表1所示.
表1 不同降噪算法的比較
從表1可以看出,改進(jìn)算法的PSNR值和SSIM值均高于傳統(tǒng)算法,兩幅圖片的PSNR平均值相比傳統(tǒng)算法提升0.804 1~2.096 2 dB,SSIM平均值相比傳統(tǒng)算法提升0.031 2~0.065 4,說明改進(jìn)算法的性能更好,適用性更強(qiáng),能有效去除腦核磁共振圖像中的高斯噪聲.
為了證明所提算法對(duì)于腦腫瘤分割的有效性,從數(shù)據(jù)庫中隨機(jī)挑選某位患者的腦腫瘤MRI圖像,抽選3個(gè)切片,并且向圖像中添加方差為0.01的高斯噪聲來模擬噪聲環(huán)境,經(jīng)改進(jìn)加權(quán)引導(dǎo)濾波算法去噪后,分別使用FCM算法[17]、KFCM算法[18]、FLICM算法[19]和本文算法進(jìn)行圖像分割,分割結(jié)果如圖5~圖7所示.
圖5 切片1分割結(jié)果對(duì)比Fig.5 Comparison of segmentation results in section 1
圖6 切片2分割結(jié)果對(duì)比Fig.6 Comparison of segmentation results in section 2
圖7 切片3分割結(jié)果對(duì)比Fig.7 Comparison of segmentation results in section 3
從圖5、圖6和圖7可以看出,本文算法能夠較好地處理腦腫瘤的邊緣部分,且分割精度更高,分割效果最接近于真值圖像.為了避免主觀評(píng)價(jià)因素的影響,本文選擇 Dice系數(shù)、Jaccard系數(shù)、Sensitivity系數(shù)、Accuracy系數(shù)四個(gè)評(píng)價(jià)指標(biāo)對(duì)算法進(jìn)行評(píng)估,其定義分別如下:
(23)
(24)
(25)
(26)
其中:TP為正確分割的目標(biāo)區(qū)域;FP為誤分割的目標(biāo)區(qū)域;FN為沒有檢測(cè)到的目標(biāo)區(qū)域;TN為一定的非目標(biāo)區(qū)域.
仍然以上述三幅腦腫瘤MRI圖像為例,將本文算法同F(xiàn)CM算法、KFCM算法和FLICM算法進(jìn)行比較,比較結(jié)果如表2所示.
表2 本文算法與其他算法的分割性能對(duì)比
由表2可以看出,本文算法在四個(gè)指標(biāo)上均有較高的數(shù)值,相較傳統(tǒng)聚類算法有一定提升,Dice平均值高出傳統(tǒng)算法1.69%~2.45%,Jaccard平均值高出傳統(tǒng)算法3.03%~4.35%,Sensitivity平均值高出傳統(tǒng)算法4.65%~5.87%,Accuracy平均值高出傳統(tǒng)算法0.16%~0.23%.在腫瘤分割實(shí)驗(yàn)中,通常以Dice和Jaccard作為主要評(píng)價(jià)指標(biāo),其值越接近1,代表分割結(jié)果越精確.在本次實(shí)驗(yàn)中,本文算法的Dice平均值達(dá)到0.955 1,Jaccard平均值達(dá)到0.914 1.
本文提出的腦腫瘤圖像分割方法,利用改進(jìn)的加權(quán)引導(dǎo)濾波算法進(jìn)行圖像預(yù)處理,具有良好的保邊抗噪性;同時(shí),通過粒子群算法確定初始聚類中心,提高了算法的穩(wěn)定性;并且利用多個(gè)高斯核函數(shù)進(jìn)行數(shù)據(jù)映射,可以優(yōu)化算法的分割性能,最后通過馬爾科夫隨機(jī)場(chǎng)先驗(yàn)概率對(duì)算法的目標(biāo)函數(shù)進(jìn)行修正,進(jìn)一步增強(qiáng)了算法的抗噪性.實(shí)驗(yàn)結(jié)果表明,所提方法能在有效去噪的同時(shí),較大程度上保留圖像的邊緣信息,與FCM算法、KFCM算法和FLICM算法相比,分割精度更高,能將腦腫瘤從腦MRI圖像中準(zhǔn)確地分割出來,是一種有效的腦腫瘤圖像分割算法.但是,在高斯核的參數(shù)調(diào)節(jié)上,需要一定的先驗(yàn)知識(shí),限制了算法的靈活性,實(shí)現(xiàn)高斯核參數(shù)的自動(dòng)調(diào)節(jié),將是下一步研究的重點(diǎn).