呂曉琪 LV Xiaoqi

李 娜 LI Na

張寶華 ZHANG Baohua

谷 宇 GU Yu

賈東征 JIA Dongzheng

內(nèi)蒙古科技大學(xué)信息工程學(xué)院 內(nèi)蒙古包頭 014010

基于體素相似性的圖像配準(zhǔn)[1,2]，是用兩圖像對應(yīng)體素對間的幾何相似性搜索全局最優(yōu)變換參數(shù)，避免了分割和特征提取的精度損失，并自動(dòng)完成配準(zhǔn)。目前，基于體素相似性的配準(zhǔn)主要有3種：①基于灰度差異，速度快，但只適用于單模下灰度差異較小的圖像間配準(zhǔn)；②基于相關(guān)系數(shù)，要求對應(yīng)體素的灰度值存在線性相關(guān)性，在多模配準(zhǔn)中會(huì)出現(xiàn)誤配；③基于信息論，以互信息為相似性測度，不依賴于圖像數(shù)據(jù)關(guān)系的假設(shè)，對多模圖像間的灰度關(guān)系也不用特殊假設(shè)，幾乎可用于任何多模圖像間配準(zhǔn)[3]。因此，基于體素相似性的互信息算法廣泛應(yīng)用于多模圖像配準(zhǔn)中[4-7]。

常用的配準(zhǔn)優(yōu)化方法有單純形法、Powell法、梯度下降法、模擬退火法、粒子群優(yōu)化法、遺傳算法等，其中前3種是局部優(yōu)化算法，收斂速度快，但易陷入局部極值；后3種是全局優(yōu)化算法，減少了陷入局部極值的概率，但計(jì)算量大，優(yōu)化速度慢。

本研究用Mattes互信息作為相似性測度以實(shí)現(xiàn)三維多模腦圖像配準(zhǔn)。此測度函數(shù)通常是不光滑的，容易受局部極值的影響，導(dǎo)致誤配準(zhǔn)。這就要求一種策略能利用測度函數(shù)連續(xù)和可微的特性，得到一個(gè)接近全局最優(yōu)的變換參數(shù)，并且計(jì)算速度快，還不影響配準(zhǔn)的準(zhǔn)確性。本研究詳細(xì)分析了互信息算法和多分辨率金字塔算法的計(jì)算原理，并將其有效結(jié)合，通過平滑濾波在一定程度上避免了局部極值，提高了配準(zhǔn)精度和魯棒性，使用下采樣減少了數(shù)據(jù)量，提高了配準(zhǔn)速度。

1 三維醫(yī)學(xué)圖像的配準(zhǔn)流程

1.1 三維配準(zhǔn)體數(shù)據(jù)的常用格式和獲得方法 三維配準(zhǔn)體數(shù)據(jù)的常用格式有Analyze文件、NifTI-1文件和MetaImage文件，可以通過兩種方式獲得體數(shù)據(jù)：讀取一系列二維切片構(gòu)成一個(gè)三維體數(shù)據(jù)，最后將體數(shù)據(jù)寫到指定文件中；利用醫(yī)學(xué)圖像轉(zhuǎn)換應(yīng)用程序包MRIConvert把二維 DICOM切片系列文件集轉(zhuǎn)換成需要的三維體數(shù)據(jù)。

本研究選用的三維體數(shù)據(jù)格式是MetaImage文件，由raw二進(jìn)制數(shù)據(jù)文件和mhd頭文件構(gòu)成，對三維圖像的存取是通過頭文件訪問對應(yīng)的數(shù)據(jù)文件。

1.2 三維配準(zhǔn)的一般流程 實(shí)現(xiàn)三維醫(yī)學(xué)圖像配準(zhǔn)需要以下基本模塊：空間變換、配準(zhǔn)測度、插值器和優(yōu)化器。其中，配準(zhǔn)測度是配準(zhǔn)框架中最關(guān)鍵的部分，用于評(píng)估配準(zhǔn)效果；插值器影響搜索空間的平滑度和整個(gè)配準(zhǔn)的執(zhí)行時(shí)間。

配準(zhǔn)組件用來協(xié)調(diào)整體，以確保在傳遞給優(yōu)化器之前準(zhǔn)備工作都已完成。Observer/Command類對象用來監(jiān)測優(yōu)化器并跟蹤配準(zhǔn)迭代過程，判斷優(yōu)化器是否正常工作，振幅長度是否合理，不必等優(yōu)化器自動(dòng)停止即可中斷配準(zhǔn)程序調(diào)整初始化參數(shù)，這樣就完成了整個(gè)配準(zhǔn)算法過程（圖1）。

