基于自適應(yīng)聚類中心的腦血管分割方法

王 喆， 趙世鳳， 田 沄， 王學(xué)松， 周明全

(1.北京師范大學(xué) 信息科學(xué)與技術(shù)學(xué)院，北京 100875； 2.文化遺產(chǎn)數(shù)字化保護(hù)與虛擬現(xiàn)實(shí)北京市重點(diǎn)實(shí)驗(yàn)室，北京 100875)

0 引言

腦血管疾病是因顱內(nèi)血液循環(huán)障礙而造成腦組織損害的一組疾病，以急性發(fā)病居多，多表現(xiàn)為半身不遂、言語障礙等.每年有大量的病人忍受著腦血管問題所帶來的常見疾病，例如頸動(dòng)脈狹窄和動(dòng)脈瘤等疾病.腦血管疾病如不能被早期診斷，將會(huì)有引發(fā)腦中風(fēng)的危險(xiǎn).組織結(jié)構(gòu)的切分對(duì)輔助診斷、治療和手術(shù)規(guī)劃具有重要的意義.因此，腦血管的精確檢測與分割，尤其是低密度細(xì)小血管的提取，有助于腦血管疾病的輔助診斷與治療.

血管造影可以進(jìn)行體內(nèi)血管的無創(chuàng)檢查.基于MR的血管造影是顱內(nèi)血管評(píng)價(jià)中最常用的非侵入性成像技術(shù).磁共振血管成像(MRA)圖像提供了血管與背景之間相對(duì)出色的對(duì)比度和較高的空間分辨率，這為血管結(jié)構(gòu)的提取創(chuàng)造了較為豐富的前提條件[1].

從復(fù)雜數(shù)據(jù)集準(zhǔn)確的可視化和量化的病情狀況來看，分割是一個(gè)基本的步驟.國內(nèi)外學(xué)者從未間斷對(duì)腦血管提取算法的探索.到目前為止，已有許多2D或3D的血管分割方法.但絕大多數(shù)2D算法無法直接用于3D數(shù)據(jù)的分割.3D血管分割技術(shù)大致可分為兩類[2]:基于活動(dòng)輪廓模型的方法和基于概率統(tǒng)計(jì)的方法.基于活動(dòng)輪廓模型的方法的基本思想是在圖像力的作用下驅(qū)動(dòng)曲線到達(dá)物體的邊界.活動(dòng)輪廓模型通常分為參數(shù)活動(dòng)輪廓模型[3]和幾何活動(dòng)輪廓模型[4].幾何活動(dòng)輪廓模型可進(jìn)一步分為基于邊界[5-6]的模型、基于區(qū)域的模型[7-8]和混合模型[9-12].分割曲線的演化通常由最小化某種變分能量泛函來實(shí)現(xiàn).常用的數(shù)值計(jì)算此類活動(dòng)輪廓演化的方法是Osher等[13]提出的水平集方法.但此類方法需針對(duì)不同的圖像設(shè)置不同的參數(shù)才能獲取較好的分割效果，也即此類算法的性能依賴于所設(shè)計(jì)的能量函數(shù).

除了活動(dòng)輪廓模型，利用統(tǒng)計(jì)信息進(jìn)行圖像分割的方法一直是一個(gè)活躍的研究領(lǐng)域.統(tǒng)計(jì)模型將圖像數(shù)據(jù)中的體素分成血管區(qū)域和非血管區(qū)域[14].從MRA數(shù)據(jù)中提取血管的自適應(yīng)統(tǒng)計(jì)方法由Wong等[15]提出.Hassouna 等[16]也提出從時(shí)間飛躍法(TOF)中采用隨機(jī)模型方法進(jìn)行分割.Roy等[17]提出了利用基于期望最大化的瑞利混合模型的瑞利分類器進(jìn)行分割.實(shí)驗(yàn)結(jié)果表明，基于統(tǒng)計(jì)信息的概率模型在信噪比較低的圖像中提高了分割方法的魯棒性，但其準(zhǔn)確性依賴于底層概率模型的設(shè)計(jì).

