基于Markov隨機場的腦部三維磁共振血管造影數(shù)據(jù)的分割

周壽軍 賈富倉 胡慶茂 謝耀欽 辜 嘉 尚 鵬

(中國科學(xué)院深圳先進技術(shù)研究院深圳518055)

基于Markov隨機場的腦部三維磁共振血管造影數(shù)據(jù)的分割

周壽軍 賈富倉 胡慶茂 謝耀欽 辜 嘉 尚 鵬

(中國科學(xué)院深圳先進技術(shù)研究院深圳518055)

文章提出了腦部核磁共振血管造影(Magnetic Resonance Angiography,MRA)的全自動分割方法，該方法有效增強了現(xiàn)有的基于Markov隨機場(Markov Random Field, MRF)的分割技術(shù)?，F(xiàn)有的三維Markov分割模型通常面臨的挑戰(zhàn)是：(1)低級MRF模型參數(shù)初始化不夠準(zhǔn)確；(2)普通的MRF 鄰域系統(tǒng)無法探測精細(xì)的血管結(jié)構(gòu)。針對這兩類問題，分別提出了基于多尺度濾波響應(yīng)閾值分析和多模式鄰域系統(tǒng)進行解決，使得MRF模型的血管分辨率提高到2個體素的細(xì)小血管。實驗中，低級模型參數(shù)的精確估計采用了最大期望算法，高階MRF參數(shù)的估計采用最大偽似然估計方法；通過三維仿真數(shù)據(jù)和實際腦部MRA數(shù)據(jù)進行驗證，分割結(jié)果顯示了較小的全局誤差。

血管分割；磁共振血管造影；馬爾可夫隨機場；鄰域系統(tǒng)

1 引 言

每年有大量患者受到由血栓、出血、動脈畸形、狹窄導(dǎo)致的中風(fēng)等腦血管疾病的影響，因此，精確地分割腦血管對于診斷圖像引導(dǎo)的神經(jīng)外科手術(shù)、術(shù)前計劃等有著重要的臨床意義。由于腦血管具有復(fù)雜的三維形態(tài)和較細(xì)的結(jié)構(gòu)，所以其圖像分割很有挑戰(zhàn)性[1]。

核磁共振血管造影(Magnetic Resonance Angiography, MRA)是一種非侵入式的流成像和三維體數(shù)據(jù)生成技術(shù)。有很多方法可以用來進行MRA的分割，這些方法可分為基于骨架和非骨架兩類方法。前者分割和重建血管通過先探測血管中心線，而后者則直接提取三維血管。最近，基于非骨架的方法受到大量關(guān)注，并可進一步劃分為尺度分析法[2,3]、形變模型法[4,5]、統(tǒng)計模型法[6-10]以及混合模型法[7]。依據(jù)某種MRA模態(tài)中各類解剖結(jié)構(gòu)的特性，統(tǒng)計模型法將一些恰當(dāng)?shù)倪吘壏植己瘮?shù)進行線性合并，建模基本要素為：(1)各種混合分布的概率函數(shù)構(gòu)造；(2)估計統(tǒng)計學(xué)模型的參數(shù)；(3)分別獲取目標(biāo)和背景的概率分布；(4)利用最大后驗概率(Maximum A Posteriori Probability, MAP)分類技術(shù)獲取全局閾值函數(shù)。本文將結(jié)合統(tǒng)計模型的基本算法進行分析，并克服其存在的不足。統(tǒng)計學(xué)方法將三維MRA數(shù)據(jù)集中不同解剖結(jié)構(gòu)的感興趣區(qū)看作不同類別數(shù)據(jù)，每一類服從一種特殊的邊緣概率分布模式。Wilson等人[6]提出了一種經(jīng)典的基于統(tǒng)計學(xué)模型的血管提取技術(shù)，能夠自適應(yīng)地從MRA中分割血管，其分別利用了兩個高斯的和一個均勻的邊緣概率分布的混合函數(shù)模擬了靜態(tài)腦脊液、腦組織和動脈血管的分布。為了估計混合模型的參數(shù)，采用了傳統(tǒng)的最大期望(Expectation Maximization,EM)算法。此外，Xin等人[7]還提出一種利用區(qū)域動態(tài)輪廓結(jié)合信號統(tǒng)計分布和形狀信息的方法分割 MRA 中的血管。為增強單純的統(tǒng)計學(xué)模型的分割效果，統(tǒng)計學(xué)分割方法[8,9]在空間信息約束方面的不足正逐漸為Markov隨機場(Markov Random Field,MRF)鄰域建模技術(shù)的增長所解決。

