胡 昊，張明慧，盧振泰，馮前進(jìn)，陳武凡

（南方醫(yī)科大學(xué) 生物醫(yī)學(xué)工程學(xué)院，廣東 廣州510515）

0 引 言

從臨床研究可知，海馬形態(tài)體積和功能的改變與許多神經(jīng)系統(tǒng)疾病有關(guān)，如顳葉癲癇和阿爾茲海默?。?］。目前，對(duì)于阿爾茲海默病，只能通過早期診斷進(jìn)而采取有效措施延緩疾病的進(jìn)程。對(duì)海馬體積的測(cè)量為辨別這種疾病的重要方法［2］。海馬體積測(cè)量建立在海馬分割的基礎(chǔ)上，所以準(zhǔn)確地從腦部圖像中將海馬分割出來顯得尤為重要。

人工分割海馬效率低，主觀性強(qiáng)，半自動(dòng)的海馬分割方法可以減少人工分割的時(shí)間與負(fù)擔(dān)，但在大樣本的研究中，仍不能成為理想的分割方式。自動(dòng)的圖像分割方法［3，4］方為圖像分割研究的發(fā)展趨勢(shì)。海馬形狀不規(guī)則，邊緣模糊，難以與周圍其余組織區(qū)分開來，其分割是圖像分割研究領(lǐng)域的難點(diǎn)。借助醫(yī)學(xué)圖譜提供的先驗(yàn)信息有助于更好地完成分割任務(wù)。配準(zhǔn)技術(shù)有效地將醫(yī)學(xué)先驗(yàn)知識(shí)融入分割過程，形成一種新的自動(dòng)分割技術(shù)，稱為基于配準(zhǔn)的分割［5］。本文提出的基于多圖譜配準(zhǔn)的自動(dòng)分割算法，在配準(zhǔn)時(shí)不但考慮圖像對(duì)應(yīng)點(diǎn)的灰度信息，亦融入空間位置信息，實(shí)現(xiàn)灰度與空間的同步標(biāo)準(zhǔn)化，提高配準(zhǔn)算法的有效性。配準(zhǔn)結(jié)果的精確度得到保證，融合后的分割結(jié)果更接近真實(shí)值。

1 基于多圖譜配準(zhǔn)的海馬分割方法

基于多圖譜配準(zhǔn)的分割方法可以分為基于單圖譜、基于平均形狀圖譜和基于多圖譜3種。其中，基于多圖譜配準(zhǔn)的分割方法準(zhǔn)確度更高，適用于對(duì)海馬進(jìn)行分割［6］。本文介紹的是一種基于多圖譜配準(zhǔn)的自動(dòng)分割算法。該類算法將預(yù)標(biāo)記的圖像向目標(biāo)圖像配準(zhǔn)，建立對(duì)應(yīng)關(guān)系并得到相應(yīng)空間變換后，利用這個(gè)空間變換，將與預(yù)標(biāo)記圖像對(duì)應(yīng)的分割結(jié)果即圖譜，映射到目標(biāo)圖像，再選擇合適的融合算法，生成目標(biāo)圖像的分割結(jié)果。算法流程如圖1所示。

圖1 算法流程

圖像配準(zhǔn)與圖譜融合為基于多圖譜配準(zhǔn)分割算法的兩個(gè)關(guān)鍵點(diǎn)。本文選擇一種新的基于最小化殘差復(fù)雜度（modified minimization of residual complexity）的配準(zhǔn)方法［7］，并使用STAPLE（simultaneous truth and performance level estimation）算法對(duì)圖譜進(jìn)行融合［8］，得到最終分割結(jié)果。

1.1 最小化殘差復(fù)雜度的圖像配準(zhǔn)算法

對(duì)磁共振圖像而言，由于磁場(chǎng)的不均勻性，造成圖像之間存在亮度不一致等問題，常用的配準(zhǔn)方法無法達(dá)到較好的配準(zhǔn)效果。最小化殘差復(fù)雜度算法（minimization of residual complexity）適合處理空間灰度不均勻的問題。假設(shè)I與J 是待配準(zhǔn)的兩幅圖像，則有

式中：I——參考圖像，J——浮動(dòng)圖像，T——對(duì)浮動(dòng)圖像J 進(jìn)行的形變，S——灰度校正場(chǎng)，η——零均值的高斯噪聲。經(jīng)過推導(dǎo)和演化，得能量函數(shù)如下

