圖像配準在肝臟疾病的計算機輔助診斷和手術規劃方面具有重要的臨床意義。基于深度學習的配準方法使得肝臟電子計算機斷層掃描（CT）圖像配準具有較高的實時性和準確度。然而，現有方法在配準具有大位移和大形變的圖像時，存在配準后圖像的紋理信息發生改變的問題，因而難以將其應用在后續的圖像處理與臨床診斷中。基于此，本文提出一種新穎的基于紋理濾波的無監督配準方法，實現了肝臟 CT 圖像的配準。該方法首先基于 L0 梯度最小化的紋理濾波算法消除 CT 圖像中肝臟表面的紋理信息，使得配準過程僅參考兩幅圖像的空間結構信息進行配準，從而解決紋理改變的問題。然后，基于級聯網絡配準具有大位移和大形變的圖像，循序漸進地將待配準圖像與參考圖像在空間結構上對齊。此外本文還提出一種新的衡量指標——直方圖相關系數，以衡量配準后紋理改變的程度。實驗結果表明，本文所提出的方法具有較高的配準精度，有效地改善了級聯網絡中存在紋理改變的問題，并且提升了在空間結構對應和抗折疊性能兩方面的配準效果。因此，本文所提方法或有助于提升醫學圖像配準的科學性，促進醫學圖像配準安全可靠地應用在肝臟疾病的計算機輔助診斷和手術規劃方面。
引用本文： 王鵬, 嚴赟琦, 錢黎俊, 所世騰, 郭翌, 許建榮, 汪源源. 基于紋理濾波的無監督配準方法及其在肝臟電子計算機斷層掃描中的應用. 生物醫學工程學雜志, 2021, 38(5): 819-827. doi: 10.7507/1001-5515.202102001 復制
引言
醫學圖像配準，指在任意兩幅醫學圖像之間建立一種空間對應關系，通過這種對應關系可以使一幅圖像變換到另一幅圖像所在的空間結構中。臨床診斷中，圖像配準可以幫助醫生從不同成像設備獲取的醫學圖像或者不同時間獲取的來自相同成像設備的醫學圖像獲取信息并進行融合互補，提高疾病診斷精度[1-3]。對于計算機輔助診斷而言，醫學圖像配準可以輔助器官和腫瘤分割、多模態圖像融合以及疾病分類等研究[4-6]。
近年來，傳統配準方法例如彈性體（Elastix）、阿恩特（Ant）已經在各種任務上實現了較好的配準效果[7-8]，但是這些方法所需的配準時間較長，不適用于實時的臨床應用。最近，基于深度學習的醫學圖像配準方法呈現出巨大的優勢，不僅在配準效果上可與眾多傳統方法媲美，而且消耗的時間遠遠小于傳統配準方法。然而，基于深度學習的配準方法大多集中于腦部圖像配準[9-11]。由于腦部圖像產生的位移較小，因此這些方法遷移到位移和形變均較大的肝臟圖像中并不能取得較好的效果。一些論文提出了通過多級級聯網絡配準具有大位移和大形變醫學圖像的方法[12-14]，在肝臟電子計算機斷層掃描（computed tomography，CT）圖像中應用取得了良好的配準效果。但是通過多級級聯網絡產生的空間變換可能使得配準后的圖像與參考圖像極為相似，不僅空間結構一致，紋理方面也趨于相同。而在肝臟 CT 圖像中，紋理特點直接反映了腫瘤或病變組織的類型，如圖 1 所示，綠色箭頭標注了肝臟病灶的位置，血管瘤和肝癌的病灶具有不同的紋理信息，這些紋理信息是臨床診斷和計算機輔助診斷過程中的重要依據，因此紋理信息的改變將會嚴重影響疾病診斷的精確度。由此可見，迫切需要一種不改變紋理信息的肝臟 CT 圖像配準方法。此外，現有的配準衡量指標一般采用戴斯（Dice）系數和特征點距離等，這些指標均針對于配準后解剖結構的對應程度而提出，未有衡量配準圖像與原圖像在紋理信息上的對應程度的指標。
圖1
肝臟 CT 圖像中的紋理特征
Figure1.
Texture features in liver CT images
針對上述問題，本文提出了一種新型的三維無監督深度學習配準框架用于肝臟 CT 圖像中。首先對肝臟 CT 圖像進行基于 L0 梯度最小化的紋理濾波處理，該方法可以消除三維醫學圖像中的大部分復雜的紋理信息和圖像噪聲，保留顯著的空間結構信息；然后采用級聯網絡對輸入圖像進行循序漸進的配準，以解決大位移和大形變問題；最后提出一種直方圖相關系數衡量配準后圖像的紋理信息與原圖像的對應程度，以關注配準前后肝臟 CT 圖像紋理信息變化。實驗證明，本文所提出的方法不僅可以實現高精度配準，而且有效地保留了肝臟影像中的紋理信息。因此，本文的方法對于提升醫學圖像配準的科學性，使得醫學圖像配準安全可靠地應用在肝臟疾病的計算機輔助診斷和手術規劃方面具有重要意義。
1 方法
本文提出的配準框架如圖 2 所示。假設 If 為配準過程中的參考圖像，Im 為待配準圖像，配準后圖像為 Ir。配準的過程即將待配準圖像 Im 配準到參考圖像 If 所在的空間結構之中，從而得到配準圖像 Ir。參考圖像 If 和待配準圖像 Im 均為三維 CT 圖像。首先，Im 與 If 經過紋理濾波后分別得到 Sm 與 Sf，然后將 Sm 與 Sf 輸入至級聯網絡進行配準。由于濾波后的兩幅圖像 Sm 與 Sf 不包含紋理信息，因此級聯網絡僅按照兩幅圖像的空間結構對應進行配準，并預測出用于指導配準的變形場 
。最后利用 
對 Im 進行插值變換得到最終的配準圖像 Ir，該過程可被記為 Ir = Im
。
圖2
基于紋理濾波的無監督級聯式配準模型
Figure2.
Unsupervised cascaded network registration model based on texture filtering
1.1 基于 L0 梯度最小化的紋理濾波
L0 梯度最小化算法最早由 Xu 等[15]提出，主要用于解決二維自然圖像的平滑問題。本文將其應用在醫學圖像中，并且在此基礎上實現了針對三維醫學圖像的紋理濾波。基于 L0 梯度最小化的紋理濾波以 L0 梯度數量為正則項，通過限制非零梯度的數量完成紋理濾波。
假設 I 和 S 分別表示輸入圖像和紋理濾波后的圖像。紋理濾波的目標函數如式（1）所示：
 |
