徐云龍　鄭元杰　鄧　翔　李　寧　湯煜春　尹義龍

(1.山東大學(xué)計算機(jī)科學(xué)與技術(shù)學(xué)院，濟(jì)南，250101； 2.山東師范大學(xué)信息科學(xué)與工程學(xué)院，濟(jì)南，250014； 3.山東省立醫(yī)院，濟(jì)南，250021； 4.山東大學(xué)醫(yī)學(xué)院，濟(jì)南，250012)

引　　言

主動脈是人體內(nèi)最重要的動脈主干，源源不斷地將血液由心臟輸送至人體各重要組織器官。隨著不健康生活方式越來越常見，罹患主動脈相關(guān)疾病的患者數(shù)量逐年增加，患者及其家庭遭受極大的痛苦。根據(jù)相關(guān)數(shù)據(jù)[1]，在發(fā)達(dá)國家，主動脈瓣膜相關(guān)疾病已成為最常見的瓣膜性疾病。主動脈相關(guān)疾病發(fā)展迅速且前期很少表現(xiàn)出明顯癥狀，這給患者的早期診斷和治療帶來極大的困難。因此，對主動脈疾病的早期精確診斷和治療具有極其重要的現(xiàn)實意義。

過去幾十年間，醫(yī)學(xué)成像技術(shù)取得了長足的發(fā)展，成像設(shè)備的推陳出新、成像方法的快速進(jìn)步、成像類別的不斷細(xì)分，使得醫(yī)學(xué)圖像在臨床醫(yī)學(xué)中發(fā)揮著越來越重要的作用。傳統(tǒng)的醫(yī)學(xué)圖像手工分割方法繁瑣且效率低下，因此不適用于所要研究的主動脈圖像分割問題。特別是對于3D的主動脈圖像，手工分割方法的劣勢進(jìn)一步凸顯。同時，由于手工分割方法的主觀性和不可重復(fù)性，加之手工分割方法嚴(yán)重依賴于標(biāo)記專家的經(jīng)驗，導(dǎo)致手工分割方法無法滿足現(xiàn)實需求。因此，計算機(jī)輔助的、高效的、自動的主動脈醫(yī)學(xué)圖像分割技術(shù)應(yīng)運而生。

按照分類原理的不同，主動脈分割方法可以大致分為如下兩類。第一類是基于最優(yōu)統(tǒng)計理論的分割方法，此類方法包括圖割[2]、馬爾科夫隨機(jī)場[3]等方法，通過優(yōu)化分割代價函數(shù)，得到圖像的最優(yōu)分割。此類分割方法理論簡單，但對于某些區(qū)分性不強(qiáng)的邊緣區(qū)域，需要進(jìn)行較多的預(yù)處理過程，提取較多的復(fù)雜特征，提高分割的維度以達(dá)到精度要求。第二類方法是基于可變模型理論的分割方法，此類方法包括水平集[4-6]、Snakes[7]等，需要設(shè)定大量參數(shù)，且需要根據(jù)具體的圖像對輪廓進(jìn)行適當(dāng)?shù)某跏蓟?，同時伴隨著漫長的計算時間。

圖1　基于圖譜的醫(yī)學(xué)圖像分割方法流程圖Fig.1　　　　Flowchart of atlas based medical image 　　　segmentation

為了解決上述兩種方法存在的問題，本文使用了一種基于圖譜的醫(yī)學(xué)圖像分割方法以實現(xiàn)對3D主動脈CT圖像的自動分割?；趫D譜的分割方法是一類重要的計算機(jī)輔助分割技術(shù)。與其他醫(yī)學(xué)圖像分割方法相比，基于圖譜的方法對多種模態(tài)的醫(yī)學(xué)圖像都有比較魯棒的分割能力[8-9]，對3D的醫(yī)學(xué)圖像也可以實現(xiàn)自動分割?；趫D譜的分割方法總是伴隨著圖像配準(zhǔn)的相關(guān)操作，將多個人工手工標(biāo)記的圖像(即圖譜)映射到待分割的目標(biāo)圖像的圖像空間中，每個圖譜會得到目標(biāo)圖像的一個相應(yīng)候選分割。對于候選分割，可以選擇不同的標(biāo)簽融合策略，將每個候選分割的標(biāo)簽進(jìn)行融合，以此得到目標(biāo)圖像最終的分割結(jié)果。該過程如圖1所示。

