基于深度卷積特征光流的形變醫學圖像配準算法

張家崗，李達平，楊曉東，鄒茂揚，2，吳 錫，胡金蓉*

（1.成都信息工程大學計算機學院，成都 610225；2.中國科學院成都計算機應用研究所，成都 610041）

（?通信作者電子郵箱hjr@cuit.edu.cn）

0 引言

醫學圖像配準是指對于一幅醫學圖像（浮動圖像）尋求一種或一系列的空間變換，使它與另一幅醫學圖像（固定圖像）上的對應點達到空間位置和解剖位置的完全一致［1］。在許多重要的臨床應用中，醫學圖像配準對醫生診斷病人病情有著十分重要的意義。如在神經外科的立體定向放療中，為提高定位精確性，需要將患者圖像和標準圖譜進行配準；在人群大腦形狀和功能研究中，需要將不同患者的圖像進行配準；在腹部以及胸腔臟器的診斷過程中，需要用到圖像配準來彌補臟器的生理形變［2］。

醫學圖像配準尤其是形變醫學圖像配準一直是醫學圖像分析領域內的一個研究熱點，國內外眾多研究人員對這一問題進行了研究并取得了大量研究成果。其中，基于光流的配準方法是一種重要的形變配準技術［3-5］。光流的概念首先在計算機視覺中提出，其模型可以用來表示物體在觀察成像平面上像素運動的瞬時速度，圖像中每一個像素的運動速度和運動方向組成的速度場就是光流場。基于配準所求形變場與光流場模型所求速度場的相似性，Thirion 等［6］首次將光流場模型引入到了圖像配準問題的求解中，提出了Demons 算法，該算法判斷出浮動圖像上各個像素點的運動方向，通過對各個像素點的移動來實現形變配準。但Demons 僅僅依靠灰度值來驅動配準，特征表示能力不足，所以Liu 等［7-8］將尺度不變特征變換（Scale Invariant Feature Transform，SIFT）特征引入到光流場的估算中，基于圖像間像素的SIFT 特征差異來計算光流，提出稱為SIFT Flow 的光流場計算方法，該方法能夠得到比Demons等傳統光流法更加精確的光流場模型。

在形變醫學圖像配準中，由于人體生理活動，或者個體與標準圖譜的生理差異所導致的大形變問題是不可避免的。Demons 和SIFT Flow 算法只適用于小形變的情況，在出現較大形變時這類傳統光流法將不能保持圖像的拓撲結構，并產生物理上的不合理形變。這是由于Demons 僅僅依靠灰度驅動配準，而忽略了像素之間的空間結構關系；SIFT Flow 雖然將SIFT 特征用于估算光流場，但是受限于手工特征描述符的表征能力不足，尤其是對于高層次抽象特征的提取能力不足，在較大形變的醫學圖像下容易導致誤匹配，從而影響配準精度。所以，如何進行更為準確的特征提取和刻畫不同特征之間的差異成為提高配準精度的關鍵，也是解決形變配準的關鍵。近年來，基于深度學習的方法在計算機視覺和人工智能各領域都取得了很好的效果。所以，本文設計了一種孿生結構的深度卷積神經網絡［9］來提取圖像稠密的深度卷積特征（Deep Convolution Feature，DCF），然后基于圖像像素間的DCF構建能量損失函數計算光流場。這種方法被稱之為基于深度卷積特征光流（Deep Convolution Feature based Optical Flow，DCFOF）的形變醫學圖像形變配準方法。

1 DCFOF方法

本文所提基于深度卷積特征光流的形變配準算法DCFOF的配準框架如圖1所示，主要包含三個部分：深度卷積特征DCF 提取、光流場（形變場）估算和變換插值，其中，DCF提取是本算法的關鍵步驟。算法先逐像素地提取固定圖像和浮動圖像的圖像塊對，并利用孿生深度卷積神經網絡提取圖像塊的特征；然后由提取到的特征構建能量損失函數，通過最小化能量損失函數，求解得形變場，即光流場；最后，根據形變場，對浮動圖像逐像素地對齊，變換插值過程中采用三次樣條插值作為插值方案，得到與固定圖像極為接近的配準后圖像。

1.1 DCF提取

