基于生成對抗網絡的多模態可變形醫學圖像配準

宋 梟，朱家明，王 瑩

(揚州大學 信息工程學院，江蘇 揚州 225000)

0 引言

圖像配準是圖像處理中的一個重要領域，旨在找尋多幅圖像之間最優空間映射，使得2幅或多幅圖像的相關像素點唯一對應[1]。醫學圖像是疾病篩查及診斷的重要依據，不同模態的醫學圖像反映人體組織或器官的側重點不同。醫學圖像配準是醫學圖像融合的前提，可以有效幫助醫護人員進行病灶定位、療效評估、術前規劃和術中導航[2]。

傳統的基于互信息和結構特征的醫學圖像配準目的在于通過長時間的迭代最大化圖像強度之間的統計相關性以及最大化基于結構表示的結構度量，不斷減小圖像映射之間的差異[3]。調整線性變換的相關參數很容易減小線性差異對配準結果的影響，然而非線性的局部形變難以配準，因此需要通過網絡充分學習圖像之間的映射以解決該問題[4]。Guha等人[5]提出了Voxelmorph配準模型，利用Unet網絡學習到的參數直接計算函數獲取配準場，大幅縮短了配準時間； Zhao等人[6]提出了級聯遞歸配準模型，將每個子網絡的輸出作為后一個子網絡的輸入，漸進生成形變場，有效降低了網絡學習的難度；Mahapatra等人[7]提出了基于生成對抗網絡(Generative Adversarial Networks,GANs)的形變配準模型，將醫學圖像的灰度作為概率測度，用概率來表示配準效果，避免了耗時的迭代，直接生成形變場與配準圖像。為了提高GANs形變配準能力，本文提出了一種結合GANs、級聯網絡和Unet的醫學圖像配準模型，考慮到標準卷積不能提取形變特征，將形變卷積嵌入到Unet的下采樣過程中，并在上采樣過程中融入形變特征。

1 形變配準模型

1.1 形變配準

由于計算機體層成像(Computed Tomography,CT)、核磁共振(Magnetic Resonance Image,MRI)的層間距、層厚不一致，導致成像結果不完全對齊，表現在人體組織或器官的非線性形變，如圖1所示。第一行中腹部分割結果展示了圖像不對齊導致的器官非線性形變。即使在同一種成像技術下，不同成像參數也會造成一定的非線性形變；第二行中紅色和藍色標注的MRI成像過程中不同弛豫時間腦部組織或病灶的非線性形變。

圖1 非線性局部形變Fig.1 Nonlinear local deformation

假設浮動圖像lm、固定圖像lf定義在二維空間域Ω?R中，以θ為參數的形變場是一個映射φθ:Ω→Ω。對于形變配準而言，旨在構建一個形變場預測函數Fθ(lm,lf)=φθ，浮動圖像lm在形變場預測函數的作用下得到與固定圖像空間對齊的形變圖像lw=lm°φθ，其中°表示形變場作用在浮動圖像上。

1.2 級聯模型

圖2 級聯模型Fig.2 Cascade model

(1)

(2)

(3)

(4)

1.3 生成器結構

在卷積神經網絡(Convolutional Neural Networks,CNNs)中同一層的激活單元具有相同的感受野，但不同位置對應著不同尺度或形變的物體，因此自適應調整感受野是精確提取形變特征的關鍵[8]。Dai等人[9]提出的可形變卷積網絡(Deformable Convolution Networks,DCNs)有效提高了CNN的形變建模能力。DCNs在特征圖中增加了額外的偏移量，從目標任務中學習到的特征矩陣相較于標準卷積增加了一個偏移矩陣[10]。Wang等人[11]將DCNs引入圖像配準，將2N個偏移輸出修改為2個偏移輸出，將變形直接應用于輸入圖像，其結構如圖3所示。

圖3 形變卷積Fig.3 Deformable convolution

可變形卷積靈活的感受野，增加了網絡空間形變的適應性。可變形卷積對特征圖的每個位置學習一個偏移量，用2個偏移場表示X軸和Y軸的偏移。對于一個3×3的卷積，每次輸出都要采樣9個位置，傳統的卷積輸出為：

(5)

式中，R={(-1,-1),(-1,0),…,(0,1),(1,1)}表示卷積核中的9個位置；w(pn)表示加權采樣；x(p0) 表示每次作用在圖像上卷積核的中心位置。

可變形卷積在傳統卷積上增加了一個偏移量Δpn，此時卷積輸出為：

(6)

本文生成結構以Unet為基礎，在每次下采樣過程中添加偏移卷積層，使空間形變在每次下采樣過程中逐層編碼，并將形變的特征信息跳躍連接至上采樣階段，使形變特征逐級還原，其結構如圖4所示。生成結構的具體網絡參數如表1所示，表中k3n64s1表示64個大小為3 pixel×3 pixel,步長為1的卷積核。

