宋 梟，朱家明，徐婷宜

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

0 引言

圖像配準(zhǔn)是圖像研究領(lǐng)域的一個(gè)典型問題和技術(shù)難點(diǎn)，其目的在于比較或融合針對(duì)同一對(duì)象在不同條件下獲取的圖像。隨著醫(yī)學(xué)成像設(shè)備的不斷進(jìn)步，可以采集更加精準(zhǔn)的器官解剖圖像，但僅憑醫(yī)務(wù)工作者的空間想象和主觀經(jīng)驗(yàn)觀察病情的發(fā)展十分困難。精確快速的醫(yī)學(xué)圖像配準(zhǔn)是醫(yī)學(xué)深入分析的基礎(chǔ)，是分析病灶、觀察器官變化、引導(dǎo)臨床手術(shù)的關(guān)鍵。為此，廣大學(xué)者提出了許多醫(yī)學(xué)圖像配準(zhǔn)方法，從傳統(tǒng)的基于特征和灰度的方法到深度學(xué)習(xí)方法，圖像配準(zhǔn)的性能得到了質(zhì)的飛躍。

傳統(tǒng)的配準(zhǔn)大多數(shù)在變分框架中實(shí)現(xiàn)[1]，即迭代最小化變形空間中的能量問題，計(jì)算量大且受局部極值影響，魯棒性差。深度學(xué)習(xí)方法直接通過深層網(wǎng)絡(luò)提取不同層次的灰度、結(jié)構(gòu)特征，最大化配準(zhǔn)圖像之間的相似度。Maes等人[2]提出的最大化互信息準(zhǔn)則的配準(zhǔn)方法，比較灰度概率分布來度量圖像配準(zhǔn)程度，利用powell算法進(jìn)行參數(shù)優(yōu)化，當(dāng)圖像對(duì)齊時(shí)互信息最大，該方法易收斂至局部極值，且采樣及插值算法對(duì)配準(zhǔn)結(jié)果影響較大；Panigrahi等人[3]提出的基于仿射變換的剛性配準(zhǔn)，通過優(yōu)化多個(gè)空間形變多項(xiàng)式的參數(shù)，消除了因成像角度、時(shí)間造成的不規(guī)則失真，但對(duì)于局部形變，不能達(dá)到良好的配準(zhǔn)效果；Fan等人[4]提出的基于生成對(duì)抗網(wǎng)(Generative Adversarial Network，GAN)的圖像配準(zhǔn)模型，通過最小化對(duì)抗損失，直接輸出配準(zhǔn)圖像，并未對(duì)配準(zhǔn)圖像進(jìn)行空間變換以及相似性度量，僅用概率表示配準(zhǔn)效果，難以達(dá)到精準(zhǔn)的配準(zhǔn)效果。為了解決上述問題，本文將B樣條、仿射變換與GAN結(jié)合，提出了二次判別機(jī)制，對(duì)B樣條仿射變換結(jié)果進(jìn)行初次判別，以弱化全局及局部形變影響，利用圖像小波變換進(jìn)行插值得到粗配準(zhǔn)結(jié)果，最后利用GAN最大化相似性度量輸出配準(zhǔn)結(jié)果。

1 B樣條仿射變換模型

1.1 B樣條仿射變換

仿射變換(AF)源于縮放、旋轉(zhuǎn)、剪切和平移變換的組合的全局線性變換，能保證方向、比例的一致性[5]。但仿射變換的全局性不能解決圖像局部形變帶來的配準(zhǔn)問題，為此引入B樣條，對(duì)形變圖像進(jìn)行粗到精配準(zhǔn)。

B樣條仿射變換的目標(biāo)是找到全局形變場(chǎng)Tglobal以及局部形變場(chǎng)Tlocal，將浮動(dòng)的圖像Ifloat映射到固定(參考)圖像Ifixed，表達(dá)為：

Ifixed=Tlocal(TglobalIfloat)。

(1)

由于B樣條的引入，仿射變換的粗配準(zhǔn)只需6個(gè)參數(shù)來完成圖像的平移、旋轉(zhuǎn)和縮放，不需要額外多的控制點(diǎn)加入。引入齊次坐標(biāo)，增廣維度，全局形變場(chǎng)Tglobal可寫成：