圖1 三維配準(zhǔn)基本流程

2 基于互信息的配準(zhǔn)方法

2.1 以互信息為相似性測度 互信息的主要優(yōu)點(diǎn)是不需要指定其依靠的實(shí)際形式。因此，待配準(zhǔn)的兩幅圖像之間的復(fù)雜映射能夠模型化，這種靈活性使互信息很好地適應(yīng)了多模配準(zhǔn)的標(biāo)準(zhǔn)。本研究選擇了Mattes等[8]提出的互信息方法作為相似性度量。

醫(yī)學(xué)圖像配準(zhǔn)本質(zhì)上是一個(gè)函數(shù)優(yōu)化問題，即一系列變換參數(shù)μ最大化一個(gè)圖像相似度函數(shù)S：

本研究用互信息作為圖像的相似度函數(shù)，假設(shè)最大化相似度函數(shù)的這一組變換參數(shù){μopt}能使變換后的浮動(dòng)圖像很好地對齊到參考圖像上。公式（1）呈現(xiàn)的是最大化問題，但是Mattes互信息計(jì)算出負(fù)的交互信息，因此，實(shí)際上要最小化這個(gè)負(fù)的相似度函數(shù)S。假設(shè)LF和LR分別是浮動(dòng)和參考圖像的離散灰度集。參考和浮動(dòng)圖像之間的負(fù)互信息S表示為變換參數(shù)μ的函數(shù)如下：

其中，P、PF和PR分別是聯(lián)合概率分布、浮動(dòng)圖像邊緣概率分布和參考圖像邊緣概率分布，后面將會(huì)導(dǎo)出。

對于高維的空間變換參數(shù)，采用梯度準(zhǔn)則有助于尋找其最大值?；バ畔⒌奶荻裙饺缦拢?/p>

單個(gè)的梯度通過公式（2）對變換參數(shù)μ求微分得到：

用于計(jì)算互信息的概率分布是基于參考和浮動(dòng)圖像的邊緣和聯(lián)合直方圖。Parzen窗用來形成連續(xù)的基本圖像直方圖估計(jì)，也減少了插值量化和二進(jìn)制數(shù)據(jù)離散化的影響。因此，聯(lián)合概率分布是可微函數(shù)。用β（3）作為3次樣條Parzen窗口，β（0）作為0次樣條Parzen窗口，聯(lián)合離散型概率分布如下：

浮動(dòng)圖像的邊緣離散概率分布通過聯(lián)合分布計(jì)算如下：

參考圖像的邊緣離散型概率分布可以獨(dú)立于變換參數(shù)來計(jì)算，通過B樣條Parzen窗來滿足分區(qū)的聯(lián)合約束。計(jì)算公式如下：

2.2 初始變換參數(shù)的估計(jì) 三維多模醫(yī)學(xué)圖像配準(zhǔn)是使兩幅待配準(zhǔn)圖像在同一空間坐標(biāo)系中達(dá)到空間位置上的對齊，用Vr和Vf分別表示參考體和浮動(dòng)體，兩點(diǎn)（xr∈Vr，xf∈Vf）間的幾何映射可以表示為：

其中，xr和xf分別表示來自體數(shù)據(jù)vr和vf的空間坐標(biāo)的體素點(diǎn)，M是兩點(diǎn)間的幾何變換矩陣。本研究中的所有實(shí)驗(yàn)圖像都來自腦部，頭骨的剛性特性使得腦部被認(rèn)為是剛性的，因此，剛性變換矩陣M由3個(gè)平移和3個(gè)旋轉(zhuǎn)來定義，可以寫成一個(gè)平移向量T(t)和一個(gè)旋轉(zhuǎn)矩陣R連接的齊次坐標(biāo)形式：

因此，矢量表示的剛體變換公式（8）可以重寫為：

一個(gè)常見的描述旋轉(zhuǎn)矩陣的方法是以分解形式表示的一系列矩陣：

旋轉(zhuǎn)矩陣可以分解為R=RγRβRα。定義一個(gè)解剖坐標(biāo)系，以體圖像的中心為原點(diǎn)，X軸方向定義從右到左，Y軸方向從前到后，Z軸方向從顱低到顱頂。Rα、Rβ和Rγ分別表示圖像繞X軸、Y軸和Z軸正發(fā)生旋轉(zhuǎn)。