最經(jīng)典的標(biāo)簽融合策略是多數(shù)投票法[10]此類方法由于完全忽略了質(zhì)量較差的圖譜所攜帶的信息，因此融合效果不佳；目前使用最多的融合方法是權(quán)重投票策略[11-12]根據(jù)各個圖譜與目標(biāo)圖像的相似度，分配給各個候選分割不同的權(quán)重進(jìn)行融合。此類方法只考慮了目標(biāo)圖像與圖譜之間的關(guān)系而忽略了圖譜之間的相關(guān)性(即不同圖譜產(chǎn)生的錯誤可能具有相關(guān)性)。為了全面考慮圖譜間以及圖譜與目標(biāo)圖像間的關(guān)系，文獻(xiàn)[13]提出了一種新穎的聯(lián)合標(biāo)簽融合策略，并且成功應(yīng)用于腦圖像分割中，取得了相較于以上兩種融合策略更好的分割結(jié)果。

本文將基于多圖譜的醫(yī)學(xué)圖像分割技術(shù)與聯(lián)合標(biāo)簽融合策略相結(jié)合，并首次應(yīng)用于3D主動脈CT圖像的自動分割問題中。為了應(yīng)對分割過程中因圖譜數(shù)量不足導(dǎo)致的精度下降問題，本文提出了一種圖譜更新算法，在提高分割精度的同時，保持了較低的計算復(fù)雜度。

1　基于多圖譜與聯(lián)合標(biāo)簽融合的主動脈CT圖像分割技術(shù)

1.1圖譜選擇

基于圖譜的醫(yī)學(xué)圖像分割方法通常需要消耗大量的計算時間。在進(jìn)行基于多圖譜的分類時，每個圖譜均需要經(jīng)過配準(zhǔn)操作，從原有圖像空間映射到目標(biāo)圖像所在的圖像空間中。這一過程須要大量的計算時間以及內(nèi)存開銷。除此之外，圖譜的質(zhì)量有優(yōu)劣之別，在對候選分割進(jìn)行最終的融合時，并不是每個候選分割都有相同的貢獻(xiàn)度，與目標(biāo)圖像相似度很低的圖譜通常會攜帶誤導(dǎo)性的信息，導(dǎo)致標(biāo)簽融合精度的下降?；谝陨蟽煞矫嬖蚩芍?，將所有圖譜都用于分割過程是不切實際的。文獻(xiàn)[14]亦指出，基于多圖譜的分割方法的精度高度依賴于所選用的圖譜質(zhì)量，與所有圖譜均參與分割相比，挑選出的少量的、高質(zhì)量的圖譜能帶來更好的分割結(jié)果。本文使用了互信息作為圖譜選擇的度量指標(biāo)，兩幅圖像之間的互信息越大，說明圖像之間的相似性越高。由于主動脈部位在整個CT圖像中所占的體積較小，因此計算了整個CT圖像之間的互信息。

圖2　圖譜選擇流程圖Fig.2　Flowchart of atlas selection

圖譜選擇具體過程如下：(1)經(jīng)過簡單的快速剛體變換和仿射配準(zhǔn)將圖譜與目標(biāo)圖像映射到同一圖像空間以確定二者之間的解剖學(xué)對應(yīng)關(guān)系；(2)計算各圖譜與目標(biāo)圖像之間的全局互信息；(3)按照互信息值的大小對圖譜進(jìn)行排序；(4)選出互信息值最大的K個圖譜。該過程的流程如圖2所示。

本文通過實驗確定了所選圖譜數(shù)量K的取值。所選圖譜數(shù)量與分割精度之間的關(guān)系如圖3所示，其中圖3(b)中的兩條線分別為Ground truth和分割結(jié)果，兩線重合度越高，則分割精度越高。圖3(b)顯示，隨著圖譜數(shù)量由5增大到10，分割精度逐漸提高并達(dá)到最高水平；隨著圖譜數(shù)量繼續(xù)增大，分割精度反而呈現(xiàn)出下降趨勢。在分割精度與計算復(fù)雜度之間達(dá)到平衡，本文選用K=10個圖譜。