圖4 生成網絡結構Fig.4 The structure of generative network

表1 生成網絡參數Tab.1 The parameters of generative network

1.4 判別器結構

判別器的目的在于鑒定形變圖像與參考的固定圖像之間的配準程度。判別器通過閾值函數將輸出映射到0和1之間，當輸出越接近1時代表配準程度越高，越接近0時配準程度越低[12]。本文判別器由8個卷積層和2個全連接層組成，結構如圖5所示。圖5中8個卷積層的卷積核大小均為3 pixel×3 pixel，步長為1或2，當步長為1時圖像尺寸不變，當步長為2時圖像尺寸減小為原來的一半，卷積核數量由64個成倍增加到512個，最終提取到512個16 pixel×16 pixel的特征圖。接著第一個全連接層用1 024個神經元將二維特征圖轉換成一維數組，第二個全連接層用一個神經元經閾值函數完成預測輸出。判別網絡的具體參數如表2所示，表中k3n64s1表示64個大小為3 pixel×3 pixel,步長為1的卷積核。

圖5 判別網絡結構Fig.5 The structure of discriminative network

表2 判別器網絡參數Tab.2 The parameters of discriminator network

2 構造損失函數

2.1 生成器損失

大部分學者提出的基于信息論的配準方法諸如NMI、CCRE等是基于香農熵構造的，而Antolin等[13]指出香農熵的可加性并未考慮2個獨立隨機變量之間的相互作用。為此，李碧草等人[14]根據具有偽性的Arimoto熵提出了可用來量化隨機變量概率分布之間距離的詹森Arimoto散度(Jensen Arimoto Divergence，JAD)。

假設一個概率分布為p=(p1,p2,…,pk)的隨機變量x(x1,x2,…,xk)，則x的Arimoto熵為：

(7)

式中，α是衡量偽可加程度的參數，α>0 且α≠1，當α趨于1時Arimoto熵等于香農熵；Aα(·) 表示Arimoto熵。概率分布P的JAD定義為：

(8)

式中，ωi表示加權因子，ωi≥0且∑ωi=1。

GANs類似極大似然估計，通過模擬數據概率分布，使得概率分布與觀測數據的概率統計分布一致或者盡可能接近[15]。對于生成的形變圖像lw=lf°φθ和固定圖像lf，將其灰度值看作隨機變量，則它們的聯合熵為：

(9)

詹森Arimoto散度為：

(10)

式中，α>0 且α≠1；fi，wj表示固定圖像和生成形變圖像的灰度級；p(fi)，p(wj)，p(wj|fi)分別表示固定圖像和形變圖像不同灰度級的概率分布以及2幅圖像的條件概率分布。

為了保證生成的形變圖像與參考的固定圖像之間結構信息相似，引入局部梯度項，那么相似性損失為：

(11)

式中，‖表示L2距離，此時生成器損失函數表示為：

(12)

(13)

式中，x，y表示圖像X、Y軸方向上的尺寸。

當lw與lf完全配準時，相似度最大，此時生成器的參數θ最優化過程為：

(14)

2.2 判別器損失

對于有限訓練樣本xi∈X,yi∈Y，GANs的目的在于訓練2個生成器G:X→Y,F:Y→X，和2個判別器DX,DY，DX用于區分樣本xi和生成數據F(y)，DY用于區分樣本yj和生成數據G(x) 。對抗損失為：

Ladv(G,DY,X,Y)=Ey∈p(y)[lgDY(y)]+

Ex∈p(x)[lg(1-DY(G(x)))]。

(15)

對于浮動圖像lf和固定圖像lm，對抗損失為Ladv(G,Dm,lf,lm)。

生成對抗網絡能夠任意變換輸入圖像以匹配目標域的分布，循環一致損失能夠保持形變的微分同胚性即可以平滑地進行正向和反向的形變[16]，表示為：

Lcyc(G,F)=Ex‖F(G(x))-x‖1+Ey‖G(F(y))-y‖1。

(16)

因此，總的目標函數為：

L(G,F,Dlm,Dlf)=Ladv(G,Dlf,lm,lf)+Ladv(F,Dlm,lf,lm)+

ηLcyc(G,F)，

(17)

式中，η是循環一致損失的加權系數。

3 驗證分析

3.1 實驗環境

本文實驗以Pytorch為框架，使用Adam優化算法，學習率設置為0.000 1，每更新2次判別器后更新生成器。由于硬件環境限制，僅將訓練批次設置為1，級聯次數設置為3。用基于ITK開發的Elastix和同樣基于Unet網絡構造的VoxelMorph模型(簡稱Vm模型)對比本文模型，為了提升運行速度，將Elastic配準工具配置在3D-slicer的GPU環境中。