(2)

式中，a11，a12，a13，a21，a22，a23是仿射變換的約束參數(shù)，仿射變換的過程就是最優(yōu)化上述參數(shù)的過程。

為了引入局部控制，將浮動(dòng)圖像網(wǎng)格化，利用控制點(diǎn)的位置更新實(shí)現(xiàn)局部配準(zhǔn)。假設(shè)二維圖像(x,y)，利用m×n網(wǎng)格控制像素點(diǎn)在x軸和y軸上的移動(dòng)，可以表達(dá)為：

(3)

式中，i，j表示控制點(diǎn)編號(hào)；u，v表示與x，y的相對(duì)位置。3次B樣條B0，B1，B2，B3表達(dá)式為：

(4)

1.2 模型代價(jià)函數(shù)

C=Csimilarity(T(Ifloat),Ifixed)+λCsmooth(T)，

(5)

式中，λCsmooth為平滑約束項(xiàng)；Csimilarity為相似性度量，該度量由歸一化誤差均值以及內(nèi)點(diǎn)比構(gòu)成，表示為：

Csimilarity(T(Ifloat),Ifixed)=E(T(Ifloat),Ifixed)-

ninliers/(ninliers+noutliers)。

(6)

1.3 插值實(shí)現(xiàn)

變換后的像素不一定落在整數(shù)上，需對(duì)非整數(shù)坐標(biāo)上的點(diǎn)進(jìn)行灰度插值，通過內(nèi)插點(diǎn)附近的像素的灰度級(jí)來估計(jì)插值灰度[8]。由于線性插值在插值過程中采用插值內(nèi)核而不考慮像素點(diǎn)所處位置，會(huì)產(chǎn)生方塊效應(yīng)從而導(dǎo)致邊緣模糊、紋理特征損失，顯然不符合醫(yī)學(xué)圖像配準(zhǔn)的要求，故本文采用基于小波變換(DWT)的非線性插值方法[9]。該方法往往具有多分辨率分析和逐漸局部細(xì)化等優(yōu)點(diǎn)，其基本思想是將信號(hào)分解到不同的尺度或者分辨率層上，因此可以在不同的尺度上獨(dú)立對(duì)信號(hào)進(jìn)行研究。

張建勛等人[10]對(duì)原始圖像進(jìn)行小波分解，用鄰接子帶的相似性計(jì)算高頻分量，最后對(duì)原始圖像低頻分量做插值運(yùn)算。小波變換的插值實(shí)現(xiàn)過程主要由小波分解和小波逆變換組成，正交小波分解不僅可將圖像的高低頻信息很好地分離，而且分解后各層子帶之間具有相似性[11]。分解后的低頻子帶信號(hào)ML中包含了圖像的絕大部分能量；高頻子帶信號(hào)MH、MV、MD則對(duì)應(yīng)圖像的邊緣信息。通過小波變換，將圖像的高低頻信息分離后，可以單獨(dú)對(duì)高頻信息進(jìn)行處理。低分辨率中的高頻分量MH2、MV2、MD2相似于高分辨率中的MH1、MV1、MD1。利用重構(gòu)理論，將得到的高頻與原有的低頻相疊加，離散小波逆變換就可以得到分辨率更高的圖像[12]，如圖1所示。

圖1 小波變換Fig.1 Wavelet transform

2 構(gòu)建配準(zhǔn)模型

2.1 生成對(duì)抗網(wǎng)絡(luò)

GANs基于兩人零和對(duì)策模型，單獨(dú)交替迭代訓(xùn)練生成模型和判別模型，自動(dòng)學(xué)習(xí)原始真實(shí)樣本集的數(shù)據(jù)分布，假樣本在訓(xùn)練過程中真假變換，以達(dá)到納什平衡[13]。

