朱文杰 丁 輝 黃 青 彭 江 王廣志#*

1(清華大學(xué)生物醫(yī)學(xué)工程系，北京 100084)2(解放軍總醫(yī)院全軍骨科研究所，北京 100853)

?

基于體素互信息的多模式和多時(shí)相圖像中股骨近端配準(zhǔn)

朱文杰1丁 輝1黃 青1彭 江2王廣志1#*

1(清華大學(xué)生物醫(yī)學(xué)工程系，北京 100084)2(解放軍總醫(yī)院全軍骨科研究所，北京 100853)

在股骨近端骨質(zhì)疏松進(jìn)程以及股骨頭壞死狀況評(píng)估方法中，圖像分析是常用的工具，通過(guò)不同時(shí)相以及不同模式的多組影像可以對(duì)病人病情進(jìn)行更全面的綜合評(píng)估。然而，在綜合評(píng)估過(guò)程中，由于病人多次在不同系統(tǒng)中成像，體位的差異使不同圖像組之間的解剖點(diǎn)位置無(wú)法一一對(duì)應(yīng)，因此分析之前需要將多組圖像對(duì)齊，才能觀(guān)察同一感興趣區(qū)在不同模式或不同時(shí)間骨組織狀況的差異。針對(duì)這個(gè)問(wèn)題，設(shè)計(jì)一種多模式、多時(shí)相圖像配準(zhǔn)的解決方案，通過(guò)圖像的前處理、雙閾值分類(lèi)并結(jié)合貝葉斯分類(lèi)的股骨分割得到股骨體素，然后通過(guò)基于歸一化互信息的圖像配準(zhǔn)獲得各組圖像中股骨之間的三維空間剛性變換矩陣，其中CT與MR圖像的配準(zhǔn)誤差在4 mm以下，CT與CT圖像的配準(zhǔn)誤差在2 mm以下。利用矩陣傳遞關(guān)系，以CT-CT多時(shí)相的配準(zhǔn)矩陣為基礎(chǔ)，可獲取任何兩組圖像間的變換矩陣。在此基礎(chǔ)上，再進(jìn)行任意兩組圖像的融合、點(diǎn)對(duì)點(diǎn)的分析以及骨質(zhì)狀況和血供狀況的定量評(píng)估。通過(guò)該方案，可以對(duì)多時(shí)相、多模式圖像分析中相同感興趣的區(qū)域進(jìn)行對(duì)比。

股骨近端；多時(shí)相；多模式；配準(zhǔn)；剛性變換矩陣

引言

股骨近端骨折、壞死是臨床常見(jiàn)的疾病。其中，股骨近端骨折包括了股骨頸、轉(zhuǎn)子間、粗隆骨折等，壞死一般指股骨頭局部壞死，表現(xiàn)為局部骨小梁的骨折、局部血供異常等。股骨近端骨折的誘因分外力壓迫以及病理轉(zhuǎn)移，而壞死通常是病理性原因造成。為研究局部骨質(zhì)變化，跟蹤病情的發(fā)展，以及評(píng)估生長(zhǎng)因子注射等新的微介入治療的效果，常采用多模式圖像對(duì)局部骨折、小梁塌陷、血供異常等狀況進(jìn)行評(píng)估。

在圖像分析評(píng)估中，CT圖像能清晰地顯示局部骨折的區(qū)域、骨質(zhì)變化以及骨質(zhì)缺損情況，而不同模式的磁共振(MR)圖像能清晰顯示出局部血供以及松質(zhì)骨內(nèi)組織代謝情況。因此，通過(guò)同一時(shí)相、不同模式的圖像聯(lián)合分析，可以對(duì)該時(shí)間點(diǎn)病人股骨近端的狀況進(jìn)行更全面的評(píng)估。而通過(guò)不同時(shí)相、同一模式的圖像，也可以在不同時(shí)間點(diǎn)對(duì)病人股骨近端的狀況進(jìn)行跟蹤評(píng)估。然而，由于多次成像的圖像組的尺寸、空間分辨率、空間坐標(biāo)參考系以及病人的體位都可能不完全一樣，造成了不同組圖像中股骨的空間位置不是一一對(duì)應(yīng)的，無(wú)法直接對(duì)不同組圖像中的同一感興趣區(qū)進(jìn)行對(duì)比，因此需要對(duì)不同時(shí)相、不同模式的多組三維圖像數(shù)據(jù)集進(jìn)行配準(zhǔn)。