血管形態(tài)結(jié)構(gòu)的復(fù)雜性使得單一的算法很難獲得滿意的分割效果，目前較為有效的方法是混合分割方法，這也是血管分割的一個(gè)研究趨勢.Gao等[18]提出了一種基于統(tǒng)計(jì)模型的快速活動(dòng)輪廓模型對(duì)腦血管進(jìn)行分割，同時(shí)為提取MRA數(shù)據(jù)集的三維腦血管提高了曲線演化的效率.Tian等[19]結(jié)合統(tǒng)計(jì)信息及血管形狀模型的向量場設(shè)計(jì)了適合血管的活動(dòng)輪廓模型進(jìn)行血管分割.混合模型算法實(shí)驗(yàn)結(jié)果是令人滿意的，但對(duì)于細(xì)小血管的提取依然有提升的空間.

筆者在傳統(tǒng)的統(tǒng)計(jì)方法基礎(chǔ)上，即采用有限混合模型對(duì)體數(shù)據(jù)密度進(jìn)行建模，采用期望最大進(jìn)行參數(shù)估計(jì)的方法，以后驗(yàn)概率的方式提取血管的主體部分，而對(duì)于血管的邊緣區(qū)域和細(xì)小的血管則首先計(jì)算其梯度值，采用基于改進(jìn)的K均值聚類方法對(duì)剩余體數(shù)據(jù)中的血管進(jìn)行進(jìn)一步提取.

1 本研究算法

筆者所提出的算法主要由兩個(gè)部分構(gòu)成.第一步由瑞利高斯混合模型提取血管的主體部分；第二步由改進(jìn)的K均值聚類算法進(jìn)一步提取血管邊界及細(xì)小血管部分.算法過程圖如圖1所示.

圖1 算法過程圖Fig.1 Process of the algorithm

1.1 基于有限混合模型的腦血管主體提取

傳統(tǒng)的基于統(tǒng)計(jì)的方法對(duì)時(shí)間飛躍法磁共振血管成像(TOF MRA)數(shù)據(jù)進(jìn)行血管分割時(shí)，首先根據(jù)TOF MRA數(shù)據(jù)密度分布特征將體素劃分為3個(gè)區(qū)域.低密度區(qū)域?qū)?yīng)于腦脊液、骨骼和背景；中密度區(qū)域?qū)?yīng)于腦組織，包括灰質(zhì)和白質(zhì)；高密度區(qū)域主要包括腦血管.采用高斯分布、瑞利分布或是統(tǒng)一分布(一般用于Phase-Contrast MRA數(shù)據(jù))分別對(duì)3種密度區(qū)域進(jìn)行建模，然后采用期望最大算法對(duì)模型進(jìn)行參數(shù)估計(jì)，最后將每一體素點(diǎn)的密度值帶入混合模型中，確定該點(diǎn)是否屬于血管部分的概率，從而實(shí)現(xiàn)對(duì)血管的分割.

筆者對(duì)血管主體部分的分割亦是采用傳統(tǒng)方式進(jìn)行.在這里，我們首先用瑞利高斯混合模型對(duì)體素進(jìn)行建模，如式(1)所示.

式中：wB、wT、wV分別為3類分布的權(quán)重，其中下標(biāo)B代表低密度區(qū)域，包括腦脊液、骨骼和背景(background)；T代表中密度區(qū)域，包括腦組織(tissue mixture)；V代表高密度區(qū)域，主要包括腦血管(vessel distribution).

R為瑞利分布：

(2)

Gk(k=1, 2)為高斯分布：

(3)

