基于B樣條仿射和GANs的醫學圖像配準

宋 梟，朱家明，徐婷宜

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

0 引言

圖像配準是圖像研究領域的一個典型問題和技術難點，其目的在于比較或融合針對同一對象在不同條件下獲取的圖像。隨著醫學成像設備的不斷進步，可以采集更加精準的器官解剖圖像，但僅憑醫務工作者的空間想象和主觀經驗觀察病情的發展十分困難。精確快速的醫學圖像配準是醫學深入分析的基礎，是分析病灶、觀察器官變化、引導臨床手術的關鍵。為此，廣大學者提出了許多醫學圖像配準方法，從傳統的基于特征和灰度的方法到深度學習方法，圖像配準的性能得到了質的飛躍。

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

1 B樣條仿射變換模型

1.1 B樣條仿射變換

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

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

Ifixed=Tlocal(TglobalIfloat)。

(1)

由于B樣條的引入，仿射變換的粗配準只需6個參數來完成圖像的平移、旋轉和縮放，不需要額外多的控制點加入。引入齊次坐標，增廣維度，全局形變場Tglobal可寫成：

(2)

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

為了引入局部控制，將浮動圖像網格化，利用控制點的位置更新實現局部配準。假設二維圖像(x,y)，利用m×n網格控制像素點在x軸和y軸上的移動，可以表達為：

(3)

式中，i，j表示控制點編號；u，v表示與x，y的相對位置。3次B樣條B0，B1，B2，B3表達式為：

(4)

1.2 模型代價函數

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

(5)

式中，λCsmooth為平滑約束項；Csimilarity為相似性度量，該度量由歸一化誤差均值以及內點比構成，表示為：

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

ninliers/(ninliers+noutliers)。

(6)

1.3 插值實現

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

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

圖1 小波變換Fig.1 Wavelet transform

2 構建配準模型

2.1 生成對抗網絡

GANs基于兩人零和對策模型，單獨交替迭代訓練生成模型和判別模型，自動學習原始真實樣本集的數據分布，假樣本在訓練過程中真假變換，以達到納什平衡[13]。

本文引入二次判別機制，由仿射變換生成變形場，經DWT插值得到中間配準圖像，由一次判別器判別，根據RANSAC算法有限迭代修正參數，直到一次判別為真，再輸入到生成模型，直到二次判別輸出為真。整個網絡由B樣條仿射變換和GAN兩部分組成，B樣條仿射變換利用RANSAC算法迭代優化數學模型參數生成粗變形場，再由GAN網絡完成學習。模型結構如圖2所示。

圖2 本文配準模型Fig.2 Proposed registration model

2.2 構造生成和判別網絡

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

圖3 生成網絡結構Fig.3 Structure of generative network

生成器輸出作為判別器輸入，經7個卷積層、2個全連層，最終由分類器輸出結果。其中卷積層由卷積核、歸一化以及激活函數組成，第1個卷積層使用32個步長為1的3×3卷積核，第2個和第3個使用64個步長為1的3×3卷積核，第4個和第5個使用128個步長為1的3×3卷積核，第6個和第7個使用256個步長為1的3×3卷積核。判別網絡結構如圖4所示。

圖4 判別網絡結構Fig.4 Structure of discriminate network

2.3 構造損失函數

GAN的判別器可以認為是一種自主學習的損失函數[15]，區別于其他固定的損失函數(如MAE、MSE等)，判別器自適應地度量Itrans，Ifixed之間連續的概率分布。判別網絡從真實數據分布pdata、生成數據分布pG中采樣，用于判別網絡的損失計算：

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

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

(7)

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

2.3.1 判別器損失函數構造

本文二次判別機制需要構造2個判別損失函數LD1，LD2，第一判別器主要目的是優化仿射變換的6個參數以估計待配準圖像對之間的形變大小，即最大化inliers數量，寫成歸一化形式：

(8)

由于判別器目標是最大化真實數據pdata的期望lg(D(x))，最小化生成數據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 生成器損失函數構造

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

(10)

式中，損失LG包含內容損失Lcon以及對抗損失Ladv，保證了Itrans與Ifixed具有相似的內容和結構特征。

互信息NMI匹配2個圖像之間的聯合灰度概率分布，廣泛用作多模態可變形配準的代價函數[17]。結構相似性指標SSIM基于灰度值，能夠準確量化不同圖像之間的對應特征關系。因此內容損失函數表達為：

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

(11)

對抗損失約束了配準圖像之間的灰度和結構信息，包括正向、反向映射轉換的對抗損失，直接引用張桂梅等在文獻[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表示前向映射中浮動圖像與參考圖像之間相似的概率；prev表示反向映射中浮動圖像與參考圖像之間相似的概率；E表示相應的期望值。

3 實驗分析

3.1 實驗數據及配置

Vanderbilt大學主持的醫學圖像配準評估項目RIRE作為目前較完善的配準項目，提供了大量的病人實例數據，包含了19位患有腦部腫瘤、顱內出血、肺部腫瘤的患者的醫學數字影像。

每位病人實例數據中包括一套CT及6套MR數據，本文挑選了480對結構紋理清晰、灰度均勻的CT-MR作為數據集。用該項目提供的測試集進行測試實驗硬件及軟件環境如表1所示。

表1 實驗環境Tab.1 Experimental environment

3.2 實驗分析

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

圖5 損失函數趨勢Fig.5 Tendency chart of loss function

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

醫學圖像包括CT、MR等通過不同灰度等級反映人體組織的密度，因此醫學圖像的灰度值可以看成是概率測度[18]，當配準效果越好時，對應像素的灰度互信息(NMI)最大，結構相似性(SSIM)最大，偏差(RMSE)越小，表達式分別為：

(13)

(14)

(15)

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

表2 不同模型圖像配準指標Tab.2 Registration indexes of different models

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

(a) 參考圖像

由圖6中可以很明顯地看出，基于Powell模型的效果圖中高亮部分導致參考CT圖像中灰度信息的缺失，而且邊緣信息損失嚴重。cycleGANs和wGAN模型效果圖灰度分布不均，軟組織以及骨骼組織信息丟失嚴重，其中cycleGANs頭骨邊緣模糊，軟組織邊緣呈鋸齒狀。而本文模型相較于其他模型不僅在輪廓及拓撲保持上有巨大的優勢，而且灰度分布更加均勻，配準更加精準。

為了驗證本文模型對形變圖像的配準能力，利用Matlab將浮動圖像進行不同程度的扭曲變形后輸入配準模型，配準效果如圖7所示。

(a) 形變圖1

由圖7可以看出，形變圖1形變量小，效果圖清晰，細節豐富；形變圖2形變量中等，畫面模糊，對比度低；形變圖3因形變量大亮部信息扭曲嚴重，配準效果圖畫面模糊嚴重。由于醫學圖像成像設備姿態穩定，信道抗干擾強，成像過程中很難發生大形變，因此本文模型在小形變醫學圖像配準中能夠取得較好的配準效果。

4 結束語

本文構建了基于仿射變換和生成對抗網絡的新型配準模型，將傳統配準方法與深度學習結合，采用基于隨機采樣一致算法的仿射變換模型，迭代估計模型參數，不斷優化變形場。并在生成對抗網絡原有的判別機制上引入了二次判別機制，實現了配準變形場的最優化。通過實驗，對比多種配準模型的性能指標，證明了該模型能夠取得較好的配準效果。