多模態(tài)醫(yī)學(xué)影像的非剛體配準(zhǔn)與多分辨率融合方法

張欣怡，張富利，王秋生

1.北京航空航天大學(xué)自動化科學(xué)與電氣工程學(xué)院，北京100191 2.中國人民解放軍陸軍總醫(yī)院，北京100700

醫(yī)學(xué)影像處理是多模態(tài)圖像處理中的一個重要分支，在醫(yī)學(xué)成像技術(shù)中，不同成像設(shè)備下的醫(yī)學(xué)影像呈現(xiàn)出不同模態(tài)，如計算機(jī)斷層成像（CT）、磁共振成像（MRI）、正電子發(fā)射斷層成像（PET）等[1]。不同模態(tài)提供了不同的人體組織信息，單一模態(tài)的醫(yī)學(xué)影像無法為放療醫(yī)生提供足夠的病灶信息，因此往往需要將不同模態(tài)的影像信息進(jìn)行綜合處理。醫(yī)學(xué)影像配準(zhǔn)是醫(yī)學(xué)影像融合的先決條件，由于醫(yī)學(xué)影像受成像時間、成像設(shè)備、患者姿勢等因素的影響，多模態(tài)影像之間難以在空間位置上配準(zhǔn)[2]。在配準(zhǔn)的基礎(chǔ)上，為了能夠更加全面、可靠地反映人體各組織器官以及病灶部位的信息，需要對多模態(tài)醫(yī)學(xué)影像進(jìn)行高效的融合處理。

在配準(zhǔn)過程中，用作參考的影像為固定影像，進(jìn)行變換的影像為浮動影像。醫(yī)學(xué)影像配準(zhǔn)的過程就是為浮動影像尋找一個最佳變換，使之與固定影像達(dá)到空間上的匹配。人體內(nèi)部組織結(jié)構(gòu)復(fù)雜且具有時變特性，例如肺部掃描影像中的組織器官會隨著患者的呼吸而位移，因此剛體配準(zhǔn)無法滿足要求，需要選擇合適的變換對醫(yī)學(xué)影像進(jìn)行復(fù)雜的非剛體配準(zhǔn)；同時，不同模態(tài)的醫(yī)學(xué)影像難以提取共同特征，因此無法基于特征來評價影像之間的配準(zhǔn)程度。醫(yī)學(xué)影像的融合不僅僅是影像之間的簡單疊加，而需要考慮不同影像間冗余信息的處理以及互補(bǔ)信息的最大化保留，并通過一系列指標(biāo)來評價融合影像的信息豐富程度。在多模態(tài)醫(yī)學(xué)影像融合方法的選擇上，多分辨率的融合方法能分別處理不同尺度下的細(xì)節(jié)信息，但需要分別選擇側(cè)重影像概貌或細(xì)節(jié)信息的融合規(guī)則。

本文以互信息為相似性測度，搭建仿射變換與B樣條變換相結(jié)合的多層次配準(zhǔn)框架，實(shí)現(xiàn)多模態(tài)醫(yī)學(xué)影像的非剛體配準(zhǔn)。同時，對配準(zhǔn)后的影像按不同融合規(guī)則進(jìn)行基于小波變換的多分辨率融合，實(shí)現(xiàn)多模態(tài)醫(yī)學(xué)影像重要信息的高效集成。

1 多模態(tài)醫(yī)學(xué)影像非剛體配準(zhǔn)方法

1.1 醫(yī)學(xué)影像配準(zhǔn)框架

通常醫(yī)學(xué)影像配準(zhǔn)主要包括了空間變換、插值方法、優(yōu)化算法與相似性測度4個部分。其中相似性測度用來描述固定影像與浮動影像之間的配準(zhǔn)程度；優(yōu)化算法根據(jù)相似性測度對變換參數(shù)進(jìn)行優(yōu)化；空間變換部分將優(yōu)化后的變換參數(shù)作用于浮動影像，并通過特定的插值方法形成新的影像[3]。通過不斷迭代直到相似性測度達(dá)到配準(zhǔn)要求。通常由于醫(yī)學(xué)影像的特殊性，醫(yī)學(xué)影像配準(zhǔn)是非常耗時的過程，也是制約醫(yī)學(xué)影像應(yīng)用的重要因素之一。提高醫(yī)學(xué)影像配準(zhǔn)精度始終是醫(yī)學(xué)影像處理的重要內(nèi)容。