圖2 展示了本文所設計的孿生深度卷積神經網絡的結構圖，該網絡由兩個權重相同的深度卷積神經網絡組成。在每個卷積神經網絡中，網絡的輸入是一個以像素為中心的圖像塊，輸出為128 維的特征向量，即所提取的DCF。網絡所有的卷積層都使用3× 3的卷積核，步長為1，使用LeakyReLU 作為激活函數。不同于卷積層的是，密集層使用ReLU 作為激活函數。此外，在Flatten 層之后添加一個Dropout 層，以避免過擬合。網絡中各層的參數也如圖2所示。

圖2 孿生深度卷積神經網絡結構Fig.2 Architecture of siamese convolutional neural network

對于輸入的樣本對，本文將沒有形變的圖像塊對稱為正樣本對，有形變的圖像塊對稱為負樣本對。據此，本文采用了對比損失（Contrastive Loss）作為網絡的損失函數，該函數能夠同時最小化正樣本對或負樣本對之間的距離，達到網絡訓練的目的，其定義如下：

其中：N是輸入樣本的個數；為每對樣本的歐氏距離；margin是一個為1的常數［11］；yi是輸入樣本對的標簽，1 為正樣本，0 為負樣本，當為正樣本對時，則優化di盡量減小，當為負樣本對時，則是讓max(margin-di，0)2盡量減小。

通過對比損失函數學習到的映射關系，可以使得在高維空間中相同類別但距離較遠的點在特征維空間中距離更近；而不同類別但距離較近的點在特征維空間中距離更遠。因此，相較于SIFT 特征以及一般的深度學習特征如VGG（Visual Geometry Group）［12］等，基于對比損失函數訓練的孿生深度卷積神經網絡提取的DCF 的區分度更高、更穩定，更適宜于差異計算和得到更加準確的計算結果。網絡訓練的數據集及具體實現見第2章。

1.2 光流場估計

本文受SIFT Flow 算法啟發，通過估計密集的形變場來進行配準。該方法的目標是最小化由DCF 構建的能量損失函數，從而獲得像素的位移。能量損失函數定義如下：

其中：S1和S2分別是固定圖像和浮動圖像的DCF向量；w(p)=(u(p)，v(p))是點p=(x，y)的流向量，u(p)和v(p)分別是流向量在水平和垂直方向上的分量；ε是點p在一定范圍內的空間鄰域；參數η和α分別是正則項系數。該函數分為三個部分：第一部分為數據項，用于計算固定圖像點p和浮動圖像中匹配點的絕對值誤差和，t為截斷閾值，用于剔除不合理的流向量；第二部分為正則項，用于約束流向量w(p)在合理范圍內；第三部分同為正則項，用于約束點p鄰域ε內一點q的流向量趨同于點p的流向量。

求解流向量的目標是最小化能量損失函數，本文采用了雙層置信度傳播算法［13］來求解光流場，同時為了加速運算，采用了從粗到精的特征匹配策略［14］。求得光流場后，對浮動圖像采用三次樣條插值算法［15］做變換插值即可得到配準圖像。

從上述能量函數的定義可以看出，所求解光流場的精確度直接與所提取像素點的特征質量相關，DCF 不僅具有更穩定和區分度更高等深度學習特征優點，還具有易于差異計算的特點。因此，本文提出的基于圖像間像素的DCF 差異計算所得的光流場更加精確和魯棒。

2 實驗和結果分析

2.1 數據集和標簽

本文使用的訓練數據主要是EMPIRE10（Evaluation of Methods for Pulmonary Image Registration 2010）數據集［16］和ACDC（Automated Cardiac Diagnosis challenge）數據集［17］。其中，EMPIRE10是一個包含30對臨床胸腔CT 掃描的肺部數據集，同時含有Van Rikxoort 等［18］用肺部分割方法制作的肺部Mask 標簽；ACDC 是由150 名患者的MRI 掃描組成的心臟數據集，其中，固定圖像和浮動圖像分別為不同時刻掃描所得，固定圖像和浮動圖像之間的差異為心臟活動產生的真實形變，并且患者的數據包含專家標注的心臟Mask標簽。