本文引入二次判別機(jī)制，由仿射變換生成變形場(chǎng)，經(jīng)DWT插值得到中間配準(zhǔn)圖像，由一次判別器判別，根據(jù)RANSAC算法有限迭代修正參數(shù)，直到一次判別為真，再輸入到生成模型，直到二次判別輸出為真。整個(gè)網(wǎng)絡(luò)由B樣條仿射變換和GAN兩部分組成，B樣條仿射變換利用RANSAC算法迭代優(yōu)化數(shù)學(xué)模型參數(shù)生成粗變形場(chǎng)，再由GAN網(wǎng)絡(luò)完成學(xué)習(xí)。模型結(jié)構(gòu)如圖2所示。

圖2 本文配準(zhǔn)模型Fig.2 Proposed registration model

2.2 構(gòu)造生成和判別網(wǎng)絡(luò)

生成器由卷積層、密集殘差塊和反卷積層組成。張桂梅等人[14]在GAN中加入5層密集殘差塊提取深層特征，提高了配準(zhǔn)精度。由于仿射變換的加入，采用一層殘差塊以減小網(wǎng)絡(luò)的學(xué)習(xí)成本。由4層卷積層組成的殘差塊進(jìn)行下采樣，每層均由64個(gè)3×3步長(zhǎng)為1的卷積核、激活函數(shù)和歸一化層組成。由6組反卷積層進(jìn)行上采樣，最后經(jīng)卷積層輸出。生成網(wǎng)絡(luò)結(jié)構(gòu)如圖3所示。

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

生成器輸出作為判別器輸入，經(jīng)7個(gè)卷積層、2個(gè)全連層，最終由分類器輸出結(jié)果。其中卷積層由卷積核、歸一化以及激活函數(shù)組成，第1個(gè)卷積層使用32個(gè)步長(zhǎng)為1的3×3卷積核，第2個(gè)和第3個(gè)使用64個(gè)步長(zhǎng)為1的3×3卷積核，第4個(gè)和第5個(gè)使用128個(gè)步長(zhǎng)為1的3×3卷積核，第6個(gè)和第7個(gè)使用256個(gè)步長(zhǎng)為1的3×3卷積核。判別網(wǎng)絡(luò)結(jié)構(gòu)如圖4所示。

圖4 判別網(wǎng)絡(luò)結(jié)構(gòu)Fig.4 Structure of discriminate network

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

GAN的判別器可以認(rèn)為是一種自主學(xué)習(xí)的損失函數(shù)[15]，區(qū)別于其他固定的損失函數(shù)(如MAE、MSE等)，判別器自適應(yīng)地度量Itrans，Ifixed之間連續(xù)的概率分布。判別網(wǎng)絡(luò)從真實(shí)數(shù)據(jù)分布pdata、生成數(shù)據(jù)分布pG中采樣，用于判別網(wǎng)絡(luò)的損失計(jì)算：

θG,θD=argminmaxV(G,D)=Ex～pdata[lg(D(x))]+

Eg～pG[lg(1-D(G(z)))]，

(7)

式中，G，D分別表示生成網(wǎng)絡(luò)和判別網(wǎng)絡(luò)；x，g分別為真實(shí)數(shù)據(jù)與生成數(shù)據(jù)；pdata為真實(shí)數(shù)據(jù)分布；pG為生成數(shù)據(jù)分布，且x=G(z);z表示生成網(wǎng)絡(luò)輸入。對(duì)于GAN，生成網(wǎng)絡(luò)的期望是最小化損失函數(shù)，而判別網(wǎng)絡(luò)的期望是最大化損失函數(shù)。

2.3.1 判別器損失函數(shù)構(gòu)造

本文二次判別機(jī)制需要構(gòu)造2個(gè)判別損失函數(shù)LD1，LD2，第一判別器主要目的是優(yōu)化仿射變換的6個(gè)參數(shù)以估計(jì)待配準(zhǔn)圖像對(duì)之間的形變大小，即最大化inliers數(shù)量，寫成歸一化形式：

(8)

由于判別器目標(biāo)是最大化真實(shí)數(shù)據(jù)pdata的期望lg(D(x))，最小化生成數(shù)據(jù)PG的期望lg(D(G(z)))，因此第二判別器的期望是最大化損失：

LD2=argmaxV(G,D)=Ex～pdata[lg(D(x))]+

Eg～pG[lg(1-D(G(z)))]。

(9)

2.3.2 生成器損失函數(shù)構(gòu)造

