宋 梟，朱家明，王 瑩

(揚(yáng)州大學(xué) 信息工程學(xué)院，江蘇 揚(yáng)州 225000)

0 引言

圖像配準(zhǔn)是圖像處理中的一個(gè)重要領(lǐng)域，旨在找尋多幅圖像之間最優(yōu)空間映射，使得2幅或多幅圖像的相關(guān)像素點(diǎn)唯一對(duì)應(yīng)[1]。醫(yī)學(xué)圖像是疾病篩查及診斷的重要依據(jù)，不同模態(tài)的醫(yī)學(xué)圖像反映人體組織或器官的側(cè)重點(diǎn)不同。醫(yī)學(xué)圖像配準(zhǔn)是醫(yī)學(xué)圖像融合的前提，可以有效幫助醫(yī)護(hù)人員進(jìn)行病灶定位、療效評(píng)估、術(shù)前規(guī)劃和術(shù)中導(dǎo)航[2]。

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

1 形變配準(zhǔn)模型

1.1 形變配準(zhǔn)

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

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

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

1.2 級(jí)聯(lián)模型

圖2 級(jí)聯(lián)模型Fig.2 Cascade model

(1)

(2)

(3)

(4)

1.3 生成器結(jié)構(gòu)

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

圖3 形變卷積Fig.3 Deformable convolution

可變形卷積靈活的感受野，增加了網(wǎng)絡(luò)空間形變的適應(yīng)性。可變形卷積對(duì)特征圖的每個(gè)位置學(xué)習(xí)一個(gè)偏移量，用2個(gè)偏移場(chǎng)表示X軸和Y軸的偏移。對(duì)于一個(gè)3×3的卷積，每次輸出都要采樣9個(gè)位置，傳統(tǒng)的卷積輸出為：

(5)

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

可變形卷積在傳統(tǒng)卷積上增加了一個(gè)偏移量Δpn，此時(shí)卷積輸出為：

(6)

本文生成結(jié)構(gòu)以Unet為基礎(chǔ)，在每次下采樣過程中添加偏移卷積層，使空間形變?cè)诿看蜗虏蓸舆^程中逐層編碼，并將形變的特征信息跳躍連接至上采樣階段，使形變特征逐級(jí)還原，其結(jié)構(gòu)如圖4所示。生成結(jié)構(gòu)的具體網(wǎng)絡(luò)參數(shù)如表1所示，表中k3n64s1表示64個(gè)大小為3 pixel×3 pixel,步長(zhǎng)為1的卷積核。

圖4 生成網(wǎng)絡(luò)結(jié)構(gòu)Fig.4 The structure of generative network

表1 生成網(wǎng)絡(luò)參數(shù)Tab.1 The parameters of generative network

1.4 判別器結(jié)構(gòu)

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

圖5 判別網(wǎng)絡(luò)結(jié)構(gòu)Fig.5 The structure of discriminative network

表2 判別器網(wǎng)絡(luò)參數(shù)Tab.2 The parameters of discriminator network

2 構(gòu)造損失函數(shù)

2.1 生成器損失

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

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

(7)

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

(8)

式中，ωi表示加權(quán)因子，ωi≥0且∑ωi=1。

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

(9)

詹森Arimoto散度為：

(10)

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

為了保證生成的形變圖像與參考的固定圖像之間結(jié)構(gòu)信息相似，引入局部梯度項(xiàng)，那么相似性損失為：

(11)

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

(12)

(13)

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

當(dāng)lw與lf完全配準(zhǔn)時(shí)，相似度最大，此時(shí)生成器的參數(shù)θ最優(yōu)化過程為：

(14)

2.2 判別器損失

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

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

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

(15)

對(duì)于浮動(dòng)圖像lf和固定圖像lm，對(duì)抗損失為L(zhǎng)adv(G,Dm,lf,lm)。

生成對(duì)抗網(wǎng)絡(luò)能夠任意變換輸入圖像以匹配目標(biāo)域的分布，循環(huán)一致?lián)p失能夠保持形變的微分同胚性即可以平滑地進(jìn)行正向和反向的形變[16]，表示為：

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

(16)

因此，總的目標(biāo)函數(shù)為：

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

ηLcyc(G,F)，

(17)

