呂曉琪 LV Xiaoqi

李 娜 LI Na

張寶華 ZHANG Baohua

谷 宇 GU Yu

賈東征 JIA Dongzheng

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

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

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

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

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

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

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

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

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

圖1 三維配準基本流程

2 基于互信息的配準方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.2 多模三維配準 選取美國Vanderbilt大學(xué)圖像數(shù)據(jù)庫中的Practice組腦圖像數(shù)據(jù)進行多模三維配準，圖像包括一個患者的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ù)被研究人員用來完成算法的初步評估（表2）。這組圖像提供了CT到MR、PET到MR配準的標準剛體變換數(shù)據(jù)，該變換數(shù)據(jù)給出了浮動體圖像8個頂點的原始物理空間位置坐標和經(jīng)過配準變換后的位置坐標。因此依據(jù)這組數(shù)據(jù)，采用8點評估標準對本研究提出算法的精確性進行評價。

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

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

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

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

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

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

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

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

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

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

4 結(jié)束語

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

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

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

[6] 劉晴, 郭希娟, 許慎洋. 基于互信息的N維多模醫(yī)學(xué)圖像配準. 中國圖象圖形學(xué)報, 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] 李喬亮, 汪國有, 劉建國, 等. 基于樣條金字塔和互信息的快速圖像配準. 計算機應(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樣條變形形變模型的配準. 中國醫(yī)學(xué)影像學(xué)雜志, 2011, 19(8): 616-619.