關(guān)于股骨近端的配準(zhǔn)，國(guó)內(nèi)外學(xué)者已經(jīng)開(kāi)展了研究。Janet Blumenfeld和Colin Studholme等人為分析股骨近端松質(zhì)骨的各種參數(shù)，基于互信息的配準(zhǔn)方法，對(duì)三維(3D)的高分辨率MR圖像進(jìn)行配準(zhǔn)[1]；Wenjun Li、Irina Kezele、D. Louis Collins等人為分析股骨近端骨質(zhì)在藥物治療后，外部壓力作用后以及隨年齡增長(zhǎng)而發(fā)生的改變，基于區(qū)域生長(zhǎng)分割以及互信息配準(zhǔn)的方法對(duì)股骨近端進(jìn)行了配準(zhǔn)[2-3]；Ferdinando等人采用有界的最近點(diǎn)迭代方法，采用4自由度(4-DOF)變換矩陣對(duì)CT圖像進(jìn)行配準(zhǔn)[4]；張建國(guó)和楊慶銘為解決股骨頭缺血性壞死時(shí)病灶的準(zhǔn)確定位問(wèn)題，采用MC算法構(gòu)造了特征曲面，再利用特征點(diǎn)與特征面進(jìn)行配準(zhǔn)[5]。

在前人的研究基礎(chǔ)上，筆者設(shè)計(jì)了一套股骨近端多時(shí)相、多模式的圖像配準(zhǔn)方案。首先利用圖像分割技術(shù)提取股骨近端的主體信息，再基于歸一化互信息的方法實(shí)現(xiàn)多時(shí)相或者多模式圖像中股骨近端的配準(zhǔn)，接著將得到的變換矩陣通過(guò)矩陣變換傳遞通路(下文會(huì)提到)計(jì)算得到任意2組圖像之間的變換矩陣?；谶@些矩陣，可將各組圖像的股骨近端對(duì)齊，以便定量觀(guān)察不同圖像組之間相同的解剖位置的骨組織狀況，實(shí)現(xiàn)圖像組的融合分析。

1 方法

1.1 變換矩陣傳遞通路

由于股骨是剛性結(jié)構(gòu)，因此采用6自由度(6-DOF)剛性變換矩陣構(gòu)建兩組圖像之間的空間變換關(guān)系。特別需要注意的是，一組圖像中有左股骨與右股骨，由于成像中體位變化，人體姿態(tài)不可能完全一致。因此，需要2個(gè)剛性變換矩陣來(lái)分別表示兩組圖像之間左、右股骨的空間變換關(guān)系，同時(shí)需要將圖像數(shù)據(jù)分為左、右側(cè)分別配準(zhǔn)。

關(guān)于圖像配準(zhǔn)時(shí)參考基準(zhǔn)圖像的選擇問(wèn)題，筆者采取下述方案：方案1，當(dāng)不同時(shí)相的CT圖之間進(jìn)行配準(zhǔn)時(shí)，以最早時(shí)間點(diǎn)的CT圖作為參考基準(zhǔn)；方案2，當(dāng)同一時(shí)相的CT與MR圖進(jìn)行配準(zhǔn)時(shí)，由于股骨在CT圖中的結(jié)構(gòu)最為清晰，因此以CT圖像為參考基準(zhǔn)；方案3，當(dāng)不同時(shí)相的MR圖之間進(jìn)行配準(zhǔn)以及不同時(shí)相的MR圖與CT圖進(jìn)行配準(zhǔn)時(shí)，通過(guò)方案1、2得到的變換矩陣可計(jì)算出方案3實(shí)際的變換矩陣，該計(jì)算的過(guò)程需要變換矩陣傳遞通路的輔助。兩個(gè)時(shí)相下兩種模式圖像的變換矩陣傳遞通路如圖1所示(以左股骨為例)。

圖1 變換矩陣傳遞通路Fig.1 The transfer pathways of transformation

為得到任意兩組圖像之間的配準(zhǔn)關(guān)系，需要求出圖1中所有的變換矩陣。在配準(zhǔn)情況1與2下，可以得到變換矩陣TMR1→CT1、TCT2→CT1和TMR2→CT2，對(duì)于情況3的其他變換矩陣，可以按照矩陣傳遞通路由前3個(gè)變換矩陣獲得。注意到6DOF-剛性變換矩陣是可逆的，可以得到

(1)

TMR2→CT1=TMR2→CT2×TCT2→CT1

(2)

(3)