與基本的統(tǒng)計學(xué)方法類似，本文將采用統(tǒng)計學(xué)模型和隨機場理論解決腦血管MRA數(shù)據(jù)集的分割問題。為了精確量化三維血管結(jié)構(gòu)并減少計算時間，我們提出的MRF分割模型與現(xiàn)有方法的主要區(qū)別在于：(1)低級MRF模型的血管參數(shù)估計采用了新的多尺度濾波和閾值分析技術(shù)；(2)采用了一種二階多模式鄰域(Multi-Pattern-NBS,MPN)系統(tǒng)構(gòu)造高級MRF模型。同時，在算法執(zhí)行當(dāng)中利用了并行技術(shù)和血管試探空間法極大地減少了計算量和節(jié)省了計算時間。本文采用的方法框架如圖1所示。

圖1 基于3D-MRF的腦部MRA分割框架Fig.1. The segmentation framework of brain MRA dataset based on 3D-MRF

本文的主要方法介紹如下：首先，采用多尺度血管濾波處理原始圖像數(shù)據(jù)，以增強血管和抑制噪聲，期間依據(jù)血管和頭部容積的經(jīng)驗比值確定一個恰當(dāng)?shù)难荛撝郸?；其次，通過直方圖自動分析確定多類邊緣分布的初始參數(shù)，并利用 EM算法估計該參數(shù)的準(zhǔn)確值；然后，利用一個更低的閾值γ2處理多尺度血管濾波響應(yīng)后，得到原始數(shù)據(jù)中血管的候選空間Sγ1，在此空間利用MRF低級模型和高級模型迭代計算時，隨機場的初始化采用極大似然(Maximum Likelihood,ML)估計的結(jié)果，結(jié)合MPN系統(tǒng)和對應(yīng)的能量函數(shù)，便可以探測血管目標(biāo)。最終，采用MRF迭代條件模式進行分割的迭代計算。上述各步驟的計算是全自動的，而且采用并行 ML 估計、期望最大化(EM)估計、血管試探空間以減少時間消耗。

2 基于MRF的統(tǒng)計學(xué)分割模型

就我們所知，最初利用MRF模型從MRA數(shù)據(jù)中有效分割出血管的研究者是Hassouna等[8]。MRF實現(xiàn)空間一點和其鄰域的交互作用，其局部相關(guān)性提供了模仿各種圖像特性的機制，其觀測數(shù)據(jù)集分別由低級和高級兩個隨機過程模擬，以此刻畫了一種基于灰度和鄰域相關(guān)性的統(tǒng)計學(xué)特性。這兩個隨機過程都定義在空間域S(即MRA數(shù)據(jù)空間)之中。令代表一套觀測隨機變量，Ys為體素s處的灰度。假使為一個Markov隨機場，Xs從標(biāo)記集中取一值LV或LB(分別代表血管和背景類標(biāo)記)。因此，已知一套觀測特征矢量Y(如果認(rèn)為三維數(shù)據(jù)集中每一體素的灰度值是條件獨立的)，并且鄰域標(biāo)記信息可以由Markov先驗概率P(X)表示，則MRF分割模型可以由最大后驗概率形式表示為：

方程右邊兩項代表低級和高級MRF過程：前者可由多類邊緣分布的混合概率表示；后者可由Hammersley-Clifford定理定義的Gibbs 分布表示為：

