在基于深度學習的圖像配準中,圖像中具有復雜解剖結構的形變區域是影響網絡配準精度的重要因素,然而現有方法很難關注到圖像的復雜解剖區域。同時,卷積神經網絡的感受野受其卷積核大小的限制,難以學習空間位置距離較遠的體素之間的關系,使其難以處理較大區域形變問題。針對以上兩個問題,本文提出了一種基于視覺變換器(Transformer)的級聯多階層配準網絡模型,并配備了一種基于均方誤差的困難形變感知機。困難形變感知機使用滑動窗口和浮動窗口技術在配準圖像中進行檢索,得到每個體素的困難形變系數,識別出配準效果最差的區域。本研究中,級聯多階層配準網絡模型采用困難形變感知機進行階層連接,在基礎配準網絡中憑借自注意力機制提取全局特征,對不同尺度的配準結果進行優化。實驗結果證明,本文提出的方法可以對復雜形變區域進行漸進配準,從而優化腦部醫學影像的配準結果,對醫生的臨床診斷工作有良好的輔助作用。
引用本文: 潘英杰, 程遠志, 劉豪, 史操. 基于視覺變換器的級聯多階層醫學影像配準方法. 生物醫學工程學雜志, 2022, 39(5): 876-886. doi: 10.7507/1001-5515.202204011 復制
版權信息: ?四川大學華西醫院華西期刊社《生物醫學工程學雜志》版權所有,未經授權不得轉載、改編
引言
可形變醫學圖像配準作為醫學影像處理和分析的一項基礎性任務,其目標是找到參考圖像和待配準圖像間體素的位移關系,識別并對齊圖像中相同或者相似的解剖結構,精確的醫學影像配準是一項具有挑戰性的工作。
基于深度學習的可形變醫學圖像配準,按照網絡的訓練策略大體分為兩類:監督學習方法和無監督學習方法。基于監督學習的配準方法,需要在訓練網絡時提供待配準圖像對的標準位移形變場,但標準位移形變場需首先通過傳統配準方法獲得,這使監督學習網絡的配準精度難以超越傳統方法[1]。為了解決監督學習中標簽限制問題,大量研究人員開始研究無監督學習的配準方法。Jaderberg等[2]提出了一種空間變換網絡(spatial transformer networks,STN),該網絡支持神經網絡的反向傳播,直接使用形變場扭曲待配準圖像,STN的發布啟發了許多無監督圖像配準方法。Balakrishnan等[3]通過結合U型網絡[4]和STN開發了體素變形網絡(voxelmorph,VM),以無監督學習的方式對核磁共振成像(magnetic resonance imaging,MRI)腦圖譜進行配準,該方法僅通過優化自定義的損失函數就可以實現圖像配準。圖像中小形變區域通常通過單次配準即可實現對齊,而大形變區域往往需要多次配準才能對齊。在Zhao等[5]提出的遞歸級聯網絡中,級聯網絡的思想被證明可以優化圖像配準的結果,但級聯網絡的每層網絡都需要輸入參考圖像和待配準圖像,這造成已經對齊的較小形變區域再次參與網絡運算,步驟冗余沒有意義。為解決此問題,Kim等[6]提出循環體素變形網絡(cyclemorph,CM),利用循環一致性進行多尺度配準,逐步優化形變場。Huang等[7]嘗試將網絡的中間特征重采樣以識別感興趣區域,這種依賴中間特征結果的識別策略在網絡訓練早期極易出現誤判和漏判感興趣區域的問題。基于以上研究結果,本文方法采用級聯多階層網絡進行多尺度配準,通過選定的策略識別對齊最差的區域,保證識別結果的穩定性。
卷積神經網絡(convolutional neural network,CNN)被廣泛應用于圖像配準領域,但卷積運算的感受野受卷積核尺寸的限制,難以學習圖像的全局特征。Li等[8]在其研究中發現,隨著卷積層的加深,距離較遠的體素點之間的相互影響會迅速衰減,這使得CNN很難學習圖像中的全局特征關系。擁有自注意力機制的視覺變換器(Transformer)的出現有效解決了CNN無法有效提取圖像全局特征的問題[9],例如Liu等[10]提出擁有移位窗口的層次化Transformer,其計算量只與窗口數量呈線性關系,改善了計算代價高昂的問題。Chen等[11]首次應用Transformer進行無監督配準研究,與V型分割網絡結合[12],并在后續任務中擴展了當前模型,提出變形網絡(transmorph,TM)捕獲待配準圖像對之間的語義關系,在定量結果上證明了其提出的架構的有效性和先進性[13]。但將Transformer應用于醫學圖像配準的研究目前仍處于起步階段。
針對上述問題,本文提出了一種基于Transformer的級聯多階層配準網絡模型進行醫學影像配準。該模型包括:① 對原始MRI腦圖譜進行標準預處理操作。② 構建CNN提取局部特征和Transformer提取全局特征的基礎配準網絡。③ 使用困難形變感知機提取復雜形變區域,采用多階層方法級聯多個基礎配準網絡,漸進優化不同尺度的配準結果,解決圖像配準中的復雜形變問題。綜上所述,課題組期望本文提出的方法可以漸進優化配準結果,提升圖像的配準精度,今后能夠幫助醫生在臨床診斷中做出更加準確的判斷。
1 方法
1.1 總體架構
局部是整體的一部分,在計算機視覺領域,整體圖像可以看作由多個局部圖像組成,局部圖像質量的提升可以提高整體圖像的質量。對齊較好的區域受參考圖像的限制,即使再優化也很難帶來整體配準性能的顯著提升,相反,對齊較差的區域可優化空間較大,進行再優化可以顯著提高整體配準性能。基于Transformer的級聯多階層配準網絡模型采用一種基于困難形變系數(difficult deformation coefficient,DDC)的困難形變感知機篩選圖像中對齊較差的區域,以分階層的方式優化不同尺度的對齊較差區域。此模型中對齊較差的區域是相對于圖像中其他局部區域而言的,在本文中也被稱為困難形變區域,當整體配準結果較好時,多階層配準網絡依然可以選擇圖像的最困難形變區域進行優化。該模型的總體架構如圖1所示。
 圖1
				級聯多階層配準網絡模型
			
												
				Figure1.
				Cascaded multi-level registration network model
						
				圖1
				級聯多階層配準網絡模型
			
												
				Figure1.
				Cascaded multi-level registration network model
			
								級聯多階層配準網絡模型擁有三個階層子網絡和一個形變場融合模塊。每個階層子網絡的基礎配準網絡的結構相同,輸入為待配準圖像對,輸出為形變場。在前兩個階層子網絡中,基礎配準網絡生成的形變場首先會輸入空間轉換網絡生成配準后的圖像,隨該圖像和參考圖像一起輸入困難形變感知機篩選出困難形變區域,該區域將輸入下一階層子網絡進行精細配準。最終,多個階層子網絡生成的不同尺度的形變場經嵌入融合后生成最終的形變場φ。級聯多階層配準網絡模型的初始輸入是圖像尺度為 的待配準圖像對,按照整體到局部的原則,第一階層子網絡應篩選大尺寸的困難形變區域,其圖像尺度設置為
的待配準圖像對,按照整體到局部的原則,第一階層子網絡應篩選大尺寸的困難形變區域,其圖像尺度設置為 ,第二階層子網絡應篩選小尺寸的困難形變區域,便于第三階層子網絡對圖像細節進行調整,所以該階層篩選的圖像尺度設置為
,第二階層子網絡應篩選小尺寸的困難形變區域,便于第三階層子網絡對圖像細節進行調整,所以該階層篩選的圖像尺度設置為 。
。
在形變場融合模塊,多個階層子網絡生成的形變場根據困難形變感知機記錄的空間位置進行嵌入融合。為了將全局和局部形變場合成為一個形變場,嵌入融合使用組合的方法,首先從所有生成的形變場中選取尺度較小的兩個形變場,兩個形變場中尺度大的形變場作為本次計算的全局形變場,再使用復合惡魔(Demons)算法扭曲全局和局部形變場得到計算結果,并根據記錄的空間位置將該結果和局部形變場相加[14]得到計算結果,該結果將和其他形變場以迭代的方式執行上述操作得到φ。
1.2 基礎配準網絡
Transformer可以有效提取全局信息,但其計算成本高昂,如果其直接應用于高分辨率的三維(three-dimensional,3D)醫學圖像,會存在顯存爆炸問題導致模型無法訓練,目前主流的方法是使用重采樣減小圖像尺寸從而降低模型整體計算成本。雖然重采樣方法降低了模型計算成本,但這一方法存在局部信息丟失問題,而CNN擁有很強的局部信息提取能力并且計算成本較低,可以應用于高分辨率圖像,很好地彌補了Transformer的缺陷。基礎配準網絡將CNN應用于高分辨率特征圖提取圖像的局部特征,Transformer應用于低分辨率特征圖提取圖像的全局特征,并在形變場生成前,額外添加參考圖像的特征,增強網絡對圖像間差異的理解能力。
基礎配準網絡分為三個模塊:CNN編碼器-解碼器、Transformer編碼器-解碼器和特征指導模塊,如圖2所示。
 圖2
				基礎配準網絡
			
												
				Figure2.
				Basic registration network
						
				圖2
				基礎配準網絡
			
												
				Figure2.
				Basic registration network
			
								1.2.1 CNN編碼器-解碼器