有限混合模型中包含一個(gè)瑞利、兩個(gè)高斯分布，分別對(duì)應(yīng)了低、中、高3個(gè)密度分布區(qū)域.高密度區(qū)域包含了目標(biāo)血管部分.由此可以將從體數(shù)據(jù)中分割血管的問題轉(zhuǎn)換為從有限混合模型判斷哪些體素屬于高密度區(qū)域所對(duì)應(yīng)高斯分布的問題.由于整體分布由幾個(gè)函數(shù)分布混合而成，無法直接得知每個(gè)體素所屬的分布，故而采用期望最大算法(EM)對(duì)分布函數(shù)的參數(shù)進(jìn)行估計(jì).

1.2 基于自適應(yīng)聚類中心的腦血管邊緣提取

完成血管主體部分的初步分割，我們發(fā)現(xiàn)剩余數(shù)據(jù)中含有部分未分割出的血管邊界和細(xì)小血管.由于剩余數(shù)據(jù)中包含血管的數(shù)據(jù)量較小，若再利用密度值進(jìn)行分割，存在較大困難.盡管剩余部分血管的數(shù)據(jù)量相對(duì)于體數(shù)據(jù)來說較小，但是目標(biāo)體素多出現(xiàn)于血管的邊界，因此將該部分?jǐn)?shù)據(jù)進(jìn)行梯度處理，以提高剩余部分血管邊界的對(duì)比度，從而進(jìn)行下一步分割.

首先將數(shù)據(jù)按垂直方向分層，分別在每一層計(jì)算數(shù)據(jù)梯度，以獲取梯度值數(shù)據(jù).因處于邊緣處血管的梯度值較大，為了提取出梯度數(shù)據(jù)中較亮的血管部分，筆者采用改進(jìn)的自適應(yīng)聚類的K均值方法，且是對(duì)梯度數(shù)據(jù)而非密度數(shù)據(jù)進(jìn)行聚類，從而找出目標(biāo)血管.

觀察到血管邊緣被包含于梯度值較高的部分，因此筆者考慮找到一個(gè)梯度值的分水嶺，比該分水嶺大的部分包含盡可能多的血管，同時(shí)比該分水嶺小的部分包含盡可能少的血管.筆者通過梯度值表示體素，并通過該梯度值的歐氏距離度量體素之間距離.由于隨著聚類中心數(shù)目的上升，分水嶺存在于K均值聚類中心值最大的類別中(盡管聚類中心數(shù)目過大會(huì)使得分水嶺從中心值最大的類別中脫離，但是在較少的類別時(shí)，分水嶺仍存在于中心值最大的類別中).因此，筆者通過自動(dòng)調(diào)節(jié)聚類數(shù)目使得聚類中心值最大的類別包含盡可能多的血管部分，同時(shí)包含盡可能少的非血管部分.

隨著聚類中心數(shù)目逐漸上升，聚類中心值最大的類別包含的非血管部分越來越少，同時(shí)在整體上，聚類中心的最大值逐漸增大.觀察到在上述過程中，當(dāng)聚類中心的最大值出現(xiàn)明顯上升時(shí)，擁有最大值的聚類中心表示的類別包含了較多的血管邊緣以及較少的非血管部分.因此，筆者通過調(diào)節(jié)聚類中心數(shù)目不斷增大并觀察聚類中心最大值的變化，從而識(shí)別聚類中心最大值的明顯上升.此時(shí)，該聚類中心的最大值可以作為上述分水嶺對(duì)血管邊緣進(jìn)行分割.

1.3 K均值聚類中心數(shù)目的確定

觀察數(shù)據(jù)后發(fā)現(xiàn)，剩余數(shù)據(jù)中血管部分的數(shù)值高于其他部分的數(shù)值，并且血管部分與非血管部分的數(shù)值大小差異相對(duì)于非血管部分的內(nèi)部差異更為明顯.上述差異給予提示可通過調(diào)節(jié)K均值聚類中心的數(shù)目，找到合適的聚類中心數(shù)使得最大的聚類中心作為分水嶺取得較好的分割效果.