以上述思想為基礎(chǔ)，本文提出了粗配準(zhǔn)與精配準(zhǔn)相結(jié)合的配準(zhǔn)方法。其中粗配準(zhǔn)通過仿射變換使固定影像與浮動影像在較少的迭代次數(shù)中達(dá)到大小和位置上的近似匹配；精配準(zhǔn)通過B樣條變換在粗配準(zhǔn)的基礎(chǔ)上對浮動影像進(jìn)行復(fù)雜形變，使2幅影像達(dá)到細(xì)節(jié)上的高精度匹配。粗配準(zhǔn)與精配準(zhǔn)均選擇基于影像自身灰度信息的互信息測度作為相似性測度，因此不需要考慮特征提取。圖1是本文的配準(zhǔn)流程圖。

圖1 配準(zhǔn)流程

1.2 粗配準(zhǔn)方法

1.2.1 多模態(tài)影像配準(zhǔn)下相似性測度方法

相似性測度是醫(yī)學(xué)影像配準(zhǔn)中的關(guān)鍵部分。在多模態(tài)醫(yī)學(xué)影像配準(zhǔn)問題中往往更側(cè)重于考慮影像本身灰度信息的關(guān)聯(lián)，而不考慮影像之間的特征提取問題。互信息是基本信息度量方法，用于度量2個事件集合之間的相關(guān)性，即一個事件集合中包含另一個事件集合信息的程度[4]?；バ畔y度是一種基于影像像素灰度的測度方法，適用于多模態(tài)下影像的配準(zhǔn)，是目前得到廣泛研究和應(yīng)用的方法[5]。

若將2幅影像看作為2個事件集合，影像之間的配準(zhǔn)程度越好，2幅影像的相關(guān)性越大，互信息值越大。假設(shè)固定影像與浮動影像分別為A和B，則2幅圖的互信息I（A，B）計算方法如式為：

式中： H(A)、 H(B)分別為影像A和影像B的信息熵； H(B|A)是 兩幅影像的條件信息熵； H(A,B)為聯(lián)合信息熵，它們定義為：

式中： pA(a)表示影像A中灰度值為a的概率；pB(b)表 示影像B中灰度值為b的概率； pA,B(a,b)表示影像A和影像B在對應(yīng)位置灰度值分別為a和b的概率值。

1.2.2 粗配準(zhǔn)變換方法

當(dāng)固定影像與浮動影像的大小、位置相差較大時，直接對影像進(jìn)行復(fù)雜形變的非剛體配準(zhǔn)的變換操作時不僅配準(zhǔn)效果不理想，而且會消耗大量迭代時間。仿射變換主要作用是對影像進(jìn)行旋轉(zhuǎn)、平移與縮放，它可以分解為線性（矩陣）變換和平移變換[6]。通過仿射變換對影像進(jìn)行粗配準(zhǔn)，可以使影像的大小與位置達(dá)到近似匹配。二維空間里的仿射變換數(shù)學(xué)公式定義為：

1.2.3 粗配準(zhǔn)結(jié)果分析

為了檢驗(yàn)粗配準(zhǔn)的有效性，本文使用的CT、MRI實(shí)驗(yàn)數(shù)據(jù)來自北京陸軍總醫(yī)院放射治療科，多模態(tài)醫(yī)學(xué)影像配準(zhǔn)部分在Visual Studio2015平臺下利用ITK工具包實(shí)現(xiàn)，選取CT影像為固定影像，MRI影像為浮動影像，其中固定影像CT與浮動影像MRI的分辨率均為512×512。如圖2所示，可看出2幅影像灰度分布不同，且腹部區(qū)域的大小和形態(tài)均有較大差別。

圖2 原始影像