上式能量函數(shù)表示為U(x)，標(biāo)記x在特定配置(0或1)下能量越高則出現(xiàn)的概率越低。我們將可能的標(biāo)記表示為X。上式中分子Z稱做隔離函數(shù)，是一個歸一化的常數(shù)，其值為分子在各種配置下的和數(shù)。下面將分別討論低級和高級MRF分割模型的構(gòu)造原理。

2.1 MRF低級模型

對于腦部MRA圖像分割來說，有限混合(Finite Mixture,FM)模型[11,12]使用較普遍，其中所有數(shù)據(jù)點被認(rèn)為群體中的獨立樣本，而并不涉及空間鄰域信息。因此FM模型通常僅適用于較低噪聲條件下的圖像應(yīng)用，盡管如此，Wilson等[6]采用局部和全局直方圖分析方法成功彌補了FM模型的不足，其FM模型參數(shù)利用了EM算法得到準(zhǔn)確的估計結(jié)果。進行MRF低級模型的構(gòu)建時，將MRA數(shù)據(jù)直方圖分為低、中、高三個區(qū)域[8]；低灰度區(qū)間包括了腦脊液、骨、背景空氣；中等灰度區(qū)間包括了腦組織的灰/白質(zhì)和部分眼組織；高灰度區(qū)間則包括了動脈血管和皮下脂肪組織。通過測試各種概率密度函數(shù)在臨床數(shù)據(jù)集中的效果，發(fā)現(xiàn)更為精確的直方圖擬合分布函數(shù)的組成是：利用一個瑞利分布和兩個高斯分布模擬中低灰度區(qū)間的背景類；而在高灰度區(qū)間用一個高斯分布模擬血管類。因此，混合概率密度函數(shù)表示為：

圖2 利用混合邊緣概率分布(1個瑞利分布和3個高斯分布)進行 MRA 直方圖曲線擬合Fig.2. Histogram curve- fi tting of TOF-MRA using mixture distributions (one Rayleigh and three Gaussian distributions)

在應(yīng)用MAP算法估計時，公式(3)和(4)共有11個參數(shù)需要估計。一種自動的估計方法是基于直方圖技術(shù)分析這些混合分布的初始參數(shù)，在此基礎(chǔ)上Hassouna等[8]提出了利用EM算法進行迭代更新的精確估計方法。估計前后的曲線擬合效果比較情況如圖2所示。由圖2(a)可見，利用初始的參數(shù)經(jīng)擬合之后不能與原始直方圖重合，在中等強度區(qū)間存在較大的差異；將初始參數(shù)作為EM迭代更新方程的初值，優(yōu)化估計后的參數(shù)能夠擬合原始直方圖，如圖2(b)。

2.2 MRF高級模型

由于MRF低級模型(即FM模型)并未考慮空間鄰域信息，它作為似然估計過程在良好的低噪聲圖像條件下可以得到較好結(jié)果，而對于MRA 這一類背景復(fù)雜和噪聲多的圖像分割是無能為力的。為克服這一問題，必須采用一種較好的空間鄰域系統(tǒng)構(gòu)造優(yōu)化的高級MRF模型。傳統(tǒng)的空間鄰域系統(tǒng)是一種基于MRF多級邏輯模型(MLL)的對點交互式的二階鄰域系統(tǒng)[8,9]，它建立在局部3×3×3大小的立方體上，其勢團能量由周圍26鄰點的能量對之和表示?；贛LL的高級MRF能量模型表示為：

其中函數(shù)為點s的高級先驗勢能：對應(yīng)于鄰域空間中血管和背景類的對點交互的能量，該能量等價于鄰域中體素標(biāo)記相同的對點數(shù)的參數(shù)加權(quán)和，與此點的先驗概率成反比，可以表示為：

參數(shù)表示鄰域中點對的交互作用強度，通常保持其在所有點對中相同。

3 基于多模式鄰域系統(tǒng)的MRF分割模型