GAN的生成器目的在于使輸出無限接近真實(shí)數(shù)據(jù)，用p(x;θG)表示生成器輸出G(z)與真實(shí)數(shù)據(jù)分布pdata相似的概率，因此需要p(x;θG)越大越好，即G的損失越小越好：

(10)

式中，損失LG包含內(nèi)容損失Lcon以及對(duì)抗損失Ladv，保證了Itrans與Ifixed具有相似的內(nèi)容和結(jié)構(gòu)特征。

互信息NMI匹配2個(gè)圖像之間的聯(lián)合灰度概率分布，廣泛用作多模態(tài)可變形配準(zhǔn)的代價(jià)函數(shù)[17]。結(jié)構(gòu)相似性指標(biāo)SSIM基于灰度值，能夠準(zhǔn)確量化不同圖像之間的對(duì)應(yīng)特征關(guān)系。因此內(nèi)容損失函數(shù)表達(dá)為：

Lcon=-NMI(Itrans,Ifixed)-SSIM(Itrans,Ifixed)。

(11)

對(duì)抗損失約束了配準(zhǔn)圖像之間的灰度和結(jié)構(gòu)信息，包括正向、反向映射轉(zhuǎn)換的對(duì)抗損失，直接引用張桂梅等在文獻(xiàn)[11]中的公式：

Ladv=EIfloat～pdata(pfor(Ifloat,Ifixed))2+

EIfixed～pdata(1-pfor(Ifloat,Gfor(Ifixed)))2+

EIfixed～pdata(prev(Ifloat,Ifixed))2+

EIfloat～pdata(1-prev(Ifixed,Grev(Ifloat)))2，

(12)

式中，pfor表示前向映射中浮動(dòng)圖像與參考圖像之間相似的概率；prev表示反向映射中浮動(dòng)圖像與參考圖像之間相似的概率；E表示相應(yīng)的期望值。

3 實(shí)驗(yàn)分析

3.1 實(shí)驗(yàn)數(shù)據(jù)及配置

Vanderbilt大學(xué)主持的醫(yī)學(xué)圖像配準(zhǔn)評(píng)估項(xiàng)目RIRE作為目前較完善的配準(zhǔn)項(xiàng)目，提供了大量的病人實(shí)例數(shù)據(jù)，包含了19位患有腦部腫瘤、顱內(nèi)出血、肺部腫瘤的患者的醫(yī)學(xué)數(shù)字影像。

每位病人實(shí)例數(shù)據(jù)中包括一套CT及6套MR數(shù)據(jù)，本文挑選了480對(duì)結(jié)構(gòu)紋理清晰、灰度均勻的CT-MR作為數(shù)據(jù)集。用該項(xiàng)目提供的測(cè)試集進(jìn)行測(cè)試實(shí)驗(yàn)硬件及軟件環(huán)境如表1所示。

表1 實(shí)驗(yàn)環(huán)境Tab.1 Experimental environment

3.2 實(shí)驗(yàn)分析

訓(xùn)練前，將所有樣本均預(yù)處理為256 pixel×256 pixel的圖像。用Adam算法，對(duì)判別器和生成器進(jìn)行交替訓(xùn)練，利用RIRE項(xiàng)目的測(cè)試集評(píng)估模型性能。設(shè)置初始學(xué)習(xí)率為0.000 1，batch-size設(shè)置為16，隨著迭代次數(shù)的增加，生成器損失函數(shù)及判別器損失函數(shù)趨勢(shì)如圖5所示。

圖5 損失函數(shù)趨勢(shì)Fig.5 Tendency chart of loss function

由圖5可以看出，生成損失隨著迭代次數(shù)的不斷增加而減小，最終收斂到0.25左右，并且在整個(gè)迭代優(yōu)化的過程中沒有明顯的突變，而對(duì)抗損失穩(wěn)定上升，說明本文模型在整個(gè)訓(xùn)練過程中穩(wěn)定可靠。