實驗中隨機選取了60 張EMPIRE10 圖像和60 張ACDC圖像作為訓練集制作所需的圖像塊對，另隨機選取了和訓練圖像不同的20 張EMPIRE10 圖像和20 張ACDC 圖像制作驗證集所需的圖像塊對。本文采用了Simard 等［19］的方法，對采集的圖像進行隨機形變，形變程度由形變強度參數σ控制，σ取值范圍為50～250，σ越大形變程度越大。將形變前的圖像作為固定圖像，形變后的圖像作為浮動圖像，然后以固定圖像和浮動圖像中的每個像素點為中心取15× 15大小的圖像塊，最后將這些圖像塊用作訓練數據集和驗證數據集。其中，來自同一固定圖像、同一位置的圖像塊作為正樣本對，來自固定圖像和對應待配準形變圖像（浮動圖像）上同一位置的圖像塊作為負樣本對。

如圖3 所示，每個縱列由上下兩個15× 15 圖像塊組成。正樣本對是完全相同的圖像塊，而負樣本對之間存在著一定程度的形變。

圖3 網絡訓練的正負樣本示例Fig.3 Positive and negative samples of network training

2.2 實驗設置

本文算法將和Demons 算法、SIFT Flow 算法和專業軟件Elastix 進行實驗對比。其中，Elastix 是一個用于剛性和非剛性圖像配準的軟件，包含一系列通常用于解決醫學圖像配準問題的算法，常作為不同配準算法配準效果的比較標準。

實驗利用2.1 節所描述的方法制作訓練數據集，共制作得到8 000 000 和1 600 000 對圖像塊分別作為孿生深度卷積神經網絡的訓練數據和驗證數據，其中正負樣本比例均為1∶1。

網絡框架采用Keras（Tensorflow 后端），網絡使用RMSprop 優化算法，batch size 設置為1 000，epochs 設置為50，每個epoch 的steps 設置為100，初始學習率設置為0.001，動量因子設置為0.9。采用1塊NVIDIA Tesla K40C GPU 作為訓練GPU，每個訓練輪次大概需要100 s。

此外，實驗中對比算法的主要參數設置如下：Demons 算法中，直方圖級別數設置為1 024，迭代次數設置為50；SIFT Flow 光流法中，正則項系數η設置為0.005，α設置為2，迭代次數設置為200；Elastix 算法中，變換類型設置為“BSPLINE”，迭代次數設置為500。

2.3 算法評測指標

2.3.1 顏色疊加

顏色疊加是可視化評估配準效果的最常用技術之一。通過將兩張圖像按照一定方式映射到RGB 色彩空間中，以可視化兩張圖片的差異。首先，將固定圖像填充到R通道中，將浮動圖像填充到G 通道中，最后將兩幅圖像的均值填充到B 通道中。本文直接使用SimpleITK 工具包［20］獲得顏色疊加的可視化，其中，準確配準的區域顯示為圖像原色，而紅色或綠色區域則表明沒有準確配準。圖4給出了顏色疊加的示例。

圖4 固定圖像和浮動圖像顏色疊加Fig.4 Color overlapping of fixed image and floating image

2.3.2 均方根差

均方根差（Root Mean Squared Difference，RMSD）可以反映兩幅圖之間的差異，其計算式如下：

其中：I1和I2是需要計算RMSD 的兩幅圖像；I1(xi)和I2(xi)是兩幅圖中相同位置點的灰度值；ΩI為I1和I2的圖像域；|ΩI|為I1的像素個數（和I2中像素個數相同）。式（3）通過計算灰度值的均方根誤差來衡量兩幅圖的相似性。固定圖像和配準后圖像的RMSD 值越小，表示兩幅圖像中相同位置點的灰度值越近似，整幅圖也越相似，配準效果越好。

2.3.3 DICE系數

DICE 系數用來評估感興趣區域的配準精度，通過計算感興趣區域（Region of Interest，ROI）的重合度來評估配準效果［21］，其表達式如下：

其中，X、Y分別是兩幅圖ROI的面積，X∩Y是重疊部分面積。固定圖像和配準后圖像的ROI所求DICE 值越大，表示ROI重合度越高，配準效果越好。對于EMPIRE10 肺部數據集，ROI為肺部Mask。圖5給出了肺部圖像和肺部Mask的示例。

圖5 肺部圖像與肺部Mask Fig.5 Lung image and lung Mask