用這種矩陣傳遞的思路獲取任何兩組圖像之間的變換關(guān)系，主要基于以下考慮：首先，股骨在MR圖像中的對(duì)比度比在CT圖像中的對(duì)比度要弱，同時(shí)MR圖像在Z軸上(層片方向)的空間分辨率一般都比較低，因此MR圖像直接與MR圖像配準(zhǔn)的精度較低。此外，同一時(shí)相的股骨形態(tài)變化很小，因此多模圖像配準(zhǔn)時(shí)選取同一時(shí)相的圖像進(jìn)行配準(zhǔn)效果會(huì)更好。而對(duì)于不同時(shí)相的圖像序列來(lái)說(shuō)，由于CT圖像的骨結(jié)構(gòu)很清晰，用不同時(shí)相的CT圖像進(jìn)行配準(zhǔn)可以得到更準(zhǔn)確的結(jié)果。因此，對(duì)于不同時(shí)相、不同模式的圖像組，應(yīng)以不同時(shí)相的CT圖之間的變換矩陣作為橋梁推導(dǎo)得出。

1.2 變換矩陣的求解方案

上述的思路需要求出N-1個(gè)相互獨(dú)立的變換矩陣，而求解變換矩陣則采用下述思路：首先，通過(guò)圖像的前處理以及分割，獲得左股骨(或右股骨)的主體結(jié)構(gòu)信息，然后基于圖像體素的歸一化互信息實(shí)現(xiàn)配準(zhǔn)，得到相應(yīng)的變換矩陣，其具體流程如圖2所示。

圖2 變換矩陣求解流程Fig.2 Flow chart of solving transformation matrix

1.2.1 圖像前處理

首先是規(guī)范化坐標(biāo)。從不同成像設(shè)備獲得的圖像坐標(biāo)參考系以及圖像原點(diǎn)可能不同，筆者將坐標(biāo)統(tǒng)一為L(zhǎng)PS坐標(biāo)系統(tǒng)。

接著是圖像的重采樣。采用最近鄰的方法進(jìn)行體素的重采樣，將圖像統(tǒng)一到一致的分辨率。Ires表示重采樣的圖，Iori表示原圖，Sres表示重采樣空間分辨率，Sori表示原圖空間分辨率，Round是最近鄰函數(shù)。其中，Pres∈Ires和Pori∈Iori分別表示重采樣后的點(diǎn)以及原圖的點(diǎn)。那么，重采樣后點(diǎn)坐標(biāo)為

(4)

然后，為了初步去除大部分干擾信息，需要對(duì)圖像進(jìn)行剪裁，分別截取包含左、右股骨和部分髖骨的圖像數(shù)據(jù)子集。在圖像剪裁后，圖像的原點(diǎn)已經(jīng)發(fā)生改變，假設(shè)S表示空間分辨率，剪裁的起點(diǎn)坐標(biāo)為(x0,y0,z0)，那么相應(yīng)的變換矩陣為

(5)

最后，根據(jù)圖像的噪聲類(lèi)型進(jìn)行圖像去噪處理。采用高斯平滑去除高斯噪聲，用中值濾波去除椒鹽噪聲。三維高斯平滑的模板函數(shù)和中值濾波函數(shù)分別為

(6)

(7)

式中，I(x,y,z)表示圖像在(x,y,z)處的灰度值，d是濾波半徑，M是求中值函數(shù)。

1.2.2 圖像中股骨的分割

股骨在CT圖像中表現(xiàn)為高亮度信號(hào)，在MR的T2圖像中表現(xiàn)為低亮度信號(hào)，采用雙閾值的方法能粗略獲取股骨的大部分結(jié)構(gòu)信息，基本上排除其他的不相關(guān)組織。在雙閾值分割后，部分結(jié)構(gòu)噪聲仍存留在股骨結(jié)構(gòu)中，采用數(shù)學(xué)形態(tài)學(xué)[6]的開(kāi)運(yùn)算去除散點(diǎn)，再用閉運(yùn)算連接股骨邊緣缺口，有

(8)

式中，Imath是數(shù)學(xué)形態(tài)學(xué)運(yùn)算后圖像，Ithre是雙閾值分割后的結(jié)果，SE是運(yùn)算模板。

股骨與髖骨在生理解剖上是不連接的，中間隔了一層髖關(guān)節(jié)軟骨。但CT圖像分割時(shí)在雙閾值寬度較大的情況下，股骨、髖關(guān)節(jié)和髖骨經(jīng)常會(huì)連到一起。由于股骨是多運(yùn)動(dòng)自由度的結(jié)構(gòu)，成像時(shí)其相對(duì)髖骨的位置往往又是變化不定的。因此，如果保留髖骨部分，會(huì)對(duì)后面的股骨配準(zhǔn)帶來(lái)很大的誤差，髖骨的去除顯得尤為重要，而去除髖骨的關(guān)鍵就是去除分割出來(lái)的髖關(guān)節(jié)部分。本研究采用Cheng等提出的貝葉斯分類(lèi)的方法[7]解決了這個(gè)問(wèn)題，很好地去除了股骨與髖骨的連接部分，得到獨(dú)立的左右股骨。