式中：N——體素個(gè)數(shù)，α——可調(diào)節(jié)的參數(shù)，qnT——對(duì)灰度校正場(chǎng)的約束演化而來的矩陣的特征向量。詳細(xì)推導(dǎo)過程請(qǐng)參考文獻(xiàn)［7］。

該方法對(duì)亮度不一致、明暗不均勻的圖像配準(zhǔn)有很好的效果。然而，由于磁共振掃描時(shí)間較長(zhǎng)，頭部運(yùn)動(dòng)增加，掃描結(jié)果易出現(xiàn)噪聲和失真。此外，磁共振磁場(chǎng)本身存在不均勻性，不同儀器、不同時(shí)間、不同病人的掃描，相互之間存在較多偶然因素。噪聲及這些偶然因素造成的突變值被融入殘差r中，從而影響式（1）的能量函數(shù)，導(dǎo)致配準(zhǔn)精度下降，魯棒性較差。對(duì)于基于多圖譜配準(zhǔn)的分割方法，只有配準(zhǔn)精度更高，最終的分割結(jié)果才能更精確。因此，本文對(duì)該方法進(jìn)行改進(jìn)，構(gòu)造一個(gè)新的配準(zhǔn)算法。

1.2 改進(jìn)后的最小化殘差復(fù)雜度的圖像配準(zhǔn)算法

式（1）中的r為參考圖像與浮動(dòng)圖像的灰度值差，這種完全依賴圖像灰度信息的配準(zhǔn)方法沒有最大化地利用圖像提供的所有信息。改進(jìn)后的新算法在考慮圖像對(duì)應(yīng)點(diǎn)的灰度信息的同時(shí)，融入空間位置信息。對(duì)不同像素點(diǎn)，根據(jù)其方向與空間信息的不同，賦予不同權(quán)重，進(jìn)行空間標(biāo)準(zhǔn)化。權(quán)重參數(shù)計(jì)算步驟如圖2所示。

圖2 權(quán)重參數(shù)計(jì)算流程

對(duì)二維圖像，鄰域選擇為8 鄰域系統(tǒng)，如圖3 所示。三維圖像則為26鄰域。以三維圖像為例描述權(quán)重參數(shù)求解方法如下。

圖3 8鄰域系統(tǒng)

梯度圖可以在一定程度上反映原圖像灰度、位置的變化，如圖4（b）所示，在各器官組織邊緣部分，梯度值較大，而灰度均勻的同一組織內(nèi)，梯度值則很小。由于磁共振掃描的特性，在兩幅圖像中，可能出現(xiàn)同樣的組織有不同的灰度值，而相同的灰度值代表不同組織的情況。對(duì)梯度圖而言，即使相同組織灰度值不同，梯度信息同樣可以準(zhǔn)確地反映該組織的位置結(jié)構(gòu)信息。因此，權(quán)重參數(shù)ω的提出，為r融入了圖像的空間位置信息，完成配準(zhǔn)圖像灰度與空間的同步標(biāo)準(zhǔn)化。這種方法可以有效降低噪聲或其它偶然因素造成的圖像灰度值突變的影響，具有更高的魯棒性與配準(zhǔn)精度。

圖4 腦部MR 圖像及對(duì)應(yīng)梯度

兩種配準(zhǔn)方法均采用自由形變模型（free form transformation）對(duì)圖像進(jìn)行形變。自由形變模型的方法因?yàn)槠浔普娴纳斫Y(jié)構(gòu)模擬特性，被廣泛地用于圖像配準(zhǔn)中［9］。此外，為避免陷入局部極值，均采用為由粗到細(xì)的多分辨率網(wǎng)格技術(shù)。將多分辨率網(wǎng)格與自由形變模型相結(jié)合，可以提高計(jì)算效率，保證配準(zhǔn)結(jié)果的準(zhǔn)確度。對(duì)形變場(chǎng)迭代求解時(shí)則選擇梯度下降法為優(yōu)化算法［10］。該方法工作量小，需要存儲(chǔ)的變量少，對(duì)初始點(diǎn)的要求較低。

1.3 STAPLE融合算法