其中，∑p（Sp ? Ip）2保證了圖像 I 和 S 在每個像素點 p 的相似度，λ 為非負參數，作用是控制紋理濾波的程度，λ 越大，紋理濾波程度越高。C（S）代表圖像 S 的 L0 梯度的統計，其定義如式（2）所示:
 |
其中，
xSp、
ySp 和 
zSp 分別表示圖像 S 的任意像素點 p 沿著 x、y 和 z 方向的偏導數。#{}的功能是計算|
xSp| + |
ySp| + |
zSp|的值不為 0 的像素點數量之和，|
xSp| + |
ySp| + |
zSp|代表了圖像 S 在像素點 p 的梯度能量。由于 C（S）對全局不連續性進行了統計建模，故無法通過傳統的梯度下降法或其他離散優化方法優化目標函數，因此，本文采用交替最小化算法優化上述目標函數。
為了使用交替最小化算法，式（2）改寫為如式（3）所示：
 |
其中，h、v、k 為引入的輔助變量，hp、vp、kp 分別表示 h、v、k 在像素點 p 的元素，C（h，v，k）= #{p||hp| + |vp| + |kp| ≠ 0}，β 是通過迭代方式更新的參數，用于控制變量（hp，vp，kp）及其對應的梯度（
xSp，
ySp，
zSp）之間的相似度。當 β 足夠大時，式（1）與式（3）接近。
交替最小化算法將式（3）分解為兩個子問題加以解決，分別為：
（1）已知 hp、vp 和 kp，計算 S 以最小化式（3）。式（3）可以另被寫為如式（4）所示：
 |
為了采用快速傅里葉變換估計式（4）中的 S，式（4）改寫為如式（5）所示：
 |
其中，F（·）代表快速傅里葉變換，F?1（·）代表快速傅里葉的反變換，F*（·）表示在快速傅里葉變換基礎上取其共軛復數形式，F（1）表示階躍函數 δ（1）的傅里葉變換。
（2）已知 S，計算 hp、vp 和 kp 以最小化式（3）。式（3）可以另被寫為如式（6）所示：
 |
式（6）可以在空間上分解，分別對每個 hp、vp 和 kp 進行估算。分解后的式（6）如式（7）所示：
 |
其中，關于 H（|hp| + |vp| + |kp|）的計算如式（8）所示：
 |
對于每一個像素點 p，式（7）可以被寫為如式（9）所示：
 |
其中，Ep 代表像素點 p 的目標函數值，為了使 Ep 達到最小值，需要滿足如式（10）所示的條件：
 |
因此，結合式（8）和（9），可以求解出使 Ep 達到最小值的（hp，vp，kp）。
利用兩個子問題進行交替最小化的算法框圖參見算法 1——L0 梯度最小化算法。參數 β 從最小值 β0 開始自動進行迭代，每次都乘以速率 κ，直到 β ≥ βmax。
算法 1：
STN1：輸入，圖像 I，紋理濾波權重 λ，參數 β0，βmax，κ；
步驟 2：初始化，S←I，β←β0；
步驟 3：已知 S，根據式（10）求解 hp、vp 和 kp；
步驟 4：已知 hp、vp 和 kp，根據式（5）更新 S；
步驟 5：β←κβ；
步驟 6：重復步驟 3、4、5，直至 β ≥ βmax；
步驟 7：輸出，處理后的圖像 S。
1.2 基于級聯網絡的無監督配準模型
參考圖像和待配準圖像經過紋理濾波后，輸入到基于級聯網絡的無監督配準模型中進行配準。如圖 3 所示，級聯網絡的基本框架由多個子網絡構成，通過多個子網絡，可以循序漸進地將待配準圖像配準到參考圖像所在的空間結構中。
圖3
級聯網絡框架
Figure3.
The framework of cascaded network
1.2.1 級聯網絡
由于肝臟 CT 圖像的配準任務具有大位移和大形變的特點，這里采用級聯網絡進行配準。級聯網絡由多個子網絡組成，如圖 3 所示。不失一般性，第 k 個子網絡以參考圖像 Sf 和待配準圖像 
（k = 1，2，
）為輸入，預測一個變形場 
。對于第一個子網絡，
= Sm。然后，每個子網絡預測的變形場 
需要和前級子網絡預測的變形場 
進行合成，形成在本級最終起作用的變形場 
，變形場合成的迭代操作記為 
，其中“
”表示由 Jaderberg 等[16]提出的空間變換器網絡（spatial transformer network，STN）。最后，在每級子網絡后，待配準圖像 
利用本級起作用的變形場 
完成一次配準，得到下一級的待配準圖像 
，該操作記為 
= 
。
和 Sf 將輸入至下一個級聯的子網絡中繼續配準。通過多次配準使得 Sf 和 
之間的相似度不斷增大，從而循序漸進地將 Sm 配準到 Sf 所在的空間結構中。由于 STN 是一個可微分的模塊，滿足神經網絡反向傳播梯度的需要，因此可以實現對級聯網絡的端到端訓練。
1.2.2 子網絡
級聯網絡中的子網絡采用的是卷間網絡（volume tween network，VTN）[12]，如圖 4 所示，VTN 網絡遵循編碼器—解碼器的結構，以參考圖像和待配準圖像沿著通道所在維度拼接而成的多維向量作為輸入，為了表示方便，四維的特征層在圖 4 中以二維表示，面積較小的特征層意味著更低的空間分辨率，通道數量標注在每一層的上方，彎曲的箭頭代表箭頭始端的特征層與箭頭末端的特征層沿著通道所在維度進行疊加。在編碼器部分，VTN 網絡采用帶步長的三維卷積層將每層的空間分辨率減少一半以及將每層的通道數量增加一倍，但是從第四個卷積層開始，每降低一次分辨率都將在該三維卷積層后增加一個不改變空間分辨率和通道數量的三維卷積層；在解碼器部分，VTN 網絡使用三維反卷積層（或稱為轉置卷積層）逐步恢復至原來空間分辨率和通道數量。類似于 U-net 網絡[17]，在三維卷積層和三維反卷積層之間加入了跳躍連接，可以將在編碼器學習到的特征直接傳播到解碼器的反卷積特征層，以幫助對變形場 
的預測。
圖4
VTN 網絡結構
Figure4.
Structure of VTN
1.2.3 損失函數
無監督配準模型不需要任何數據標簽，直接以參考圖像與配準后的圖像的差異指導模型的訓練，為了以無監督的方式訓練配準模型，本文對每級子網絡產生的變形場 
計算總變差損失值 LTV 以保證產生的變形場具有平滑性，以互相關損失函數計算參考圖像和待配準圖像之間的相似度損失值 Lsim。在訓練過程中，直到所有損失值計算完成后，再將所有損失值相加進行梯度反向傳播以更新網絡的參數，從而實現對所有子網絡的端到端訓練。兩個損失函數的詳細定義如下所示：
（1）互相關系數損失函數：
互相關系數損失函數 Lsim 的作用是衡量兩幅圖像之間線性相關的程度，其定義如式（11）和式（12）所示：
 |
 |