接著通過(guò)交互可以將股骨提取出來(lái)，并對(duì)提取的股骨進(jìn)行填充，以減少配準(zhǔn)時(shí)邊緣的權(quán)重。

通過(guò)以上處理，可以從一組三維圖像數(shù)據(jù)集中得到左側(cè)股骨和右側(cè)股骨的兩組獨(dú)立的體素?cái)?shù)據(jù)，用于進(jìn)一步的配準(zhǔn)。

1.2.3 股骨的圖像配準(zhǔn)

1.2.3.1 平移粗配

在統(tǒng)一坐標(biāo)系、重采樣以及圖像剪裁等操作后，需要配準(zhǔn)的股骨與固定的股骨位置可能相距還比較遠(yuǎn)，為了防止配準(zhǔn)時(shí)陷入局部最優(yōu)以及提高配準(zhǔn)效率，把兩組股骨體素?cái)?shù)據(jù)通過(guò)平移實(shí)現(xiàn)重心的重合，在此基礎(chǔ)上再進(jìn)行精細(xì)配準(zhǔn)。假設(shè)基準(zhǔn)圖像重心坐標(biāo)為(x1,y1,z1)，配準(zhǔn)圖像的重心坐標(biāo)為(x2,y2,z2)，空間分辨率為S，那么平移粗配的變換矩陣Tregist-roughly為

(9)

1.2.3.2 降采樣配準(zhǔn)

圖3 配準(zhǔn)框圖Fig.3 chart of registration

在本研究中，精細(xì)配準(zhǔn)的相似性測(cè)度采用歸一化互信息，空間搜索策略采用Powell算法，收斂?jī)?yōu)化采用降采樣策略，最后的配準(zhǔn)效果與迭代次數(shù)有關(guān)，處理的流程如圖3所示。

互信息可用于描述兩個(gè)系統(tǒng)間的相關(guān)性，通常用“熵”表示，常作為多模圖像配準(zhǔn)的相似性測(cè)度[8]。H(A)表示A的熵值，H(B)表示B的熵值，H(A,B)表示聯(lián)合熵，那么歸一化互信息(NMI)[9]為

(10)

Powell算法[10]是一種空間迭代搜索算法，核心思想是在每個(gè)方向依次找到極值點(diǎn)，然后進(jìn)行迭代，直至滿(mǎn)足迭代停止條件。

為優(yōu)化配準(zhǔn)中的優(yōu)化速度，采用降采樣策略進(jìn)行處理，先對(duì)降采樣的低分辨率圖像粗配，然后逐級(jí)提高分辨率達(dá)到精細(xì)的配準(zhǔn)。假定Ifix,D和Imove,D分別表示基準(zhǔn)圖像和配準(zhǔn)圖像經(jīng)過(guò)降采D倍后得到的圖像，Tregist,D表示在圖像在降采D倍后經(jīng)過(guò)配準(zhǔn)得到的變換矩陣。那么，D倍降采迭代得到的配準(zhǔn)后變換矩陣為

(11)

在新一輪降采迭代配準(zhǔn)中，降采系數(shù)降到D/2，那么初始的變換矩陣Tregist,D/2為

(12)

當(dāng)降采系數(shù)降到1時(shí)，得到最終的精細(xì)配準(zhǔn)矩陣Tregist-fine為

(13)

綜合配準(zhǔn)過(guò)程中的所有變換矩陣，兩組圖像之間的實(shí)際變換矩陣Tregist為

(14)

2 結(jié)果

2.1 配準(zhǔn)結(jié)果

所采用的圖像數(shù)據(jù)來(lái)自臨床，其中CT圖像由SIEMENSEmotion16成像系統(tǒng)得到。多序列的磁共振圖像由GEDiscoveryMR750 3.0TMRI成像系統(tǒng)得到，采集了T1、T2、動(dòng)態(tài)增強(qiáng)的掩模圖像，以及隨后按一定時(shí)間間隔采集的對(duì)比圖像等多個(gè)序列的三維圖像數(shù)據(jù)集，分別用于骨結(jié)構(gòu)與組織代謝功能的分析。其中，注射造影劑后的動(dòng)態(tài)增強(qiáng)圖像，共有1 200張圖片，分30個(gè)時(shí)間點(diǎn)采集，反映股骨近端血供的情況。通過(guò)配準(zhǔn)可以對(duì)CT圖像中感興趣區(qū)域?qū)?yīng)的血供狀況進(jìn)行跟蹤評(píng)估。各組圖像的成像時(shí)間、成像范圍和空間分辨率信息如表1所示。

