基于多通道稀疏編碼的非剛性多模態醫學圖像配準

王麗芳,成 茜,秦品樂,高 媛

0 引言

多模態醫學圖像配準是對從不同成像設備上獲取的圖像進行配準,使兩幅圖像上的同一解剖點或者具有診斷意義的關鍵點達到空間上的一致[1]。當兩幅圖像配準完成后,可以對它們進行比較和分析。多模態醫學圖像配準作為醫學圖像處理中的一個重要組成部分,對手術定位、患者術前和追蹤掃描[2]、多模態圖像融合[3]等具有重要意義。

醫學圖像配準方法主要有剛性配準和非剛性配準兩大類。剛性配準僅適用于不存在形變的配準,目前已經不能滿足人們的需求;當浮動圖像與參考圖像存在較大形變時, 需要進行非剛性配準。非剛性多模態醫學圖像配準主要分為三類:1)基于特征的方法。雖在計算量上要比基于圖像像素的方法小很多,但該方法配準的精度依賴于所提取的特征[4],在多數情況下,特別是在不同模態醫學圖像之間,當圖像特征點不突出、位置復雜的情況下,很難進行準確的特征提取。因此,基于特征的方法對多模態醫學圖像配準并不理想。2)模態減少方法。該方法是將兩種不同的模態轉換成具有直接可比性的單一模態,這種模態可以是待配準圖像中其中一種模態,也可以轉換成一個全新的模態或灰度空間。它的優點是轉換后,多模態配準問題變為單模態的問題,這樣可以使用單模態相似性測度及其相關的優化和變形配準算法。雖然不同的模態可以以不同的方式呈現圖像信息,但由于在某些情況下,一些圖像信息僅在單一一種模態下出現,當轉換成一種全新的模態,此過程中可能造成信息丟失,之后便不能再使用這些圖像信息來進行配準。3)信息理論的方法。互信息(Mutual Information, MI)源于信息論,被廣泛應用于進行多模態醫學圖像配準,是一種基于圖像體素的相似性測度;由于它只定義了對應像素的灰度關系,未考慮像素之間的空間依賴性,通常對含有灰度偏移場的圖像魯棒性差,且在優化過程中易陷入局部極值,最終不能產生滿意的配準結果。為了解決互信息等傳統相似性測度對于灰度不均勻性圖像魯棒性較差的問題,Myronenko等[5]提出在殘差復雜性(Residual Complexity, RC)的分析上解決灰度級偏移場,盧振泰等[6]提出了局部方差和殘差復雜性(Local Variance and Residual Complexity, LVRC)的醫學圖像配準方法,Ghaffari等[7-8]提出了稀疏誘導的相似性測度(Sparse-Induced Similarity Measure, SISM)和Rank Induced相似性測度(Rank Induced Similarity Measure, RISM)。這些方法可以產生較準確的配準結果,但適用范圍較小,多只適用于單模態醫學圖像配準。Cheng等[9]提出使用深度相似性學習方法訓練一個二元分類器來學習兩個圖像塊的對應關系；Simonovsky等[10]提出使用有卷積神經網絡(Convolutional Neural Network, CNN)的相似性測度,來估計來自不同模態兩個圖像塊之間的相似性；但文獻[9-10]中關于深度學習的方法訓練成本高,需要大量的訓練數據才能達到滿意的效果,而要獲得足夠的未配準的醫學圖像數據集并為這些圖像進行人工標注需要耗費大量的工作量[11]。

為了使多模態醫學圖像配準對灰度偏移場具有較好的魯棒性,本文提出了基于多通道稀疏編碼的非剛性多模態醫學圖像配準方法。該方法將多通道與基于學習的稀疏編碼的L1范數相似性測度(L1norm Similarily Measure， L1SM)相結合,使多模態配準問題變為一個多通道配準問題,劃分通道后每個通道中圖像的模態相同。由于最初缺乏在這個多通道配準框架中可以使用的對應模式，為了解決這個問題,本文使用文獻[12]中圖像合成方法創建的代理圖像來彌補由于轉換而可能丟失的信息,使多通道配準框架在浮動圖像和參考圖像中都產生相同模態的圖像,從而配準來自不同模態的圖像。將多通道配準框架與只適用于單模態配準的基于稀疏編碼的相似性測度相結合,不僅克服了醫學圖像中存在灰度不均勻性造成的灰度偏移場,而且實現了基于稀疏表示的多模態醫學圖像配準。