粗配準(zhǔn)框架下配準(zhǔn)結(jié)果如圖3所示。從配準(zhǔn)結(jié)果可看出，經(jīng)過剛體配準(zhǔn)后，MRI影像的大小與位置已經(jīng)與CT影像達(dá)到近似配準(zhǔn)，但2幅影像的輪廓與內(nèi)部組織依然存在一定差異，圖3（b）棋盤圖像中右下角部分能看到明顯的錯位。

圖3 粗配準(zhǔn)結(jié)果

圖4是粗配準(zhǔn)過程中的互信息變化曲線。在前10次迭代中，由于2幅影像差異很大，經(jīng)仿射變換可看出互信息明顯上升；11~17次迭代中不斷尋找最佳變換參數(shù)，互信息曲線上下波動；隨后的迭代中，迭代結(jié)果不斷向最優(yōu)配準(zhǔn)逼近，互信息曲線趨于穩(wěn)定。經(jīng)過40次迭代后，CT影像與MRI影像的互信息由0.10上升到0.38，說明2幅影像之間的相關(guān)性增強(qiáng)，與未配準(zhǔn)時相比，重疊程度明顯上升，影像配準(zhǔn)程度變好。

圖4 粗配準(zhǔn)互信息變化曲線

1.3 非剛體配準(zhǔn)框架的方法

1.3.1 精配準(zhǔn)框架中的變換方法

由于式（1）所示的仿射變換無法對影像進(jìn)行復(fù)雜變換，浮動影像的外部輪廓和內(nèi)部組織結(jié)構(gòu)仍然無法實(shí)現(xiàn)精確配準(zhǔn)，這就需要在粗配準(zhǔn)的基礎(chǔ)上對浮動影像進(jìn)行非線性形變的非剛體精配準(zhǔn)。

在非剛體配準(zhǔn)問題中，由于B樣條空間變換[7?8]方法能夠模擬生理結(jié)構(gòu)，已經(jīng)成為非剛體配準(zhǔn)中的重要變換模型。B樣條變換的主要過程包括設(shè)置網(wǎng)格空間、嵌入浮動影像、對網(wǎng)格點(diǎn)變形以及重建變形后的浮動影像。在B樣條變換中，變形場可以描述為

B樣條模型的優(yōu)點(diǎn)在于它對影像操作是局部的，即控制點(diǎn)位置變化只影響控制點(diǎn)鄰域的變化，這樣配準(zhǔn)問題就等價于三次B樣條的局部控制問題，從而使B樣條計算效率更高，精度也能達(dá)到亞像素級[9]。

1.3.2 精配準(zhǔn)結(jié)果分析

精配準(zhǔn)框架下配準(zhǔn)結(jié)果如圖5所示。精配準(zhǔn)后MRI影像的外部輪廓與內(nèi)部組織結(jié)構(gòu)與CT影像一一匹配，在棋盤差異圖下幾乎沒有邊緣錯位，配準(zhǔn)效果較好。圖6是精配準(zhǔn)過程中的互信息變化曲線。

圖5 精配準(zhǔn)結(jié)果

圖6 精配準(zhǔn)互信息變化曲線

由圖6可知，前50次迭代中，由于2幅影像差異較明顯，配準(zhǔn)程度提高較快，互信息曲線明顯上升；在隨后的迭代中，迭代結(jié)果不斷向最優(yōu)配準(zhǔn)逼近，互信息曲線趨于穩(wěn)定。在粗配準(zhǔn)基礎(chǔ)上，精配準(zhǔn)經(jīng)過200次迭代后互信息達(dá)到0.411，說明精配準(zhǔn)進(jìn)一步提高了影像間的相關(guān)性，使配準(zhǔn)結(jié)果與固定影像達(dá)到細(xì)節(jié)上的匹配，為影像融合打下良好基礎(chǔ)。

2 多模態(tài)醫(yī)學(xué)影像多分辨率融合方法

將經(jīng)過預(yù)處理后的CT影像與精配準(zhǔn)后的MRI影像作為影像融合部分的輸入影像。為了利用不同來源醫(yī)學(xué)影像的各自特征，需要對影像進(jìn)行多源信息融合。本文第1章論述的影像配準(zhǔn)為影像融合打下了良好基礎(chǔ)，在此基礎(chǔ)下，采用基于小波變換的融合方法。