本文采用STAPLE算法實(shí)現(xiàn)圖譜融合。STAPLE 算法的目標(biāo)為利用函數(shù)f（D，Tr｜p，q）的似然函數(shù)估計(jì)性能水平參數(shù)（p，q）。f（D，Tr｜p，q）的最大化就實(shí)現(xiàn)了參數(shù)（p，q）的最優(yōu)值。公式表示為

式中：D——N×R 的矩陣，表示圖譜的每個(gè)像素在每個(gè)分割中的結(jié)果，R——圖譜數(shù)量，N——每個(gè)圖譜的體素個(gè)數(shù)，Tr——目標(biāo)圖像的真實(shí)分割結(jié)果。對(duì)每個(gè)像素用二進(jìn)制0或1表示，其中1指示分割對(duì)象，0指示非分割對(duì)象。f（D，Tr｜p，q）表示概率群分布函數(shù)。其中pj表示敏感度，qj表示特異性

式中：pj和qj的取值范圍為［0，1］。假設(shè)真實(shí)的分割結(jié)果、pj和qj、圖譜的分割結(jié)果是條件獨(dú)立的，即不同分割結(jié)果之間相互獨(dú)立，沒有影響，那么分割質(zhì)量就可以用pj和qj來評(píng)價(jià)。

對(duì)式（3）可采用最大期望算法估計(jì)其值，詳細(xì)步驟可參照文獻(xiàn)［11，12］。

2 配準(zhǔn)及分割結(jié)果的評(píng)價(jià)方法

選擇合適的評(píng)價(jià)方法，可以幫助定量、科學(xué)、有效地對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行分析，并對(duì)不同實(shí)驗(yàn)方法進(jìn)行深入比較，通過正確整合和嚴(yán)格評(píng)價(jià)數(shù)據(jù)來提取有用數(shù)據(jù)支持結(jié)論。

2.1 配準(zhǔn)結(jié)果的評(píng)價(jià)

本文選擇互信息作為配準(zhǔn)結(jié)果的衡量標(biāo)準(zhǔn)?；バ畔⑹且环N度量?jī)蓚€(gè)系統(tǒng)的統(tǒng)計(jì)相關(guān)性的量，其計(jì)算復(fù)雜度低，魯棒性好。計(jì)算公式如下式

式中：PS、PA’——參考圖像S和形變后的浮動(dòng)圖像A’的邊緣概率密度分布函數(shù)，PA’S——聯(lián)合概率密度分布函數(shù)。計(jì)算的互信息值越大，表示兩幅圖像的空間位置越一致，配準(zhǔn)效果越好［13］。

假設(shè)用NRC方法配準(zhǔn)后得到的圖像計(jì)算的互信息值為mrcMI，用RC方法配準(zhǔn)后得到的圖像計(jì)算的互信息值為rcMI，由下式計(jì)算比較兩種方法得到的互信息值

2.2 分割結(jié)果的評(píng)價(jià)

相似性測(cè)度Dice為分割結(jié)果的常用評(píng)價(jià)方法。相似性測(cè)度是對(duì)目標(biāo)圖像的真實(shí)分割結(jié)果和融合得到的分割結(jié)果的重疊率進(jìn)行評(píng)價(jià)，計(jì)算公式如下

式中：Sseg——目標(biāo)圖像的真實(shí)分割結(jié)果，Sfus——融合得到的分割結(jié)果，V——體積大小。相似性測(cè)度Dice值越大，說明目標(biāo)圖像的真實(shí)分割結(jié)果和融合得到的分割結(jié)果的重疊率越高，結(jié)果越精確［14］。

3 實(shí)驗(yàn)結(jié)果及分析

磁共振圖像具有非損傷性的優(yōu)點(diǎn)，對(duì)軟組織有極好的分辨力，被廣泛用于臨床醫(yī)學(xué)研究。本文選擇的18 個(gè)T1加權(quán)的磁共振腦部圖像均來自美國的Henry Ford醫(yī)院。18個(gè)圖像分別由兩臺(tái)不同場(chǎng)強(qiáng)的磁共振掃描儀獲得，部分病人患有顳葉癲癇癥，出現(xiàn)海馬體萎縮的現(xiàn)象，其余海馬體的形態(tài)功能正常。