1 基于稀疏編碼的L1SM

假設兩個圖像之間的灰度關系如下:

R=F(T)+S

(1)

其中:R表示參考圖像;F表示浮動圖像;S表示灰度不均勻性造成的灰度偏移場;T是空間變換。

在分析模型中,用式(2)來定義一幅圖像的稀疏表示:

α=ΩY

(2)

其中:Y是稀疏系數矩陣;Ω是分析字典;α是圖像塊的稀疏表示。

根據式(2),得到參考圖像和浮動圖像的稀疏表示分別為:

αR=ΩR

(3)

αF=ΩF

(4)

分析式(1)和它的稀疏表示之間的關系,在式(1)兩邊都乘以Ω,在稀疏域中得到如下的線性關系:

Ω×R=Ω×(F(T)+S)?αR=αF+αS

(5)

其中:αS=ΩS,是灰度校正場的稀疏表示。這里將L1范數作為相似性測度。

通過最大后驗概率(Maximum A Posteriori, MAP)的方法來近似估計T,最大化以下的可能性:

P(T|αR,αF)∝P(αR,αF|T)P(T)=

PαS(αR-αF)P(T)

(6)

T=arg max(ln(P(αR,αF|T)P(T)))=

arg min(‖αS‖1)=arg min(‖αR-αF‖1)=

arg min(‖Ω(R-F)‖1)

因此有如下的L1范數相似性測度:

L1SM(R,F)=‖Ω(R-F)‖1

(7)

其中:‖·‖1表示L1范數。

2 多通道配準框架

多通道配準是將不同模態的圖像劃分為多個通道,每個模態在一個單獨的通道下運行[13-14]。在單通道配準中,假設R(x)表示參考圖像,F(x)表示浮動圖像,T為空間變換。通過T來最小化能量函數H:

H=Csim(F°T,R)

(8)

其中：Csim是相似性測度。

在多通道配準中,將單個通道的能量函數的加權和作為最終的相似性測度函數HMC。HMC的計算公式為:

(9)

其中：M表示模態的總數；(F1,F2,…,FM)為浮動圖像集；(R1,R2,…,RM)為參考圖像集；wm為每個模態m對最終能量函數貢獻的權重；Csim-m是與每個模態相關的通道所使用的相似性測度。

假設浮動圖像為第一模態F→F1,參考圖像為第二模態R→R2。在多通道配準框架中,需使用圖像合成方法,能量函數為:

(10)

圖1 多通道配準框架

3 本文配準方法

3.1 基于多通道稀疏編碼的相似性測度

由于腦影像中都存在緩慢變化的灰度偏移場,而基于稀疏編碼的L1范數相似性測度(L1SM)對含有灰度不均勻性圖像魯棒性較好,因此分別在多通道配準框架的每個通道中使用L1SM作為相似性測度,總的基于多通道稀疏編碼的相似性測度(Multi-Channel L1SM, MCL1SM)函數為:

其中:w1、w2分別為T1加權、T2加權圖像通道對HMCL1SM的貢獻權重;Ω1、Ω2分別為兩個通道中的分析字典;α為對應圖像塊的稀疏表示;‖·‖1表示L1范數。

3.1.1 圖像合成

3.1.2 圖像正則化

3.1.3 劃分通道和圖像塊

3.1.4 稀疏編碼和字典更新

使用K奇異值分解(K-means based Singular Value Decomposition,K-SVD)算法分別訓練兩個通道中兩幅圖像的圖像塊得到分析字典Ω1、Ω2,計算每個圖像塊的稀疏系數Yi。

本過程主要分兩個階段：

1)稀疏編碼的求解過程階段。需要計算圖像塊

在字典Ω1、Ω2上的稀疏系數Yi。

(11)

s.t. ‖α-ΩY‖≤ε