2.1 小波變換融合方法

基于小波變換的影像融合方法是在不同頻率分量下進(jìn)行不同規(guī)則的融合處理，可以獲得與人眼視覺更接近的融合效果。基于小波變換的影像融合方法[10]首先對2幅影像分別進(jìn)行小波分解，將分解后的高頻小波系數(shù)與低頻小波系數(shù)分別按一定的規(guī)則進(jìn)行融合，得到融合后的小波系數(shù)，并進(jìn)行小波逆變換得重構(gòu)后的融合影像[11]。設(shè)I1和I2為待融合的2幅影像，則小波融合過程如圖7所示。

圖7 小波融合過程示意

小波融合過程中的重點(diǎn)在于小波分解與分解后融合規(guī)則的選取上[12]。小波分解形式包括小波基與分解層數(shù)的選擇，為了提高小波分解效果，本文采用了db3小波并采用二次分解方法；融合規(guī)則包括點(diǎn)對點(diǎn)的像素融合以及窗口對窗口的區(qū)域融合方法選擇，融合規(guī)則直接影響融合影像的質(zhì)量。

2.2 二維小波分解與重構(gòu)

利用連續(xù)小波變換進(jìn)行計算，計算量非常大，因此要考慮小波變換的快速算法，目前所采用的技術(shù)是沿用Mallat快速算法。二維小波分解Mallat快速算法的分解公式為[13]

由式（2）所示的小波系數(shù)可以對小波變換進(jìn)行重構(gòu)，小波逆變換算法為：

2.3 小波變換融合規(guī)則

CT和MRI影像分解后不同頻帶系數(shù)反映的影像信息不同，低頻系數(shù)融合規(guī)則側(cè)重于反映多模影像的共同輪廓特征，高頻系數(shù)融合規(guī)則側(cè)重于反映盡可能多的細(xì)節(jié)信息。因此對不同頻帶小波系數(shù)需要選取不同的融合規(guī)則。

2.3.1 低頻系數(shù)融合規(guī)則

本文采用加權(quán)平均方法確定低頻子帶系數(shù)的融合規(guī)則，即將影像I1、I2分解后的低頻小波系數(shù)各乘上一個權(quán)重因子，按式（3）求出融合后的系數(shù)，是基于像素的融合規(guī)則。

2.3.2 高頻系數(shù)融合規(guī)則

1）絕對值最大法

絕對值最大法是基于像素的融合規(guī)則。高頻系數(shù)在零值上下波動，絕對值越大說明灰度變化越明顯。因此如式（4）所示選取對應(yīng)高頻系數(shù)中各對應(yīng)點(diǎn)的最大值，可以得到更多細(xì)節(jié)信息。

該方法原理清晰且便于理解，計算效率較高，然而絕對值最大法僅對各系數(shù)單獨(dú)處理，對系數(shù)間的關(guān)聯(lián)性考慮不夠全面，同時融合結(jié)果在對比度上可能有所降低。

2）區(qū)域方差最大法

區(qū)域方差最大法是基于區(qū)域的融合規(guī)則[14]。首先計算高頻系數(shù)中對應(yīng) M×N （本文選擇 3×3）區(qū)域下各中心點(diǎn)的區(qū)域方差，將方差大的系數(shù)點(diǎn)選入新的系數(shù)矩陣F中，最后比較F矩陣各區(qū)域內(nèi)來自不同影像高頻系數(shù)的點(diǎn)的個數(shù)，得出最后的融合結(jié)果。區(qū)域方差最大法具體流程如圖8所示。

圖8 區(qū)域方差最大法

該方法可以反應(yīng)影像灰度分布的離散程度，在灰度反差明顯的部分能很好地保留源影像的特征信息；然而在灰度分布平緩的區(qū)域，區(qū)域方差法難以最大化保留源影像的信息，同時該方法計算量較大。

3)區(qū)域能量融合法