因此，參考體vr和浮動(dòng)體vf之間的相關(guān)位置由這樣一組參數(shù)決定，其中，tx、ty和tz是平移量，α、β和γ分別是浮動(dòng)體沿著相關(guān)的參考體三維軸旋轉(zhuǎn)的角度。從而，為了使浮動(dòng)體更好地對齊到參考體上，圖像配準(zhǔn)就轉(zhuǎn)化為尋找最優(yōu)參數(shù)組P的過程。

2.3 圖像多分辨率策略 多分辨率配準(zhǔn)框架需要一對圖像金字塔來平滑和下采樣圖像[9-11]。高斯金字塔、小波變換、拉普拉斯金字塔和steerable金字塔等常用于醫(yī)學(xué)圖像配準(zhǔn)。本研究使用的是高斯金字塔，其底層金字塔T0是原始圖像，通過低通濾波和下采樣獲取下一個(gè)金字塔層次T1，然后T1以同樣的方式濾波和下采樣獲得T2，重復(fù)以上步驟，生成剩下的金字塔層次結(jié)構(gòu)。顯然，這個(gè)金字塔結(jié)構(gòu)是通過迭代得到的，迭代公式如下：

其中，Tk(i, j)表示第k層金字塔圖像，i是圖像的列數(shù)，j是圖像的行數(shù)，k是金字塔的層數(shù)，p(m, n)是一個(gè)5×5的高斯模板窗口函數(shù)。

2.4 結(jié)合多分辨率的互信息配準(zhǔn) 互信息配準(zhǔn)函數(shù)通常是不光滑的，有多個(gè)局部最大值。局部最大值可能由以下原因引起：一部分局部最大值來自兩幅待配準(zhǔn)圖像的局部對齊；另一部分來自配準(zhǔn)算法執(zhí)行過程，如插值算法或者圖像重疊部分的變化等，這一部分是不合理的。面對龐大的三維體數(shù)據(jù)集，本研究把Mattes等提出的負(fù)互信息度量作為相似性測度的方法和多分辨率金字塔算法相結(jié)合。這種混合算法先將金字塔分解兩幅待配準(zhǔn)圖像，得到分層次多分辨率的圖像；接著對金字塔的粗糙尺度層進(jìn)行參數(shù)的最優(yōu)搜索，由于減小了圖像尺寸，從而加快了收斂速度；然后將粗糙水平的搜索結(jié)果用來初始化下一層較細(xì)尺度的圖像執(zhí)行配準(zhǔn)，該尺度層有了初始化估計(jì)，也加快了變換參數(shù)的搜索；最后重復(fù)以上過程，直到達(dá)到最細(xì)尺度層，得到最優(yōu)的配準(zhǔn)結(jié)果。此混合配準(zhǔn)算法的大部分迭代過程是在粗尺度層上進(jìn)行，并給出了初始化估計(jì)，因此，在細(xì)尺度層的配準(zhǔn)迭代減少，加快了配準(zhǔn)的收斂速度。

多分辨率策略是結(jié)合下采樣和平滑操作的一種圖像金字塔結(jié)構(gòu)，通過由粗糙到平滑的方式解決配準(zhǔn)問題。金字塔分解使平滑濾波后的圖像盡量保留所需要的全局有用信息，在一定程度上減少了局部極值的影響，提高了配準(zhǔn)的成功率、速度和魯棒性。

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

為了驗(yàn)證本研究提出的算法對多模三維體圖像的有效性，分別進(jìn)行了三維單模和多模配準(zhǔn)。采用3種算法進(jìn)行對比分析：算法①用均方差[12]作為度量，算法②用Mattes互信息作為度量，算法③用加入了多分辨率金字塔算法的Mattes互信息作為度量，即本研究提出的算法。

本研究實(shí)驗(yàn)程序用VC++2005實(shí)現(xiàn)，配準(zhǔn)程序中主要參數(shù)設(shè)置：收斂步長為0.0001，空間取樣數(shù)量用大概像素的1%，平均信息數(shù)的二進(jìn)制位數(shù)為64，最大迭代次數(shù)為500。

3.1 單模三維配準(zhǔn) 選取來自BrainWeb網(wǎng)站的圖像為實(shí)驗(yàn)對象。參考圖像為質(zhì)子密度（PD）加權(quán)的腦部MR體圖像，體數(shù)據(jù)為181×217×180，沿每個(gè)方向像素之間的間隔均為1.0 mm；浮動(dòng)圖像在參考圖像的基礎(chǔ)上繞原點(diǎn)旋轉(zhuǎn)20°，并在X軸方向上移動(dòng)20 mm。