在噪聲、偽影、強度不均勻的MRA數(shù)據(jù)集中，對于微小的3×3×3空間來說，傳統(tǒng)的基于6鄰域和26鄰域的MRF三維模型對于目標(biāo)管體、管壁、末梢等微觀結(jié)構(gòu)傾向于產(chǎn)生相同的高級能量，為此，傳統(tǒng)模型很難辨別上述微觀結(jié)構(gòu)。本節(jié)下述內(nèi)容首先提出新的MRF分割目標(biāo)函數(shù)，而后定義多種模式的勢團結(jié)構(gòu)，在此基礎(chǔ)上給出試探法參數(shù)估計方案。由此提高模型的分辨率，尤其是微細(xì)結(jié)構(gòu)探測能力。

3.1 MRF目標(biāo)函數(shù)

通過低級MRF分割模型進行極大似然估計，可以得到初始標(biāo)記場，體素s的標(biāo)記值xs可以利用以下Baysian方程估計：

分別利用 LV和LB表示血管和背景標(biāo)記，上述MRF可以由Ising模型定義，血管和背景的后驗概率可以分別表示為：

參照Baysian準(zhǔn)則，如果一個體素xs滿足便可判斷其為血管點，反之則為背景點。條件概率和的寫為：

3.2 預(yù)處理和血管候選空間的獲取

由于三維MRF隨機場的計算非常耗時，我們采用多尺度血管濾波增強技術(shù)對原始數(shù)據(jù)集進行預(yù)處理，進而獲得血管候選空間(VCS)——該空間位于直方圖高強度區(qū)域，通常僅包括血管和皮下脂肪組織。由于VCS遠(yuǎn)小于頭顱容積，因此可以極大減少計算量。方法如下：

由于基于Hessian矩陣的多尺度濾波技術(shù)能夠增強管狀目標(biāo)特征[2,3]，它已被廣泛應(yīng)用于醫(yī)學(xué)圖像處理領(lǐng)域；我們只需在不同尺度用一個二階高斯導(dǎo)數(shù)核G卷積三維原始數(shù)據(jù)Iimg得到各點的Hessian特征值，再利用Law等[3]提出的多尺度血管函數(shù)進行濾波，便可以獲得多尺度血管增強結(jié)果Rf。后續(xù)采用閾值分析手段便可獲取血管候選空間：首先，利用閾值處理多尺度濾波響應(yīng)Rf，可以得到不同的偽血管空間Ωn；其次，計算各閾值對應(yīng)空間的均值和方差并與血管經(jīng)驗參數(shù)進行比較，當(dāng)滿足下式時，可以認(rèn)為找到了臨界點nt和閾值 γt。

3.3 空間鄰域系統(tǒng)

大多數(shù)基于MRF分類和分割的應(yīng)用中，利用各向同性和均勻性構(gòu)造MRF模型與Gibbs分布函數(shù)[8,9]。比如當(dāng)獨立于勢團c在中的相對位置時，我們認(rèn)為MRF是均勻的；當(dāng)獨立于勢團c的方向時，我們認(rèn)為MRF是各向同性的。勢團的空間布局模式需依據(jù)鄰域系統(tǒng)的定義，比如：圖像分析中，許多均勻的區(qū)域表現(xiàn)為集中大量相同標(biāo)記的勢團模式；而非均勻區(qū)域傾向于集中大量不同標(biāo)記的勢團模式。本文中，我們在3×3×3空間上定義一個二階MPN系統(tǒng)，以此提供一種均勻的各向同性的Markov先驗?zāi)Ｐ?。該二階MPN系統(tǒng)包括3類勢團模式(CP)，如圖3(a)～3(c)的CP，分別顯示了點s周圍鄰點數(shù)Nm＝6,9,7的三種代表情況：每種CP模式可以在鄰點數(shù)不便的條件下，以立方體的中心和中位面派生出布局對稱的同類模式。圖3(a)表示一點s位于血管體內(nèi)部；圖3(b)表示一點s位于血管管壁的個例；圖3(c)表示一點s位于血管末梢的個例。以上個例對稱拓展后將可以表示血管的各種局部形態(tài)。因此，MPN系統(tǒng)提供了完整的20種相互形態(tài)獨立的CP模式。我們還可以定義其它類型的CP，考慮到噪聲干擾和CP間可能的重復(fù)性，更多的CP模式并非能夠達(dá)到最佳效果。

