基于光流計算與增強網絡的頸椎CT 圖像層間插值算法研究

王 明，吳曉紅，滕奇志，楊 毅，任 超*

（1.四川大學電子信息學院，成都 610065；2.四川大學華西醫院骨科，成都 610041）

0 引言

人體頸椎CT 序列圖像通過X 射線對人體頸椎作連續的斷面掃描得到。通過觀察患者的頸椎CT圖像，醫療人員可以診斷出患者頸椎是否發生病變。但由于X 射線對人體組織有損傷[1]，為節省拍攝時間，CT 圖像的采樣間隔較大，導致CT 圖像的縱向點長度遠大于片層內的橫向點長度。因此，直接由人體頸椎CT 序列圖像的橫截面生成的矢狀面和冠狀面圖像會產生明顯變形（如圖1 所示），與頸椎的真實情況不符。變形的矢狀面和冠狀面圖像會對后續頸椎椎骨的終板形態分型研究帶來影響。為解決此問題，一般先對人體頸椎CT 序列圖像中相鄰層的CT 圖像進行層間插值，縮短CT 序列圖像的縱向點長度，再生成矢狀面和冠狀面圖像，其形態結構就會與頸椎真實情況相符，更有利于后續對人體頸椎模型進行分析。

圖1 原始頸椎CT 圖像的矢狀面和冠狀面

基于配準的插值方法[2-4]在醫學領域得到了廣泛應用，這種插值方法先通過配準得到2 張CT 圖像的形變場，由形變場進行擬合插值。張杰等[3]提出了基于曲面擬合與配準的插值算法，針對非剛性配準后的B 樣條插值算法對異常點不敏感的問題，引進了薄板曲面擬合的方法，有效地去除了異常點的影響，并且更進一步利用雙向形變場將2 個形變場融合后再結合樣條插值與曲面擬合來得到插值圖像。這種做法一定程度上提高了插值圖像的精度，但是插值圖像在細節的局部放大中仍然有毛刺、偽影現象。張耀等[4]利用非下采樣輪廓波變換（nonsubsampled contourlet transform，NSCT）對圖像進行分解，對NSCT 分解后的高頻和低頻子圖分別采用不同插值算法生成新插值圖像的高頻和低頻子圖，再將二者逆NSCT 生成新插值圖像。這種算法生成的圖像仍然有模糊現象，并且由于迭代運算和非下采樣技術，運算十分耗時，在速度上需要進一步優化。

基于光流信息配準后再進行插值的研究中，Ehrhardt 等[5]基于光流對序列圖像插值，不僅提高了插值圖像的精度還盡可能保留了目標原始結構，但是邊緣輪廓依然不清晰。林毓秀等[6]先利用聯合全局與局部（combined local and global，CLG）算法生成光流模型，通過此光流模型進行配準與生成插值圖像，再利用序列圖像的幀間非局部自相似性來解決插值后圖像中的像素缺失問題，從而能夠生成高質量的中間層圖像。但是，當圖像層間距較大時，GLG 算法生成的中間圖像準確性不高。

視頻插幀[7-13]和頸椎CT 圖像層間插值具有一定的相似性，可以借鑒視頻插幀算法完成CT 圖像的層間插值任務。視頻插幀算法是給定2 個連續幀生成中間幀，以形成空間和時間上一致的視頻序列。文獻[7]使用卷積神經網絡（convolutional neural network，CNN）直接合成中間幀，由于網絡沒有學習相鄰幀之間光學變化信息導致生成的圖像十分模糊。文獻[8]使用空間自適應卷積核來估計每個像素值，當卷積核變大時，運算成本很高。文獻[9]先對相鄰幀的光流進行估計，再對遮擋進行處理。文獻[10]引入了上下文紋理特征進行幀插值，減少了插值圖像的偽影現象。文獻[11]綜合了光流估計、深度圖和上下文紋理特征對輸入幀進行處理，有效增強了中間幀的質量。文獻[12]提出了一種基于雙邊運動估計的視頻插幀算法，首先使用近似雙向運動網絡來預測雙邊運動，然后利用2 個輸入幀和雙向運動生成中間幀，最后使用動態混合過濾器網絡過濾中間幀來生成輸出幀。文獻[13]提出了自適應協作流網絡（adaptive collaboration of flows，AdaCoF）來估計每個像素點的核權重和偏移向量以合成輸出幀，其可以處理復雜的運動且能生成更真實的輸出圖像。