分別用上述3種算法對兩幅待配準(zhǔn)圖像進(jìn)行配準(zhǔn)，由于從體圖像顯示上難以看出這3種算法的細(xì)微差別，所以僅對其中1種算法的配準(zhǔn)輸出進(jìn)行顯示來說明配準(zhǔn)問題。結(jié)果發(fā)現(xiàn)，本研究提出的算法比前2種算法有更高的精確度，同時(shí)運(yùn)行速度也得到了很大的提高（表1）。采用體三視圖對比顯示配準(zhǔn)輸出的結(jié)果發(fā)現(xiàn)，浮動(dòng)體圖像發(fā)生了很大的偏移和旋轉(zhuǎn)，浮動(dòng)體圖像和參考體圖像得到了很好的配準(zhǔn)（圖2）。

3.2 多模三維配準(zhǔn) 選取美國Vanderbilt大學(xué)圖像數(shù)據(jù)庫中的Practice組腦圖像數(shù)據(jù)進(jìn)行多模三維配準(zhǔn)，圖像包括一個(gè)患者的CT圖像數(shù)據(jù)、PET圖像數(shù)據(jù)和不同成像參數(shù)（T1、T2、PD加權(quán)，幾何失真校正的T1-rectif i ed、T2-rectif i ed、PD-rectif i ed）的MR圖像數(shù)據(jù)。該組數(shù)據(jù)被研究人員用來完成算法的初步評(píng)估（表2）。這組圖像提供了CT到MR、PET到MR配準(zhǔn)的標(biāo)準(zhǔn)剛體變換數(shù)據(jù)，該變換數(shù)據(jù)給出了浮動(dòng)體圖像8個(gè)頂點(diǎn)的原始物理空間位置坐標(biāo)和經(jīng)過配準(zhǔn)變換后的位置坐標(biāo)。因此依據(jù)這組數(shù)據(jù)，采用8點(diǎn)評(píng)估標(biāo)準(zhǔn)對本研究提出算法的精確性進(jìn)行評(píng)價(jià)。

表1 三種算法配準(zhǔn)結(jié)果對比

圖2 Mattes互信息算法單模配準(zhǔn)前后對比，上排為橫斷面對比圖，中排為矢狀面對比圖，下排為冠狀面對比圖。A~E分別為參考圖、原始浮動(dòng)圖、配準(zhǔn)后浮動(dòng)圖、配準(zhǔn)前差別圖、配準(zhǔn)后差別圖

表2 Practice組腦圖像數(shù)據(jù)的大小和像素間隔

經(jīng)驗(yàn)證，算法①進(jìn)行多模配準(zhǔn)時(shí)由于超過限定的迭代次數(shù)而溢出，配準(zhǔn)失敗，說明算法①只適合進(jìn)行單模配準(zhǔn)。因此，只用Mattes互信息法和改進(jìn)后的混合算法對CT和MR體圖像進(jìn)行配準(zhǔn)比較，MR作為參考體圖像，CT作為浮動(dòng)體圖像。把經(jīng)過配準(zhǔn)變換后的CT浮動(dòng)體圖像的8個(gè)頂點(diǎn)的空間坐標(biāo)qj,MR和標(biāo)準(zhǔn)的8個(gè)對應(yīng)點(diǎn)qj,ref按公式（13）求平均幾何距離△，用以評(píng)價(jià)配準(zhǔn)結(jié)果：

圖3 CT到不同成像參數(shù)的MR圖像配準(zhǔn)誤差

從圖3可以看出，使用算法②的配準(zhǔn)誤差中，只有CT-T1和CT-T2配準(zhǔn)模式誤差小于一個(gè)像素大小，其余模式由于陷入局部極值等原因，造成配準(zhǔn)結(jié)果出現(xiàn)誤差，導(dǎo)致精度較低。然而加入多分辨率金字塔的Mattes互信息算法采用多步分解配準(zhǔn)策略，在較低分辨率層次消除局部極值，從而提高了配準(zhǔn)的精度，增強(qiáng)了配準(zhǔn)的魯棒性，所有配準(zhǔn)誤差都小于一個(gè)像素大小，可認(rèn)為配準(zhǔn)結(jié)果達(dá)到了亞像素級(jí)標(biāo)準(zhǔn)。

圖4 本研究算法多模配準(zhǔn)前后對比 ，上排為橫斷面對比圖，中排為矢狀面對比圖，下排為冠狀面對比圖。A~E分別為最初CT浮動(dòng)圖、MR參考圖、配準(zhǔn)后CT圖、配準(zhǔn)前差別圖、配準(zhǔn)后差別圖