其中，CorrCoef（·）的作用是計算兩幅圖像之間的相關系數，Cov 表示兩幅圖像之間的方差，例如，Cov［
，Sf］的定義如式（13）所示：
 |
其中，Ω 表示圖像所定義的空間。互相關系數損失函數值 Lsim 越小，表示兩幅圖像之間的線性相關程度越高，配準效果越好。
（2）總變差損失函數：
總變差損失函數的作用是使得產生的變形場 
具有平滑特性，其定義如式 (14) 所示：
 |
綜上所述，總的損失函數如式（15）所示：
 |
其中，k 代表級聯的子網絡個數，μ 代表總變差損失函數的權重。
2 實驗與結果
2.1 實驗數據
本文實驗采用的數據集包括醫學十項全能分割比賽的數據集（Medical Segmentation Decathlon，MSD）（網址為：https://decathlon-10.grand-challenge.org/）、北京友誼醫院提供的公開數據集（public datasets provided by Beijing Friend Hospital，BFD）（網址為：https://drive.google.com/open?id=17IiuM74HPj1fsWwkAfq-5Rc6r5vpxUJF）以及肝臟分割挑戰賽的公開數據集（Segmentation of Liver，Sliver）（網址為：http://www.sliver07.org/index.php），其中，MSD 包含 1 064 例 CT 圖像，BFD 包含 92 例 CT 圖像，Sliver 包含 22 例 CT 圖像，MSD 和 BFD 數據集在實驗中作為訓練集，Sliver 數據集在實驗中作為測試集。
2.2 實驗設計
原始三維肝臟 CT 圖像首先經過裁剪和去除不相關的區域后，重采樣為 128 × 128 × 128 的尺寸。紋理濾波算法中的 λ 設置為 0.08，κ 設置為 2，β0 和 βmax 分別設置為 2λ 和 105。式（15）中的 μ 取值為 1.0，級聯網絡采用亞當（Adam）優化器共訓練 105 次[18]，初始學習率為 10?4，每運行 6 × 104 次學習率減半，每運行 8 × 104 次學習率再次減半。本文共訓練了三個模型，分別為 S_VTN_2、S_VTN_4、S_VTN_6，模型名字后面的數字代表級聯的子網絡的個數。
本文所需的實驗配置如下：處理器（IntelCore i7-3770 CPU@4GHz，Intel，美國）、深度學習框架（TensorFlow, Google，美國）、獨立顯卡（Tesla v100，Nvidia，美國）。
2.2.1 比較方法
（1）傳統方法：Elastix 在傳統方法中對肝臟的配準表現最好，因此本文與其進行了比較。
（2）基于單級網絡的配準模型：Balakrishnan 等[11]提出的體素變形網絡（VoxelMorph）是基于深度學習的配準方法中的基線方法，具有較強的比較意義。
（3）基于多級級聯網絡的配準模型：為了衡量本文提出的紋理濾波算法對多級級聯網絡帶來的性能提升，本文與級聯數量分別為 2、4、6 的多級級聯網絡配準模型進行比較，這些模型的名字分別記為 VTN_2、VTN_4、VTN_6。
2.2.2 衡量指標
（1）直方圖相關系數
肝臟 CT 圖像的紋理信息主要由血管和腫瘤組織的亮度變化決定，發生較大位移或者較大形變的肝臟 CT 圖像在配準前后的灰度值分布應該成線性相關，否則，即表示配準之后肝臟的紋理信息發生了變化。基于此，本文提出直方圖相關系數以衡量在紋理信息對應上的配準效果。直方圖相關系數表示為 HistCoeff，其定義如式（16）所示：
 |
其中，X 代表肝臟區域所在空間，M（·）和 T（·）分別代表兩幅圖像在肝臟區域內的灰度值的直方圖統計。直方圖相關系數越高，代表紋理保持特性越好。
（2）Dice 系數
Dice 系數可以衡量兩個區域之間的相似性，從而反映了在空間結構對應上的配準質量。Dice 系數定義如式（17）所示：
 |