2.4 實驗結果

2.4.1 EMPIRE10模擬形變圖像實驗

圖6 中（a）和（b）分別為固定圖像和浮動圖像，浮動圖像有較大形變（σ=250）。

圖6 模擬形變實驗中固定圖像和浮動圖像（σ=250）Fig.6 Fixed image and floating image in simulated deformation experiment（σ=250）

圖7 是四種對比算法的顏色疊加圖。其中，（a）～（d）分別為Demons、SIFT Flow、Elastix 和本文算法DCFOF 的配準結果和固定圖像的顏色疊加圖。特別的，圖中方框放大所展示的結果中，DCFOF 所展示結果中，沒有綠色或者紅色的偽影，表明其對肺葉纖維的配準精度遠遠高于其他三種方法。

圖7 圖6配準圖像與固定圖像顏色疊加圖（σ=250）Fig.7 Color overlapping of registered image and fixed image for Fig.6（σ=250）

圖8 和圖9 分別展示了在不同形變強度σ下（依次取50、100、150、200、250），四種算法的RMSD 值和DICE 值變化趨勢，其中，每種形變強度做30組平行實驗取平均值。

圖8 配準圖像與固定圖像的RMSD值Fig.8 RMSD values of registered image and fixed image

圖8 中，隨著形變強度的增大，DCFOF 相較其他三種方法，其RMSD 值比較穩定，在σ=250 的較大形變下，其RMSD值小于5，處理大形變能力較強。

圖9 配準圖像ROI與固定圖像ROI的DICE值Fig.9 DICE values of registered image ROI and fixed image ROI

圖9 中，DCFOF 的DICE 值高于其他幾種方法，表明DCFOF能夠更好地對齊ROI。

2.4.2 ACDC真實形變圖像實驗

為了驗證DCFOF 對真實形變醫學圖像的配準能力，選取了ACDC作為真實形變醫學圖像進行配準。圖10中（a）和（b）分別不同時刻的兩幅圖像，時間差值為15 s，（a）、（b）之間的形變為心臟生理活動產生的真實形變，將圖（a）作為固定圖像，圖（b）作為浮動圖像。

圖10 真實形變實驗中固定圖像和浮動圖像（σ=250）Fig.10 Fixed image and floating image in real deformation experiment（σ=250）

圖11 展示了四種對比算法對ROI（心臟區域）的配準能力，圖片截取了心臟區域進行對比展示。每張圖片中，紅色（深）輪廓線為固定圖像ROI的輪廓線，作為基準線，綠色（淺）為配準后圖像ROI 的輪廓線，所有輪廓線均根據ROI 的Mask輪廓勾畫出。其中，（a）～（e）分別為浮動圖像、Demons 配準結果、SIFT Flow 配準結果、Elastix 配準結果和本文方法DCFOF配準結果的ROI輪廓線與固定圖像ROI輪廓線的重疊。特別的，圖片分三列展示了不同的ROI對齊效果，第一列為左心室的輪廓，第二列為右心室的輪廓，第三列為心肌外壁的輪廓。

從圖11 可以看出，在生理形變中，DCFOF 對ROI 的配準后輪廓線（綠色（淺））與基準線（紅色（深））重合較好，表明DCFOF能夠很好地對齊ROI。

表1 展示了四種算法對ROI 配準的DICE 值，ROI 分為心臟整體、左心室、右心室和心肌壁，表中結果為30 組平行實驗的平均值。

從表1 可以看出，DCFOF 對各個ROI 的配準精度均高于其他幾種算法。

圖11 圖10 ROI配準的輪廓線圖Fig.11 Contour of ROI registration for Fig.10

表1 ROI配準的DICE值Tab.1 DICE values of ROI registration

3 分析和討論

本文針對現有光流法特征提取不足、無法精確配準的問題，提出了基于深度卷積特征光流的形變醫學圖像配準算法。通過對EMPIRE10 和ACDC 醫學圖像進行配準實驗，結果表明本文算法DCFOF 與Demons 算法、SIFT Flow 算法以及Elastix軟件相比，有著更好的配準精度。

