王 喆， 趙世鳳， 田 沄， 王學(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)為半身不遂、言語(yǔ)障礙等.每年有大量的病人忍受著腦血管問(wèn)題所帶來(lái)的常見(jiàn)疾病，例如頸動(dòng)脈狹窄和動(dòng)脈瘤等疾病.腦血管疾病如不能被早期診斷，將會(huì)有引發(fā)腦中風(fēng)的危險(xiǎn).組織結(jié)構(gòu)的切分對(duì)輔助診斷、治療和手術(shù)規(guī)劃具有重要的意義.因此，腦血管的精確檢測(cè)與分割，尤其是低密度細(xì)小血管的提取，有助于腦血管疾病的輔助診斷與治療.

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

從復(fù)雜數(shù)據(jù)集準(zhǔn)確的可視化和量化的病情狀況來(lái)看，分割是一個(gè)基本的步驟.國(guó)內(nèi)外學(xué)者從未間斷對(duì)腦血管提取算法的探索.到目前為止，已有許多2D或3D的血管分割方法.但絕大多數(shù)2D算法無(wú)法直接用于3D數(shù)據(jù)的分割.3D血管分割技術(shù)大致可分為兩類(lèi)[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].分割曲線的演化通常由最小化某種變分能量泛函來(lái)實(shí)現(xiàn).常用的數(shù)值計(jì)算此類(lèi)活動(dòng)輪廓演化的方法是Osher等[13]提出的水平集方法.但此類(lèi)方法需針對(duì)不同的圖像設(shè)置不同的參數(shù)才能獲取較好的分割效果，也即此類(lèi)算法的性能依賴(lài)于所設(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]提出了利用基于期望最大化的瑞利混合模型的瑞利分類(lèi)器進(jìn)行分割.實(shí)驗(yàn)結(jié)果表明，基于統(tǒng)計(jì)信息的概率模型在信噪比較低的圖像中提高了分割方法的魯棒性，但其準(zhǔn)確性依賴(lài)于底層概率模型的設(shè)計(jì).

血管形態(tài)結(jié)構(gòu)的復(fù)雜性使得單一的算法很難獲得滿意的分割效果，目前較為有效的方法是混合分割方法，這也是血管分割的一個(gè)研究趨勢(shì).Gao等[18]提出了一種基于統(tǒng)計(jì)模型的快速活動(dòng)輪廓模型對(duì)腦血管進(jìn)行分割，同時(shí)為提取MRA數(shù)據(jù)集的三維腦血管提高了曲線演化的效率.Tian等[19]結(jié)合統(tǒng)計(jì)信息及血管形狀模型的向量場(chǎng)設(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均值聚類(lèi)方法對(duì)剩余體數(shù)據(jù)中的血管進(jìn)行進(jìn)一步提取.

1 本研究算法

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

圖1 算法過(guò)程圖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類(lèi)分布的權(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ù)中分割血管的問(wèn)題轉(zhuǎn)換為從有限混合模型判斷哪些體素屬于高密度區(qū)域所對(duì)應(yīng)高斯分布的問(wèn)題.由于整體分布由幾個(gè)函數(shù)分布混合而成，無(wú)法直接得知每個(gè)體素所屬的分布，故而采用期望最大算法(EM)對(duì)分布函數(shù)的參數(shù)進(jìn)行估計(jì).

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

完成血管主體部分的初步分割，我們發(fā)現(xiàn)剩余數(shù)據(jù)中含有部分未分割出的血管邊界和細(xì)小血管.由于剩余數(shù)據(jù)中包含血管的數(shù)據(jù)量較小，若再利用密度值進(jìn)行分割，存在較大困難.盡管剩余部分血管的數(shù)據(jù)量相對(duì)于體數(shù)據(jù)來(lái)說(shuō)較小，但是目標(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)聚類(lèi)的K均值方法，且是對(duì)梯度數(shù)據(jù)而非密度數(shù)據(jù)進(jìn)行聚類(lèi)，從而找出目標(biāo)血管.

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

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

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

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

如圖2所示，筆者的方法是逐步增加聚類(lèi)中心點(diǎn)數(shù)目，通過(guò)觀察聚類(lèi)中心最大值變化的差異來(lái)判定合適的聚類(lèi)中心點(diǎn)數(shù)，相應(yīng)聚類(lèi)中心的最大值即作為分割的分水嶺.K均值在N簇(cluster)中心點(diǎn)下進(jìn)行聚類(lèi)時(shí)，血管部分往往包含于聚類(lèi)中心較大的簇中.比如，當(dāng)N=2時(shí)，由K均值計(jì)算，將整體數(shù)據(jù)聚為兩類(lèi)，即cluster1和cluster2，其聚類(lèi)中心點(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ù)聚成兩類(lèi)，在cluster2中仍包含了大量的非血管部分的數(shù)據(jù).因此可以將cluster2中的數(shù)據(jù)表示為nonvessel+vessel兩個(gè)部分，其中nonvessel是非血管數(shù)據(jù)，vessel是血管數(shù)據(jù).隨著聚類(lèi)數(shù)目N的增大，在最大聚類(lèi)中心所對(duì)應(yīng)的clusterN中nonvessel的比例逐漸減少，vessel的比例逐漸上升.從數(shù)據(jù)特征上看，非血管部分的數(shù)據(jù)值與血管部分?jǐn)?shù)據(jù)值存在較為明顯邊界，表現(xiàn)在數(shù)值上差異較大.因此，當(dāng)聚類(lèi)數(shù)目N上升至Ntarget，相應(yīng)聚類(lèi)中心點(diǎn)在數(shù)值上出現(xiàn)明顯的突增，最大中心對(duì)應(yīng)的簇中包含的可以將Ntarget對(duì)應(yīng)的聚類(lèi)中心的最大值作為血管邊緣與非血管部分梯度值的分水嶺.

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

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

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

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

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