醫(yī)學(xué)圖像包括CT、MR等通過不同灰度等級(jí)反映人體組織的密度，因此醫(yī)學(xué)圖像的灰度值可以看成是概率測(cè)度[18]，當(dāng)配準(zhǔn)效果越好時(shí)，對(duì)應(yīng)像素的灰度互信息(NMI)最大，結(jié)構(gòu)相似性(SSIM)最大，偏差(RMSE)越小，表達(dá)式分別為：

(13)

(14)

(15)

式中，H(Gtrans)，H(Gfixed)表示配準(zhǔn)后的圖像與參考圖像的香農(nóng)熵；H(Gtrans,Gfixed)表示2幅圖像之間的聯(lián)合熵；μItrans，μIfixed分別表示圖像Itrans，Ifixed的均值；σItrans，σIfixed分別表示圖像Itrans，Ifixed的標(biāo)準(zhǔn)差；σItrans，Ifixed表示Itrans和Ifixed的協(xié)方差。不同配準(zhǔn)模型的配準(zhǔn)結(jié)果如表2所示。

表2 不同模型圖像配準(zhǔn)指標(biāo)Tab.2 Registration indexes of different models

表2對(duì)比了基于當(dāng)下幾種常見的配準(zhǔn)算法的相關(guān)數(shù)據(jù)，從表中可以看出，本文的算法與傳統(tǒng)的Powell相比，NMI提高了30.48%，RANSAC提高了35.41%，SSIM提高了4.35%，RMSE減小了47.93%；與其他基于GANs的配準(zhǔn)模型相比NMI提高了5.72%，1.39%，RANSAC提高了2.77%，9.59%，SSIM提高了3.85%，2.28%，RMSE減小了17.57%，5.93%，因此本文所提出的配準(zhǔn)模型能夠取得較好的配準(zhǔn)效果。不同模型配準(zhǔn)效果如圖6所示。

(a) 參考圖像

由圖6中可以很明顯地看出，基于Powell模型的效果圖中高亮部分導(dǎo)致參考CT圖像中灰度信息的缺失，而且邊緣信息損失嚴(yán)重。cycleGANs和wGAN模型效果圖灰度分布不均，軟組織以及骨骼組織信息丟失嚴(yán)重，其中cycleGANs頭骨邊緣模糊，軟組織邊緣呈鋸齒狀。而本文模型相較于其他模型不僅在輪廓及拓?fù)浔３稚嫌芯薮蟮膬?yōu)勢(shì)，而且灰度分布更加均勻，配準(zhǔn)更加精準(zhǔn)。

為了驗(yàn)證本文模型對(duì)形變圖像的配準(zhǔn)能力，利用Matlab將浮動(dòng)圖像進(jìn)行不同程度的扭曲變形后輸入配準(zhǔn)模型，配準(zhǔn)效果如圖7所示。

(a) 形變圖1

由圖7可以看出，形變圖1形變量小，效果圖清晰，細(xì)節(jié)豐富；形變圖2形變量中等，畫面模糊，對(duì)比度低；形變圖3因形變量大亮部信息扭曲嚴(yán)重，配準(zhǔn)效果圖畫面模糊嚴(yán)重。由于醫(yī)學(xué)圖像成像設(shè)備姿態(tài)穩(wěn)定，信道抗干擾強(qiáng)，成像過程中很難發(fā)生大形變，因此本文模型在小形變醫(yī)學(xué)圖像配準(zhǔn)中能夠取得較好的配準(zhǔn)效果。

4 結(jié)束語(yǔ)

本文構(gòu)建了基于仿射變換和生成對(duì)抗網(wǎng)絡(luò)的新型配準(zhǔn)模型，將傳統(tǒng)配準(zhǔn)方法與深度學(xué)習(xí)結(jié)合，采用基于隨機(jī)采樣一致算法的仿射變換模型，迭代估計(jì)模型參數(shù)，不斷優(yōu)化變形場(chǎng)。并在生成對(duì)抗網(wǎng)絡(luò)原有的判別機(jī)制上引入了二次判別機(jī)制，實(shí)現(xiàn)了配準(zhǔn)變形場(chǎng)的最優(yōu)化。通過實(shí)驗(yàn)，對(duì)比多種配準(zhǔn)模型的性能指標(biāo)，證明了該模型能夠取得較好的配準(zhǔn)效果。