其中(a)～(c)定義了3類CP，每一類參照立方體的中心和中位面分別可以衍生出其它的CP。其中(a)類僅此一種情況；(b)類可對稱拓展至6種情況；(c)類可對稱拓展至8種情況。

其中，為勢團Cm的總模式數(shù)，有；；。上式反映了一點s的高級MRF能量，兩個求和符號分別作用于類和模式。為了實現(xiàn)各向同性血管分割，我們采用；對于每一勢團且，能量方程遵循公式(12)：如果多有的CP中xs＝xr＝LV，式(14)提供高級 MRF 模型的血管能量；反之xs＝xr＝LB時，式(14)提供高級MRF 模型的背景能量；當(dāng) xs≠xr時，點s的高級MRF能量為。總之，以上CP定義了一種均勻的各向同性的Ising-MRF 勢團模型，可以分辨小至 2 個體素的目標(biāo)直徑。

圖3 為三維MRF高級模型定義二階MPN系統(tǒng)Fig.3. The second-order MPN system is def i ned for the high level model of 3D-MRF

3.4 高級MRF的參數(shù)估計

現(xiàn)有的各種高級模型的參數(shù)大多采用試探誤差分析方法，本文利用最大偽似然算法進行試探誤差分析，并采用空間編碼模式將分割空間拆分為 k＝1，2 兩部分，并分別計算高階模型參數(shù)后求平均值。依據(jù)上述Ising-MRF模型，其似然概率為用來初步判斷血管，即：參照Bayesian規(guī)則和MRF理論，我們可以在 S(k)空間構(gòu)造偽似然能量函數(shù)的條件概率，并用乘積形式表示為：

3.5 迭代求解

我們采用最普遍使用的迭代條件模式(ICM)進行上述Ising-MRF分割模型求解。ICM采用貪婪策略最大化局部后驗概率，主要過程為：算法第n步，在已知觀測數(shù)據(jù)ys和標(biāo)記場的條件下，算法利用最大后驗概率順序從更新，初始標(biāo)記場的設(shè)置采用最大似然估計，當(dāng)?shù)螖?shù)等于指定數(shù)值N或者所有標(biāo)記不在變化時算法終止，如表1所示。

4 實驗與討論

實驗采用多套臨床MRA數(shù)據(jù)進行腦血管的分割提取，其效果可與三維MIP投影數(shù)據(jù)進行血管形態(tài)對比，然而具體的量化評估依然需要仿真數(shù)據(jù)實驗，用手工分割臨床MRA數(shù)據(jù)獲得的標(biāo)準(zhǔn)腦血管樹進行算法驗證是不現(xiàn)實的，原因取決于于血管形態(tài)結(jié)構(gòu)的復(fù)雜性、主觀判斷的差異性以及人力和經(jīng)驗的匱乏性。因此實驗采用合成的3D仿真數(shù)據(jù)進行實驗確認(rèn)。類似于Hassouna等[8]所用的方法，仿真體膜中包含了目標(biāo)和背景：目標(biāo)的形狀可以是彎曲的圓管、圓柱或方柱；目標(biāo)和背景的混合比例可以按照實際腦血管和頭顱容積率之比進行調(diào)整，而背景數(shù)據(jù)混合的比例按照實際MRA中的背景分布調(diào)整?？傊畼?gòu)造的仿真體包含了4種分布模式(強度分別服從1個瑞利分布和3個高斯分布)，混合比例和分布參數(shù)可以參照臨床數(shù)據(jù)直方圖人工指定。其二維切層圖像如圖4。

表1 迭代算法Table1. Iteration algorith

圖4 仿真體數(shù)據(jù)截層圖像Fig.4. The cross section of phantom volume data

4.1 仿真實驗