表1 圖像組信息

圖4列出了沒(méi)有對(duì)齊的原圖像(以第一次成像為例)、不同時(shí)相CT圖像的配準(zhǔn)結(jié)果、同一時(shí)相不同模式圖像的配準(zhǔn)結(jié)果，以及通過(guò)變換矩陣傳遞通路后其他任意2組圖像的配準(zhǔn)結(jié)果。其中，配準(zhǔn)的結(jié)果采用棋盤(pán)格的方法展示兩圖像的對(duì)比。

在圖像采集時(shí)，CT1圖像與其對(duì)應(yīng)的MR圖像時(shí)間上相差1天，可以認(rèn)為是同一時(shí)相的圖像組，其股骨形態(tài)變化很小。從圖4可以明顯看到，不同時(shí)相的CT圖像之間的配準(zhǔn)程度最好，其次是同一時(shí)相CT與MR的配準(zhǔn)。CT與MR圖像的多模配準(zhǔn)準(zhǔn)確度取決于MR圖像的質(zhì)量，比如MR12與MR11圖像中股骨的對(duì)比度高，配準(zhǔn)精度較高。對(duì)于MR13與MR14，股骨在圖像中相對(duì)周?chē)M織的對(duì)比度較低，而且空間分辨率低，Z軸方向的空間體素點(diǎn)間隔達(dá)6.0 mm，因此直接進(jìn)行配準(zhǔn)的誤差會(huì)比較大。在實(shí)際中，由于所有磁共振圖像是在一次掃描中得到的，在掃描中要求病人保持體位不變，因此可以通過(guò)一系列磁共振圖像中最容易配準(zhǔn)的一個(gè)序列的配準(zhǔn)，得到所有序列的空間變換矩陣。

圖4(c)～(f)是時(shí)相1的各組圖像經(jīng)過(guò)矩陣傳遞通路后得到的配準(zhǔn)結(jié)果，都能達(dá)到了配準(zhǔn)的要求，由此驗(yàn)證了矩陣通路的思路是可行的。在此基礎(chǔ)上可以做更多時(shí)相以及多模式的股骨近端病情評(píng)估。

為了定量衡量配準(zhǔn)結(jié)果，按照文獻(xiàn)[3]的方法，選取CT1中比較有特征的解剖點(diǎn)(見(jiàn)圖5)，從左到右分別為方形結(jié)節(jié)(quadrate tubercle)、轉(zhuǎn)子窩(trochanteric fossa)、股骨頭中央凹(fovea capitis)。然后在CT2、MR11、MR12、MR13、MR14圖中選擇相應(yīng)的解剖點(diǎn)，并比較通過(guò)配準(zhǔn)后解剖點(diǎn)的位置誤差，誤差分析見(jiàn)表2?？梢钥闯觯錅?zhǔn)前的平均誤差都在45 mm以上，而配準(zhǔn)后的平均誤差在4 mm以下，其中同一時(shí)相的CT圖之間的配準(zhǔn)誤差最小，在2 mm以下。

圖4 各圖像組配準(zhǔn)前后比較。(a)CT1、CT2、MR11、MR12、MR13、MR14沒(méi)有對(duì)齊時(shí)的原圖；(b)CT1與 CT2、MR11、MR12、MR13、MR14的配準(zhǔn)結(jié)果；(c)不同時(shí)相CT2與MR11、MR12、MR13、MR14在經(jīng)過(guò)矩陣傳遞通路后的配準(zhǔn)結(jié)果；(d)同一時(shí)相MR12與MR11、MR13、MR14在經(jīng)過(guò)矩陣傳遞通路后的配準(zhǔn)結(jié)果；(e)同一時(shí)相MR11和MR13、MR14在經(jīng)過(guò)矩陣傳遞通路后的配準(zhǔn)結(jié)果；(f)相同時(shí)相M1R3與MR14在經(jīng)過(guò)矩陣傳遞通路后的配準(zhǔn)結(jié)果Fig.4 The comparison of images before and after registration.(a)original images of CT1,CT2,MR11,MR12,MR13 and MR14 without alignment (b) registration results of CT1 and CT2,MR11,MR12,MR13,MR14. (c) registration results of CT2 and MR11,MR12,MR13,MR14 with different phases through matrix transfer pathways,(d)registration results of MR12 and MR11,MR13,MR14 with the same phase through matrix transfer pathways,(e)registration results of MR11 and MR13,MR14 with the same phase through matrix transfer pathways,(f)registration result of MR13 and MR14 with the same phase through matrix transfer pathways