如圖2所示，筆者的方法是逐步增加聚類中心點(diǎn)數(shù)目，通過觀察聚類中心最大值變化的差異來判定合適的聚類中心點(diǎn)數(shù)，相應(yīng)聚類中心的最大值即作為分割的分水嶺.K均值在N簇(cluster)中心點(diǎn)下進(jìn)行聚類時(shí)，血管部分往往包含于聚類中心較大的簇中.比如，當(dāng)N=2時(shí)，由K均值計(jì)算，將整體數(shù)據(jù)聚為兩類，即cluster1和cluster2，其聚類中心點(diǎn)分別是center1和center2.其中center1

圖2 基于K均值血管邊緣提取示意圖Fig.2 Boundary extraction of vessel with K-Means

由于血管數(shù)據(jù)在整體梯度數(shù)據(jù)中是較亮的部分，因此cluster2中包含血管數(shù)據(jù).但是由于僅僅將整體數(shù)據(jù)聚成兩類，在cluster2中仍包含了大量的非血管部分的數(shù)據(jù).因此可以將cluster2中的數(shù)據(jù)表示為nonvessel+vessel兩個(gè)部分，其中nonvessel是非血管數(shù)據(jù)，vessel是血管數(shù)據(jù).隨著聚類數(shù)目N的增大，在最大聚類中心所對(duì)應(yīng)的clusterN中nonvessel的比例逐漸減少，vessel的比例逐漸上升.從數(shù)據(jù)特征上看，非血管部分的數(shù)據(jù)值與血管部分?jǐn)?shù)據(jù)值存在較為明顯邊界，表現(xiàn)在數(shù)值上差異較大.因此，當(dāng)聚類數(shù)目N上升至Ntarget，相應(yīng)聚類中心點(diǎn)在數(shù)值上出現(xiàn)明顯的突增，最大中心對(duì)應(yīng)的簇中包含的可以將Ntarget對(duì)應(yīng)的聚類中心的最大值作為血管邊緣與非血管部分梯度值的分水嶺.

因此，尋找目標(biāo)值的整體算法步驟如下：nonvessel比例相對(duì)較小，而vessel比例相對(duì)較大，可以將Ntarget對(duì)應(yīng)的聚類中心的最大值作為血管邊緣與非血管部分梯度值的分水嶺.因此尋找目標(biāo)值的整體算法步驟如下.

輸入：密度梯度數(shù)據(jù).

輸出：血管邊緣與非血管部分梯度值分水嶺.

Step1：對(duì)梯度數(shù)據(jù)依次聚為2類，3類，…，20類，記錄每次聚類中心點(diǎn)的最大值的有序序列.

Step2：對(duì)Step1中序列的每個(gè)數(shù)據(jù)依次進(jìn)行兩兩求差，并記錄差值中的最大值.

Step3：Step2中的差值d=centerj-centeri,其中centeri和centerj屬于步驟Step1中聚類中心點(diǎn)最大值的集合且j=i+1.當(dāng)d為最大值時(shí)，centerj標(biāo)值.

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

實(shí)驗(yàn)所用數(shù)據(jù)是從海軍總醫(yī)院獲取的TOF MRA數(shù)據(jù)，大小為512 mm×512 mm×136 mm，空間分辨率為0.69 mm×0.69 mm×1.4 mm,TR=18 ms,TE=2.4 ms.測試平臺(tái)為普通PC機(jī)，Windows 7,64位操作系統(tǒng).

這里，選擇瑞利-高斯模型對(duì)體數(shù)據(jù)密度分布進(jìn)行建模分析.瑞利-高斯模型的估計(jì)參數(shù)如表1所示.圖3為提取結(jié)果.從圖3(a)可以看出，瑞利-高斯混合模型的誤差為0.179 4.統(tǒng)計(jì)方法獲取的血管主體部分結(jié)果如圖3(b)所示.圖3(c)是在血管主體部分基礎(chǔ)上，進(jìn)一步采用改進(jìn)自適應(yīng)聚類中心的K均值方法進(jìn)行血管邊緣提取后的效果圖.從圖3可以看出，傳統(tǒng)的統(tǒng)計(jì)方法所提取的血管主體部分相對(duì)較完整，但在血管邊緣和細(xì)小血管部分有不連續(xù)的分支存在.而改進(jìn)的方法在提取血管邊緣和細(xì)小血管分支的完整性及連續(xù)性方面則有較為明顯的提高.