如果將視頻插幀算法應用到層間插值任務上，還需要考慮CT 圖像本身每一像素點的橫向點長度和縱向點長度信息，不能在2 張CT 圖像間任意插入圖像。并且由于人體CT 圖像采樣間隔較大，CT 圖像的變化比較劇烈，所以相鄰上下層圖像對中間層圖像有更重要的參考性。

因此，為了解決人體頸椎CT 序列圖像橫向點長度和縱向點長度不匹配帶來的矢狀面和冠狀面結構變形問題，本文基于殘差密集網絡和UNet 網絡，針對人體頸椎CT 圖像的特點制作專用數據集，通過CT圖像雙向光流計算與增強網絡來進行人體頸椎CT圖像的層間插值。

1 人體頸椎CT 圖像層間插值算法

1.1 算法整體概述

CT 圖像層間插值算法利用相鄰的上下層圖像之間的變化生成中間層圖像，而視頻插幀算法是在時間序列上的兩幀之間根據前后幀圖像之間的變化和時間信息生成中間幀，二者具有較強的相似性，因此，借鑒視頻插幀算法[9]，將相鄰的3 張序列圖像分為上層圖像Iu，需要插值生成的中間層圖像Im和下層圖像Id。先計算Iu和Id之間的雙向光流Fu→d、Fd→u，由此雙向光流通過公式（1）和公式（2）計算出Im到Iu的近似光流以及Im到Id的近似光流F?m→d：

式中，t=0.5，表示Im到Iu以及Im到Id的光流應為Iu和Id之間光流的一半。

式中，g（·，·）根據光流和光流的一端圖像進行雙線性插值反向生成光流另一端圖像；參數a 的值取決于兩邊圖像對中間圖像的貢獻。根據光流信息，部分像素塊可能在另一邊圖像上消失或者存在灰度變化，因此本研究引入可視圖Vm←u和Vm←d表示像素塊從圖像Iu到Id的灰度變化。由此公式（3）可以轉換為公式（4）：

1.2 CT 圖像雙向光流計算網絡

雙向光流計算網絡基于殘差密集網絡[14]，利用3個殘差密集塊（residual dense blocks，RDB）來提取特征，網絡框架如圖3 所示。其由淺層特征提取部分、由3 個RDB 構成的殘差密集特征提取部分和多級特征融合部分組成。在淺層特征提取部分，使用2 層卷積核為3 的卷積網絡來提取初始光流特征。接著通過3 個RDB 將不同RDB 的輸出融合后輸入到多級特征融合網絡中。RDB 既有殘差網絡的特點也有密集連接的卷積網絡的特點，其模型結構如圖4 所示。RDB 將每一級的卷積輸出經過ReLU 連接起來形成一個密集的殘差結構。多級特征融合部分通過1 個1×1 卷積和1 個3×3 卷積提取多級光流特征，使用殘差連接，再經過3×3 卷積后輸出最終的雙向光流。

圖2 算法整體框圖

圖3 CT 圖像雙向光流計算網絡框架圖

圖4 RDB 模型結構圖

1.3 CT 圖像雙向光流增強網絡

通過CT 圖像雙向光流計算網絡得到上下層圖像之間的光流后，由公式（1）和（2）可以得到中間圖像到上下層圖像的雙向近似光流。因為此近似光流是由上下層圖像之間的光流得到的，仍然需要考慮上下層圖像間目標結構變化的問題，所以光流精度還需加強。本節基于UNet 設計了一個CT 圖像雙向光流增強網絡。