Step3：Step2中的差值d=centerj-centeri,其中centeri和centerj屬于步驟Step1中聚類(lèi)中心點(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.測(cè)試平臺(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)聚類(lèi)中心的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均值梯度聚類(lèi)的血管邊緣提取方法相比于基于參數(shù)估計(jì)的統(tǒng)計(jì)方法，對(duì)細(xì)小血管部分的提取效果在視覺(jué)上有明顯改善.

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

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

(4)

式中：M和G分別為算法結(jié)果與真實(shí)結(jié)果；N為血管體素個(gè)數(shù),當(dāng)M和G中的元素都相同時(shí)為1時(shí)，若M和G沒(méi)有任何相同元素則為0.同時(shí)，放射科專(zhuān)家提供了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)則很難將此部分提取出來(lái).而活動(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ù)的密度值而并沒(méi)有考慮像一階導(dǎo)數(shù)等其他因素，且依賴(lài)于混合模型的選擇.活動(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)聚類(lèi)中心的K均值梯度聚類(lèi)方法.該方法首先使用瑞利-高斯混合模型的方法對(duì)時(shí)間飛躍法磁共振血管成像數(shù)據(jù)建模，并采用期望最大參數(shù)估計(jì)的方法獲取混合模型的參數(shù)，進(jìn)而利用后驗(yàn)概率提取出血管的主體部分；隨后對(duì)剩余體素?cái)?shù)據(jù)進(jìn)行梯度化處理，并采用改進(jìn)的K均值聚類(lèi)方法在剩余梯度數(shù)據(jù)中提取出血管的邊緣以及細(xì)小血管部分，將兩部分結(jié)果合并進(jìn)行可視化處理便可得到血管的三維結(jié)構(gòu).實(shí)驗(yàn)結(jié)果表明，筆者進(jìn)一步采用基于改進(jìn)自適應(yīng)聚類(lèi)中心的梯度聚類(lèi)方法，不管從視覺(jué)上還是精確度上，分割結(jié)果在血管的細(xì)節(jié)部分都有較為明顯的改善.在下一步的研究工作中，筆者將融合更多血管幾何形態(tài)特征，進(jìn)一步提高血管分割的精度.