圖3　不同圖譜數(shù)量的分割精度表現(xiàn)Fig.3　Performance of different numbers of atlases

1.2　圖譜更新

圖譜是人手工標(biāo)記的、描述各解剖學(xué)結(jié)構(gòu)位置、形狀及空間關(guān)系的原型描述信息[15]。圖譜可以通過整合若干已分割圖像的信息來生成，但大多數(shù)情況下，圖譜是通過專家的手工分割來獲得。手工分割生成圖譜是一項耗時費力的工作，同時需要極強(qiáng)的領(lǐng)域知識。以主動脈圖譜生成為例，一個包含100切片的主動脈圖譜，需要一個接受過專業(yè)訓(xùn)練的標(biāo)記人員進(jìn)行至少60 min的工作?；谶@個原因，很多情況下沒有足夠多的圖譜成為多圖譜分割的巨大挑戰(zhàn)。此外，有限數(shù)量的圖譜所具有的泛化能力不高。對于一些新的待分割圖像，其解剖學(xué)信息可能與現(xiàn)有圖譜差異很大，導(dǎo)致其與各圖譜的相似性都很低，進(jìn)而導(dǎo)致在確定最終標(biāo)簽的過程中現(xiàn)有圖譜只能提供很少有價值的指導(dǎo)信息。

為了解決上述問題，本文提出了一種全新的圖譜更新策略。該圖譜更新策略工作原理如下。每當(dāng)有一個新的待分割圖像IT，進(jìn)行剛體變換和仿射配準(zhǔn)并計算其與各圖譜Ai的相似度Sim(IT,Ai) 。如果對于i=1,2,…,n，其中n為圖譜的數(shù)量，Sim(IT,Ai) 均小于某一給定閾值ε(實驗中該值設(shè)定為0.6)，則可以認(rèn)定現(xiàn)有圖譜無法對目標(biāo)圖像的分割提供足夠的指導(dǎo)信息，需要對現(xiàn)有圖譜進(jìn)行更新操作。圖譜更新流程如下：(1)選出與目標(biāo)圖像相似度最小的圖譜；(2)利用該圖譜與目標(biāo)圖像，構(gòu)造一個新的模板作為一個新的圖譜；(3)新產(chǎn)生的圖譜替換掉相似度最小的圖譜。第2步中的模板構(gòu)造可以通過ANTs[16]軟件包實現(xiàn)。經(jīng)過上述操作，新產(chǎn)生的圖譜包含了目標(biāo)圖像和原圖譜的信息，能夠為分割過程提供更多的指導(dǎo)。同時，圖譜數(shù)量保持不變，分割過程的計算復(fù)雜度并未增加。

1.3　圖像配準(zhǔn)

本文中，兩個不同階段均用到了配準(zhǔn)操作。在圖譜選擇階段，快速剛體變換和仿射配準(zhǔn)被用來將各圖譜與目標(biāo)圖像映射到同一圖像空間以計算二者之間的相似度。在這一階段，由于所有圖譜均需要與目標(biāo)圖像進(jìn)行一次配準(zhǔn)操作，所以空間映射的速度成為這一步的關(guān)鍵。在該過程中，兩種配準(zhǔn)操作均采用4層變換，每層執(zhí)行100,50,25和10次迭代，平滑系數(shù)σ分別為3，2，1和0 Voxel。Convergence閾值設(shè)定為1×10-6，Convergence 窗口大小為10。

圖譜選擇階段完成后，K個高質(zhì)量的圖譜分別與目標(biāo)圖像進(jìn)行可變模型配準(zhǔn)操作。本文所用的可變模型配準(zhǔn)算法是ANTs的Symmetric image normalization method (Syn)[17]。Syn在微分同胚映射空間中最大化了互相關(guān)，以此來提高配準(zhǔn)的精度。在本文中，使用Syn的默認(rèn)參數(shù)組合，即執(zhí)行5層變換操作，每層變換的迭代次數(shù)為100，100，70，50和20，平滑系數(shù)σ分別為5，3，2，1和0 Voxel，Syn收縮系數(shù)為10,6,4,2和1。 Convergence閾值設(shè)定為1×10-6，Convergence 窗口大小設(shè)置為10。