Ronneberger 等[15]提出的UNet 網絡是由對稱的編碼器和解碼器組成的卷積神經網絡。編碼器對圖像進行下采樣提取特征，解碼器對圖像進行上采樣恢復圖像，編碼器和解碼器之間有跳層連接。在編碼器中，共有4 個下采樣層，每個下采樣層先進行池化再經過2 個卷積層。下采樣時每經過一層池化層，特征圖尺寸就縮小一半，卷積核數目增加為原來的2倍。在解碼器中，共有4 個上采樣層，每個上采樣層先對特征圖進行上采樣，并和編碼器對應同層數據做尺度融合后再經過2 次卷積操作。上采樣時先對特征圖使用雙線性插值拓展圖像尺寸為原來的2倍，經過跳層連接后再經過2 次卷積操作。

UNet 使用卷積對空間局部特征進行提取，缺少了對全局信息的感知能力。傳統的注意力機制通過獨立的鍵值對（key-value pair）生成注意力矩陣，通過注意力矩陣對權重進行分配，從而增強關鍵特征對輸出的影響，但是此種做法忽略了相鄰鍵與鍵之間的關聯信息。而CoTBlock 注意力結構[16]充分利用了輸入鍵之間上下文信息來指導動態注意力矩陣的學習，從而提升了特征表達能力。其模塊結構圖如圖5 所示。

圖5 CoTBlock 模塊結構圖

對于給定的H×W×C 的二維特征圖X，其鍵矩陣定義為K=X，查詢矩陣定義為Q=X，值矩陣定義為V=XWv，其中Wv是利用1×1 卷積實現的嵌入矩陣。通過k×k 組卷積對輸入特征X 進行編碼，得到輸入特征X 的靜態局部關聯性權重K1，其反映了局部鍵之間的靜態上下文信息。將K1和X 進行拼接后通過2 個連續的1×1 卷積得到動態的多頭注意力矩陣A：

式中，Wθ和Wδ是分別利用2 個1×1 卷積實現的嵌入矩陣。注意力矩陣A 具有查詢特征和上下文鍵特征，將其與V 相乘后得到特征圖K2：

式中，特征圖K2代表輸入之間的動態交互特征，其反映了局部鍵之間的動態上下文信息。將靜態上下文K1和動態上下文K2通過注意力機制融合得到輸出。經過CoTBlock 層，根據靜態上下文信息指導注意力矩陣的生成，再通過注意力矩陣增強關鍵特征的權重可以有效提取光流特征。

CT 圖像雙向光流增強網絡在UNet 架構基礎上改造了編碼階段的結構，如圖6 所示。首先，分別使用7×7 卷積層和5×5 卷積層增加感受野。其次，每級特征圖經過池化層進行下采樣后，先通過1 個卷積層，再通過CoTBlock 層得到下一級特征圖，其中每經過一次操作都需要額外再通過1 個Leakly ReLU層。CoTBlock層充分利用了輸入鍵之間上下文信息來指導動態注意力矩陣的學習，從而提升了特征表達能力。最后，在編碼階段的末尾額外增加了一次下采樣操作，該下采樣操作只對特征圖進行下采樣而不增加通道數。

圖6 CT 圖像雙向光流增強網絡框架圖

1.4 損失函數

對于圖像層間插值任務，可以利用損失函數計算生成的中間層圖像與真實圖像之間的差異。本文首先用圖像損失函數lr來計算目標圖像Im與預測圖像之間的誤差：

同時，引入高維特征損失函數lp來度量目標圖像與預測圖像高維度特征圖之間的誤差：

式中，φ 是ImageNet 預訓練的VGG16[17]模型中conv4_3 輸出的特征。