為增加實(shí)驗(yàn)結(jié)果的可信度，避免出現(xiàn)偶然性，隨機(jī)抽取10個(gè)圖像分別作為待分割圖像，編號(hào)為01到10號(hào)，剩余8個(gè)圖像編號(hào)為11 至18 號(hào)。用Matlab7.12.0在PC 機(jī)上進(jìn)行了實(shí)驗(yàn)，并分別用最小化殘差復(fù)雜度和改進(jìn)后的最小化殘差復(fù)雜度兩種不同的配準(zhǔn)方法進(jìn)行對(duì)比。兩種不同配準(zhǔn)方法最終分別得到10個(gè)融合結(jié)果。

3.1 圖像預(yù)處理

18個(gè)圖像分別由兩臺(tái)不同場(chǎng)強(qiáng)的磁共振掃描儀獲得，圖像的分辨率存在較大差異，配準(zhǔn)前需進(jìn)行預(yù)處理。處理方法為以圖01為基準(zhǔn)，將其余17個(gè)圖像歸一化到基準(zhǔn)圖像的空間上，使它們有一致的原點(diǎn)與尺寸。此外，圖像的對(duì)比度也存在一定差別，需對(duì)圖像的灰度值進(jìn)行調(diào)節(jié)，使得18個(gè)圖像的灰度值和對(duì)比度處于同一水平。

本算法是對(duì)目標(biāo)圖像的海馬進(jìn)行分割，配準(zhǔn)部位不必涵蓋整個(gè)腦部。因此，為提高配準(zhǔn)效率，降低算法所需時(shí)間，以海馬為中心進(jìn)行ROI提取。處理方法為減小圖像的體積，提取一個(gè)以海馬為中心的長(zhǎng)方體。提取后的圖像大小均為64×128×64。

圖5為01 號(hào)圖像分別在預(yù)處理前及預(yù)處理后，橫斷面、冠狀面和矢狀面的二維圖像，其中第一行為預(yù)處理前，第二行為預(yù)處理后。箭頭指示區(qū)域?yàn)閳D像對(duì)應(yīng)的海馬。

圖5 預(yù)處理前及預(yù)處理后的01圖像

經(jīng)過預(yù)處理后圖像以海馬為中心，去除與分割目標(biāo)無關(guān)的大部分組織，能更好地對(duì)海馬區(qū)域進(jìn)行配準(zhǔn)，大大降低算法所需時(shí)間。

3.2 基于最小化殘差復(fù)雜度算法的配準(zhǔn)

為進(jìn)行對(duì)照實(shí)驗(yàn)，采用了最小化殘差復(fù)雜度和改進(jìn)后的最小化殘差復(fù)雜度兩種方法進(jìn)行配準(zhǔn)。其中，兩種配準(zhǔn)方法的參數(shù)選擇如下：網(wǎng)格大小為8，參數(shù)α為：0.01，參數(shù)λ為：0.001，其中λ為式（1）中約束項(xiàng)qnT的參數(shù)。

圖6為式（4）計(jì)算結(jié)果的條形圖，橫坐標(biāo)表示不同目標(biāo)圖像的十組配準(zhǔn)結(jié)果，每組分別計(jì)算得到17個(gè)不同值，縱坐標(biāo)表示MIdif 值。

圖6 MIdif 值條形圖

由式（5）可知，MIdif 值若為正數(shù)，表示由式（4）計(jì)算的mrcMI值大于rcMI 值。從圖6可觀察到，MIdif 值多為正數(shù)，說明MRC方法配準(zhǔn)得到的圖像計(jì)算的互信息值比RC方法配準(zhǔn)得到的圖像計(jì)算的互信息值要大，配準(zhǔn)效果更好。

配準(zhǔn)完成后，每組分別得到17個(gè)形變場(chǎng)，利用該形變場(chǎng)將已分割的海馬映射到目標(biāo)圖像上，即可得到17個(gè)形變后的海馬。圖8為07號(hào)圖像作為待分割目標(biāo)圖像時(shí)，映射后得到的目標(biāo)圖像的17個(gè)分割結(jié)果。其中第一幅為07號(hào)圖像對(duì)應(yīng)的真實(shí)海馬，后17幅為映射形變后的海馬。圖7只顯示橫斷面。

圖7 17號(hào)圖像的真實(shí)分割圖及映射后的各海馬