我們首先采用仿真體數(shù)據(jù)驗證MPN系統(tǒng)優(yōu)于傳統(tǒng)6鄰域和26鄰域系統(tǒng)。每一種體數(shù)據(jù)具有任意類型的噪聲、紋理和混合參數(shù)。分割誤差表示為誤分類像素占三維仿真體數(shù)據(jù)容積的百分比。仿真體數(shù)據(jù)中，最小直徑1～3像素，數(shù)據(jù)截層如圖4所示。

實驗中，對于不同的目標(biāo)直徑和曲率，所提出的MPN-MRF分割模型的計算時間較傳統(tǒng)26鄰域和6鄰域系統(tǒng)有所增加，主要是因為統(tǒng)計預(yù)處理、MLE、ICM過程的耗時，最大偽似然參數(shù)估計時間消耗與此過程相當(dāng)。結(jié)果表明：利用MPN-MRF方法分割兩套仿真數(shù)據(jù)時的誤差率分別為0.19%和0.20%，實現(xiàn)了更低的分割誤差，較6鄰域MRF分割模型的平均誤差減少近1.5倍；較26鄰域MRF分割模型的平均誤差減少10多倍；相對MLE的平均誤差20.75%～28.71%來說，誤差率降低了約100倍。各仿真體數(shù)據(jù)的誤差統(tǒng)計和參數(shù)估計結(jié)果如表2所示。

圖5(a)和(b)分別為真實目標(biāo)和和基于MPN鄰域系統(tǒng)的分割結(jié)果；(c)～(e)分別為對應(yīng)某一切層的MLE結(jié)果、真實目標(biāo)和MPN-MRF分割結(jié)果。

表2 各種鄰域模型的參數(shù)估計和分割誤差統(tǒng)計Table2. Statistical results of the parameter estimation and segmentation error with respect to various NBS

圖5 MRF分割模型效果比較Fig.5. Comparisons of the segmentation results

4.2 MRA臨床數(shù)據(jù)分割

在實際臨床數(shù)據(jù)的分割方面，我們采用的是MIDAS平臺數(shù)據(jù)庫中提供的10套患者和受試者MRA數(shù)據(jù)，空間分辨率0.5×0.5×0.8 mm3；數(shù)據(jù)的像素空間大小為448×448×128。算法編程采用MATLAB2010b，主機硬件：Intel(R)Core(TM)i5、CPU3.2GHz、8G-DDR-RAM和Win7-64Bit操作系統(tǒng)。圖6比較顯示了4套分割結(jié)果，圖中的血管三維表面渲染采用Matlab的isosurface.m和patch.m函數(shù)。

圖6 臨床MRA數(shù)據(jù)的分割結(jié)果比較Fig.6. Comparison of segmentation results of clinical MRA dataset

實驗對每一套數(shù)據(jù)采用6鄰域、26鄰域和本文提出的MPN鄰域系統(tǒng)進行比較，對大量個體化MRA數(shù)據(jù)能夠無一例外地準(zhǔn)確估計低級模型參數(shù)，高級模型參數(shù)統(tǒng)一采用最大偽似然方程(16)進行計算，三種MRF鄰域模型的高級能量統(tǒng)一采用公式(14)進行計算。圖6顯示了4套MRA數(shù)據(jù)的分割結(jié)果。

圖6中，a1～d1為原始MRA的軸位MIP；a2～d2多尺度血管濾波后的軸位MIP；a3～d3為血管候選空間；a4～d4為26鄰域MRF模型的分割結(jié)果(＝0.40)；a5～d5為6鄰域MRF模型的分割結(jié)果(＝0.20)；a6～d6為MPN鄰域MRF模型的分割結(jié)果(＝0.19)。