1.4　聯(lián)合標(biāo)簽融合

在引言中提到，目前絕大多數(shù)融合策略獨立地對每個圖譜進(jìn)行權(quán)重分配，在此過程中只考慮了各個圖譜與待分割圖像之間的相關(guān)性。這會導(dǎo)致一個嚴(yán)重的問題：當(dāng)新添加的圖譜中包含較多冗余信息時，由于其對最終分割的確定幾乎無貢獻(xiàn)，采用現(xiàn)有融合策略會導(dǎo)致分割精度的下降。為了解決上述問題，文獻(xiàn)[13]提出了一種新的聯(lián)合標(biāo)簽融合策略(Joint label fusion，JLF)。聯(lián)合融合策略在對每個圖譜進(jìn)行權(quán)重分配時考慮了圖譜之間的相關(guān)性。有關(guān)人腦海馬體分割及海馬體子區(qū)分割的實驗表明，與當(dāng)前最新的幾種融合策略(STAPLE[18]，LWGaussian[11]，LWInverse[12])相比，聯(lián)合標(biāo)簽融合策略能夠取得更高的分割精度。聯(lián)合標(biāo)簽融合策略的具體原理如下。

為了方便描述，用IT表示待分割的目標(biāo)圖像，ST(x) 為目標(biāo)圖像IT的真實分割結(jié)果。A={Ai|i=1,2,…,n}代表n個圖譜，每個圖譜表示為Ai=(Ii,Si)的形式，其中，Ii和Si分別為第i個圖譜配準(zhǔn)后的圖像數(shù)據(jù)和對應(yīng)的分割結(jié)果。假設(shè)分割問題是二分類問題，即目標(biāo)圖像中的每個Voxel被標(biāo)記為0或1。

在該二分類問題中，分割產(chǎn)生的誤差表示為

δi(x)=Si(x)-ST(x)

(1)

式中：δi(x)為第i個圖譜的分割結(jié)果與目標(biāo)圖像真實分割在x位置的差異。經(jīng)過對每個圖譜進(jìn)行權(quán)重分配并融合后，x位置的最終分割結(jié)果可表示成

(2)

(3)

式中：wx為各圖譜權(quán)重組成的列向量；Mx為n*n矩陣，每一元素Mx(i,j) 表示第i個與第j個圖譜同時在x位置產(chǎn)生錯誤的幾率。Mx反映了圖譜兩兩之間的相關(guān)性，其具體形式為

Mx(i,j)=Eδi(x),δj(x)[δi(x),δj(x)|IT,I1,…,In]=p(δi(x)δj(x)=1|IT,I1,…,In)

(4)

矩陣Mx元素的取值可根據(jù)灰度值的相似度進(jìn)行估計

Mx(i,j)～[∑y∈N(x)|IT(y)-Ii(y)||IT(y)-Ij(y)|]β

(5)

式中：N(x)為x的鄰域，β為設(shè)定的模型參數(shù)。

(6)

式(6)可用拉普拉斯乘子法解決，即

(7)

(8)

式中α為0.01。該問題可以很容易地用二次規(guī)劃進(jìn)行解決。

2　實驗結(jié)果與分析

2.1實驗數(shù)據(jù)與預(yù)處理

本文選用了15例3D主動脈CT掃描圖像進(jìn)行實驗。15例CT掃描主動脈圖像均采集自Siemens SOMATOM definition flash scanner(100 KVp)。每例圖像數(shù)據(jù)包含的切片數(shù)目介于230與363之間，切片尺寸為512×512，切片厚度為1.5 mm。每例掃描數(shù)據(jù)的位深度為16 bit/pixel，體素空間大小為0.7×0.7×1.2。所有數(shù)據(jù)均經(jīng)過專家按照事先制定的標(biāo)記協(xié)議進(jìn)行標(biāo)記，整個標(biāo)記過程在ITK-SNAP軟件下進(jìn)行。