圖5 定量衡量配準(zhǔn)結(jié)果所選取的解剖點(diǎn)。 (a)方形結(jié)節(jié)；(b)轉(zhuǎn)子窩；(c)股骨頭中央凹Fig.5 Anatomical landmarks for quantitative measurement of registration results. (a)Quadrate tubercle; (b)Trochanteric fossa; (c)fovea capitis

配準(zhǔn)前后圖像組方形結(jié)節(jié)誤差/mm轉(zhuǎn)子窩誤差/mm股骨頭中央凹/mm平均誤差/mm配準(zhǔn)前配準(zhǔn)后CT2911900792868MR111876193020771961MR12506559648571MR132257223023402276MR14388415604469CT227240820MR11374026835MR1224192723MR1318455439MR1423202422

2.2 圖像融合結(jié)果

通過(guò)變換矩陣傳遞通路，即可得到任意兩組圖像之間股骨的變換關(guān)系，借分割得到股骨的體素?cái)?shù)據(jù)集并以此作為掩模，可以實(shí)現(xiàn)任何兩組圖像組之間股骨的融合。圖6展示了CT2數(shù)據(jù)與對(duì)應(yīng)的MR12圖像的融合結(jié)果，從橫斷面和冠狀面進(jìn)行觀(guān)察，其中(a)～(d)從橫斷面進(jìn)行觀(guān)察，(e)～(j)從冠狀面進(jìn)行觀(guān)察，各圖中左邊高亮的部分來(lái)自CT2的股骨，其他部分來(lái)自MR T2圖像。通過(guò)融合圖像，可以看到股骨與周?chē)M織邊緣銜接性很好，從多角度、多層面看是對(duì)齊的，在此基礎(chǔ)上可以對(duì)股骨骨質(zhì)及其周?chē)慕M織結(jié)構(gòu)做綜合評(píng)估。

圖6 圖像融合結(jié)果。 (a)～(d)橫斷面下的融合結(jié)果； (e)～(j)冠狀面下的融合結(jié)果Fig.6 The results of image fusion. (a)～(d) Fusion results in axial view; (e)～(j) are fusion results in coronal view

2.3 對(duì)應(yīng)點(diǎn)配準(zhǔn)結(jié)果

在得到所有圖組之間的變換矩陣后，以其中一組圖像作為基準(zhǔn)，其他所有圖像組的點(diǎn)通過(guò)變換矩陣與基準(zhǔn)圖像組的點(diǎn)一一對(duì)應(yīng)。以時(shí)相1的CT、MR圖像以及時(shí)相2的CT圖像配準(zhǔn)為例，選取了4個(gè)特征明顯的點(diǎn)的位置，從橫斷面和冠狀面對(duì)配準(zhǔn)后的點(diǎn)的位置來(lái)觀(guān)察對(duì)應(yīng)點(diǎn)的配準(zhǔn)結(jié)果，其對(duì)應(yīng)結(jié)果如圖7所示。圖中第1行是從橫斷面進(jìn)行觀(guān)察，第2行是從冠狀面進(jìn)行觀(guān)察，從左到右每列依次為CT1、CT2、MR11、MR12、MR13和MR14。雖然各組圖像在尺寸、空間分辨率以及病人體位都不一致，但經(jīng)過(guò)上述多時(shí)相和多模式的配準(zhǔn)處理后，各組圖像之間的生理解剖點(diǎn)能夠?qū)?yīng)到一起，這是后續(xù)定量評(píng)估的基礎(chǔ)。

2.4 觀(guān)察血供狀況

在實(shí)現(xiàn)配準(zhǔn)以后，采用CT1與MR12、MR22(動(dòng)態(tài)增強(qiáng)像)圖像對(duì)兩側(cè)股骨進(jìn)行血供情況分析。其中，MR12和MR22圖像均為30個(gè)時(shí)相、共2 400張圖，該30個(gè)時(shí)相可以表征股骨近端血供隨時(shí)間的變化。CT1與MR12的點(diǎn)配準(zhǔn)結(jié)果見(jiàn)圖8(a)，CT1與MR22的點(diǎn)配準(zhǔn)結(jié)果見(jiàn)圖8(b)，相應(yīng)點(diǎn)的增強(qiáng)曲線(xiàn)見(jiàn)圖8(c)、(d)。其中，標(biāo)號(hào)為11、15的點(diǎn)是CT圖像中骨質(zhì)缺損明顯的解剖點(diǎn)。