其中: ‖Y‖0表示稀疏系數向量Y中非零元素的個數；ε表示允許偏差的精度。

(12)

(13)

(14)

最后，K-SVD是一種迭代算法,通過交替地執行稀疏編碼和字典更新這兩步驟得到分析字典Ω1、Ω2。

3.2 算法流程

本文提出的基于多通道稀疏編碼的非剛性多模態醫學圖像配準方法的算法流程如圖2所示。

圖2 本文算法流程

1)讀入兩幅分辨率相同的待配準的醫學圖像,分別記為參考圖像R和浮動圖像F。

2)初始化多層P樣條自由形式變換(Free-Form Deformation, FFD)模型參數和優化模塊中優化的迭代次數和迭代的最小誤差,參考圖像和浮動圖像之間的配準誤差為10-4,梯度下降法的最大迭代次數為30。

4)使用K-SVD算法分別訓練每個通道中的圖像塊,得到分析字典Ω1、Ω2和每個圖像塊的稀疏表示。在每個通道分別使用L1SM,計算每個通道L1SM的初始值:

并計算兩個通道相似性測度的能量函數的加權和:

5)將多模態相似性測度函數

作為優化模塊目標函數的第一部分HMCL1SM,另一部分加入相關的平滑變換,使待配準的兩幅圖像之間保持權衡,并且變換光滑。

C(R,F)=HMCL1SM+λCsmooth

(15)

其中:A是區域的面積;λ是權重參數。

6)本文利用梯度下降法迭代地優化目標函數,進而更新變形場,本文HMCL1SM的梯度表達式計算過程如下:

(16)

7)使用多層P樣條FFD[15]對浮動圖像的相應區域進行更新來擬合圖像形變。變形過程是通過重復迭代,直到得到最優變形場。

在二維圖像中,控制點位置為φij,間距為Δ的δx×δy網格,其多層P樣條函數的自由形式變換可描述為:

(17)

總的懲罰P為:

P=λ0P1?Ik0+λ1Ik1?P0+λ2P3?Ik2+

λ3Ik3?P2

其中:P0、P1、P2、P3分別代表邊緣基B0、B1、B2、B3的懲罰；λ0、λ1、λ2、λ3為平滑系數；Ik0、Ik1、Ik2、Ik3為單位矩陣。

8)輸出配準后圖像。

4 實驗結果與分析

本文的實驗環境為: CPU為Intel CoreTM i5-7500,內存8 GB,操作系統為Windows 10,編程平臺為Matlab2013a。

為了驗證MCL1SM在非剛性多模態醫學圖像配準中具有對灰度偏移場魯棒性好的特點,本文進行了大量的仿真實驗,并與局部互信息(Regional Mutual Information, RMI)[6]、LVRC[6]、SISM[7]和RISM[8]在多通道配準框架下的多通道局部方差和殘差復雜性(Multi-Channel Local Variance and Residual Complexity, MCLVRC)、多通道稀疏誘導的相似性測度(Multi-Channel Sparse-Induced Similarity Measure, MCSISM)、多通道Rank Induced相似性測度(Multi-Channel Rank Induced Similarity Measure, MCRISM)進行了比較。

實驗首先選用BrainWeb database[17]的腦部MR數據,這里對其進行人工形變并施加不同程度的灰度不均勻場,圖3顯示了待配準的切片厚度為1 mm的三組圖像,均是加入40%或20%灰度偏移場的非均勻腦部圖像。分別使用局部互信息(RMI)、MCLVRC、MCSISM、MCRISM與本文MCL1SM共五種測度進行配準對比實驗,配準結果分別如圖4所示。

其次,選用哈佛大學醫學院[18]的兩組腦部數據,對其進行人工形變并施加不同程度的灰度不均勻場,圖5為待配準的兩組圖像。分別使用局部互信息(RMI)、MCLVRC、MCSISM、MCRISM與本文MCL1SM共5種測度進行配準對比實驗,配準結果分別如圖6所示。

圖4、6對應五組待配準圖像分別使用局部互信息(RMI)、MCLVRC、MCSISM、MCRISM和本文MCL1SM配準后所得圖像和經差值運算得到的差值圖像。