由于每例圖像數(shù)據(jù)的切片數(shù)目差異較大，且有大量切片不具有分割價值，因此用ImageJ對每例圖像數(shù)據(jù)進(jìn)行切片截取。每例圖像截取100片，范圍從心臟下方延伸到腹部。這一區(qū)域的主動脈呈現(xiàn)出長直管狀結(jié)構(gòu)，但由于該區(qū)域存在大量周邊組織和器官，且主動脈疾病會導(dǎo)致主動脈嚴(yán)重變形，因此具有較大的分割難度。3D CT圖像中的圖像方向及體素空間大小是配準(zhǔn)過程中兩個重要參數(shù)，為了解決原始圖像中兩參數(shù)不一致的問題，本文對圖像頭文件進(jìn)行了歸一化操作。

2.2　實驗實現(xiàn)及評估

本文采用了Leave-one-out的交叉驗證方式，一例圖像作為待分割的目標(biāo)圖像，其余14例圖像作為圖譜。文中應(yīng)用到的所有配準(zhǔn)操作及聯(lián)合標(biāo)簽融合策略均基于ANTs實現(xiàn)。所有實驗均在配備12核心至強(qiáng)E5處理器以及64 GB內(nèi)存的計算環(huán)境下運行。

圖4展示了本文方法對3D主動脈CT圖像的分割結(jié)果。與圖4(a) Ground truth相比，圖4(b)的分割結(jié)果總體上保持了主動脈的結(jié)構(gòu)，除主動脈兩端邊緣會出現(xiàn)輕微的模糊及扭曲外，總體表現(xiàn)出了很好的分割能力。對于分割難度較高的彎曲主動脈(如第2行)，本文方法亦可進(jìn)行有效分割。

為了證明聯(lián)合標(biāo)簽融合策略的有效性，將基于多圖譜的主動脈分割方法與多數(shù)投票法(Majority voting，MV)、 STAPLE(Simultaneous truth and performance level estimation)、局部相似度權(quán)重融合法(Local weighted fusion，LW)以及JLF結(jié)合，進(jìn)行分割精度對比。為了對各分割算法的分割精度進(jìn)行衡量，本文選用Dice similarity coefficient(DSC)方法。DSC計算方法如下

(9)

式中：S為算法的分割結(jié)果，G為專家手工分割的結(jié)果。

圖4　3D主動脈CT圖像分割結(jié)果Fig.4　3D segmentation results with the proposed method

表1展示了4種方法對相同3D主動脈CT圖像的分割精度。本文方法的平均分割精度達(dá)到了0.904±0.018，高于STAPLE，MV和LW所取得的精度，算法穩(wěn)定性也優(yōu)于3種對比方法。

表1　4種標(biāo)簽融合策略的分割精度

在進(jìn)行圖像6和圖像8的分割時，發(fā)現(xiàn)這兩幅圖像與其他14個圖譜的相似度都低于閾值，因而需要進(jìn)行圖譜更新以提高圖譜攜帶的有用信息量。更新前后，兩幅圖像的分割精度如表2所示。表2結(jié)果表明，經(jīng)過圖譜更新操作后，兩幅圖像的分割精度均有4%左右的提升，在一定程度上表明了圖譜更新方法的有效性。

3　結(jié)束語

本文采用了基于多圖譜的醫(yī)學(xué)圖像分割方法，并首次將之與JLF策略相結(jié)合應(yīng)用于3D主動脈CT圖像的自動分割問題中。聯(lián)合標(biāo)簽融合策略在為每個圖譜分配權(quán)重時考慮了圖譜之間的相關(guān)性，對圖譜中干擾融合的冗余信息有更強(qiáng)的抑制作用，因而與傳統(tǒng)的標(biāo)簽融合策略相比有更高的融合準(zhǔn)確性。為了解決圖譜數(shù)量不足的問題，本文提出了一種圖譜更新策略以豐富圖譜包含的信息，為融合過程提供更多指導(dǎo)。實驗結(jié)果表明：本文的方法能夠?qū)?D主動脈CT圖像進(jìn)行有效的自動分割，與基于3種傳統(tǒng)融合策略的多圖譜分割方法相比，本文取得了更高的分割精度和算法穩(wěn)定性。在未來的工作中，計劃從3D主動脈圖像中提取更高層次的、區(qū)分性更強(qiáng)的特征，以將算法擴(kuò)展到解決有嚴(yán)重疾病的主動脈圖像自動分割問題中。

參考文獻(xiàn):