區(qū)域能量融合法是基于像素的融合規(guī)則[15]，但引入了區(qū)域能量的概念。首先計算高頻系數(shù)中對應(yīng)M×N （本文選擇 3×3）區(qū)域下各中心點(diǎn)的區(qū)域能量，并由此計算區(qū)域匹配度，通過區(qū)域匹配度與閾值計算區(qū)域融合法的權(quán)值，最終計算出權(quán)值下各點(diǎn)融合結(jié)果。區(qū)域能量融合法具體流程如圖9所示。

圖9 區(qū)域能量融合法

區(qū)域方差法結(jié)合了前兩種融合規(guī)則的特點(diǎn)，既應(yīng)用了加權(quán)平均的思想，又在確定權(quán)重系數(shù)時充分考慮了像素之間的關(guān)聯(lián)性；然而該方法容易忽略灰度較低的細(xì)節(jié)部分，使最終結(jié)果中出現(xiàn)失真現(xiàn)象，且計算量大，處理過程耗時長。

2.4 影像融合結(jié)果分析

在MATLAB R2014b平臺下實(shí)現(xiàn)影像小波分解的Mallat算法，從緊支性、正交性、對稱性、消失矩4個方面考慮，選擇db3小波作為小波分解與重構(gòu)的小波基。影像的2層小波分解結(jié)果如圖10所示，低頻系數(shù)反映輪廓信息，高頻系數(shù)反映細(xì)節(jié)信息。

圖10 小波分解結(jié)果

對分解后CT影像與MRI影像的小波系數(shù)進(jìn)行融合處理，在融合規(guī)則選取上，對低頻系數(shù)進(jìn)行加權(quán)平均融合，對高頻系數(shù)采用上述絕對值最大、區(qū)域方差最大、區(qū)域能量融合3種方法進(jìn)行融合。圖11是多源影像融合結(jié)果，從3幅影像中可看出，融合后的胸腹腔不僅顯示出了CT影像中的骨骼結(jié)構(gòu)，也顯示出了MRI影像中的血管、軟組織等細(xì)節(jié)，最關(guān)鍵的是腫瘤區(qū)域輪廓明顯、亮度增加，更利于放療醫(yī)生進(jìn)行病灶識別。

圖11 影像融合結(jié)果

為了客觀地比較不同方法的融合結(jié)果，本文利用均值、信息熵、標(biāo)準(zhǔn)差、清晰度、空間頻率[16]5個指標(biāo)進(jìn)行評價，結(jié)果如表1所示?？梢钥闯觯煌u價方法在不同指標(biāo)下所得的評價值有各自的優(yōu)勢，有經(jīng)驗(yàn)的醫(yī)師可以根據(jù)不同的評價指標(biāo)選擇最合適的融合方法。本文所提融合算法已在多組CT與MRI斷層掃描影像序列下進(jìn)行了有效性驗(yàn)證。

表1 不同融合方法的評價指標(biāo)

3 結(jié)論

本文主要討論了多模態(tài)醫(yī)學(xué)影像的非剛體配準(zhǔn)與多分辨率融合方法。

1)在影像配準(zhǔn)部分采用基于仿射變換的粗配準(zhǔn)與基于B樣條變換的精配準(zhǔn)相結(jié)合的非剛體配準(zhǔn)方法，并用互信息測度評價配準(zhǔn)程度，提高了配準(zhǔn)的精度與計算速度。

2)在影像融合部分，提出在不同頻帶使用不同融合規(guī)則的基于小波變換的多分辨率融合方法，對低頻系數(shù)進(jìn)行加權(quán)平均融合，對高頻系數(shù)分別采用絕對值最大法、區(qū)域方差最大法、區(qū)域能量融合法進(jìn)行融合，融合結(jié)果最大限度地保留了多模態(tài)醫(yī)學(xué)影像的有效信息。

本文使用腹部CT影像與MRI影像驗(yàn)證了所提配準(zhǔn)與融合方法的有效性。本文工作是未來實(shí)現(xiàn)多模態(tài)醫(yī)學(xué)影像的三維重構(gòu)、三維配準(zhǔn)融合以及三維靶區(qū)勾畫的重要基礎(chǔ)。上述方法還可以廣泛應(yīng)用于其他多源圖像配準(zhǔn)與融合，例如多源遙感圖像。