3.2 實驗數據

本文使用MICCAI BraTS 2018分割挑戰賽和Learn2Reg 2021配準挑戰賽中的數據作為實驗數據，其中BraTS 2018收集了大量腦部腫瘤MRI數據，訓練集共有285個病例，每個病例包括4種MRI序列：T1、T1ce、T2和Flair，尺寸為240 pixel×240 pixel×155。Learn2Reg 2021收錄了來自TCIA、BCV、CHAOS的122對經過裁剪、重采樣等預處理的腹部MRI和CT成像。

3.3 評估指標

均方差誤差(Root Mean Square Error,RMSE)是反映估計量和被估計量之間的差異程度的一種度量，本文用來度量配準前后圖像灰度概率分布之間的差異[18]。其值越小，模型擬合程度越高，生成的形變圖像越接近參考的固定圖像 ，定義為：

(18)

式中，t表示像素點的位置；P表示所有像素點的位置。

峰值信噪比(Peak Signal to Noise Ratio,PSNR)是峰值信號能量與噪聲平均能量之間的比值，本文用于評估配準圖像的質量，其值越小失真越大。當PSNR值低于20 dB時，配準圖像質量將嚴重失真[17]，定義為：

(19)

式中，MAX為圖片最大像素值。

結構相似性(Structural Similarity,SSIM)，從圖像組成的角度將結構信息定義為獨立于亮度、對比度的反映場景中物體結構的屬性，并將失真建模為亮度、對比度和結構3個不同因素的組合[19]，其值越大結構越相似，定義為：

(20)

式中，μlw為形變圖像的均值；μlf為固定圖像的均值；σlw為形變圖像的方差；σlf為固定圖像的方差；σlwlf為形變圖像和固定圖像的協方差；c1,c2為常數。

3.4 實驗結果

對核磁共振T2和Flair圖像進行雙向配準，配準結果如圖6所示。圖中第一行配準結果的固定圖像是Flair，浮動圖像是T2，第二行配準結果的固定圖像是T2，浮動圖像是Flair。T2圖像能夠清晰反應腦部腫瘤及腫瘤團塊周圍的水腫帶，而Flair圖像對腫瘤團塊反映不明顯，且2幅圖像在水腫帶邊沿存在明顯的形變。

T2

為了更加直觀地對比配準效果，對配準結果進行偽彩色處理，效果如圖 7所示。用藍色、白色線框分別標注腦部腫瘤以及側腦室三角區。

Elastix(T2-Flair)

由圖7可以看出，本文模型在腦腫瘤配準結果中水腫帶邊沿紅綠偽影較少，且在側腦三角室部分重合度高于其他2種模型，配準效果明顯優于其他2種模型。

Flair和T2圖像配準的指標結果如表3所示。

表3 Flair，T2配準數據Tab.3 Registration data for Flair and T2

從表中可以看出，本文模型雙向配準的PSNR、SSIM均高于Elastix和Vm模型，RMSE、TIME均低于Elastix和Vm模型，說明本文的配準精度和配準時間優于其他2種模型。此外，Elastix模型雙向配準結果差異率的絕對值分別為3.82%,12.50%,6.78%和1.65%，Vm模型雙向配準的差異率的絕對值分別為1.39%,10.42%,4.33%和22.22%，本文模型雙向配準的差異率的絕對值分別為1.43%,2.50%,4.63%和12.5%。上述數據中本文模型的雙向配準差異率與Vm模型相當，但前三者明顯低于Elastix模型，說明基于本文模型具有較好的雙向配準能力。

為了對比不同成像技術下的形變配準效果，本文對腹部的CT和MR圖像進行配準，配準結果如圖8所示。用紅色箭頭標注局部形變配準欠缺的地方，可以看出Elastix模型中箭頭數量較多，Vm模型其次，本文模型最少，說明本文模型在CT和MR形變配準中具有明顯的優勢。

圖8 MR,CT配準案例Fig.8 Example results for MR and CT registration

MR和CT圖像配準指標結果如表4所示。可以看出，本文模型的PSNR、SSIM、TIME均優于Elastix和Vm模型，盡管RMSE略高于Elastix，但本文模型的形變圖像與參考的固定圖形具有較高的相似度且失真較小。

表4 MR、 CT配準數據Tab.4 Registration data for MR and CT

4 結束語

本文提出了一種基于生成對抗網絡的可變形醫學圖像配準模型。該模型將可變形卷積和級聯結構引入生成模塊，提取形變特征，從而使模型具有形變配準的能力； 用JAD散度代替原始GAN網絡中的JS散度，考慮了2個隨機概率分布之間的相關性，提高了配準的精度；級聯了3個生成器，使形變場逐級優化，進一步提高了配準的精度。基于不同評估指標的實驗表明，本文模型具有較好的配準精度和形變配準能力。