[1]Kochanek K D, Xu J, Murphy S L, et al. Deaths: Final data for 2009[J]. National Vital Statistics Reports: From the Centers for Disease Control and Prevention, National Center for Health Statistics, National Vital Statistics System, 2011, 60(3): 1-116.

[2]Ayed I B, Wang M, Miles B, et al. TRIC: Trust region for invariant compactness and its application to abdominal aorta segmentation[M]//Medical Image Computing and Computer-Assisted Intervention-MICCAI 2014. [S.l.]: Springer International Publishing, 2014: 381-388.

[3]Rueckert D, Burger P, Forbat S M, et al. Automatic tracking of the aorta in cardiovascular MR images using deformable models[J]. Medical Imaging, IEEE Transactions on, 1997, 16(5): 581-590.

[4]Zhuge F, Rubin G D, Sun S, et al. An abdominal aortic aneurysm segmentation method: Level set with region and statistical information [J]. Medical Physics, 2006,33(5):1440-1453.

[5]呂曉琪, 石靜, 任曉穎, 等. 一種基于水平集的三維肝臟磁共振圖像混合分割方法[J]. 數(shù)據(jù)采集與處理, 2015, 30(2): 350-358.

Lu Xiaoqi, Shi Jing, Ren Xiaoying, et al. Hybrid segmentation for 3D liver magnetic resonance imaging based on level set method[J]. Journal of Data Acquisition and Processing, 2015, 30(2): 350-358.

[6]唐利明, 黃大榮, 李可人. 基于變分水平集的圖像分割模型[J]. 數(shù)據(jù)采集與處理, 2014, 29(5): 704-712.

Tang Liming, Huang Darong, Li Keren. New model based on variational level set for image segmentation[J]. Journal of Data Acquisition and Processing, 2014, 29(5):704-712.

[8]I?gum I, Staring M, Rutten A, et al. Multi-atlas-based segmentation with local decision fusion-application to cardiac and aortic segmentation in CT scans[J]. Medical Imaging, IEEE Transactions on, 2009,28(7):1000-1010.

[9]何小海, 梁子飛, 唐曉穎, 等. 圖譜法腦部 MRI 圖像自動分割技術(shù)發(fā)展及應(yīng)用[J].數(shù)據(jù)采集與處理, 2015,30(5): 956-964.

He Xiaohai, Liang Zifei, Tang Xiaoying, et al. Development and application for atlas-based brain MRI image segmentation technology[J]. Journal of Data Acquisition and Processing, 2015,30(5):956-964.

[10] Collins D L, Pruessner J C. Towards accurate, automatic segmentation of the hippocampus and amygdala from MRI by augmenting ANIMAL with a template library and label fusion[J]. Neuroimage, 2010,52(4):1355-1366.

[11] Artaechevarria X, Munoz-Barrutia A, Ortiz-de-Solórzano C. Combination strategies in multi-atlas image segmentation: Application to brain MR data[J]. Medical Imaging, IEEE Transactions on, 2009,28(8):1266-1277.

[12] Heckemann R A, Hajnal J V, Aljabar P, et al. Automatic anatomical brain MRI segmentation combining label propagation and decision fusion[J]. NeuroImage, 2006,33(1):115-126.

[13] Wang H, Suh J W, Das S R, et al. Multi-atlas segmentation with joint label fusion[J]. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 2013,35(3):611-623.

[14] Sanroma G, Wu G, Gao Y, et al. Learning to rank atlases for multiple-atlas segmentation[J]. Medical Imaging,IEEE Transactions on, 2014,33(10):1939-1953.

[15] Elnakib A, Gimel′ Farb G, Suri J S, et al. Medical image segmentation: A brief survey[M]. New York: Springer, 2011:1-39.

[16] Avants B B, Tustison N, Song G. Advanced normalization tools (ANTS)[J]. Insight J, 2009,2:1-35.

[17] Avants B B, Epstein C L, Grossman M, et al. Symmetric diffeomorphic image registration with cross-correlation: Evaluating automated labeling of elderly and neurodegenerative brain[J]. Medical Image Analysis, 2008,12(1):26-41.

[18] Warfield S K, Zou K H, Wells W M. Simultaneous truth and performance level estimation (STAPLE): An algorithm for the validation of image segmentation[J]. Medical Imaging, IEEE Transactions on, 2004,23(7):903-921.