多模態醫學影像的非剛體配準與多分辨率融合方法

張欣怡，張富利，王秋生

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

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

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

本文以互信息為相似性測度，搭建仿射變換與B樣條變換相結合的多層次配準框架，實現多模態醫學影像的非剛體配準。同時，對配準后的影像按不同融合規則進行基于小波變換的多分辨率融合，實現多模態醫學影像重要信息的高效集成。

1 多模態醫學影像非剛體配準方法

1.1 醫學影像配準框架

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

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

圖1 配準流程

1.2 粗配準方法

1.2.1 多模態影像配準下相似性測度方法

相似性測度是醫學影像配準中的關鍵部分。在多模態醫學影像配準問題中往往更側重于考慮影像本身灰度信息的關聯，而不考慮影像之間的特征提取問題。互信息是基本信息度量方法，用于度量2個事件集合之間的相關性，即一個事件集合中包含另一個事件集合信息的程度[4]。互信息測度是一種基于影像像素灰度的測度方法，適用于多模態下影像的配準，是目前得到廣泛研究和應用的方法[5]。

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

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

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

1.2.2 粗配準變換方法

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

1.2.3 粗配準結果分析

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

圖2 原始影像

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

圖3 粗配準結果

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

圖4 粗配準互信息變化曲線

1.3 非剛體配準框架的方法

1.3.1 精配準框架中的變換方法

由于式（1）所示的仿射變換無法對影像進行復雜變換，浮動影像的外部輪廓和內部組織結構仍然無法實現精確配準，這就需要在粗配準的基礎上對浮動影像進行非線性形變的非剛體精配準。

在非剛體配準問題中，由于B樣條空間變換[7?8]方法能夠模擬生理結構，已經成為非剛體配準中的重要變換模型。B樣條變換的主要過程包括設置網格空間、嵌入浮動影像、對網格點變形以及重建變形后的浮動影像。在B樣條變換中，變形場可以描述為

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

1.3.2 精配準結果分析

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

圖5 精配準結果

圖6 精配準互信息變化曲線

由圖6可知，前50次迭代中，由于2幅影像差異較明顯，配準程度提高較快，互信息曲線明顯上升；在隨后的迭代中，迭代結果不斷向最優配準逼近，互信息曲線趨于穩定。在粗配準基礎上，精配準經過200次迭代后互信息達到0.411，說明精配準進一步提高了影像間的相關性，使配準結果與固定影像達到細節上的匹配，為影像融合打下良好基礎。

2 多模態醫學影像多分辨率融合方法

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

2.1 小波變換融合方法

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

圖7 小波融合過程示意

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

2.2 二維小波分解與重構

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

由式（2）所示的小波系數可以對小波變換進行重構，小波逆變換算法為：

2.3 小波變換融合規則

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

2.3.1 低頻系數融合規則

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

2.3.2 高頻系數融合規則

1）絕對值最大法

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

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

2）區域方差最大法

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

圖8 區域方差最大法

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

3)區域能量融合法

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

圖9 區域能量融合法

區域方差法結合了前兩種融合規則的特點，既應用了加權平均的思想，又在確定權重系數時充分考慮了像素之間的關聯性；然而該方法容易忽略灰度較低的細節部分，使最終結果中出現失真現象，且計算量大，處理過程耗時長。

2.4 影像融合結果分析

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

圖10 小波分解結果

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

圖11 影像融合結果

為了客觀地比較不同方法的融合結果，本文利用均值、信息熵、標準差、清晰度、空間頻率[16]5個指標進行評價，結果如表1所示。可以看出，不同評價方法在不同指標下所得的評價值有各自的優勢，有經驗的醫師可以根據不同的評價指標選擇最合適的融合方法。本文所提融合算法已在多組CT與MRI斷層掃描影像序列下進行了有效性驗證。

表1 不同融合方法的評價指標

3 結論

本文主要討論了多模態醫學影像的非剛體配準與多分辨率融合方法。

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

2)在影像融合部分，提出在不同頻帶使用不同融合規則的基于小波變換的多分辨率融合方法，對低頻系數進行加權平均融合，對高頻系數分別采用絕對值最大法、區域方差最大法、區域能量融合法進行融合，融合結果最大限度地保留了多模態醫學影像的有效信息。

本文使用腹部CT影像與MRI影像驗證了所提配準與融合方法的有效性。本文工作是未來實現多模態醫學影像的三維重構、三維配準融合以及三維靶區勾畫的重要基礎。上述方法還可以廣泛應用于其他多源圖像配準與融合，例如多源遙感圖像。