通過(guò)多模配準(zhǔn)，可以觀(guān)察股骨近端在CT圖像中感興趣區(qū)域的血供狀況。其中，點(diǎn)11、15顯示出了該部位明顯增強(qiáng)所反映的血供異常。通過(guò)MR22與MR12的曲線(xiàn)比較，可以發(fā)現(xiàn)15號(hào)解剖點(diǎn)在時(shí)相2時(shí)相對(duì)時(shí)相1的血供量減少。

通過(guò)多模配準(zhǔn)，可以觀(guān)察股骨近端在CT圖像中感興趣區(qū)域的血供狀況。其中，點(diǎn)11、15顯示出了該部位明顯增強(qiáng)所反映的血供異常。

圖8 CT1與MR12、MR22的配準(zhǔn)結(jié)果和感興趣區(qū)域的血供分析。 (a)CT1與MR12配準(zhǔn)后18個(gè)點(diǎn)的對(duì)應(yīng)關(guān)系； (b)CT1與MR22配準(zhǔn)后18個(gè)點(diǎn)的對(duì)應(yīng)關(guān)系； (c)MR12的18個(gè)點(diǎn)的血供曲線(xiàn)； (d)MR22的18個(gè)點(diǎn)的血供曲線(xiàn)Fig.8 Registration result of CT1，MR12 and MR22 and analysis of blood supply from interesting places. (a)Correspondence of 18 points after registration of CT1 and MR12; (b) Correspondence of 18 points after registration of CT1 and MR22; (c)Blood supply curve of 18 points from MR12; (d) blood supply curve of 18 points from MR22

3 討論

本研究的核心思想是：首先通過(guò)對(duì)三維圖像數(shù)據(jù)的分割，提取股骨體素的信息，進(jìn)一步通過(guò)體素互信息計(jì)算將股骨進(jìn)行配準(zhǔn)，得到不同模式或不同時(shí)相數(shù)據(jù)集中股骨間的變換矩陣，接著通過(guò)變換矩陣的傳遞通路得到所有圖像之間的空間變換關(guān)系，在此基礎(chǔ)上再進(jìn)行后續(xù)的圖像融合、點(diǎn)對(duì)點(diǎn)分析以及定量評(píng)估。

本研究的分割采用雙閾值的分割方法而沒(méi)有采用基于梯度或能量信息的分割方法(如Canny[11]、GVF、Snake[12]分割算法)，主要基于以下考慮：首先，本研究涉及多模圖像，在一些MR圖像中，股骨內(nèi)部以及邊緣信息與圖中其他部分組織結(jié)構(gòu)的邊緣信息相當(dāng)豐富，股骨邊緣確實(shí)檢測(cè)出來(lái)，但其他組織的邊緣也同時(shí)檢測(cè)出來(lái)，在處理中就成為了干擾信息。其次，本解決方案應(yīng)該盡量減少需要調(diào)整的參量，對(duì)于Snake等一些基于能量的算法，在圖像分割領(lǐng)域中的分割效果比較好，但需要調(diào)整彈性系數(shù)、剛性系數(shù)、迭代次數(shù)等參數(shù)，分割的靈活性不能令人滿(mǎn)意。此外，本研究最后的目的是配準(zhǔn)，雖然利用雙閾值分割的結(jié)果會(huì)引入一些難以剔除的誤差和干擾信息，特別是MR圖，但由于配準(zhǔn)是采用基于歸一化互信息對(duì)分割得到的三維體素間的關(guān)系進(jìn)行計(jì)算，所以在股骨內(nèi)部以及邊緣信息足夠的情況下，配準(zhǔn)的結(jié)果仍能達(dá)到要求。

4 結(jié)論

在股骨近端多模式、多時(shí)相圖像組分析時(shí)，股骨空間位置不對(duì)應(yīng)，難以進(jìn)行關(guān)聯(lián)分析。針對(duì)該問(wèn)題，筆者提出了一種解決方案，借此能獲得任意兩組圖像之間的左、右股骨配準(zhǔn)的變換矩陣，從而使股骨對(duì)齊，便于進(jìn)行后續(xù)的分析工作，如圖像融合、不同圖像組之間同一解剖點(diǎn)的觀(guān)察對(duì)比，以及通過(guò)體素強(qiáng)度值定量分析不同時(shí)相圖組中感興趣區(qū)的體素與組織結(jié)構(gòu)的變化。