針對每套數(shù)據(jù)的分割，計算過程實際在圖6(a3～d3)得到的血管候選空間進行，三種鄰域模型得到三個分割結(jié)果分別顯示在圖6(a4～d4)、(a5～d5)、(a6～d6)中。從MRA數(shù)據(jù)分割的一般效果來看存在部分非血管組織被當(dāng)做血管分割開來，原因在于：(1)腦部皮下脂肪組織和靜脈組織(矢狀竇、乙狀竇、枕下竇等)在灰度分布區(qū)間上位于高灰度腦部血管的中低端；(2)這些結(jié)構(gòu)在空間上相距較近(最近時往往間隔幾個體素)，容易被血管目標(biāo)函數(shù)擬合。為此，圖6的分割結(jié)果只顯示了最大連通區(qū)域——血管，以盡可能完全地剔除其它非血管組織。然而，仍有部分粘連的靜脈組織被當(dāng)做血管分割出來，如圖6(a4～a6)中存在的矢狀竇，以及圖6(d4～d6)中存在的乙狀竇。不妨從總體上來看，圖6(a6~d6)的分割結(jié)果更加接近其原始數(shù)據(jù)和多尺度濾波在MIP中的血管結(jié)構(gòu)形態(tài)。原因在于空間單位立方體中，MPN系統(tǒng)能夠較大程度避免26鄰域系統(tǒng)存在的多對一問題(其實質(zhì)表現(xiàn)為一種MRF模型高級能量V2的多對一問題)，而且避免了6鄰域系統(tǒng)存在的血管連續(xù)性缺失問題(其實質(zhì)表現(xiàn)為血管收到噪聲破壞而斷開的問題)。因此，基于MPN的MRF模型能夠分割出更加完整的血管網(wǎng)絡(luò)。

計算耗時方面：實驗中采用鄰域大小為26、6、MNP的MRF分割模型對各套MRA數(shù)據(jù)的血管候選空間進行計算，消耗的時間平均為2.4分鐘、2.12分鐘、3.82分鐘；我們還測試了算法作用于整個MRA數(shù)據(jù)空間，分割結(jié)果完全相同，計算量分別增加到2小時30分、2小時5分、3小時20分。由此可見，本方法在血管分割的精度和速度方面可以為臨床應(yīng)用提供重要的理論依據(jù)。

需要特別注意的是：本算法依據(jù)的MRF分割模型建立在數(shù)據(jù)空間的各向同性和均勻性基礎(chǔ)上。為此，分割前需要對原始MRA數(shù)據(jù)進行空間均勻性重采樣和消除磁場非均勻性。

5 結(jié) 論

本文提供了一種全自動的基于三維MRF模型的腦部MRA數(shù)據(jù)分割方法。該方法與傳統(tǒng)的分割算法相比，主要優(yōu)勢在于：提出了一種新穎的MRF高級模型的鄰域構(gòu)造方法，其分割效果優(yōu)于傳統(tǒng)鄰域模型；其次，提出了一種基于血管候選空間的快速計算策略，在結(jié)果一致的前提下極大地滿足了臨床實時性的要求。此外，在極大似然估計、偽似然高級參數(shù)估計、EM低級模型參數(shù)估計方面，利用了大型矩陣邏輯運算這一并行處理的方式，也較大程度減少了計算量。總之，本方法不僅對臨床腦血管快速、自動、精確分割提供了有效方法，也為基于三維MRF的醫(yī)學(xué)圖像分割提供了很好的支持。隨著當(dāng)前計算機硬件的快速發(fā)展，以海量計算為特點的統(tǒng)計學(xué)建模方法將有著更大的應(yīng)用前景。

[1] Bogunovic H, Pozo JM, Villa-Uriol MC, et al. Automated segmentation of cerebral vasculaturewith aneurysms in 3DRA and TOF-MRA using geodesic active regions: an evaluation study [J]. Medical Physics, 2011, 38(1), 210-222.

[2] Frangi AF, Niessen WJ, Vincken KL, et al. Multiscale vessel enhancement filtering [C] // Proceedings of the International Conference on Medical Image Computing and Computer-assisted Interventation - MICCAI’98 Lecture Notes in Computer Science, 1998, 1496: 130-137.

[3] Law MW, Chung AC. Eff i cient implementation for spherical flux computation and its application to vascular segmentation [J]. IEEE Transactions on Image Processing, 2009, 18(3): 596-612.

[4] Yan PK, Kassim AA. Segmentation of volumetric MRA images by using capillary active contour [J]. Medical Image Analysis, 2006, 10(3): 317-329.