采用橫斷面、矢狀面、冠狀面三視圖分別對比顯示CT到MR配準(zhǔn)的結(jié)果。利用本研究算法可得到很好的配準(zhǔn)結(jié)果（圖4）。在相同實(shí)驗(yàn)條件下，由于算法②存在局部極值使配準(zhǔn)精度比本研究提出的算法精度要低得多，僅用橫斷面顯示算法②的配準(zhǔn)結(jié)果（圖5）。

通過圖4可知，參考體與浮動(dòng)體原始空間位置相差較大。對比橫斷面圖4C和圖5A、圖4E和圖5B可以得出，加入多分辨率金字塔的Mattes互信息算法使兩幅體圖像在空間位置上得到了很好的對準(zhǔn)。

圖5 Mattes互信息算法多模配準(zhǔn)橫斷面輸出圖。A. 配準(zhǔn)后CT圖；B. 配準(zhǔn)后差別圖

4 結(jié)束語

基于體素相似性的Mattes互信息算法是利用不同模態(tài)下的圖像灰度信息的統(tǒng)計(jì)特性來執(zhí)行配準(zhǔn)，不需要分割和特征提取，從而提高了配準(zhǔn)的精度和魯棒性，但是互信息配準(zhǔn)函數(shù)通常是不光滑的，在插值或者圖像重疊部分變化的過程中，容易陷入局部最優(yōu)，導(dǎo)致配準(zhǔn)的精度降低，甚至出現(xiàn)很大的誤配。多分辨率金字塔算法采用分層的、由粗到精的遞進(jìn)搜索，可以配合互信息算法有效地避免在收斂過程中陷入局部極值的情況。實(shí)驗(yàn)證明，用這種混合配準(zhǔn)算法對三維CT、MR圖像進(jìn)行多模配準(zhǔn)，兩幅待配準(zhǔn)圖像在空間位置上得到了很好的對齊，配準(zhǔn)的精度達(dá)到了亞像素級(jí)。

[1] Nyúl LG, Udupa JK, Saha PK. Incorporating a measure of local scale in voxel-based 3-D image registration. IEEE Trans Med Imaging, 2003, 22(2): 228-237.

[2] Loeckx D, Maes F, Vandermeulen D, et al. Voxel based nonrigid image registration using local and partial volume similarity measures. Proceedings 7th IEEE international symposium on biomedical imaging, 2010: 348-351.

[3] Maes F, Collignon A, Vandermeulen D, et al. Multimodality image registration by maximization of mutual information.IEEE Trans Med Imaging, 1997, 16(2): 187-198.

[4] 余慧婷, 張杰, 潘萌. 噪聲對三維圖像歸一化互信息配準(zhǔn)的影響. 中國醫(yī)學(xué)影像學(xué)雜志, 2011, 19(11): 844-849.

[5] 胡凱, 王衛(wèi)東, 邱本勝, 等. 顱頜面CT與MR圖像的配準(zhǔn).中國醫(yī)學(xué)影像學(xué)雜志, 2002, 10(2): 123-125.

[6] 劉晴, 郭希娟, 許慎洋. 基于互信息的N維多模醫(yī)學(xué)圖像配準(zhǔn). 中國圖象圖形學(xué)報(bào), 2009, 14(10): 2061-2068.

[7] Lee D, Hofmann M, Steinke F, et al. Learning similarity measure for multi-modal 3D image registration. IEEE Conference on CVPR. Tubingen: Germany, 2009: 186-193.

[8] Mattes D, Haynor DR, Vesselle H, et al. PET-CT image registration in the chest using free-form deformations. IEEE Trans Med Imaging, 2003, 22(1): 120-128.

[9] Thévenaz P, Unser M. Optimization of mutual information for multiresolution image registration. IEEE Trans Image Process,2000, 9(12): 2083-2099.

[10] 李喬亮, 汪國有, 劉建國, 等. 基于樣條金字塔和互信息的快速圖像配準(zhǔn). 計(jì)算機(jī)應(yīng)用研究, 2009, 26(5): 1949-1950, 1960.

[11] Bunting P, Labrosse F, Lucas R. A multi-resolution area-based technique for automatic multi-modal image registration. Image Vis Comput, 2010, 28(8): 1203-1219.

[12] 謝小輝, 楊光, 那奇, 等. 肝臟介入治療中基于B樣條變形形變模型的配準(zhǔn). 中國醫(yī)學(xué)影像學(xué)雜志, 2011, 19(8): 616-619.