表1 瑞利-高斯混合模型擬合參數(shù)Tab.1 Fitting parameters of Rayleigh-Gauss mixed model

在圖3中也給出了部分改善區(qū)域的標(biāo)記(藍(lán)色).從標(biāo)記的部分處，即在細(xì)小血管處，可以看出，基于改進(jìn)K均值梯度聚類的血管邊緣提取方法相比于基于參數(shù)估計(jì)的統(tǒng)計(jì)方法，對(duì)細(xì)小血管部分的提取效果在視覺上有明顯改善.

圖3 瑞利-高斯混合分布血管提取結(jié)果Fig.3 Results of Rayleigh-Gauss mixed model

圖4是兩種方法(瑞利-高斯模型方法與改進(jìn)的方法)的血管提取結(jié)果密度統(tǒng)計(jì)分布圖.其中藍(lán)色部分是第一步，即由基于參數(shù)估計(jì)的傳統(tǒng)統(tǒng)計(jì)方法獲取血管主體部分的分布情況.紅色部分是第二步，即在第一步的基礎(chǔ)上采用改進(jìn)K均值的血管邊緣提取出的結(jié)果分布情況.從圖4可以看出，傳統(tǒng)統(tǒng)計(jì)的方法最終的結(jié)果是找一個(gè)相對(duì)合適的全局目標(biāo)值，使得血管能夠盡可能地提取出來.但這個(gè)目標(biāo)值大多是基于數(shù)據(jù)的密度值的，因而最終的分割結(jié)果在密度分布圖上可以看到明顯的切分點(diǎn).

圖4 兩步血管提取密度統(tǒng)計(jì)分布圖Fig.4 Intensity distribution of two stages of vessel extraction

但因血管的密度范圍分布較廣，單純靠一個(gè)基于密度的目標(biāo)值很難獲得完整的血管結(jié)構(gòu).因此采用了基于梯度的方式來彌補(bǔ)由密度值分割所帶來的缺陷.因血管在邊緣處才會(huì)有較為明顯的梯度值，且經(jīng)過第一步提取后，血管主體部分已獲取，剩余部分是密度值較小且與周圍組織較為接近的邊緣及細(xì)小血管部分，只需對(duì)此部分重點(diǎn)分析提取即可.因此并非在整個(gè)體數(shù)據(jù)中計(jì)算梯度，而是在經(jīng)過第一步分割之后，對(duì)剩余部分的數(shù)據(jù)進(jìn)行梯度計(jì)算.從圖4中可以看出，經(jīng)過K均值梯度聚類后的結(jié)果，是對(duì)藍(lán)色部分的截面處進(jìn)行明顯的補(bǔ)充，而且此部分(紅色曲線)也更符合血管的密度值分布(范圍廣)，也使得最終結(jié)果的密度分布更加完整.

為了更好地評(píng)估算法分割效果的準(zhǔn)確率，采用公認(rèn)的相似性系數(shù)DSC來評(píng)價(jià)分割算法.

(4)

式中：M和G分別為算法結(jié)果與真實(shí)結(jié)果；N為血管體素個(gè)數(shù),當(dāng)M和G中的元素都相同時(shí)為1時(shí)，若M和G沒有任何相同元素則為0.同時(shí)，放射科專家提供了13套手動(dòng)切割的結(jié)果作為真實(shí)結(jié)果,數(shù)據(jù)如表2所示.

表2 時(shí)間飛躍法磁共振血管成像數(shù)據(jù)參數(shù)Tab.2 TOF MRA parameter