觀察圖7可發(fā)現(xiàn)，17個(gè)初步分割結(jié)果在形狀位置上都與真實(shí)值有一定差異，若選擇其中一個(gè)圖像進(jìn)行單圖譜配準(zhǔn)分割，得到的結(jié)果必不能作為有效的海馬分割圖。因此，本文采用多圖譜配準(zhǔn)分割方法，對(duì)17個(gè)初步分割結(jié)果進(jìn)行融合，生成一個(gè)與真實(shí)值更為接近、更精確的海馬圖。

3.3 基于STAPLE算法的圖譜融合

本文采用的融合方法為STAPLE 算法。參數(shù)選擇如下：收斂因子為0.00001，初始化類型為0，插值類型選擇鄰域插值法。其中初始化類型的選擇與最大期望算法的初始估計(jì)值有關(guān)。對(duì)兩種不同配準(zhǔn)方法分別得到的十組結(jié)果融合后，最終分割的海馬如圖8所示。其中第一行與第四行為目標(biāo)圖像的真實(shí)分割結(jié)果，第二行與第五行為用RC方法進(jìn)行配準(zhǔn)融合后得到的結(jié)果，第三行與第六行為用MRC方法進(jìn)行配準(zhǔn)融合后得到的結(jié)果。第一列第一行至第三行為第一組分割結(jié)果對(duì)比，第二列第一行至第三行為第二組分割結(jié)果對(duì)比，以此類推。

由圖8可觀察到，MRC配準(zhǔn)方法融合得到的海馬，在形狀、大小上都與真實(shí)分割結(jié)果更為接近。說明基于該配準(zhǔn)方法將圖譜映射到目標(biāo)圖像上得到的結(jié)果與RC 配準(zhǔn)方法相比更精確，更接近目標(biāo)圖像。

通過圖8可對(duì)分割結(jié)果進(jìn)行主觀定性分析，但定量分析的方法更為科學(xué)，在對(duì)實(shí)驗(yàn)結(jié)果的評(píng)價(jià)中更為重要。通過式（6）計(jì)算兩種不同配準(zhǔn)方法最后融合的結(jié)果與待分割的目標(biāo)圖像的真實(shí)分割數(shù)據(jù)的Dice值，計(jì)算結(jié)果見表1。

圖8 目標(biāo)圖像的真實(shí)分割海馬與兩種不同配準(zhǔn)方法融合結(jié)果

表1 兩種不同配準(zhǔn)方法融合結(jié)果Dice值

對(duì)比表1數(shù)據(jù)可以發(fā)現(xiàn)，用RC方法配準(zhǔn)后融合的結(jié)果的Dice值有兩組大于0.8，而用MRC方法配準(zhǔn)后融合的結(jié)果的Dice值均在0.8以上。改進(jìn)后的配準(zhǔn)方法明顯地提高了海馬分割的精度。

為更直觀地對(duì)表1數(shù)據(jù)進(jìn)行對(duì)比分析，分別計(jì)算兩種方法Dice值的均值并作折線，如圖9所示。

圖9 兩種不同配準(zhǔn)方法融合結(jié)果Dice值折線

觀察圖9可以看出，MRC配準(zhǔn)方法融合得到的十組分割結(jié)果Dice均值遠(yuǎn)大于RC 配準(zhǔn)方法融合得到的結(jié)果，實(shí)心圓點(diǎn)折線位于空心圓點(diǎn)折線之上，且更平緩，即分割結(jié)果更精確、更穩(wěn)定。由此說明MRC 配準(zhǔn)方法的精確度更高，本文選用的基于MRC的配準(zhǔn)方法的多圖譜配準(zhǔn)分割算法是一種有效的海馬分割算法。

4 結(jié)束語

本文提出一種基于多圖譜配準(zhǔn)的海馬分割算法，該算法在配準(zhǔn)時(shí)融入圖像的空間與方向信息，增加配準(zhǔn)算法的魯棒性與有效性，從而實(shí)現(xiàn)更精確的海馬分割。從實(shí)驗(yàn)結(jié)果可以觀察到，該分割算法更穩(wěn)定，對(duì)海馬分割的精度更高。由此表明本文提出的基于NRC配準(zhǔn)方法的多圖譜配準(zhǔn)分割算法是一種有效的海馬分割算法。分析實(shí)驗(yàn)結(jié)果Dice值可知，分割結(jié)果仍有一定的提高空間。后續(xù)工作中，可增加圖譜數(shù)量并進(jìn)行擇優(yōu)篩選，或采用多種不同的融合算法，進(jìn)行對(duì)比研究，找到最適合進(jìn)行海馬分割的基于多圖譜配準(zhǔn)的分割方法。