[5] Chen J, Amini AA. Quantifying 3-D vascular structures in MRA images using hybrid PDE and geometric deformable models [J]. IEEE Transactions on Medical Imaging, 2004, 23(10): 1251-1262.

[6] Wilson DL, Noble JA. An adaptive segmentation algorithm for time-of-flight MRA data [J]. IEEE Transactions on Medical Imaging, 1999, 18(10): 938-945.

[7] Gao X, Uchiyama Y, Zhou X, et al. A fast and fully automatic method for cerebrovascular segmentation on time-of-f l ight (TOF) MRA image [J]. Journal of Digital Imaging, 2011, 24(4): 609-625.

[8] Hassouna MS, Farag AA, Hushek S, et al. Cerebrovascular segmentation from TOF using stochastic models [J]. Medical Image Analysis, 2006, 10: 2-18.

[9] Hao JT, Li ML, Tang FL, et al. Adaptive segmentation of cerebrovascular tree in time-of-flight magnetic resonance angiography [J]. Medical & Biological Engineering & Computing , 2008, 46: 75-83.

[10] Ayman EB, Aly AF, Gimelfarb G, et al. Automatic cerebrovascular segmentation by accurate probabilistic modeling of MRA images [C] // Medical Image Computing and Computer-assisted Intervention - Lecture Notes in Computer Science, 2005, 8: 34-42.

[11] Jasra A, Holmes CC, Stephens DA, et al. Markov chain Monte Carlo methods and the label switching problem in Bayesian mixture modeling [J]. Statistical Science, 2005, 20(1): 50-67.

[12] Lin TI, Lee JC, Yin SY. Finite mixture modelling using the skew normal distribution [J]. Statistica Sinica, 2007, 17: 909-927.

Segmentation of Three-Dimensional Data of Brain Magnetic Resonance Angiography Based on Markov Random Field

ZHOU Shoujun JIA Fucang HU Qingmao XIE Yaoqing GU Jia SHANG Peng
( Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, Shenzhen 518055, China )

In this paper, a full automatic method was proposed for the segmentation of brain magnetic resonance angiography (MRA) dataset, which improved the technologies of Markov random fi eld (MRF). Existing 3D-MRF models generally faced some challenges including: (1) The parameter initialization of low level MRF model is not accurate; (2) The ordinary neighborhood system cannot deal with local fi ne vessel structure. Aiming to solve the two problems, the multiscale fi ltering with threshold analysis and a multi-pattern neighborhood system were proposed, respectively. Such method enabled the MRF model delineating vessels to be as small as two voxels in diameters. In the experiments, the parameters of the low level MRF model were estimated using the expectation maximization algorithm, while the parameters of the high level MRF models were estimated based on the maximum pseudo likelihood algorithm. A set of phantoms and some MRA clinical datasets were used to validate the algorithms, to yield smaller segmentation errors.

vessel segmentation; magnetic resonance angiography; Markov random fi eld; neighborhood system

R 445

A

2013-11-20

廣東省引進科研創(chuàng)新團隊項目(2011S013)；國家自然科學(xué)基金項目(61179020) 。

周壽軍(通訊作者)，博士，副研究員，研究方向為醫(yī)學(xué)圖像處理與計算機輔助診斷，E-mail：sj.zhou@siat.ac.cn；賈富倉，博士，高級工程師，研究方向為計算機輔助醫(yī)學(xué)、醫(yī)學(xué)數(shù)據(jù)分析處理等；胡慶茂，博士，研究員，研究方向為生物醫(yī)學(xué)圖像分析、重大疾病的計算機輔助診療；謝耀欽，博士，研究員，研究方向為影像引導(dǎo)治療、多模影像配準(zhǔn)等；辜 嘉，博士，研究員，研究方向為各種仿生型醫(yī)療機器人的研制；尚鵬，研究方向為人工關(guān)節(jié)、人體肌骨系統(tǒng)運動學(xué)與動力學(xué)、計算機成像輔助外科手術(shù)系統(tǒng)及軟組織生物力學(xué)。