在高分辨率特征處理階段,卷積操作可以更好地學習圖像的局部特征關系。CNN編碼器-解碼器左側為編碼器,右側為解碼器,假設輸入圖像尺寸為 ,編碼器和解碼器均使用4個卷積核尺寸為3的3D卷積層提取圖像局部特征,改變特征圖尺寸和通道數。解碼器中經上采樣的特征圖通過跳躍連接與編碼器輸出的特征圖級聯,跳躍連接會跳躍神經網絡中的某些層,并將一層的輸出作為下一層的輸入,能夠幫助網絡獲得更加精細的細節特征,同時可以解決模型訓練過程中的梯度消失問題。每個卷積層之后都有一個整流線性單元(rectified linear units,ReLU)激活函數,ReLU 是一個分段線性函數,能夠幫助網絡更好地挖掘特征關系,擬合訓練數據。
,編碼器和解碼器均使用4個卷積核尺寸為3的3D卷積層提取圖像局部特征,改變特征圖尺寸和通道數。解碼器中經上采樣的特征圖通過跳躍連接與編碼器輸出的特征圖級聯,跳躍連接會跳躍神經網絡中的某些層,并將一層的輸出作為下一層的輸入,能夠幫助網絡獲得更加精細的細節特征,同時可以解決模型訓練過程中的梯度消失問題。每個卷積層之后都有一個整流線性單元(rectified linear units,ReLU)激活函數,ReLU 是一個分段線性函數,能夠幫助網絡更好地挖掘特征關系,擬合訓練數據。
1.2.2 Transformer編碼器-解碼器
在低分辨率特征處理階段,網絡使用Transformer學習網絡學習圖像的全局特征關系。Transformer編碼器-解碼器首先使用一個卷積層將特征圖尺寸變為 (其中C為設置的通道數),隨后該特征圖輸入由Transformer模塊、圖像塊合并層、轉置卷積層串聯的V型結構,并將相同分辨率級別的特征圖跳躍連接。Transformer模塊不改變特征圖的尺寸,使用基于體積的多頭自注意力機制(volume-based multi-headedself-attention,V-MSA)進行特征提取,V-MSA通過計算查詢-鍵-值向量學習輸入序列間的相關性 [15-16]。圖像塊合并層用于圖像下采樣操作,它連接8個相鄰位置的圖像塊向量,并應用一個線性層產生2倍通道數的輸出,經過上述操作特征圖尺寸減半,通道數翻倍[10]。轉置卷積層的步長為2,用于將低分辨率特征圖上采樣映射成高分辨率特征圖,相比于向上插值采樣 + 卷積,轉置卷積擁有更高的執行效率[17]。
(其中C為設置的通道數),隨后該特征圖輸入由Transformer模塊、圖像塊合并層、轉置卷積層串聯的V型結構,并將相同分辨率級別的特征圖跳躍連接。Transformer模塊不改變特征圖的尺寸,使用基于體積的多頭自注意力機制(volume-based multi-headedself-attention,V-MSA)進行特征提取,V-MSA通過計算查詢-鍵-值向量學習輸入序列間的相關性 [15-16]。圖像塊合并層用于圖像下采樣操作,它連接8個相鄰位置的圖像塊向量,并應用一個線性層產生2倍通道數的輸出,經過上述操作特征圖尺寸減半,通道數翻倍[10]。轉置卷積層的步長為2,用于將低分辨率特征圖上采樣映射成高分辨率特征圖,相比于向上插值采樣 + 卷積,轉置卷積擁有更高的執行效率[17]。
1.2.3 特征指導模塊
在特征指導模塊部分,參考圖像的特征可以更好地指導復雜網絡學習圖像映射關系,首先使用一個步長為1的轉置卷積層用來消除轉置卷積核大小為奇數可能帶來的網格棋盤效應,并用一個卷積層單獨提取參考圖像特征,隨后使用一個卷積層融合來自參考圖像和CNN解碼器的特征圖,最后使用一個卷積層將特征映射為形變場。
1.3 困難形變感知機
圖像中困難形變區域是影響整體配準結果的主要原因,為了篩選圖像中困難形變區域,本文提出了一種基于DDC的困難形變感知機,將其作為多階層方法的橋接裝置。由于無法通過形變場直接衡量困難形變程度,該模塊采用滑動窗口和浮動窗口掃描圖像,計算并比較DDC的大小以篩選困難形變區域。圖像中相似度較低的區域往往也是對齊效果較差的區域,DDC的計算基于局部圖像塊的均方誤差(mean square error,MSE)評分[18],該系數[式(1)中,以符號DDC(·)表示]衡量了圖像中所有體素點的平均困難形變程度,其數值越大,則感知機認為該區域圖像差異越大、越難形變對齊。具有較大解剖差異的圖像塊很難對齊,應盡量選取器官結構信息豐富的圖像塊,為了使感知機聚焦有更多前景體素的有價值區域,滑動窗口會統計當前窗口的前景體素數,并除以當前窗口中的體素總數得到前景比重,用常量加上該比重作為MSE評分[式(2)中,以符號MSE(·)表示]的權重系數,其定義如式(1)和式(2)所示:
|  | 
|  | 
其中, 代表FL和ML所有體素點的DDC,
代表FL和ML所有體素點的DDC, 代表FL和ML所有體素點的MSE評分,FL為參考圖像的局部圖像塊,ML為已配準圖像的局部圖像塊,N為圖像體素總數,NF為前景體素數,p為體素位置,
代表FL和ML所有體素點的MSE評分,FL為參考圖像的局部圖像塊,ML為已配準圖像的局部圖像塊,N為圖像體素總數,NF為前景體素數,p為體素位置, 為圖像域。
為圖像域。
感知機會自動記錄已掃描過圖像塊的DCC最大值以及其在圖像中的位置信息,當窗口掃描到其他圖像塊時,會比較其系數大小,若當前圖像塊的DDC大于已記錄的最大值,則會更新該最大值和窗口位置信息,最終選取該系數值最大的區域作當前圖像的最困難形變區域輸入下一階層子網絡。窗口的尺寸是預先設定的,當圖像中存在多個不同區域大小的復雜形變區域時,相較于大區域,小區域對窗口DDC的影響較弱,同時,DDC代表了局部區域的困難形變程度,感知機會根據該系數選擇當前圖像中的困難形變程度最大的區域。若使用步長為1滑動窗口遍歷整幅圖像,其計算代價較為高昂,為了減小計算量,感知機采用滑動窗口和浮動窗口結合的方式,分兩個階段對圖像進行掃描檢查。困難形變感知機工作原理如圖3所示。
 圖3
				困難形變感知機工作原理圖
			
												
				Figure3.
				The working principle of the difficult deformation perceptron
						
				圖3
				困難形變感知機工作原理圖
			
												
				Figure3.
				The working principle of the difficult deformation perceptron
			
								第一階段為滑動遍歷階段,該階段目標是選取配準對齊較差的區域,通過一個3D滑動窗口遍歷整幅圖并計算DDC,遍歷完成后,記錄第一階段該系數的最大值以及圖像塊的位置信息。本階段滑動窗口的步長設置為5,圖3中藍色實線立方體表示本階段遍歷檢測的滑動窗口示例,紅色實線立方體表示本階段的選中窗口。