其中，A 和 B 分別代表配準圖像和參考圖像的肝臟分割掩模。
（3）特征點距離
特征點距離可以衡量空間結構對應方面的配準質量，其定義是配準圖像的解剖學特征點坐標與真實解剖學特征點坐標之間的平均歐式距離。Sliver 數據集包含由三位專家標注的解剖學特征點，可以據此計算各個特征點之間的平均距離。特征點距離越小，配準后的空間對應特性越好。
（4）折疊比例
圖像折疊使得產生的變形場不可逆，違反了解剖學意義上的對應。de Vos 等[14]提出了一種可以衡量變形場引起圖像折疊程度的方法，該方法確定了變形場中每個點 p（i，j，k）的雅可比行列式，如式（18）所示：
 |
如果式（18）的值為 1 則表示在 p 點沒有發生體積變化，大于 1 則表示膨脹，在 0～1 之間則表示收縮，小于等于 0 則表示該點為奇點，奇點會引起圖像折疊。通過計算引起圖像折疊的點在整個變形場中所占的比例，就可以計算出圖像折疊的比例，以此衡量變形場的拓撲保留特性。
2.3 結果分析與討論
2.3.1 可視化分析
基于 L0 梯度最小化的紋理濾波方法對圖像進行濾波后的效果如圖 5 所示，肝臟表面具有一塊非常明顯的病變區域，當紋理濾波權重 λ 由 0.02 增大至 0.08 時，紋理濾除程度逐漸增高，當紋理濾波權重 λ 為 0.08 時，紋理濾波后的肝臟表面已經幾乎沒有任何病變區域的紋理特征，并且較好地保留了原圖的空間結構，因此本文在后續實驗中固定參數 λ = 0.08。該結果證明本文提出的紋理濾波方法可以有效實現對三維醫學圖像的紋理濾波，滿足后續采用無監督方法進行配準的要求。
圖5
紋理濾波算法效果展示
Figure5.
The results of texture filtering algorithm
不同方法的配準效果如圖 6 所示，為了方便觀察肝臟中的腫瘤，本文選取了三維 CT 圖像中最具代表性的切片并以彩色形式展示，圖中紅色圈標注了具有腫瘤的肝臟表面，待配準圖像的肝臟表面存在一個非常明顯的腫瘤邊緣。經過配準后，VoxelMorph、VTN_4 和 VTN_6 配準后的圖像可以較好地對齊到參考圖像，但是其腫瘤的邊緣經過配準后變得模糊，而通過本文方法（S_VTN_4、S_VTN_6）配準后的圖像的腫瘤邊緣維持了待配準圖像中腫瘤邊緣的特點。
圖6
不同方法的配準結果比較
Figure6.
Comparison of registration performance of different methods
2.3.2 量化分析
本文在 Sliver 數據集上量化了基于紋理濾波的級聯網絡對配準性能產生的影響，并與傳統方法 Elastix、單級網絡 VoxelMorph 以及無紋理濾波的級聯網絡進行了對比，結果如表 1 所示。根據表 1 中的結果，傳統方法 Elastix 的配準結果在 Dice 系數上僅能達到 0.83，在特征點距離上高達 16.81，VoxelMorph 的配準結果在 Dice 系數上可以達到 0.918，在特征點距離上達到 14.19，而級聯網絡（VTN_2、VTN_4、VTN_6、S_VTN_2、S_VTN_4、S_VTN_6）在 Dice 系數上均可以達到 0.92 以上，在特征點距離上小于 13.75。因此，從配準精度上看，級聯網絡 > VoxelMorph > Elastix。隨著級聯網絡級數的增加，無紋理濾波的級聯網絡（VTN_2、VTN_4、VTN_6）的 Dice 系數雖然逐漸提升到 0.954，但是其直方圖相關系數下降至 0.959，而本文提出的基于紋理濾波的級聯網絡（S_VTN_2、S_VTN_4、S_VTN_6）在 Dice 系數上和普通的級聯網絡表現相當，并且在直方圖相關系數上其結果一直維持在 0.962 以上。值得注意的是，本文提出的基于紋理濾波的級聯網絡在反映空間對應程度的特征點距離指標上也優于無紋理濾波的級聯網絡，例如，VTN_6 的特征點距離為 11.52，而 S_VTN_6 的特征點距離為 11.47；并且本文方法也減少了配準過程中發生的折疊現象，例如，VTN_4 中的折疊比例為 7.96 × 10?3，而 S_VTN_4 中的折疊比例為 7.92 × 10?3。基于無監督學習的配準方法是直接比較參考圖像和待配準圖像之間的相似度訓練神經網絡，而神經網絡強大的擬合能力容易使得配準過度的現象發生，而且參考圖像和待配準圖像中的噪聲也會影響神經網絡預測的變形場，由于配準是一個無金標準的任務，因此基于深度學習的配準方法難以消除這些噪聲帶來的影響；本文方法通過紋理濾波使得級聯網絡僅僅關注兩幅圖像之間的空間信息，使得級聯網絡在配準精度和配準科學性上均取得優異的效果。
表1
不同方法的性能比較（
）
Table1.
Comparison of registration performance of different methods (
)
綜上所述，實驗結果證明了本文提出的基于紋理濾波的級聯網絡不僅解決了配準肝臟 CT 圖像中存在的大位移、大形變的問題，同時也解決了紋理改變的問題，并在一定程度上增加了空間對應程度，減少了折疊現象的發生。
雖然上述實驗結果可以證明本文所提出的方法在配準精度和配準科學性上均具有優勢，但是本文也存在著兩方面不足，需要在日后的研究中進一步完善。首先，本文在實驗中固定了基于 L0 梯度最小化的紋理濾波算法的參數，在實際使用過程中，由于肝臟病變區域的紋理復雜多變，仍需要根據紋理濾波的效果適當調整參數，以獲得更為可靠的結果。除此之外，本文方法的時間消耗在 8～12 s，相比于其他基于深度學習的配準算法，增加了一定的時間消耗，但是由于依然優于傳統配準算法，因此該方法適用于肝臟疾病診斷及手術規劃等實際臨床應用。
3 結論
本文提出了一種針對于肝臟 CT 圖像的基于紋理濾波的無監督配準方法。該方法首先基于 L0 梯度最小化原理實現了一種針對于三維醫學圖像的紋理濾波方法，可以有效地濾除參考圖像和待配準圖像的紋理細節，避免后續的配準過程被紋理信息所影響。然后采用級聯網絡對紋理濾波后的參考圖像和待配準圖像進行循序漸進的配準，以解決大位移和大形變問題。實驗結果表明，本文所提出的方法具有較高的配準精度，可以有效改善配準過程中發生的紋理改變的情況，并且提升了在空間結構上的對應程度以及抗圖像折疊的配準性能，對醫學配準應用在臨床和計算機輔助診斷上具有較強的指導意義。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
引言
醫學圖像配準，指在任意兩幅醫學圖像之間建立一種空間對應關系，通過這種對應關系可以使一幅圖像變換到另一幅圖像所在的空間結構中。臨床診斷中，圖像配準可以幫助醫生從不同成像設備獲取的醫學圖像或者不同時間獲取的來自相同成像設備的醫學圖像獲取信息并進行融合互補，提高疾病診斷精度[1-3]。對于計算機輔助診斷而言，醫學圖像配準可以輔助器官和腫瘤分割、多模態圖像融合以及疾病分類等研究[4-6]。
近年來，傳統配準方法例如彈性體（Elastix）、阿恩特（Ant）已經在各種任務上實現了較好的配準效果[7-8]，但是這些方法所需的配準時間較長，不適用于實時的臨床應用。最近，基于深度學習的醫學圖像配準方法呈現出巨大的優勢，不僅在配準效果上可與眾多傳統方法媲美，而且消耗的時間遠遠小于傳統配準方法。然而，基于深度學習的配準方法大多集中于腦部圖像配準[9-11]。由于腦部圖像產生的位移較小，因此這些方法遷移到位移和形變均較大的肝臟圖像中并不能取得較好的效果。一些論文提出了通過多級級聯網絡配準具有大位移和大形變醫學圖像的方法[12-14]，在肝臟電子計算機斷層掃描（computed tomography，CT）圖像中應用取得了良好的配準效果。但是通過多級級聯網絡產生的空間變換可能使得配準后的圖像與參考圖像極為相似，不僅空間結構一致，紋理方面也趨于相同。而在肝臟 CT 圖像中，紋理特點直接反映了腫瘤或病變組織的類型，如圖 1 所示，綠色箭頭標注了肝臟病灶的位置，血管瘤和肝癌的病灶具有不同的紋理信息，這些紋理信息是臨床診斷和計算機輔助診斷過程中的重要依據，因此紋理信息的改變將會嚴重影響疾病診斷的精確度。由此可見，迫切需要一種不改變紋理信息的肝臟 CT 圖像配準方法。此外，現有的配準衡量指標一般采用戴斯（Dice）系數和特征點距離等，這些指標均針對于配準后解剖結構的對應程度而提出，未有衡量配準圖像與原圖像在紋理信息上的對應程度的指標。
圖1
肝臟 CT 圖像中的紋理特征
Figure1.
Texture features in liver CT images
針對上述問題，本文提出了一種新型的三維無監督深度學習配準框架用于肝臟 CT 圖像中。首先對肝臟 CT 圖像進行基于 L0 梯度最小化的紋理濾波處理，該方法可以消除三維醫學圖像中的大部分復雜的紋理信息和圖像噪聲，保留顯著的空間結構信息；然后采用級聯網絡對輸入圖像進行循序漸進的配準，以解決大位移和大形變問題；最后提出一種直方圖相關系數衡量配準后圖像的紋理信息與原圖像的對應程度，以關注配準前后肝臟 CT 圖像紋理信息變化。實驗證明，本文所提出的方法不僅可以實現高精度配準，而且有效地保留了肝臟影像中的紋理信息。因此，本文的方法對于提升醫學圖像配準的科學性，使得醫學圖像配準安全可靠地應用在肝臟疾病的計算機輔助診斷和手術規劃方面具有重要意義。
1 方法
本文提出的配準框架如圖 2 所示。假設 If 為配準過程中的參考圖像，Im 為待配準圖像，配準后圖像為 Ir。配準的過程即將待配準圖像 Im 配準到參考圖像 If 所在的空間結構之中，從而得到配準圖像 Ir。參考圖像 If 和待配準圖像 Im 均為三維 CT 圖像。首先，Im 與 If 經過紋理濾波后分別得到 Sm 與 Sf，然后將 Sm 與 Sf 輸入至級聯網絡進行配準。由于濾波后的兩幅圖像 Sm 與 Sf 不包含紋理信息，因此級聯網絡僅按照兩幅圖像的空間結構對應進行配準，并預測出用于指導配準的變形場 。最后利用 對 Im 進行插值變換得到最終的配準圖像 Ir，該過程可被記為 Ir = Im。
圖2
基于紋理濾波的無監督級聯式配準模型
Figure2.
Unsupervised cascaded network registration model based on texture filtering
1.1 基于 L0 梯度最小化的紋理濾波
L0 梯度最小化算法最早由 Xu 等[15]提出，主要用于解決二維自然圖像的平滑問題。本文將其應用在醫學圖像中，并且在此基礎上實現了針對三維醫學圖像的紋理濾波。基于 L0 梯度最小化的紋理濾波以 L0 梯度數量為正則項，通過限制非零梯度的數量完成紋理濾波。
假設 I 和 S 分別表示輸入圖像和紋理濾波后的圖像。紋理濾波的目標函數如式（1）所示：
|
其中，∑p（Sp ? Ip）2保證了圖像 I 和 S 在每個像素點 p 的相似度，λ 為非負參數，作用是控制紋理濾波的程度，λ 越大，紋理濾波程度越高。C（S）代表圖像 S 的 L0 梯度的統計，其定義如式（2）所示:
|
其中，xSp、ySp 和 zSp 分別表示圖像 S 的任意像素點 p 沿著 x、y 和 z 方向的偏導數。#{}的功能是計算|xSp| + |ySp| + |zSp|的值不為 0 的像素點數量之和，|xSp| + |ySp| + |zSp|代表了圖像 S 在像素點 p 的梯度能量。由于 C（S）對全局不連續性進行了統計建模，故無法通過傳統的梯度下降法或其他離散優化方法優化目標函數，因此，本文采用交替最小化算法優化上述目標函數。
為了使用交替最小化算法，式（2）改寫為如式（3）所示：
|
其中，h、v、k 為引入的輔助變量，hp、vp、kp 分別表示 h、v、k 在像素點 p 的元素，C（h，v，k）= #{p||hp| + |vp| + |kp| ≠ 0}，β 是通過迭代方式更新的參數，用于控制變量（hp，vp，kp）及其對應的梯度（xSp，ySp，zSp）之間的相似度。當 β 足夠大時，式（1）與式（3）接近。
交替最小化算法將式（3）分解為兩個子問題加以解決，分別為：
（1）已知 hp、vp 和 kp，計算 S 以最小化式（3）。式（3）可以另被寫為如式（4）所示：
|
為了采用快速傅里葉變換估計式（4）中的 S，式（4）改寫為如式（5）所示：
|
其中，F（·）代表快速傅里葉變換，F?1（·）代表快速傅里葉的反變換，F*（·）表示在快速傅里葉變換基礎上取其共軛復數形式，F（1）表示階躍函數 δ（1）的傅里葉變換。
（2）已知 S，計算 hp、vp 和 kp 以最小化式（3）。式（3）可以另被寫為如式（6）所示：
|
式（6）可以在空間上分解，分別對每個 hp、vp 和 kp 進行估算。分解后的式（6）如式（7）所示：
|
其中，關于 H（|hp| + |vp| + |kp|）的計算如式（8）所示：
|
對于每一個像素點 p，式（7）可以被寫為如式（9）所示：
|
其中，Ep 代表像素點 p 的目標函數值，為了使 Ep 達到最小值，需要滿足如式（10）所示的條件：
|
因此，結合式（8）和（9），可以求解出使 Ep 達到最小值的（hp，vp，kp）。
利用兩個子問題進行交替最小化的算法框圖參見算法 1——L0 梯度最小化算法。參數 β 從最小值 β0 開始自動進行迭代，每次都乘以速率 κ，直到 β ≥ βmax。
算法 1：
STN1：輸入，圖像 I，紋理濾波權重 λ，參數 β0，βmax，κ；
步驟 2：初始化，S←I，β←β0；
步驟 3：已知 S，根據式（10）求解 hp、vp 和 kp；
步驟 4：已知 hp、vp 和 kp，根據式（5）更新 S；
步驟 5：β←κβ；
步驟 6：重復步驟 3、4、5，直至 β ≥ βmax；
步驟 7：輸出，處理后的圖像 S。
1.2 基于級聯網絡的無監督配準模型
參考圖像和待配準圖像經過紋理濾波后，輸入到基于級聯網絡的無監督配準模型中進行配準。如圖 3 所示，級聯網絡的基本框架由多個子網絡構成，通過多個子網絡，可以循序漸進地將待配準圖像配準到參考圖像所在的空間結構中。
圖3
級聯網絡框架
Figure3.
The framework of cascaded network
1.2.1 級聯網絡
由于肝臟 CT 圖像的配準任務具有大位移和大形變的特點，這里采用級聯網絡進行配準。級聯網絡由多個子網絡組成，如圖 3 所示。不失一般性，第 k 個子網絡以參考圖像 Sf 和待配準圖像 （k = 1，2，）為輸入，預測一個變形場 。對于第一個子網絡， = Sm。然后，每個子網絡預測的變形場 需要和前級子網絡預測的變形場 進行合成，形成在本級最終起作用的變形場 ，變形場合成的迭代操作記為 ，其中“”表示由 Jaderberg 等[16]提出的空間變換器網絡（spatial transformer network，STN）。最后，在每級子網絡后，待配準圖像 利用本級起作用的變形場 完成一次配準，得到下一級的待配準圖像 ，該操作記為 = 。 和 Sf 將輸入至下一個級聯的子網絡中繼續配準。通過多次配準使得 Sf 和 之間的相似度不斷增大，從而循序漸進地將 Sm 配準到 Sf 所在的空間結構中。由于 STN 是一個可微分的模塊，滿足神經網絡反向傳播梯度的需要，因此可以實現對級聯網絡的端到端訓練。
1.2.2 子網絡
級聯網絡中的子網絡采用的是卷間網絡（volume tween network，VTN）[12]，如圖 4 所示，VTN 網絡遵循編碼器—解碼器的結構，以參考圖像和待配準圖像沿著通道所在維度拼接而成的多維向量作為輸入，為了表示方便，四維的特征層在圖 4 中以二維表示，面積較小的特征層意味著更低的空間分辨率，通道數量標注在每一層的上方，彎曲的箭頭代表箭頭始端的特征層與箭頭末端的特征層沿著通道所在維度進行疊加。在編碼器部分，VTN 網絡采用帶步長的三維卷積層將每層的空間分辨率減少一半以及將每層的通道數量增加一倍，但是從第四個卷積層開始，每降低一次分辨率都將在該三維卷積層后增加一個不改變空間分辨率和通道數量的三維卷積層；在解碼器部分，VTN 網絡使用三維反卷積層（或稱為轉置卷積層）逐步恢復至原來空間分辨率和通道數量。類似于 U-net 網絡[17]，在三維卷積層和三維反卷積層之間加入了跳躍連接，可以將在編碼器學習到的特征直接傳播到解碼器的反卷積特征層，以幫助對變形場 的預測。
圖4
VTN 網絡結構
Figure4.
Structure of VTN
1.2.3 損失函數
無監督配準模型不需要任何數據標簽，直接以參考圖像與配準后的圖像的差異指導模型的訓練，為了以無監督的方式訓練配準模型，本文對每級子網絡產生的變形場 計算總變差損失值 LTV 以保證產生的變形場具有平滑性，以互相關損失函數計算參考圖像和待配準圖像之間的相似度損失值 Lsim。在訓練過程中，直到所有損失值計算完成后，再將所有損失值相加進行梯度反向傳播以更新網絡的參數，從而實現對所有子網絡的端到端訓練。兩個損失函數的詳細定義如下所示：
（1）互相關系數損失函數：
互相關系數損失函數 Lsim 的作用是衡量兩幅圖像之間線性相關的程度，其定義如式（11）和式（12）所示：
|
|
其中，CorrCoef（·）的作用是計算兩幅圖像之間的相關系數，Cov 表示兩幅圖像之間的方差，例如，Cov［，Sf］的定義如式（13）所示：
|
其中，Ω 表示圖像所定義的空間。互相關系數損失函數值 Lsim 越小，表示兩幅圖像之間的線性相關程度越高，配準效果越好。
（2）總變差損失函數：
總變差損失函數的作用是使得產生的變形場 具有平滑特性，其定義如式 (14) 所示：
|
綜上所述，總的損失函數如式（15）所示：
|
其中，k 代表級聯的子網絡個數，μ 代表總變差損失函數的權重。
2 實驗與結果
2.1 實驗數據
本文實驗采用的數據集包括醫學十項全能分割比賽的數據集（Medical Segmentation Decathlon，MSD）（網址為：https://decathlon-10.grand-challenge.org/）、北京友誼醫院提供的公開數據集（public datasets provided by Beijing Friend Hospital，BFD）（網址為：https://drive.google.com/open?id=17IiuM74HPj1fsWwkAfq-5Rc6r5vpxUJF）以及肝臟分割挑戰賽的公開數據集（Segmentation of Liver，Sliver）（網址為：http://www.sliver07.org/index.php），其中，MSD 包含 1 064 例 CT 圖像，BFD 包含 92 例 CT 圖像，Sliver 包含 22 例 CT 圖像，MSD 和 BFD 數據集在實驗中作為訓練集，Sliver 數據集在實驗中作為測試集。
2.2 實驗設計
原始三維肝臟 CT 圖像首先經過裁剪和去除不相關的區域后，重采樣為 128 × 128 × 128 的尺寸。紋理濾波算法中的 λ 設置為 0.08，κ 設置為 2，β0 和 βmax 分別設置為 2λ 和 105。式（15）中的 μ 取值為 1.0，級聯網絡采用亞當（Adam）優化器共訓練 105 次[18]，初始學習率為 10?4，每運行 6 × 104 次學習率減半，每運行 8 × 104 次學習率再次減半。本文共訓練了三個模型，分別為 S_VTN_2、S_VTN_4、S_VTN_6，模型名字后面的數字代表級聯的子網絡的個數。
本文所需的實驗配置如下：處理器（IntelCore i7-3770 CPU@4GHz，Intel，美國）、深度學習框架（TensorFlow, Google，美國）、獨立顯卡（Tesla v100，Nvidia，美國）。
2.2.1 比較方法
（1）傳統方法：Elastix 在傳統方法中對肝臟的配準表現最好，因此本文與其進行了比較。
（2）基于單級網絡的配準模型：Balakrishnan 等[11]提出的體素變形網絡（VoxelMorph）是基于深度學習的配準方法中的基線方法，具有較強的比較意義。
（3）基于多級級聯網絡的配準模型：為了衡量本文提出的紋理濾波算法對多級級聯網絡帶來的性能提升，本文與級聯數量分別為 2、4、6 的多級級聯網絡配準模型進行比較，這些模型的名字分別記為 VTN_2、VTN_4、VTN_6。
2.2.2 衡量指標
（1）直方圖相關系數
肝臟 CT 圖像的紋理信息主要由血管和腫瘤組織的亮度變化決定，發生較大位移或者較大形變的肝臟 CT 圖像在配準前后的灰度值分布應該成線性相關，否則，即表示配準之后肝臟的紋理信息發生了變化。基于此，本文提出直方圖相關系數以衡量在紋理信息對應上的配準效果。直方圖相關系數表示為 HistCoeff，其定義如式（16）所示：
|
其中，X 代表肝臟區域所在空間，M（·）和 T（·）分別代表兩幅圖像在肝臟區域內的灰度值的直方圖統計。直方圖相關系數越高，代表紋理保持特性越好。
（2）Dice 系數
Dice 系數可以衡量兩個區域之間的相似性，從而反映了在空間結構對應上的配準質量。Dice 系數定義如式（17）所示：
|
其中，A 和 B 分別代表配準圖像和參考圖像的肝臟分割掩模。
（3）特征點距離
特征點距離可以衡量空間結構對應方面的配準質量，其定義是配準圖像的解剖學特征點坐標與真實解剖學特征點坐標之間的平均歐式距離。Sliver 數據集包含由三位專家標注的解剖學特征點，可以據此計算各個特征點之間的平均距離。特征點距離越小，配準后的空間對應特性越好。
（4）折疊比例
圖像折疊使得產生的變形場不可逆，違反了解剖學意義上的對應。de Vos 等[14]提出了一種可以衡量變形場引起圖像折疊程度的方法，該方法確定了變形場中每個點 p（i，j，k）的雅可比行列式，如式（18）所示：
|
如果式（18）的值為 1 則表示在 p 點沒有發生體積變化，大于 1 則表示膨脹，在 0～1 之間則表示收縮，小于等于 0 則表示該點為奇點，奇點會引起圖像折疊。通過計算引起圖像折疊的點在整個變形場中所占的比例，就可以計算出圖像折疊的比例，以此衡量變形場的拓撲保留特性。
2.3 結果分析與討論
2.3.1 可視化分析
基于 L0 梯度最小化的紋理濾波方法對圖像進行濾波后的效果如圖 5 所示，肝臟表面具有一塊非常明顯的病變區域，當紋理濾波權重 λ 由 0.02 增大至 0.08 時，紋理濾除程度逐漸增高，當紋理濾波權重 λ 為 0.08 時，紋理濾波后的肝臟表面已經幾乎沒有任何病變區域的紋理特征，并且較好地保留了原圖的空間結構，因此本文在后續實驗中固定參數 λ = 0.08。該結果證明本文提出的紋理濾波方法可以有效實現對三維醫學圖像的紋理濾波，滿足后續采用無監督方法進行配準的要求。
圖5
紋理濾波算法效果展示
Figure5.
The results of texture filtering algorithm
不同方法的配準效果如圖 6 所示，為了方便觀察肝臟中的腫瘤，本文選取了三維 CT 圖像中最具代表性的切片并以彩色形式展示，圖中紅色圈標注了具有腫瘤的肝臟表面，待配準圖像的肝臟表面存在一個非常明顯的腫瘤邊緣。經過配準后，VoxelMorph、VTN_4 和 VTN_6 配準后的圖像可以較好地對齊到參考圖像，但是其腫瘤的邊緣經過配準后變得模糊，而通過本文方法（S_VTN_4、S_VTN_6）配準后的圖像的腫瘤邊緣維持了待配準圖像中腫瘤邊緣的特點。
圖6
不同方法的配準結果比較
Figure6.
Comparison of registration performance of different methods
2.3.2 量化分析
本文在 Sliver 數據集上量化了基于紋理濾波的級聯網絡對配準性能產生的影響，并與傳統方法 Elastix、單級網絡 VoxelMorph 以及無紋理濾波的級聯網絡進行了對比，結果如表 1 所示。根據表 1 中的結果，傳統方法 Elastix 的配準結果在 Dice 系數上僅能達到 0.83，在特征點距離上高達 16.81，VoxelMorph 的配準結果在 Dice 系數上可以達到 0.918，在特征點距離上達到 14.19，而級聯網絡（VTN_2、VTN_4、VTN_6、S_VTN_2、S_VTN_4、S_VTN_6）在 Dice 系數上均可以達到 0.92 以上，在特征點距離上小于 13.75。因此，從配準精度上看，級聯網絡 > VoxelMorph > Elastix。隨著級聯網絡級數的增加，無紋理濾波的級聯網絡（VTN_2、VTN_4、VTN_6）的 Dice 系數雖然逐漸提升到 0.954，但是其直方圖相關系數下降至 0.959，而本文提出的基于紋理濾波的級聯網絡（S_VTN_2、S_VTN_4、S_VTN_6）在 Dice 系數上和普通的級聯網絡表現相當，并且在直方圖相關系數上其結果一直維持在 0.962 以上。值得注意的是，本文提出的基于紋理濾波的級聯網絡在反映空間對應程度的特征點距離指標上也優于無紋理濾波的級聯網絡，例如，VTN_6 的特征點距離為 11.52，而 S_VTN_6 的特征點距離為 11.47；并且本文方法也減少了配準過程中發生的折疊現象，例如，VTN_4 中的折疊比例為 7.96 × 10?3，而 S_VTN_4 中的折疊比例為 7.92 × 10?3。基于無監督學習的配準方法是直接比較參考圖像和待配準圖像之間的相似度訓練神經網絡，而神經網絡強大的擬合能力容易使得配準過度的現象發生，而且參考圖像和待配準圖像中的噪聲也會影響神經網絡預測的變形場，由于配準是一個無金標準的任務，因此基于深度學習的配準方法難以消除這些噪聲帶來的影響；本文方法通過紋理濾波使得級聯網絡僅僅關注兩幅圖像之間的空間信息，使得級聯網絡在配準精度和配準科學性上均取得優異的效果。
表1
不同方法的性能比較（）
Table1.
Comparison of registration performance of different methods ()
綜上所述，實驗結果證明了本文提出的基于紋理濾波的級聯網絡不僅解決了配準肝臟 CT 圖像中存在的大位移、大形變的問題，同時也解決了紋理改變的問題，并在一定程度上增加了空間對應程度，減少了折疊現象的發生。
雖然上述實驗結果可以證明本文所提出的方法在配準精度和配準科學性上均具有優勢，但是本文也存在著兩方面不足，需要在日后的研究中進一步完善。首先，本文在實驗中固定了基于 L0 梯度最小化的紋理濾波算法的參數，在實際使用過程中，由于肝臟病變區域的紋理復雜多變，仍需要根據紋理濾波的效果適當調整參數，以獲得更為可靠的結果。除此之外，本文方法的時間消耗在 8～12 s，相比于其他基于深度學習的配準算法，增加了一定的時間消耗，但是由于依然優于傳統配準算法，因此該方法適用于肝臟疾病診斷及手術規劃等實際臨床應用。
3 結論
本文提出了一種針對于肝臟 CT 圖像的基于紋理濾波的無監督配準方法。該方法首先基于 L0 梯度最小化原理實現了一種針對于三維醫學圖像的紋理濾波方法，可以有效地濾除參考圖像和待配準圖像的紋理細節，避免后續的配準過程被紋理信息所影響。然后采用級聯網絡對紋理濾波后的參考圖像和待配準圖像進行循序漸進的配準，以解決大位移和大形變問題。實驗結果表明，本文所提出的方法具有較高的配準精度，可以有效改善配準過程中發生的紋理改變的情況，并且提升了在空間結構上的對應程度以及抗圖像折疊的配準性能，對醫學配準應用在臨床和計算機輔助診斷上具有較強的指導意義。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