由圖4、6中五組圖像在不同相似性測度下配準后圖像差不難看出,使用本文相似性測度得到的配準圖像與參考圖像的圖像差要比RMI、MCLVRC、MCSISM、MCRISM得到的配準圖像與參考圖像的圖像差差異性要小,配準后兩幅圖像差異性越小,配準效果越好。

為客觀評價配準效果,本文主要從配準精度和耗時同已有相關文獻方法進行比較,并采用均方根誤差(Root Mean Square Error, RMSE)來衡量圖像配準精度。RMSE值越小,配準效果越好。

為定量比較, 圖3中三組圖像在各相似性測度下配準后圖像與參考圖像的RMSE值與配準運行時間如表1所示。

圖3 待配準的三組圖像(Brainweb database)

圖4 對T1/T2加權圖像用不同相似性測度得到的配準結果及差值圖像

圖5 待配準的兩組圖像(哈佛大學醫學院)

圖6 對圖5用不同方法得到的配準結果及差值圖像

組別配準方案RMSE配準時間/s第一組(矢狀面)第二組(橫斷面)第三組(冠狀面)RMI1.2148493.2MCLVRC1.0917480.3MCSISM1.1941593.1MCRISM0.9682658.4本文算法0.7900696.0RMI1.0810423.2MCLVRC0.9749455.3MCSISM1.0149573.3MCRISM0.8678637.0本文算法0.7242660.6RMI1.1508471.8MCLVRC1.0736455.7MCSISM1.0532692.3MCRISM0.9539679.4本文算法0.7737709.6

圖5中兩組圖像在各相似性測度下配準后所得圖像與參考圖像的RMSE值與配準運行時間如表2所示。

表1～2給出5組圖像在局部互信息(RMI)、MCLVRC、MCSISM、MCRISM與本文MCL1SM下配準結果的定量指標(數據是5次實驗數據的平均值)。依據對實驗結果中配準結果、差值圖像及表1～2的分析,得出基于MCL1SM相對于RMI、MCLVRC、MCSISM、MCRISM在多模態非剛性配準中對含有灰度偏移場的圖像有更好的配準表現;RMSE平均分別下降了30.86%、22.24%、26.84%和16.49%。

出現上述現象主要因為本文方法是基于圖像塊且字典是通過K-SVD算法對圖像塊訓練得到的,并在此基礎上對來自不同模態的圖像進行合成和正則化,與固定字典和基于信息理論的相似性測度相比,基于字典學習的L1范數相似性測度在多通道配準框架下能夠對存在灰度不均勻場的多模態醫學圖像有很好的魯棒性。但從配準時間上來看, MCSISM、MCRISM、MCL1SM的效率要略低于RMI、 MCLVRC。

表2 CT/MRI、PET/MRI在各相似性測度下配準均方根誤差和時間比較

5 結語

本文將基于學習的稀疏編碼和多通道配準框架相結合用于多模態非剛性醫學圖像配準,解決了多模態配準中傳統相似性測度對灰度不均勻性魯棒性較差的問題。仿真實驗結果表明,本文算法能夠較準確地對多模態醫學圖像進行配準并提升了圖像配準的魯棒性,但由于需要對圖像進行合成、正則化,并對圖像塊進行訓練,該算法相對于傳統的方法時間效率較低,因此如何提升該算法的配準效率是接下來進一步研究的重點。

參考文獻(References)

[1] FRANCISCO P M, MANUEL J O, TAVARES R S. Medical image registration: a review[J]. Computer Methods in Biomechanics and Biomedical Engineering, 2014, 17(2): 73-93.

[2] KWON D, ZENG K, BILELLO M, et al. Estimating patient specific templates for pre-operative and follow-up brain tumor registration[C]// MICCAI 2015: Proceedings of the 18th International Conference on Medical Image Computing and Computer-Assisted Intervention, LNCS 9350. Berlin: Springer, 2015: 222-229.

[3] BHATTACHARYA M, ARPITA D. Multimodality medical image registration and fusion techniques using mutual information and genetic algorithm-based approaches[J]. Software Tools and Algorithms for Biological Systems, 2011, 696(3): 441-449.