式中，η是循環(huán)一致?lián)p失的加權(quán)系數(shù)。

3 驗(yàn)證分析

3.1 實(shí)驗(yàn)環(huán)境

本文實(shí)驗(yàn)以Pytorch為框架，使用Adam優(yōu)化算法，學(xué)習(xí)率設(shè)置為0.000 1，每更新2次判別器后更新生成器。由于硬件環(huán)境限制，僅將訓(xùn)練批次設(shè)置為1，級(jí)聯(lián)次數(shù)設(shè)置為3。用基于ITK開發(fā)的Elastix和同樣基于Unet網(wǎng)絡(luò)構(gòu)造的VoxelMorph模型(簡(jiǎn)稱Vm模型)對(duì)比本文模型，為了提升運(yùn)行速度，將Elastic配準(zhǔn)工具配置在3D-slicer的GPU環(huán)境中。

3.2 實(shí)驗(yàn)數(shù)據(jù)

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

3.3 評(píng)估指標(biāo)

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

(18)

式中，t表示像素點(diǎn)的位置；P表示所有像素點(diǎn)的位置。

峰值信噪比(Peak Signal to Noise Ratio,PSNR)是峰值信號(hào)能量與噪聲平均能量之間的比值，本文用于評(píng)估配準(zhǔn)圖像的質(zhì)量，其值越小失真越大。當(dāng)PSNR值低于20 dB時(shí)，配準(zhǔn)圖像質(zhì)量將嚴(yán)重失真[17]，定義為：

(19)

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

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

(20)

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

3.4 實(shí)驗(yàn)結(jié)果

對(duì)核磁共振T2和Flair圖像進(jìn)行雙向配準(zhǔn)，配準(zhǔn)結(jié)果如圖6所示。圖中第一行配準(zhǔn)結(jié)果的固定圖像是Flair，浮動(dòng)圖像是T2，第二行配準(zhǔn)結(jié)果的固定圖像是T2，浮動(dòng)圖像是Flair。T2圖像能夠清晰反應(yīng)腦部腫瘤及腫瘤團(tuán)塊周圍的水腫帶，而Flair圖像對(duì)腫瘤團(tuán)塊反映不明顯，且2幅圖像在水腫帶邊沿存在明顯的形變。

T2

為了更加直觀地對(duì)比配準(zhǔn)效果，對(duì)配準(zhǔn)結(jié)果進(jìn)行偽彩色處理，效果如圖 7所示。用藍(lán)色、白色線框分別標(biāo)注腦部腫瘤以及側(cè)腦室三角區(qū)。

Elastix(T2-Flair)

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

Flair和T2圖像配準(zhǔn)的指標(biāo)結(jié)果如表3所示。

表3 Flair，T2配準(zhǔn)數(shù)據(jù)Tab.3 Registration data for Flair and T2

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

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

圖8 MR,CT配準(zhǔn)案例Fig.8 Example results for MR and CT registration

MR和CT圖像配準(zhǔn)指標(biāo)結(jié)果如表4所示?？梢钥闯?，本文模型的PSNR、SSIM、TIME均優(yōu)于Elastix和Vm模型，盡管RMSE略高于Elastix，但本文模型的形變圖像與參考的固定圖形具有較高的相似度且失真較小。

表4 MR、 CT配準(zhǔn)數(shù)據(jù)Tab.4 Registration data for MR and CT

4 結(jié)束語

本文提出了一種基于生成對(duì)抗網(wǎng)絡(luò)的可變形醫(yī)學(xué)圖像配準(zhǔn)模型。該模型將可變形卷積和級(jí)聯(lián)結(jié)構(gòu)引入生成模塊，提取形變特征，從而使模型具有形變配準(zhǔn)的能力； 用JAD散度代替原始GAN網(wǎng)絡(luò)中的JS散度，考慮了2個(gè)隨機(jī)概率分布之間的相關(guān)性，提高了配準(zhǔn)的精度；級(jí)聯(lián)了3個(gè)生成器，使形變場(chǎng)逐級(jí)優(yōu)化，進(jìn)一步提高了配準(zhǔn)的精度。基于不同評(píng)估指標(biāo)的實(shí)驗(yàn)表明，本文模型具有較好的配準(zhǔn)精度和形變配準(zhǔn)能力。