[1]BlumenfeldJ,StudholmeC,Carballido-GamioJ, et al.Three-dimensionalimageregistrationofMRproximalfemurimagesfortheanalysisoftrabecularboneparameters[J].MedicalPhysics, 2008, 35(10): 4630-4639.

[2]LiWenjun,SodeM,SaeedI, et al.AutomatedregistrationofhipandspineforlongitudinalQCTstudies:integrationwith3Ddensitometricandstructuralanalysis[J].Bone, 2006,38:273-279.

[3]LiWenjun,KezeleI,CollinsDL, et al.Voxel-basedmodelingandquantificationoftheproximalfemurusinginter-subjectregistrationofquantitativeCTimages[J].Bone, 2007, 41(5): 888-895.

[4]BaenaFR,HawkeT,JakopecM.Aboundediterativeclosestpointmethodforminimallyinvasiveregistrationofthefemur[J].ProceedingsoftheInstitutionofMechanicalEngineers,PartH:JournalofEngineeringinMedicine, 2013,227(10):1135-1144.

[5] 張建國(guó)，楊慶銘. 股骨頭壞死病灶清除計(jì)算機(jī)輔助系統(tǒng)配準(zhǔn)模塊和導(dǎo)航模塊的建立[J]. 中國(guó)組織工程研究, 2011, 15(26): 4768-4772.

[6]SerraJ.Introductiontomathematicalmorphology[J].Computervision,Graphics,andImageProcessing, 1986, 35(3): 283-305.

[7]ChengYuanzhi,ZhouShengjun,WangYadong, et al.AutomaticsegmentationtechniqueforacetabulumandfemoralheadinCTimages[J].PatternRecognition, 2013, 46(11): 2969-2984.

[8]MaesF,CollignonA,VandermeulenD, et al.Multimodalityimageregistrationbymaximizationofmutualinformation[J].IEEETransactionsonMedicalImaging, 1997, 16(2): 187-198.

[9]StudholmeC,HillDLG,HawkesDJ.Anoverlapinvariantentropymeasureof3Dmedicalimagealignment[J].PatternRecognition, 1999, 32(1): 71-86.

[10]PowellMJD.Onsearchdirectionsforminimizationalgorithms[J].MathematicalProgramming, 1973, 4(1): 193-201.

[11]CannyJ.Acomputationalapproachtoedgedetection[J].IEEETransactionsonPatternAnalysisandMachineIntelligence, 1986 (6): 679-698.

[12]XuChenyang,PrinceJL.Snakes,shapes,andgradientvectorflow[J].IEEETransactionsonImageProcessing, 1998, 7(3): 359-369.

Registration for Multi-Modal and Multi-Phase Images of Proximal Femur Based on Voxel Mutual Information

Zhu Wenjie1Ding Hui1Huang Qing1Peng Jiang2Wang Guangzhi1#*

1(DepartmentofBiomedicalEngineering,TsinghuaUniversity,Beijing100084,China)2(InstituteofOrthopedicsinPLAGeneralHospital,Beijing100853,China)

Image analysis is a common method for evaluation of proximal femur osteoporosis and femoral head necrosis. This method evaluates patient’s condition based on images of proximal femur with different modes and phases. However, because the images are generated in different systems and the positions of the patient relative to systems are different, the positions of anatomical points in different images are not one-to-one consistent. We need to align points of images with different modes and phases before we can research on interesting area. To solve the problem, we proposed a solution to get spatial rigid transformation through image pre-processing, voxels segmentation of femur based on dual-threshold combined with Bayes decision rule, and femurs registration based on normalized mutual information. The error of CT-MR registration and CT-CT registration were below 4 mm and 2 mm respectively. Using the matrix transfer relationship, rigid transformation between any two images based on multi-phase CT-CT registration matrix was obtained. On this basis, image fusion of any two images, point-to-point analysis and quantitative evaluation of osseous and blood supply condition were processed. The interesting area in images with different modes and phases through the solution was compared.

proximal femur; multi-phase; multi-modal; registration; rigid transformation

10.3969/j.issn.0258-8021. 2015. 05.001

2015-06-15， 錄用日期:2015-08-18

清華大學(xué)實(shí)驗(yàn)室創(chuàng)新基金；北京市科技項(xiàng)目(Z131100006413027)；國(guó)家自然科學(xué)基金(81127003，51361130032)

R318

A

0258-8021(2015) 05-0513-09

# 中國(guó)生物醫(yī)學(xué)工程學(xué)會(huì)會(huì)員(Member, Chinese Society of Biomedical Engineering)

*通信作者(Corresponding author)， E-mail: wgz-dea@mail.tsinghua.edu.cn