[4] 王觀英, 許新征, 丁世飛. 基于3D-PCNN和互信息的3D-3D醫學圖像配準方法[J]. 計算機應用, 2017, 37(增刊1): 215-219.(WANG G Y, XU X Z, DING S F. 3D-3D medical image registration method based on 3D-PCNN and mutual information [J]. Journal of Computer Applications, 2017, 37(S1): 215-219.)

[5] MYRONENKO A, SONG X. Intensity-based image registration by minimization of residual complexity[J]. IEEE Transactions on Medical Imaging, 2010, 29(11): 1882-1891.

[6] 盧振泰, 張娟, 馮前進, 等.基于局部方差與殘差復雜性的醫學圖像配準[J]. 計算機學報, 2015, 38(12): 2400-2411.(LU Z T, ZHANG J, FENG Q J, et al. Medical image registration based on local variance and residual complexity [J]. Chinese Journal of Computers, 2015, 38(12): 2400-2411.)

[7] GHAFFARI A, FATEMIZADEH E. Sparse-induced similarity measure: mono-modal image registration via sparse-induced similarity measure[J]. IET Image Processing, 2014, 8(12): 728-741.

[8] GHAFFARI A, FATEMIZADEH E. RISM: single-modal image registration via rank-induced similarity measure[J]. IEEE Transactions on Image Processing, 2015, 24(12): 5567-5580.

[9] CHENG X, ZHANG L, ZHENG Y. Deep similarity learning for multimodal medical images[J]. Computer Methods in Biomechanics and Biomedical Engineering: Imaging & Visualization, 2016, 6(4): 1-5.

[10] SIMONOVSKY M, GUTIERREZ-BECKER B, MATEUS D, et al. A deep metric for multimodal registration [C]// MICCAI 2016: Proceedings of the 19th International Conference on Medical Image Computing and Computer-Assisted Intervention, LNCS 9902. Berlin: Springer, 2016: 10-18.

[11] LITJENS G, KOOI T, BEJNORDI B E, et al. A survey on deep learning in medical image analysis[J]. Medical Image Analysis, 2017, 42(9): 60-88.

[12] JOG A, ROY S, CARASS A, et al. Magnetic resonance image synthesis through patch regression [C]// Proceedings of the 2013 IEEE 10th International Symposium on Biomedical Imaging. Piscataway, NJ: IEEE, 2013: 350-353.

[13] CHEN M, JOG A, CARASS A, et al. Using image synthesis for multi-channel registration of different image modalities[C]// Proceedings of SPIE 9413. Bellingham, WA: SPIE Press, 2015: 21-26.

[14] FORSBERG D, RATHI Y, BOUIX S, et al. Improving registration using multi-channel diffeomorphic demons combined with certainty maps[C]// MBIA 2011: Proceedings of the First International Conference on Multimodal Brain Image Analysis. Berlin: Springer, 2011: 19-26.

[15] 王麗芳,成茜,秦品樂,高媛.基于多層P樣條和稀疏編碼的非剛性醫學圖像配準方法[J/OL].計算機應用研究:1-4.[2017- 03- 27].http://kns.cnki.net/kcms/detail/51.1196.TP.2017072 7.2121.124.html.(WANG L F, CHENG X, QIN P L, et al. Non-rigid medical image registration method based on multilayer P-spline and sparse coding[J]. Applications Research of Computers:1-4.[2017- 03- 27].http://kns.cnki.net/kcms/detail/51.1196.TP.20170727.2121.124.html.)

[16] PRADHAN S, PATRA D. RMI based non-rigid image registration using BF-QPSO optimization and P-spline[J]. International Journal of Electronics and Communications, 2015, 69(3): 609-621.

[17] COCOSCO C A, KOLLOKIAN V. BrainWeb: simulated brain database [EB/OL]. [2018- 01- 23]. http: //brainweb.bic.mni.mcgill.ca/brainweb/.

[18] JOHNSON K A, BECKER J A. The whole brain Atlas[DB/OL]. [2017- 11- 15]. http: //www.med.harvard.edu/aanlib/home.html.

This work is partially supported by the Natural Science Foundation of Shanxi Province (2015011045).