第二階段為空間浮動階段,該階段目標是檢查記錄圖像塊周圍區域的困難形變程度,得到最困難形變區域及其位置信息。在第一階段記錄的空間位置處,浮動窗口以滑動遍歷階段記錄的圖像塊為基準進行浮動檢查,沿每個坐標軸有3種圖像塊偏移選擇,即可選擇沿坐標軸正方向或負方向分別偏移浮動步長的體素點位,或者選擇不偏移,3D坐標系排列共計27個方向(包括三個坐標系均不偏移的情況,圖3中僅展示4個方向)的空間浮動,空間浮動完成后,感知機會更新DDC的最大值以及圖像塊的位置信息。之后減小步長并重復上述操作,兩次浮動的步長依次設置為3和1。圖3中的深紅色實線立方體為本階段選中的最終窗口,橙色、藍色、綠色虛線立方體分別表示僅沿Y軸正方向、僅沿X軸負方向、僅沿Z軸正方向偏移浮動步長的浮動窗口示例。
在本文的方法中,滑動窗口和浮動窗口的尺寸是可變的,但由于Tansformer配準網絡在處理圖像時須最小采樣到原圖像尺寸的 ,計算過程需要線性投影對圖像進行分塊處理[9],這要求篩選的圖像尺寸必須是32的倍數。按照多階層配準網絡模型的結構設計,第一階層感知機的窗口尺寸設置為
,計算過程需要線性投影對圖像進行分塊處理[9],這要求篩選的圖像尺寸必須是32的倍數。按照多階層配準網絡模型的結構設計,第一階層感知機的窗口尺寸設置為 ,第二階層感知機的窗口尺寸設置為
,第二階層感知機的窗口尺寸設置為 。
。
1.4 空間變換網絡
為了獲取經過形變場變形的圖像,采用STN對待配準圖像進行空間變換,STN是一個支持反向傳播的可微分模塊,它根據基礎配準網絡提供的形變場φ對輸入圖像進行空間扭曲。對于待配準圖像M的每一個體素p,其扭曲后的體素p'計算公式為 ,u(p)為p的空間位移。p'的大小不一定是整數,但圖像中的體素必須在整數位置處定義,為解決此問題,本文使用適合3D圖像的三線性插值法對體素點進行重新計算。進行體素變換的公式如式(3)所示:
,u(p)為p的空間位移。p'的大小不一定是整數,但圖像中的體素必須在整數位置處定義,為解決此問題,本文使用適合3D圖像的三線性插值法對體素點進行重新計算。進行體素變換的公式如式(3)所示:
|  '/> | 
其中,q表示選定的體素點, 表示p'相鄰位置體素,d為選定的迭代方向,x、y、z分別表示3D坐標系的三個迭代方向,I(·)表示已選定的三線性插值方法。
表示p'相鄰位置體素,d為選定的迭代方向,x、y、z分別表示3D坐標系的三個迭代方向,I(·)表示已選定的三線性插值方法。
1.5 損失函數
無監督學習的方法無需在損失函數中提供標簽信息,在訓練過程中直接通過最小化損失函數即可得到最優形變場。本模型的損失函數由兩部分組成:一部分使用圖像相似損失項,減小圖像間差異;另一部分使用擴散正則化對形變場進行平滑性約束[4]。模型的總體損失函數如式(4)所示:
|  | 
其中, 為模型的總體損失函數,F表示參考圖像,M表示待配準的圖像,
為模型的總體損失函數,F表示參考圖像,M表示待配準的圖像, 使用MSE度量圖像的線性相似度,λ1是為該項的權重系數,其大小通常設置為1。
使用MSE度量圖像的線性相似度,λ1是為該項的權重系數,其大小通常設置為1。 是保持φ平滑的正則化項,λ2為該項的權重系數。總體損失函數優化的主要目標為減小圖像間差異,正則化項對總體損失函數的貢獻應顯著小于圖像相似損失項。對于腦部MRI配準,Balakrishnan等[3]的研究中報告該權重系數設置為0.02是最佳值,因此本文中正則化權重系數λ2大小設置為0.02。
是保持φ平滑的正則化項,λ2為該項的權重系數。總體損失函數優化的主要目標為減小圖像間差異,正則化項對總體損失函數的貢獻應顯著小于圖像相似損失項。對于腦部MRI配準,Balakrishnan等[3]的研究中報告該權重系數設置為0.02是最佳值,因此本文中正則化權重系數λ2大小設置為0.02。 的定義如式(5)所示:
的定義如式(5)所示:
|  | 
由于最小化 可能會導致網絡生成一個不平滑的φ,所以本文使用擴散正則化平滑φ。
可能會導致網絡生成一個不平滑的φ,所以本文使用擴散正則化平滑φ。 的定義如式(6)所示:
的定義如式(6)所示:
|  | 
其中,Δu(p)表示φ中p的空間位移的前向差分,差分結果反映了離散位移值之間的變化,形變場中可能會存異常值,相比于后向差分,前向差分更能鼓勵某一點的位移值相似于其相鄰位置的值,Δu(p)的定義如式(7)所示:
|  | 
2 實驗
2.1 數據集
本文所使用數據集,包括:神經影像學實驗室的洛尼概率腦圖譜(Loniprobabilistic brain atlas 40,LPBA40),該圖譜來自公開數據庫阿爾茨海默病神經影像學倡議(Alzheimer’s disease neuroimaging initiative,ADNI)(網址為:
使用MRI圖像處理開源軟件Freesurfer(v7.1.0,MIT Health HST,美國)對數據集進行顱骨去除、仿射對齊、圖像歸一化、圖像裁剪等標準預處理操作[21],圖像裁剪后大小為 。LPBA40數據集通過不重復的兩兩配對生成760個待配準圖像組,OASIS-1數據集隨機配對生成600個待配準圖像組。兩個數據集的訓練集、驗證集和測試集均按照7∶1∶2的比例分配。
。LPBA40數據集通過不重復的兩兩配對生成760個待配準圖像組,OASIS-1數據集隨機配對生成600個待配準圖像組。兩個數據集的訓練集、驗證集和測試集均按照7∶1∶2的比例分配。
2.2 實驗配置
對比方法采用對稱歸一化(symmetric normalization,SYN)[22]、CM、VM和TM。在14種典型非線性形變算法中,SYN是性能最好的配準算法之一[23],同時它是一種基于迭代的傳統配準算法,本文使用高級歸一化工具軟件包實現[24],實驗中每個級別的最大迭代次數分別設置為160、80、40。CM是使用循環一致性的多尺度配準模型,超參數α、β和λ分別對應其損失函數中循環損失、身份損失和形變場正則化的權重,本文使用Kim等[6]報告的最佳值,超參數α、β和λ分別設置為0.1、0.5和1。VM是一種目前較為流行的無監督可形變配準方法,使用CNN預測形變場,本文使用該方法的2號變體進行對比實驗[4]。TM是在Transformer基礎上開發的一種無監督配準網絡,對于腦部MRI配準,本文使用其推薦的默認參數設置。
級聯多階層配準網絡模型的神經元丟棄率設置為0.3,通道數C設置為64,訓練的批尺寸設置為1,初始學習率為1 × 10?4,學習率衰減策略采用衰減率為0.9的指數衰減策略,訓練迭代次數設置為250,優化器為自適應矩估計(adaptive moment estimation,Adam)優化器。本文采用數據集提供的分割圖來評估各個配準方法的效果。
本文方法訓練模型的實驗設備配置如下:深度學習框架(PyTorch 1.8.0,FAIR,美國) [25],獨立顯卡(Nvidia Geforce RTX 3090,Nvidia,美國),處理器(IntelCore i9-10900K CPU@3.70GHz,Intel,美國)。
2.3 評估指標
本文采用戴斯相似系數(Dice similarity coefficient,DSC)計算分割圖的相似度[26],DSC量化了分割圖中解剖結構之間的體積重疊程度,DSC值越高,說明圖像重疊程度越高,配準性能越好。DSC[式(8)中,以符號DSC(·)表示]的定義如式(8)所示:
|  | 
其中, 表示SF和SM的DSC,SF表示參考分割圖,SM表示配準后的分割圖,SF∩SM表示兩分割圖樣本的交集。
表示SF和SM的DSC,SF表示參考分割圖,SM表示配準后的分割圖,SF∩SM表示兩分割圖樣本的交集。
同時,使用負雅可比行列式(negative Jacobian determinant,NJD)對預測的形變場進行平滑度評價[27],NJD系數[式(9)中,以符號NJD(·)表示]越接近0時,形變場的平均折疊程度越小,形變場越平滑。對于φ中每個p,NJD系數的定義如式(9)所示:
|  | 
其中, 表示φ中p處的NJD系數, i、j、z分別表示3D坐標系中p在
表示φ中p處的NJD系數, i、j、z分別表示3D坐標系中p在 、
、 、
、 三個方向的坐標。
三個方向的坐標。
在圖像的相似性評估上,使用三個指標對圖像的相似性進行評價,這些指標分別是均方誤差平方根(root mean square error,RMSE)[18]、結構相似性(structural similarity,SSIM)[28]、互信息(mutual information,MI)[29]。若圖像的RMSE越高、SSIM和MI越低,代表圖像相似性越差。對于F和M,RMSE[式(10)中,以符號RMSE(·)表示]的定義如式(10)所示:
|  | 
其中, 表示F和M的RMSE值,用來衡量圖像的標準誤差。
表示F和M的RMSE值,用來衡量圖像的標準誤差。
SSIM[式(11)中,以符號SSIM(·)表示]的具體計算過程如式(11)~式(14)所示:
|  | 
|  | 
|  | 
|  | 
其中, 主要關注F和M的邊緣和紋理的結構相似性,以衡量圖像的相似程度,
主要關注F和M的邊緣和紋理的結構相似性,以衡量圖像的相似程度, 、
、 、
、 分別代表亮度比較、對比度比較和結構比較。μF和μM代表F和M的平均值,σF和σM代表F、M的標準差,σFM代表F與M的協方差,c1、c2、c3為非零常數,避免分母為零帶來的系統錯誤,a、b、c均設置為1。
分別代表亮度比較、對比度比較和結構比較。μF和μM代表F和M的平均值,σF和σM代表F、M的標準差,σFM代表F與M的協方差,c1、c2、c3為非零常數,避免分母為零帶來的系統錯誤,a、b、c均設置為1。
MI[式(15)中,以符號MI(·)表示]的定義如式(15)所示:
|  | 
其中, 表示圖像F和M相互包含的信息量,圖像匹配程度越高該值越大,
表示圖像F和M相互包含的信息量,圖像匹配程度越高該值越大, 為F和M的聯合概率密度函數,
為F和M的聯合概率密度函數, 和
和 為分別表示F和M的邊緣概率密度函數。
為分別表示F和M的邊緣概率密度函數。
2.4 結果分析
2.4.1 配準精度
本研究量化了實驗中所有方法在兩個數據集上的平均DSC和NJD系數,結果如表1所示。根據表1結果,本文方法在兩個數據集上的平均DSC均高于對比方法,提供了更好的圖像配準質量,同時生成形變場NJD系數在LPBA40數據集上低于對比方法,在OASIS-1數據集上略高于對比方法。如圖4所示,舉例說明了所有方法的可視化配準結果,包括圖像的細節放大圖、分割圖(包括四個解剖結構皮質、灰質、白質和腦脊液的分割結果)和分割細節,本文使用的方法生成的扭曲圖像和分割圖在外觀上更加接近于原始參考圖像及其分割圖。解剖結構的箱線圖計算了所有方法在OASIS-1數據集上的各個解剖標簽的平均DSC,其中左腦和右腦中的相同解剖結構被合并成一個結構進行計算,具體結果如圖5所示,縱坐標為DSC,橫坐標為解剖結構類別。統計結果顯示,本文提出的方法取得了最好的配準效果。
 表1
                不同方法的配準結果比較
		 	
		 			 				Table1.
    			Comparison of registration results of different methods
			
						表1
                不同方法的配準結果比較
		 	
		 			 				Table1.
    			Comparison of registration results of different methods
       		
       				 圖4
				不同方法的配準結果的可視化示例
			
												
				Figure4.
				Visual examples of registration results for different methods
						
				圖4
				不同方法的配準結果的可視化示例
			
												
				Figure4.
				Visual examples of registration results for different methods
			
								 圖5
				解剖結構分組繪制的箱線圖
			
												
				Figure5.
				Boxplots of anatomical structures grouped
						
				圖5
				解剖結構分組繪制的箱線圖
			
												
				Figure5.
				Boxplots of anatomical structures grouped
			
								2.4.2 階層對比
本文方法設置的階層數為3,理論上,階層數的設置是任意的,為了驗證多階層的方法可以實現對形變場的逐步優化,本研究嘗試使用更多階層進行實驗,但由于本文的基礎配準網絡對圖像尺寸存在限制(必須為32的倍數),同時第一階層子網絡已經篩選了大尺寸的困難形變區域,導致本研究僅能在第二階層和第三階層子網絡中間增加圖像尺度為 的階層子網絡,因此,本研究使用4階層子網絡進行實驗,對訓練好的最佳預訓練模型進行微調訓練,對比了使用不同階層數子網絡生成的形變場扭曲分割圖的平均DSC,如表2所示,可以發現隨著階層的增多,平均DSC在不斷提高,證明多階層方法可以實現對圖像的漸進優化配準,提升配準效果。如圖6所示,對不同階層子網絡關注的困難形變區域及該區域的配準結果和形變場進行可視化展示,實線框出的區域代表當前階層子網絡篩選的困難形變區域,不同階層由不同顏色的實線框進行標注。
的階層子網絡,因此,本研究使用4階層子網絡進行實驗,對訓練好的最佳預訓練模型進行微調訓練,對比了使用不同階層數子網絡生成的形變場扭曲分割圖的平均DSC,如表2所示,可以發現隨著階層的增多,平均DSC在不斷提高,證明多階層方法可以實現對圖像的漸進優化配準,提升配準效果。如圖6所示,對不同階層子網絡關注的困難形變區域及該區域的配準結果和形變場進行可視化展示,實線框出的區域代表當前階層子網絡篩選的困難形變區域,不同階層由不同顏色的實線框進行標注。
 表2
                不同階層數形變場的配準結果比較
		 	
		 			 				Table2.
    			Comparison of registration results using different levels of deformation fields
			
						表2
                不同階層數形變場的配準結果比較
		 	
		 			 				Table2.
    			Comparison of registration results using different levels of deformation fields
       		
       				 圖6
				不同階層關注的感興趣區域對比示例圖
			
												
				Figure6.
				Comparative examples of regions of interest concerned by different levels
						
				圖6
				不同階層關注的感興趣區域對比示例圖
			
												
				Figure6.
				Comparative examples of regions of interest concerned by different levels
			
								2.4.3 感興趣區域識別
困難形變感知機的識別結果直接影響了后續階層能否關注到圖像中有價值信息的區域,但形變場的評價指標沒有金標準,無法直接判斷所選區域的復雜性,因此本研究使用OASIS-1測試集進行實驗,對第一階層子網絡中感知機識別前后的圖像進行相似度評估,進而判斷感知機能否識別圖像的困難形變區域,對比結果如表3所示,其中全局圖像代表已配準圖像對,局部圖像代表經過識別后困難形變區域。參數結果顯示,局部圖像的結構相似度更差,互相包含的信息量低于全局圖像,證明該感知機可以識別困難形變區域。
 表3
                圖像相似性對比結果
		 	
		 			 				Table3.
    			Image similarity comparison results
			
						表3
                圖像相似性對比結果
		 	
		 			 				Table3.
    			Image similarity comparison results
       		
       				3 結論
本文介紹了一種基于Transformer的級聯多階層醫學影像配準方法。該方法首先基于Transformer構建了基礎配準網絡,有效地融合圖像的全局信息和局部信息,同時配備了一種基于DDC的困難形變感知機,可以識別并裁剪圖像的感興趣區域。然后以困難形變感知機為橋接方法,級聯了三個包含基礎配準網絡的階層子網絡,并且采用嵌入融合的方式生成形變場,構建了基于Transformer的級聯多階層配準網絡模型。實驗結果表明,與當前流行的無監督配準方法相比,本文方法可以漸進地優化配準結果,解決圖像的復雜形變問題,提升配準精度,幫助醫生做出更準確的臨床診斷,促進計算機輔助醫療技術的發展。但目前本文的方法只適用于單模態圖像配準,在后續工作中會考慮將該方法擴展到多模態圖像配準。
重要聲明
利益沖突聲明:本文全體作者均聲明不存在利益沖突。
作者貢獻聲明:潘英杰負責本研究的算法設計、程序編寫、實驗數據處理和論文撰寫;程遠志是本研究的負責人,指導實驗設計和論文寫作,提出修改意見;劉豪和史操負責論文修改和數據分析。
引言
可形變醫學圖像配準作為醫學影像處理和分析的一項基礎性任務,其目標是找到參考圖像和待配準圖像間體素的位移關系,識別并對齊圖像中相同或者相似的解剖結構,精確的醫學影像配準是一項具有挑戰性的工作。
基于深度學習的可形變醫學圖像配準,按照網絡的訓練策略大體分為兩類:監督學習方法和無監督學習方法。基于監督學習的配準方法,需要在訓練網絡時提供待配準圖像對的標準位移形變場,但標準位移形變場需首先通過傳統配準方法獲得,這使監督學習網絡的配準精度難以超越傳統方法[1]。為了解決監督學習中標簽限制問題,大量研究人員開始研究無監督學習的配準方法。Jaderberg等[2]提出了一種空間變換網絡(spatial transformer networks,STN),該網絡支持神經網絡的反向傳播,直接使用形變場扭曲待配準圖像,STN的發布啟發了許多無監督圖像配準方法。Balakrishnan等[3]通過結合U型網絡[4]和STN開發了體素變形網絡(voxelmorph,VM),以無監督學習的方式對核磁共振成像(magnetic resonance imaging,MRI)腦圖譜進行配準,該方法僅通過優化自定義的損失函數就可以實現圖像配準。圖像中小形變區域通常通過單次配準即可實現對齊,而大形變區域往往需要多次配準才能對齊。在Zhao等[5]提出的遞歸級聯網絡中,級聯網絡的思想被證明可以優化圖像配準的結果,但級聯網絡的每層網絡都需要輸入參考圖像和待配準圖像,這造成已經對齊的較小形變區域再次參與網絡運算,步驟冗余沒有意義。為解決此問題,Kim等[6]提出循環體素變形網絡(cyclemorph,CM),利用循環一致性進行多尺度配準,逐步優化形變場。Huang等[7]嘗試將網絡的中間特征重采樣以識別感興趣區域,這種依賴中間特征結果的識別策略在網絡訓練早期極易出現誤判和漏判感興趣區域的問題。基于以上研究結果,本文方法采用級聯多階層網絡進行多尺度配準,通過選定的策略識別對齊最差的區域,保證識別結果的穩定性。
卷積神經網絡(convolutional neural network,CNN)被廣泛應用于圖像配準領域,但卷積運算的感受野受卷積核尺寸的限制,難以學習圖像的全局特征。Li等[8]在其研究中發現,隨著卷積層的加深,距離較遠的體素點之間的相互影響會迅速衰減,這使得CNN很難學習圖像中的全局特征關系。擁有自注意力機制的視覺變換器(Transformer)的出現有效解決了CNN無法有效提取圖像全局特征的問題[9],例如Liu等[10]提出擁有移位窗口的層次化Transformer,其計算量只與窗口數量呈線性關系,改善了計算代價高昂的問題。Chen等[11]首次應用Transformer進行無監督配準研究,與V型分割網絡結合[12],并在后續任務中擴展了當前模型,提出變形網絡(transmorph,TM)捕獲待配準圖像對之間的語義關系,在定量結果上證明了其提出的架構的有效性和先進性[13]。但將Transformer應用于醫學圖像配準的研究目前仍處于起步階段。
針對上述問題,本文提出了一種基于Transformer的級聯多階層配準網絡模型進行醫學影像配準。該模型包括:① 對原始MRI腦圖譜進行標準預處理操作。② 構建CNN提取局部特征和Transformer提取全局特征的基礎配準網絡。③ 使用困難形變感知機提取復雜形變區域,采用多階層方法級聯多個基礎配準網絡,漸進優化不同尺度的配準結果,解決圖像配準中的復雜形變問題。綜上所述,課題組期望本文提出的方法可以漸進優化配準結果,提升圖像的配準精度,今后能夠幫助醫生在臨床診斷中做出更加準確的判斷。
1 方法
1.1 總體架構
局部是整體的一部分,在計算機視覺領域,整體圖像可以看作由多個局部圖像組成,局部圖像質量的提升可以提高整體圖像的質量。對齊較好的區域受參考圖像的限制,即使再優化也很難帶來整體配準性能的顯著提升,相反,對齊較差的區域可優化空間較大,進行再優化可以顯著提高整體配準性能。基于Transformer的級聯多階層配準網絡模型采用一種基于困難形變系數(difficult deformation coefficient,DDC)的困難形變感知機篩選圖像中對齊較差的區域,以分階層的方式優化不同尺度的對齊較差區域。此模型中對齊較差的區域是相對于圖像中其他局部區域而言的,在本文中也被稱為困難形變區域,當整體配準結果較好時,多階層配準網絡依然可以選擇圖像的最困難形變區域進行優化。該模型的總體架構如圖1所示。
 圖1
				級聯多階層配準網絡模型
			
												
				Figure1.
				Cascaded multi-level registration network model
						
				圖1
				級聯多階層配準網絡模型
			
												
				Figure1.
				Cascaded multi-level registration network model
			
								級聯多階層配準網絡模型擁有三個階層子網絡和一個形變場融合模塊。每個階層子網絡的基礎配準網絡的結構相同,輸入為待配準圖像對,輸出為形變場。在前兩個階層子網絡中,基礎配準網絡生成的形變場首先會輸入空間轉換網絡生成配準后的圖像,隨該圖像和參考圖像一起輸入困難形變感知機篩選出困難形變區域,該區域將輸入下一階層子網絡進行精細配準。最終,多個階層子網絡生成的不同尺度的形變場經嵌入融合后生成最終的形變場φ。級聯多階層配準網絡模型的初始輸入是圖像尺度為 的待配準圖像對,按照整體到局部的原則,第一階層子網絡應篩選大尺寸的困難形變區域,其圖像尺度設置為
的待配準圖像對,按照整體到局部的原則,第一階層子網絡應篩選大尺寸的困難形變區域,其圖像尺度設置為 ,第二階層子網絡應篩選小尺寸的困難形變區域,便于第三階層子網絡對圖像細節進行調整,所以該階層篩選的圖像尺度設置為
,第二階層子網絡應篩選小尺寸的困難形變區域,便于第三階層子網絡對圖像細節進行調整,所以該階層篩選的圖像尺度設置為 。
。
在形變場融合模塊,多個階層子網絡生成的形變場根據困難形變感知機記錄的空間位置進行嵌入融合。為了將全局和局部形變場合成為一個形變場,嵌入融合使用組合的方法,首先從所有生成的形變場中選取尺度較小的兩個形變場,兩個形變場中尺度大的形變場作為本次計算的全局形變場,再使用復合惡魔(Demons)算法扭曲全局和局部形變場得到計算結果,并根據記錄的空間位置將該結果和局部形變場相加[14]得到計算結果,該結果將和其他形變場以迭代的方式執行上述操作得到φ。
1.2 基礎配準網絡
Transformer可以有效提取全局信息,但其計算成本高昂,如果其直接應用于高分辨率的三維(three-dimensional,3D)醫學圖像,會存在顯存爆炸問題導致模型無法訓練,目前主流的方法是使用重采樣減小圖像尺寸從而降低模型整體計算成本。雖然重采樣方法降低了模型計算成本,但這一方法存在局部信息丟失問題,而CNN擁有很強的局部信息提取能力并且計算成本較低,可以應用于高分辨率圖像,很好地彌補了Transformer的缺陷。基礎配準網絡將CNN應用于高分辨率特征圖提取圖像的局部特征,Transformer應用于低分辨率特征圖提取圖像的全局特征,并在形變場生成前,額外添加參考圖像的特征,增強網絡對圖像間差異的理解能力。
基礎配準網絡分為三個模塊:CNN編碼器-解碼器、Transformer編碼器-解碼器和特征指導模塊,如圖2所示。
 圖2
				基礎配準網絡
			
												
				Figure2.
				Basic registration network
						
				圖2
				基礎配準網絡
			
												
				Figure2.
				Basic registration network
			
								1.2.1 CNN編碼器-解碼器
在高分辨率特征處理階段,卷積操作可以更好地學習圖像的局部特征關系。CNN編碼器-解碼器左側為編碼器,右側為解碼器,假設輸入圖像尺寸為 ,編碼器和解碼器均使用4個卷積核尺寸為3的3D卷積層提取圖像局部特征,改變特征圖尺寸和通道數。解碼器中經上采樣的特征圖通過跳躍連接與編碼器輸出的特征圖級聯,跳躍連接會跳躍神經網絡中的某些層,并將一層的輸出作為下一層的輸入,能夠幫助網絡獲得更加精細的細節特征,同時可以解決模型訓練過程中的梯度消失問題。每個卷積層之后都有一個整流線性單元(rectified linear units,ReLU)激活函數,ReLU 是一個分段線性函數,能夠幫助網絡更好地挖掘特征關系,擬合訓練數據。
,編碼器和解碼器均使用4個卷積核尺寸為3的3D卷積層提取圖像局部特征,改變特征圖尺寸和通道數。解碼器中經上采樣的特征圖通過跳躍連接與編碼器輸出的特征圖級聯,跳躍連接會跳躍神經網絡中的某些層,并將一層的輸出作為下一層的輸入,能夠幫助網絡獲得更加精細的細節特征,同時可以解決模型訓練過程中的梯度消失問題。每個卷積層之后都有一個整流線性單元(rectified linear units,ReLU)激活函數,ReLU 是一個分段線性函數,能夠幫助網絡更好地挖掘特征關系,擬合訓練數據。
1.2.2 Transformer編碼器-解碼器
在低分辨率特征處理階段,網絡使用Transformer學習網絡學習圖像的全局特征關系。Transformer編碼器-解碼器首先使用一個卷積層將特征圖尺寸變為 (其中C為設置的通道數),隨后該特征圖輸入由Transformer模塊、圖像塊合并層、轉置卷積層串聯的V型結構,并將相同分辨率級別的特征圖跳躍連接。Transformer模塊不改變特征圖的尺寸,使用基于體積的多頭自注意力機制(volume-based multi-headedself-attention,V-MSA)進行特征提取,V-MSA通過計算查詢-鍵-值向量學習輸入序列間的相關性 [15-16]。圖像塊合并層用于圖像下采樣操作,它連接8個相鄰位置的圖像塊向量,并應用一個線性層產生2倍通道數的輸出,經過上述操作特征圖尺寸減半,通道數翻倍[10]。轉置卷積層的步長為2,用于將低分辨率特征圖上采樣映射成高分辨率特征圖,相比于向上插值采樣 + 卷積,轉置卷積擁有更高的執行效率[17]。
(其中C為設置的通道數),隨后該特征圖輸入由Transformer模塊、圖像塊合并層、轉置卷積層串聯的V型結構,并將相同分辨率級別的特征圖跳躍連接。Transformer模塊不改變特征圖的尺寸,使用基于體積的多頭自注意力機制(volume-based multi-headedself-attention,V-MSA)進行特征提取,V-MSA通過計算查詢-鍵-值向量學習輸入序列間的相關性 [15-16]。圖像塊合并層用于圖像下采樣操作,它連接8個相鄰位置的圖像塊向量,并應用一個線性層產生2倍通道數的輸出,經過上述操作特征圖尺寸減半,通道數翻倍[10]。轉置卷積層的步長為2,用于將低分辨率特征圖上采樣映射成高分辨率特征圖,相比于向上插值采樣 + 卷積,轉置卷積擁有更高的執行效率[17]。
1.2.3 特征指導模塊
在特征指導模塊部分,參考圖像的特征可以更好地指導復雜網絡學習圖像映射關系,首先使用一個步長為1的轉置卷積層用來消除轉置卷積核大小為奇數可能帶來的網格棋盤效應,并用一個卷積層單獨提取參考圖像特征,隨后使用一個卷積層融合來自參考圖像和CNN解碼器的特征圖,最后使用一個卷積層將特征映射為形變場。
1.3 困難形變感知機
圖像中困難形變區域是影響整體配準結果的主要原因,為了篩選圖像中困難形變區域,本文提出了一種基于DDC的困難形變感知機,將其作為多階層方法的橋接裝置。由于無法通過形變場直接衡量困難形變程度,該模塊采用滑動窗口和浮動窗口掃描圖像,計算并比較DDC的大小以篩選困難形變區域。圖像中相似度較低的區域往往也是對齊效果較差的區域,DDC的計算基于局部圖像塊的均方誤差(mean square error,MSE)評分[18],該系數[式(1)中,以符號DDC(·)表示]衡量了圖像中所有體素點的平均困難形變程度,其數值越大,則感知機認為該區域圖像差異越大、越難形變對齊。具有較大解剖差異的圖像塊很難對齊,應盡量選取器官結構信息豐富的圖像塊,為了使感知機聚焦有更多前景體素的有價值區域,滑動窗口會統計當前窗口的前景體素數,并除以當前窗口中的體素總數得到前景比重,用常量加上該比重作為MSE評分[式(2)中,以符號MSE(·)表示]的權重系數,其定義如式(1)和式(2)所示:
|  | 
|  | 
其中, 代表FL和ML所有體素點的DDC,
代表FL和ML所有體素點的DDC, 代表FL和ML所有體素點的MSE評分,FL為參考圖像的局部圖像塊,ML為已配準圖像的局部圖像塊,N為圖像體素總數,NF為前景體素數,p為體素位置,
代表FL和ML所有體素點的MSE評分,FL為參考圖像的局部圖像塊,ML為已配準圖像的局部圖像塊,N為圖像體素總數,NF為前景體素數,p為體素位置, 為圖像域。
為圖像域。
感知機會自動記錄已掃描過圖像塊的DCC最大值以及其在圖像中的位置信息,當窗口掃描到其他圖像塊時,會比較其系數大小,若當前圖像塊的DDC大于已記錄的最大值,則會更新該最大值和窗口位置信息,最終選取該系數值最大的區域作當前圖像的最困難形變區域輸入下一階層子網絡。窗口的尺寸是預先設定的,當圖像中存在多個不同區域大小的復雜形變區域時,相較于大區域,小區域對窗口DDC的影響較弱,同時,DDC代表了局部區域的困難形變程度,感知機會根據該系數選擇當前圖像中的困難形變程度最大的區域。若使用步長為1滑動窗口遍歷整幅圖像,其計算代價較為高昂,為了減小計算量,感知機采用滑動窗口和浮動窗口結合的方式,分兩個階段對圖像進行掃描檢查。困難形變感知機工作原理如圖3所示。
 圖3
				困難形變感知機工作原理圖
			
												
				Figure3.
				The working principle of the difficult deformation perceptron
						
				圖3
				困難形變感知機工作原理圖
			
												
				Figure3.
				The working principle of the difficult deformation perceptron
			
								第一階段為滑動遍歷階段,該階段目標是選取配準對齊較差的區域,通過一個3D滑動窗口遍歷整幅圖并計算DDC,遍歷完成后,記錄第一階段該系數的最大值以及圖像塊的位置信息。本階段滑動窗口的步長設置為5,圖3中藍色實線立方體表示本階段遍歷檢測的滑動窗口示例,紅色實線立方體表示本階段的選中窗口。
第二階段為空間浮動階段,該階段目標是檢查記錄圖像塊周圍區域的困難形變程度,得到最困難形變區域及其位置信息。在第一階段記錄的空間位置處,浮動窗口以滑動遍歷階段記錄的圖像塊為基準進行浮動檢查,沿每個坐標軸有3種圖像塊偏移選擇,即可選擇沿坐標軸正方向或負方向分別偏移浮動步長的體素點位,或者選擇不偏移,3D坐標系排列共計27個方向(包括三個坐標系均不偏移的情況,圖3中僅展示4個方向)的空間浮動,空間浮動完成后,感知機會更新DDC的最大值以及圖像塊的位置信息。之后減小步長并重復上述操作,兩次浮動的步長依次設置為3和1。圖3中的深紅色實線立方體為本階段選中的最終窗口,橙色、藍色、綠色虛線立方體分別表示僅沿Y軸正方向、僅沿X軸負方向、僅沿Z軸正方向偏移浮動步長的浮動窗口示例。
在本文的方法中,滑動窗口和浮動窗口的尺寸是可變的,但由于Tansformer配準網絡在處理圖像時須最小采樣到原圖像尺寸的 ,計算過程需要線性投影對圖像進行分塊處理[9],這要求篩選的圖像尺寸必須是32的倍數。按照多階層配準網絡模型的結構設計,第一階層感知機的窗口尺寸設置為
,計算過程需要線性投影對圖像進行分塊處理[9],這要求篩選的圖像尺寸必須是32的倍數。按照多階層配準網絡模型的結構設計,第一階層感知機的窗口尺寸設置為 ,第二階層感知機的窗口尺寸設置為
,第二階層感知機的窗口尺寸設置為 。
。
1.4 空間變換網絡
為了獲取經過形變場變形的圖像,采用STN對待配準圖像進行空間變換,STN是一個支持反向傳播的可微分模塊,它根據基礎配準網絡提供的形變場φ對輸入圖像進行空間扭曲。對于待配準圖像M的每一個體素p,其扭曲后的體素p'計算公式為 ,u(p)為p的空間位移。p'的大小不一定是整數,但圖像中的體素必須在整數位置處定義,為解決此問題,本文使用適合3D圖像的三線性插值法對體素點進行重新計算。進行體素變換的公式如式(3)所示:
,u(p)為p的空間位移。p'的大小不一定是整數,但圖像中的體素必須在整數位置處定義,為解決此問題,本文使用適合3D圖像的三線性插值法對體素點進行重新計算。進行體素變換的公式如式(3)所示:
|  '/> | 
其中,q表示選定的體素點, 表示p'相鄰位置體素,d為選定的迭代方向,x、y、z分別表示3D坐標系的三個迭代方向,I(·)表示已選定的三線性插值方法。
表示p'相鄰位置體素,d為選定的迭代方向,x、y、z分別表示3D坐標系的三個迭代方向,I(·)表示已選定的三線性插值方法。
1.5 損失函數
無監督學習的方法無需在損失函數中提供標簽信息,在訓練過程中直接通過最小化損失函數即可得到最優形變場。本模型的損失函數由兩部分組成:一部分使用圖像相似損失項,減小圖像間差異;另一部分使用擴散正則化對形變場進行平滑性約束[4]。模型的總體損失函數如式(4)所示:
|  | 
其中, 為模型的總體損失函數,F表示參考圖像,M表示待配準的圖像,
為模型的總體損失函數,F表示參考圖像,M表示待配準的圖像, 使用MSE度量圖像的線性相似度,λ1是為該項的權重系數,其大小通常設置為1。
使用MSE度量圖像的線性相似度,λ1是為該項的權重系數,其大小通常設置為1。 是保持φ平滑的正則化項,λ2為該項的權重系數。總體損失函數優化的主要目標為減小圖像間差異,正則化項對總體損失函數的貢獻應顯著小于圖像相似損失項。對于腦部MRI配準,Balakrishnan等[3]的研究中報告該權重系數設置為0.02是最佳值,因此本文中正則化權重系數λ2大小設置為0.02。
是保持φ平滑的正則化項,λ2為該項的權重系數。總體損失函數優化的主要目標為減小圖像間差異,正則化項對總體損失函數的貢獻應顯著小于圖像相似損失項。對于腦部MRI配準,Balakrishnan等[3]的研究中報告該權重系數設置為0.02是最佳值,因此本文中正則化權重系數λ2大小設置為0.02。 的定義如式(5)所示:
的定義如式(5)所示:
|  | 
由于最小化 可能會導致網絡生成一個不平滑的φ,所以本文使用擴散正則化平滑φ。
可能會導致網絡生成一個不平滑的φ,所以本文使用擴散正則化平滑φ。 的定義如式(6)所示:
的定義如式(6)所示:
|  | 
其中,Δu(p)表示φ中p的空間位移的前向差分,差分結果反映了離散位移值之間的變化,形變場中可能會存異常值,相比于后向差分,前向差分更能鼓勵某一點的位移值相似于其相鄰位置的值,Δu(p)的定義如式(7)所示:
|  | 
2 實驗
2.1 數據集
本文所使用數據集,包括:神經影像學實驗室的洛尼概率腦圖譜(Loniprobabilistic brain atlas 40,LPBA40),該圖譜來自公開數據庫阿爾茨海默病神經影像學倡議(Alzheimer’s disease neuroimaging initiative,ADNI)(網址為:
使用MRI圖像處理開源軟件Freesurfer(v7.1.0,MIT Health HST,美國)對數據集進行顱骨去除、仿射對齊、圖像歸一化、圖像裁剪等標準預處理操作[21],圖像裁剪后大小為 。LPBA40數據集通過不重復的兩兩配對生成760個待配準圖像組,OASIS-1數據集隨機配對生成600個待配準圖像組。兩個數據集的訓練集、驗證集和測試集均按照7∶1∶2的比例分配。
。LPBA40數據集通過不重復的兩兩配對生成760個待配準圖像組,OASIS-1數據集隨機配對生成600個待配準圖像組。兩個數據集的訓練集、驗證集和測試集均按照7∶1∶2的比例分配。
2.2 實驗配置
對比方法采用對稱歸一化(symmetric normalization,SYN)[22]、CM、VM和TM。在14種典型非線性形變算法中,SYN是性能最好的配準算法之一[23],同時它是一種基于迭代的傳統配準算法,本文使用高級歸一化工具軟件包實現[24],實驗中每個級別的最大迭代次數分別設置為160、80、40。CM是使用循環一致性的多尺度配準模型,超參數α、β和λ分別對應其損失函數中循環損失、身份損失和形變場正則化的權重,本文使用Kim等[6]報告的最佳值,超參數α、β和λ分別設置為0.1、0.5和1。VM是一種目前較為流行的無監督可形變配準方法,使用CNN預測形變場,本文使用該方法的2號變體進行對比實驗[4]。TM是在Transformer基礎上開發的一種無監督配準網絡,對于腦部MRI配準,本文使用其推薦的默認參數設置。
級聯多階層配準網絡模型的神經元丟棄率設置為0.3,通道數C設置為64,訓練的批尺寸設置為1,初始學習率為1 × 10?4,學習率衰減策略采用衰減率為0.9的指數衰減策略,訓練迭代次數設置為250,優化器為自適應矩估計(adaptive moment estimation,Adam)優化器。本文采用數據集提供的分割圖來評估各個配準方法的效果。
本文方法訓練模型的實驗設備配置如下:深度學習框架(PyTorch 1.8.0,FAIR,美國) [25],獨立顯卡(Nvidia Geforce RTX 3090,Nvidia,美國),處理器(IntelCore i9-10900K CPU@3.70GHz,Intel,美國)。
2.3 評估指標
本文采用戴斯相似系數(Dice similarity coefficient,DSC)計算分割圖的相似度[26],DSC量化了分割圖中解剖結構之間的體積重疊程度,DSC值越高,說明圖像重疊程度越高,配準性能越好。DSC[式(8)中,以符號DSC(·)表示]的定義如式(8)所示:
|  | 
其中, 表示SF和SM的DSC,SF表示參考分割圖,SM表示配準后的分割圖,SF∩SM表示兩分割圖樣本的交集。
表示SF和SM的DSC,SF表示參考分割圖,SM表示配準后的分割圖,SF∩SM表示兩分割圖樣本的交集。
同時,使用負雅可比行列式(negative Jacobian determinant,NJD)對預測的形變場進行平滑度評價[27],NJD系數[式(9)中,以符號NJD(·)表示]越接近0時,形變場的平均折疊程度越小,形變場越平滑。對于φ中每個p,NJD系數的定義如式(9)所示:
|  | 
其中, 表示φ中p處的NJD系數, i、j、z分別表示3D坐標系中p在
表示φ中p處的NJD系數, i、j、z分別表示3D坐標系中p在 、
、 、
、 三個方向的坐標。
三個方向的坐標。
在圖像的相似性評估上,使用三個指標對圖像的相似性進行評價,這些指標分別是均方誤差平方根(root mean square error,RMSE)[18]、結構相似性(structural similarity,SSIM)[28]、互信息(mutual information,MI)[29]。若圖像的RMSE越高、SSIM和MI越低,代表圖像相似性越差。對于F和M,RMSE[式(10)中,以符號RMSE(·)表示]的定義如式(10)所示:
|  | 
其中, 表示F和M的RMSE值,用來衡量圖像的標準誤差。
表示F和M的RMSE值,用來衡量圖像的標準誤差。
SSIM[式(11)中,以符號SSIM(·)表示]的具體計算過程如式(11)~式(14)所示:
|  | 
|  | 
|  | 
|  | 
其中, 主要關注F和M的邊緣和紋理的結構相似性,以衡量圖像的相似程度,
主要關注F和M的邊緣和紋理的結構相似性,以衡量圖像的相似程度, 、
、 、
、 分別代表亮度比較、對比度比較和結構比較。μF和μM代表F和M的平均值,σF和σM代表F、M的標準差,σFM代表F與M的協方差,c1、c2、c3為非零常數,避免分母為零帶來的系統錯誤,a、b、c均設置為1。
分別代表亮度比較、對比度比較和結構比較。μF和μM代表F和M的平均值,σF和σM代表F、M的標準差,σFM代表F與M的協方差,c1、c2、c3為非零常數,避免分母為零帶來的系統錯誤,a、b、c均設置為1。
MI[式(15)中,以符號MI(·)表示]的定義如式(15)所示:
|  | 
其中, 表示圖像F和M相互包含的信息量,圖像匹配程度越高該值越大,
表示圖像F和M相互包含的信息量,圖像匹配程度越高該值越大, 為F和M的聯合概率密度函數,
為F和M的聯合概率密度函數, 和
和 為分別表示F和M的邊緣概率密度函數。
為分別表示F和M的邊緣概率密度函數。
2.4 結果分析
2.4.1 配準精度
本研究量化了實驗中所有方法在兩個數據集上的平均DSC和NJD系數,結果如表1所示。根據表1結果,本文方法在兩個數據集上的平均DSC均高于對比方法,提供了更好的圖像配準質量,同時生成形變場NJD系數在LPBA40數據集上低于對比方法,在OASIS-1數據集上略高于對比方法。如圖4所示,舉例說明了所有方法的可視化配準結果,包括圖像的細節放大圖、分割圖(包括四個解剖結構皮質、灰質、白質和腦脊液的分割結果)和分割細節,本文使用的方法生成的扭曲圖像和分割圖在外觀上更加接近于原始參考圖像及其分割圖。解剖結構的箱線圖計算了所有方法在OASIS-1數據集上的各個解剖標簽的平均DSC,其中左腦和右腦中的相同解剖結構被合并成一個結構進行計算,具體結果如圖5所示,縱坐標為DSC,橫坐標為解剖結構類別。統計結果顯示,本文提出的方法取得了最好的配準效果。
 表1
                不同方法的配準結果比較
		 	
		 			 				Table1.
    			Comparison of registration results of different methods
			
						表1
                不同方法的配準結果比較
		 	
		 			 				Table1.
    			Comparison of registration results of different methods
       		
       				 圖4
				不同方法的配準結果的可視化示例
			
												
				Figure4.
				Visual examples of registration results for different methods
						
				圖4
				不同方法的配準結果的可視化示例
			
												
				Figure4.
				Visual examples of registration results for different methods
			
								 圖5
				解剖結構分組繪制的箱線圖
			
												
				Figure5.
				Boxplots of anatomical structures grouped
						
				圖5
				解剖結構分組繪制的箱線圖
			
												
				Figure5.
				Boxplots of anatomical structures grouped
			
								2.4.2 階層對比
本文方法設置的階層數為3,理論上,階層數的設置是任意的,為了驗證多階層的方法可以實現對形變場的逐步優化,本研究嘗試使用更多階層進行實驗,但由于本文的基礎配準網絡對圖像尺寸存在限制(必須為32的倍數),同時第一階層子網絡已經篩選了大尺寸的困難形變區域,導致本研究僅能在第二階層和第三階層子網絡中間增加圖像尺度為 的階層子網絡,因此,本研究使用4階層子網絡進行實驗,對訓練好的最佳預訓練模型進行微調訓練,對比了使用不同階層數子網絡生成的形變場扭曲分割圖的平均DSC,如表2所示,可以發現隨著階層的增多,平均DSC在不斷提高,證明多階層方法可以實現對圖像的漸進優化配準,提升配準效果。如圖6所示,對不同階層子網絡關注的困難形變區域及該區域的配準結果和形變場進行可視化展示,實線框出的區域代表當前階層子網絡篩選的困難形變區域,不同階層由不同顏色的實線框進行標注。
的階層子網絡,因此,本研究使用4階層子網絡進行實驗,對訓練好的最佳預訓練模型進行微調訓練,對比了使用不同階層數子網絡生成的形變場扭曲分割圖的平均DSC,如表2所示,可以發現隨著階層的增多,平均DSC在不斷提高,證明多階層方法可以實現對圖像的漸進優化配準,提升配準效果。如圖6所示,對不同階層子網絡關注的困難形變區域及該區域的配準結果和形變場進行可視化展示,實線框出的區域代表當前階層子網絡篩選的困難形變區域,不同階層由不同顏色的實線框進行標注。
 表2
                不同階層數形變場的配準結果比較
		 	
		 			 				Table2.
    			Comparison of registration results using different levels of deformation fields
			
						表2
                不同階層數形變場的配準結果比較
		 	
		 			 				Table2.
    			Comparison of registration results using different levels of deformation fields
       		
       				 圖6
				不同階層關注的感興趣區域對比示例圖
			
												
				Figure6.
				Comparative examples of regions of interest concerned by different levels
						
				圖6
				不同階層關注的感興趣區域對比示例圖
			
												
				Figure6.
				Comparative examples of regions of interest concerned by different levels
			
								2.4.3 感興趣區域識別
困難形變感知機的識別結果直接影響了后續階層能否關注到圖像中有價值信息的區域,但形變場的評價指標沒有金標準,無法直接判斷所選區域的復雜性,因此本研究使用OASIS-1測試集進行實驗,對第一階層子網絡中感知機識別前后的圖像進行相似度評估,進而判斷感知機能否識別圖像的困難形變區域,對比結果如表3所示,其中全局圖像代表已配準圖像對,局部圖像代表經過識別后困難形變區域。參數結果顯示,局部圖像的結構相似度更差,互相包含的信息量低于全局圖像,證明該感知機可以識別困難形變區域。
 表3
                圖像相似性對比結果
		 	
		 			 				Table3.
    			Image similarity comparison results
			
						表3
                圖像相似性對比結果
		 	
		 			 				Table3.
    			Image similarity comparison results
       		
       				3 結論
本文介紹了一種基于Transformer的級聯多階層醫學影像配準方法。該方法首先基于Transformer構建了基礎配準網絡,有效地融合圖像的全局信息和局部信息,同時配備了一種基于DDC的困難形變感知機,可以識別并裁剪圖像的感興趣區域。然后以困難形變感知機為橋接方法,級聯了三個包含基礎配準網絡的階層子網絡,并且采用嵌入融合的方式生成形變場,構建了基于Transformer的級聯多階層配準網絡模型。實驗結果表明,與當前流行的無監督配準方法相比,本文方法可以漸進地優化配準結果,解決圖像的復雜形變問題,提升配準精度,幫助醫生做出更準確的臨床診斷,促進計算機輔助醫療技術的發展。但目前本文的方法只適用于單模態圖像配準,在后續工作中會考慮將該方法擴展到多模態圖像配準。
重要聲明
利益沖突聲明:本文全體作者均聲明不存在利益沖突。
作者貢獻聲明:潘英杰負責本研究的算法設計、程序編寫、實驗數據處理和論文撰寫;程遠志是本研究的負責人,指導實驗設計和論文寫作,提出修改意見;劉豪和史操負責論文修改和數據分析。
 
        

 
                 
				 
																   	
                                                                    
                                                                    
																	 
																   	
                                                                    
                                                                    
																	 
																   	
                                                                    
                                                                    
																	 
                                                                    
                                                                        
                                                                        
                                                                         
																   	
                                                                    
                                                                    
																	 
																   	
                                                                    
                                                                    
																	