此外，筆者所提出的算法還與傳統(tǒng)的有限混合模型方法[16]以及當(dāng)前流行的基于水平集框架的活動(dòng)輪廓模型[19]方法的分割效果進(jìn)行了對(duì)比.圖5給出了其中一套數(shù)據(jù)的3種方法的分割結(jié)果.從圖5可以看出，對(duì)于大的血管，3種方法都能得到較為滿意的結(jié)果.但對(duì)于細(xì)小血管，傳統(tǒng)的有限混合模型出現(xiàn)血管分支間斷現(xiàn)象，這是由于細(xì)小血管部分的密度值相對(duì)較低，在整個(gè)體數(shù)據(jù)中，采用唯一的判別準(zhǔn)則很難將此部分提取出來.而活動(dòng)輪廓模型，將各種因素融合設(shè)計(jì)新的能量泛函，利用水平集求解能量函數(shù)，能提取出部分細(xì)小血管的分支，如圖5(b)圓角矩形框所示.而筆者所提出的改進(jìn)的方法，在細(xì)小血管的提取方面相對(duì)于統(tǒng)計(jì)模型也有不同程度的提高，而且在血管的邊緣處，能看到更為明顯的效果，此外改進(jìn)的方法在獲取細(xì)小血管的連續(xù)性方面效果更好.

圖5 3種方法的分割效果Fig.5 Extraction results comparison among three method

圖6給出了13套數(shù)據(jù)的3種方法的DSC的對(duì)比結(jié)果.筆者所提方法的DSC平均在80% 以上，有的數(shù)據(jù)甚至達(dá)到了94%，比傳統(tǒng)統(tǒng)計(jì)模型的結(jié)果平均高出10.7%，比活動(dòng)輪廓模型的結(jié)果高出8.5%.傳統(tǒng)有限混合模型僅基于數(shù)據(jù)的密度值而并沒有考慮像一階導(dǎo)數(shù)等其他因素，且依賴于混合模型的選擇.活動(dòng)輪廓模型將密度與梯度等參數(shù)融合至統(tǒng)一的能量函數(shù)，其分割結(jié)果要比有限混合模型的效果好，尤其是在細(xì)小血管分支的連續(xù)性方面有較為明顯的提高.但活動(dòng)輪廓模型對(duì)整個(gè)體數(shù)據(jù)同時(shí)應(yīng)用密度和梯度值.而筆者的算法是在前一步密度處理結(jié)果的基礎(chǔ)上運(yùn)用梯度進(jìn)行血管邊緣的提取，因此在細(xì)小血管分支的連續(xù)性及血管的邊緣方面優(yōu)于活動(dòng)輪廓模型.

圖6 3種方法相似性系數(shù)DSC的比較Fig.6 DSC comparison among three methods

3 結(jié)論

筆者在有限混合模型的血管分割基礎(chǔ)上，提出了一種基于自適應(yīng)聚類中心的K均值梯度聚類方法.該方法首先使用瑞利-高斯混合模型的方法對(duì)時(shí)間飛躍法磁共振血管成像數(shù)據(jù)建模，并采用期望最大參數(shù)估計(jì)的方法獲取混合模型的參數(shù)，進(jìn)而利用后驗(yàn)概率提取出血管的主體部分；隨后對(duì)剩余體素?cái)?shù)據(jù)進(jìn)行梯度化處理，并采用改進(jìn)的K均值聚類方法在剩余梯度數(shù)據(jù)中提取出血管的邊緣以及細(xì)小血管部分，將兩部分結(jié)果合并進(jìn)行可視化處理便可得到血管的三維結(jié)構(gòu).實(shí)驗(yàn)結(jié)果表明，筆者進(jìn)一步采用基于改進(jìn)自適應(yīng)聚類中心的梯度聚類方法，不管從視覺上還是精確度上，分割結(jié)果在血管的細(xì)節(jié)部分都有較為明顯的改善.在下一步的研究工作中，筆者將融合更多血管幾何形態(tài)特征，進(jìn)一步提高血管分割的精度.