另外，利用光流損失函數lw來計算光流誤差，由于沒有真實光流對照值，由計算出的光流通過g（·，·）反向生成圖像與原始圖像對照來近似替代光流誤差：

本研究還添加平滑損失函數ls用來使光流預測結果更加平滑：

式中，代表哈密頓算子。

最后，結合公式（7）、（8）、（9），本研究CT 圖像光流計算與增強網絡的損失函數l 的計算公式為

式中，λr、λw、λp和λs分別為每部分損失函數的權重，由于每部分損失函數計算得到的結果相差較大，為了調節每部分反向傳播的能力，將公式中λr、λw、λp和λs的值分別設置為204、102、0.005 和1。

2 實驗方法和結果

2.1 數據集的制作

為了使CT 圖像雙向光流計算與增強網絡學習到相鄰頸椎CT 圖像光流變化的細節特征，本文構建了針對人體頸椎CT 圖像層間插值的數據集。為了研究人體頸椎椎間盤形態，四川大學電子信息學院圖像信息研究所實驗室與四川大學華西醫院合作，選取了32 名患者的人體頸椎CT 序列圖像來制作數據集。這32 名患者的頸椎CT 序列圖像形態結構良好且圖像的橫向點長度為0.23～0.26 mm，縱向點長度為0.625 mm。在每位患者的人體頸椎CT 序列圖像中，大約含有200余張CT 圖像，其中每相鄰的3 張圖像分為一小組，每位患者的CT 圖像大約能分成200 小組。每小組中，上層圖像和下層圖像作為輸入，中間圖像作為層間插值后的驗證圖像供神經網絡訓練，每張圖像大小均為512×512 像素。從32 名患者的人體頸椎CT 序列圖像中選擇25 名患者，即共5 226 組圖像作為訓練集；選擇2 名患者，共422 組圖像作為驗證集；選擇5 名患者，共1 061 組圖像作為測試集。在訓練時，正常的訓練順序是第一張、第三張作為輸入，但為了擴充數據避免過擬合，且相鄰圖像的輸入順序并不影響插值結果，本文隨機采取第一張、第三張和第三張、第一張2 種方式加載圖像來擴充數據集。

2.2 實驗設置

本實驗中，實驗計算機的系統為Ubuntu20.04系統，CPU 為Intel?CoreTMi7-9700，內存為32 GiB，圖形處理器（graphics processing unit，GPU）為NVIDIA RTX2080Ti。

本研究使用Pytorch 框架實現算法，在人體頸椎CT 圖像插值數據集上訓練網絡，網絡學習速率設為1×10-4，batchsize 設置為2，使用Adam 方法來優化網絡模型參數，訓練65 個Epoch 后網絡收斂。

2.3 實驗結果與分析

為了驗證本文網絡模型的有效性，將本文算法與基于配準的線性插值算法、Penny 插值算法[2]、基于曲面擬合的雙向配準層間插值（using surface fitting with bi-directional registration，USFBR）算法[3]、SuperSlomo（超級慢動作）幀間插值算法[9]和雙邊運動估計（bilateral motion with bilateral cost，BMBC）幀間插值算法[12]進行比較，并采用平均峰值信噪比（peak signal-tonoise ratio，PSNR）和平均結構相似性（structural similarity，SSIM）來比較生成圖像的質量與原始圖像的相似性。實驗測試了5 組共1 061 張人體頸椎CT 圖像的插值效果，結果見表1、2。同時，以病例47 為例，列出了同一病例不同算法生成的中間層的局部放大圖像，如圖7 所示。

表1 不同算法PSNR 結果 單位：dB

表2 不同算法SSIM 結果

圖7 同一病例不同算法生成的中間層圖像對比