在圖7的顏色疊加圖中，DCFOF有著最精確的配準結果。Demons、SIFT Flow、DCFOF 雖然都是基于光流場模型的算法，它們根據像素之間的差異進行配準，但是三種方法對像素的特征描述不同。Demons 僅僅依靠像素灰度值驅動配準，而不考慮像素周圍的空間信息，較大的形變使得圖像的拓撲結構發生較大變化，僅僅依靠灰度值無法記錄像素周圍原有的空間信息，所以不能精確配準。SIFT Flow 雖然將圖像的梯度信息用作記錄像素的空間信息，但是過大的形變導致像素梯度信息發生強烈變化，使得原本匹配的像素點變得不再匹配，從而導致也無法精確配準。DCFOF 提取了以像素為中心的圖像塊的高層次抽象的特征，實驗結果表明，即使發生較大形變，也能夠實現特征點的準確匹配，估算出精確的形變場，從而精確地處理形變配準。另外，Elastix 作為專業的配準軟件，最優化求解整幅固定圖像和浮動圖像之間的變換參數，而忽略了圖像的局部空間信息，因此其處理形變的能力有限。而DCFOF 逐像素地提取圖像塊的DCF 特征作為圖像的局部空間信息，因此對形變更為敏感，所以Elastix的配準精度要低于DCFOF。

在圖8 的RMSD 對比結果中，四種算法隨著形變強度增大表現出了相似的變化趨勢。在50 和100 的形變程度下，四種算法的RMSD 相差不大，表明其配準精度也相差不大；但是當形變程度增大到150 時，Demons 算法性能急劇下降，如圖8中RMSD 值已接近10，在σ=250 時，Demons 算法的RMSD 值已經超過20，此時算法已經失效，說明Demons 無法處理較大形變，這與上面的分析結果一致。SIFT Flow 的配準效果從σ=150 開始降低，同樣無法精確處理150 以上的形變。Elastix 作為專業配準軟件，在50～250 的形變程度下，表現出了不錯的配準效果，但仍然不如本文算法DCFOF，具體表現為隨著形變強度增大，Elastix 的RMSD 值和DCFOF 逐漸拉開差距，性能下降較大，而DCFOF 性能表現穩定，性能下降較小。這是由于Elastix 并沒有像DCFOF 一樣考慮了圖像的局部特征，所以在對形變的敏感程度上不如DCFOF，配準精度也不如DCFOF。

圖9 所展示的ROI 的DICE 值對比結果可以看出，DCFOF在所有形變強度下DICE 值最大，即使在σ=250 的較大形變強度下，其DICE值也在0.99以上，表明DCFOF能很好地對齊肺部ROI，配準效果好。Demons 和SIFT Flow 算法從σ=150開始，DICE 值開始迅速下降，表明Demons 和SIFT Flow 并不能很好地處理較大形變圖像的配準。Elastix 和DCFOF 的差距也越來越大，表現不如DCFOF穩定。

ACDC真實形變圖像實驗進一步驗證了DCFOF對真實形變醫學圖像的處理能力，從圖11 的ROI 配準的輪廓線圖可以看出，DCFOF 對ROI 的對齊能力較高，配準后的綠線和紅線重合度高，而Demons 并不能很好地對齊ROI，SIFT Flow 和Elastix 與DCFOF 相比，對ROI 邊緣的平滑程度不如DCFOF，在心肌外壁的對齊結果中，也能看出DCFOF 優于SIFT Flow和Elastix。表1 展示了30 組平行實驗的平均DICE 值，可以看出，DCFOF 的確能夠處理真實的形變醫學圖像，且優于專業軟件Elastix。

4 結語

本文提出了一種基于深度卷積特征光流的醫學圖像配準算法，其提取的深度卷積特征區分度高、穩定性強，基于該特征所計算得到的光流場更接近于真實形變場。通過實驗的定性、定量分析，實驗結果表明本文算法DCFOF 能夠精確地進行形變醫學圖像配準，且能處理真實的形變醫學圖像，效果要優于Demons算法、SIFT Flow算法以及Elastix軟件。

在未來的工作中，我們將進一步擴展現有的工作：1）研究基于端到端深度學習特征的光流方法，進一步擴展DCFOF，進行形變圖像配準，同時減少運行時間；2）研究基于深度學習特征的光流方法進行多模態形變圖像配準。