［1］Yushkevich Paul A，Wang Hongzhi，Pluta Jogn，et al.Nearly automatic segmentation of hippocampal subfields in vivo focal T2－weighted MRI ［J］.Neuroimage，2010，53 （4）：1208－1224.

［2］Morra Jonathan H，Tu Zhuowen，Apostolova Liana G，et al.Validation of a fully automated 3D hippocampal segmentation method using subjects with Alzheimer’s disease，mild cognitive impairment，and elderly controls ［J］.Neuroimage，2008，43（1）：59－68.

［3］Chupin Marie，Gérardin Emilie，Cuingnet Rémi，et al.Fully automatic hippocampus segmentation and classification in Alzheimer’s disease and mild cognitive impairment applied on data from ADNI［J］.Hippocampus，2009，19 （6）：579－587.

［4］Leemput Koen Van，Bakkour Akram，Benner Thomas，et al.Automated segmentation of hippocampal subfields from Ultrahigh resolution in vivo MRI ［J］.Hippocampus，2009，19（6）：549－557.

［5］Chen Wenyan，Li Shutao，Jia Fucang，et al.Segmentation of hippocampus based on ROI atlas registration ［R］.International Symposium on IT in Medicine and Education，2011：226－230.

［6］Shan Liang，Charles Cecil，Niethammer Marc.Automatic multi－atlas－based cartilage segmentation from knee MR images［R］.9th IEEE International Symposium on Biomedical Imaging，2012：1028－1031.

［7］Myronenko Andriy，Song Xubo.Image registration by minimization of residual complexity ［R］.IEEE Conference on Computer Vision and Pattern Recognition，2009：49－56.

［8］Andrew J A，Bennett A L.Non－local STAPLE：an intensitydriven multi－atlas rater model［J］.Med Image Comput Comput Assist Interv，2012，15 （Pt 3）：426－434.

［9］Rueckert Daniel，Aljabar Paul，Heckemann Rolf A，et al.Diffeomorphic registration using B－splines ［J］.Med Image Comput Comput Assist Interv，2006，9 （2）：702－709.

［10］LI Bin，OU Shanxing，TIAN Lianfang，et al.Thorax multimodal medical image registration based on adaptive free－form deformation and gradient descent［J］.Application Research of Computers，2009，26 （10）：3978－3982 （in Chinese）. ［李彬，歐陜興，田聯(lián)房，等.基于自適應(yīng)自由變形法和梯度下降法的胸部多模醫(yī)學(xué)圖像配準(zhǔn) ［J］.計(jì)算機(jī)應(yīng)用研究，2009，26 （10）：3978－3982.］

［11］Wikipedia.Expectation maximization，http：／／en.wikipedia.org／wiki／Expectation＿maximization［DB／OL］.2013－7－20（in Chinese）.［維基百科.最大期望算法，http：／／zh.wikipedia.org／wiki／最大期望算法［DB／OL］.2013－7－20.］

［12］CHEN Wenyan，LI Shutao.Segmentation of hippocampus based on ROI atlas registration ［D］.Changsha：Hunan University，2012 （in Chinese）.［陳雯艷，李樹濤.基于ROI多圖譜配準(zhǔn)的海馬磁共振圖像分割 ［D］.長(zhǎng)沙：湖南大學(xué)，2012.］

［13］SHI Yingqi，GU Lixu.Optimizied multistage medical image registration method based on mutual information ［J］.Computer Engineering，2006，32 （22）：187－190 （in Chinese）.［施穎琦，顧力栩.基于互信息多步驟優(yōu)化的醫(yī)學(xué)圖像配準(zhǔn)［J］.計(jì)算機(jī)工程，2006，32 （22）：187－190.］

［14］Han Xiao，F(xiàn)ischl Bruce.Atlas renormalization for improved brain MR image segmentation across scanner platforms ［J］.IEEE Transactions on Medical imaging，2007，26 （4）：479－486.