由表1、2 結果可知，本文算法的PSNR 和SSIM均優于其他算法。觀察圖7 可以看出，基于配準的線性插值算法生成的圖像不僅會出現信息的明顯缺失，而且圖像存在明顯的毛刺和重影現象；Penny 插值算法相較于基于配準的線性插值算法生成的圖像結構上更符合原始圖像，但仍存在模糊、重影現象；USFBR 算法相較于Penny 插值算法生成的圖像質量有了進一步提升，但邊緣輪廓依然不清晰；Super Slomo 算法生成的圖像基本消除了偽影、毛刺問題，但是邊緣結構仍與原始圖像有差異；BMBC 算法生成的圖像和原始圖像較接近，但是圖像質量還有待提高；本文算法生成的圖像不僅消除了毛刺和偽影問題，并且邊緣輪廓與原始圖像也更為相似，同時減少了模糊現象。綜上，本文算法確實提升了插值圖像的質量，并且優于其他算法。

將不同算法所需的參數量、不同算法生成圖像的運行時間進行對比，結果見表3。由于深度學習算法可以使用顯卡加速，其生成圖像的速度快于傳統算法。本文算法在降低參數量與運行時間的同時能夠得到最佳的插值圖像。

表3 不同算法參數量和運行時間對比

由圖8 可知，經過插值后的CT 序列圖像水平面生成的矢狀面和冠狀面不僅符合頸椎組織結構的變化，也更貼近人體頸椎的真實情況。

圖8 經過層間插值的頸椎CT 圖像的矢狀面和冠狀面

2.4 消融實驗

本文算法是通過搭建網絡模型實現的，為驗證本文提出的CT 圖像雙向光流計算網絡中的殘差密集網絡和CT 圖像雙向光流增強網絡中的CoTBlock模塊對頸椎CT 圖像層間插值任務的有效性，進行消融實驗比較上述模塊對結果的影響。其中，在第一階段與第二階段均使用UNet 網絡作為Baseline，上采樣均使用雙線性插值。實驗結果見表4、5，可以看出本文提出的模型在PSNR 和SSIM 方面均取得了最佳效果，殘差密集網絡和CoTBlock 模塊均能有效提升插值圖像的精度。殘差密集網絡可以較為精確地提取圖像之間雙向光流的特征，而CoTBlock 模塊利用了注意力機制和相鄰鍵之間的上下文信息加強了光流精度。

表4 不同網絡模型的PSNR 結果 單位：dB

3 結語

為解決人體頸椎CT 圖像層間插值問題，本文構建了人體頸椎CT 序列層間插值數據集，并使用此數據集訓練了一個兩階段的深度學習網絡。網絡的第一階段計算輸入2 張圖像的雙向光流場，第二階段提高近似光流場精度，并最終生成插值圖像。在第一階段中使用了CT 圖像雙向光流計算網絡，該網絡使用RDB 計算2 張圖像之間的雙向光流。第二階段使用了CT 圖像雙向光流增強網絡，該網絡使用CoTBlock 模塊改進了UNet 網絡，CoTBlock 模塊利用自注意力機制和鍵與鍵之間的上下文信息加強了光流精度并生成最終的中間層圖像。實驗表明，本文算法可以生成較高質量的中間層圖像。而且本文算法消耗的時間均在訓練網絡模型上，訓練完成后，本文算法在生成中間層圖像的速度上有更好的表現。脊柱外科專家觀察了經過插值后的頸椎CT 序列圖像，認為其矢狀面和冠狀面形態更符合真實的人體頸椎結構形態，可將插值后的頸椎CT 序列圖像用于頸椎終板形態結構的研究和分析。但是當相鄰圖像差異過大時，本文算法只能盡可能模擬層間圖像的變化，精度還有待提高。后續可以考慮用人體頸椎的矢狀面和冠狀面信息來修繕插值圖像，提升插值圖像的質量。另外，本文網絡模型是一個兩階段的網絡，網絡參數量相對較大，可以進一步研究，減少網絡參數量。

表5 不同網絡模型的SSIM 結果