基于姿態編碼器的2D/3D脊椎醫學圖像實時配準方法

徐少康，張戰成，2*，姚浩男，鄒智偉，張寶成

（1.蘇州科技大學 電子與信息工程學院，江蘇 蘇州 215009；2.蘇州科技大學 蘇州市虛擬現實智能交互及應用技術重點實驗室，江蘇 蘇州 215009；3.中部戰區總醫院 骨科，武漢 430070）

0 引言

2D/3D 醫學圖像配準是三維手術導航系統中的關鍵技術之一，在脊椎手術中，椎弓根螺釘植入是臨床脊柱外科手術中最常用的內固定方法，術中精確植入螺釘的同時，不損傷血管和神經是一項挑戰性的工作［1］。醫生術中通過X 射線圖像（2D）對病人脊椎位置進行分析，查看螺釘的植入情況。但2D 圖像的清晰度能夠反映出的信息較少，想要精確地分析出植入螺釘位置信息較為困難。而3D 圖像可以提供此類的多維信息，如果在手術過程中多次實施CT 成像會對醫生以及病人造成傷害，將術前3D 浮動圖像對照術中2D 參考圖像進行空間轉換達到坐標系相統一，從3D 浮動圖像中引導出對應位置的清晰的2D 圖像，即可更精確地引導手術，提高植入椎弓根螺釘的精確度與外科醫生的手術效率，同時減少醫生和病人的輻射傷害。

傳統的基于迭代回歸優化的2D/3D 配準方法［2-4］：從一個假設的姿態開始，通過對CT（Computed Tomography）進行模擬虛擬X 射線的衰減創建數字重建影像（Digitally Reconstructed Radiographs，DRR），評估DRR 圖像與術中X 射線圖像的相似性，并優化搜索下一個潛在的姿態；然后繼續為這個新姿態生成DRR，直到相似性達到最大時，完成配準。然而，在該傳統方法中涉及多次的DRR 成像以及大量的圖像相似度計算，計算成本都很高。因此該類方法無法滿足手術過程中對實時配準的需求。此外，在配準優化過程中，優化器容易陷入局部最優值，導致配準的成功率無法滿足臨床應用的精度要求。

近年來，基于機器學習的配準方法已經逐漸運用于醫學圖像配準任務當中［5-9］。Miao 等［10］提出基于卷積神經網絡回歸方法，分層學習姿勢估計（Pose Estimation Hierarchical Learning，PEHL），以實現具有較大捕獲范圍的實時2D/3D 配準。為了能夠提高2D/3D配準速度，Pan等［11］提出了一種基于閾值分割的方法來加速配準過程，通過使用基于Bresenham方法的加速DRR 算法，以減少生成數字重建射線照片的時間，使得光線投射的時間下降到0.937 s；但該方法僅是提高了DRR 的計算速度，并沒有真正避免迭代計算的過程，在整個配準的時間中，仍然耗費了91.301 s。Wohlhart 等［12］提出訓練卷積神經網絡（Convolutional Neural Network，CNN），從距離圖像中學習姿勢區分描述符，并使用最近鄰進行姿勢估計。雖然使用該方法可以實現全局姿態估計，但精度相對較低，即角度誤差小于5°的成功率低于60%。Liao等［13］通過用于跟蹤和三角測量的興趣點網絡（Point-Of-Interest Network for Tracking and Triangulation，POINT2）來解決多視圖2D/3D 剛性配準問題。在2D 和3D 中建立干預前和干預內數據之間的點對點對應關系，然后在匹配點之間進行形狀對齊以估計CT的位姿，從而避免了通常昂貴且不可靠的迭代姿勢搜索；但該方法需要來自至少兩個視圖的X 射線可用。Gouveia 等［14］從X 射線圖像中提取了一個手工制作的特征，并訓練了一個多層感知器（Multi Layer Perceptron，MLP）回歸器來估計三維變換參數；然而該方法所能達到的準確度遠低于使用基于優化的方法所能達到的準確度。

Sundermeyer 等［15］提出了一種基于RGB 的實時目標檢測和6D姿態估計方法，提供了一個由隱空間中的樣本定義的對象位姿的隱空間表示。數據的隱空間表示包含表示原始數據點所需的重要信息，隱空間表示用于將更復雜的原始數據形式（如圖像、視頻）轉換為更“易于處理”和分析的簡單表示［16-18］。Gu 等［19］提出將隱空間分解為兩個次隱空間：模態特定的隱空間和姿勢特定的隱空間，用于三維手部姿態估計，在基于單圖像的手部姿態估計方面的性能明顯優于目前最先進的方法。Jiang 等［20］提出了一種通用的自適應對抗隱空間框架，該框架主要由自適應隱空間發生器（Adaptive Latent Space Generator，ALSG）和基于約束的對抗自動編碼器（Constrained Antiversarial AutoEncoder，CAAE）兩部分組成。通過自適應映射器和對抗學習方法，建立ALSG，獲得真實的潛在空間分布。這些研究表明隱空間可以較好地重構姿態信息。

為了提高2D/3D醫學圖像配準的速度，達到術中實時配準的要求，受隱空間回歸啟發，將2D/3D 醫學圖像配準問題表達為基于自編碼的回歸網絡，并加入隱空間回歸約束，以捕獲X射線圖像中姿態信息，從而精確地回歸脊椎的旋轉的姿態參數。對比了傳統的優化迭代配準方法，本文方法避免了大量的DRR渲染時間，使得2D/3D圖像配準的時間滿足實時的要求。

1 相關工作

1.1 DRR投影模型

DRR 圖像生成算法有許多種：拋雪球法（Splatting）、光線投射法（Ray-Casting）、光流場法（Light-Field）等。在2D/3D 醫學圖像配準中主要使用光線投射算法。Gao 等［21］提出的投影空間變換器模塊（Projective Spatial Transformer，ProST）可將空間變換器推廣到投影幾何，從而實現可微分的體積渲染。

從三維立體場景到投影透射圖像Im的映射可以建模為Im=T(θ)V，其中T(θ)是依賴于姿態參數θ的變換矩陣，V是3D CT體積圖像。在基于強度的2D/3D配準中，尋求檢索姿勢參數θ，以便從3D CT掃描模擬的圖像盡可能趨近于采集的圖像。投影空間變換通過學習控制點G的變換來擴展標準投影幾何，每條投影光線在CT 內有K個均勻間隔的控制點。給定θ∈SE(3)，通過仿射變換矩陣T(θ)獲得一組變換后的控制點，投影平面長寬為M×N。T(θ)通常由三個平移參數tx、ty、tz和三個旋轉參數rx、ry、rz來參數化，可以寫成齊次坐標下的4×4矩陣：

其中R(θ)是控制CT 體積繞原點旋轉的旋轉矩陣。由于這些控制點位于體積內、位于體素之間，因此對控制點的值進行插值。

其中Gs∈RM×N×K，最后得到一個二維圖像I∈RM×N沿每條射線積分，通過“堆疊”(m，n)處的Gs的k維度來實現。

上述過程利用了空間變換網格，將投影操作簡化為一系列線性變換，使得2D/3D 配準的DRR 圖像生成過程變得快速精確。

1.2 三維轉換參數

一個剛體三維變換可以用一個含有6 個參數的矢量來參數化［10］。在本文方法中，通過3個平面內和3個平面外的轉換參數來參數化這個三維變換。其中，平面內的變換參數包含2個平移參數tx、ty和1個旋轉參數rz。平面內的變換參數的結果看作為二維剛體變換。平面外的變換參數包括1 個面外平移參數tz和2 個面外旋轉參數rx和ry。平面外平移和旋轉的結果分別對應著剛體圖像的縮放和形狀變化。

在脊椎手術中，為了能夠提示外科醫生在鉆探椎弓根螺釘導孔時可能引發的椎體皮質骨折，通常會在植入椎弓根螺釘前向脊椎內植入一枚探針用于探查脊椎管道情況，再選擇合適的椎弓根螺釘進行植入。在進行脊椎圖像配準前，外科醫生可以直接通過探針植入的位置準確地知道需要配準的是第幾節脊椎位置，所以對于X射線圖像的固定平移拍攝，平移的三個參數是已知的，只需要計算出旋轉姿態參數rx、ry、rz即可完成配準。

2 脊椎姿態回歸方法

2D/3D 配準的目標是搜索CT 圖像的6 個未知自由度姿態，以該姿態對CT 體積模擬投射X 射線在平面上產生DRR與X 射線圖像進行比較。因此2D/3D 配準的問題最終可以歸結到如何從X射線圖像中尋找最佳CT投影姿態的問題。

2.1 配準模型結構

訓練模型結構如圖1 所示，前半部分由Encoder 加全連接層組成，本文Encoder 采用卷積自編碼器架構。該結構由4 個卷積組成，每層使用3×3的卷積核，步長為2，在編碼過程中應用2D 卷積。在實驗中，輸入圖像大小為128×128×1，但框架不受特定大小的限制。每個卷積后面都有一個參數為0.2 的LeakyReLU 層。模型的輸入為脊椎X 射線圖像，經過Encoder將X射線圖像特征映射到隱空間并從中提取隱變量。然后通過接入128 → 64 → 32 → 16 → 3 全連接層，回歸出對應的姿態參數向量。通過投影得到DRR 圖像，即配準圖像。在訓練過程中，將配準圖像輸入前者的Encoder 模塊，從配準圖像中繼續提取隱變量，用于輔助損失的計算（具體計算方式見2.3節）。

圖1 配準網絡框架Fig.1 Framework of registration network

式（5）中：ICT為CT圖像，IX為X射線圖像，φ為Encoder+MLP姿態回歸過程模型，Proj 為DRR 投影，Ireg為配準圖像。編碼器用F表示，該網絡將原始數據X映射到隱空間F中(φ：X→F)。通過用激活函數σ傳遞的標準神經網絡函數表示，其中Z是隱變量。

2.2 隱空間生成

在隱空間特征表示中，相似的樣本圖像之間的差別特征會被作為非核心信息被剔除，只有其中核心特征信息會被保留并且學習，所以在數據特征點被映射到隱空間后，相似特征點距離會更近。在本文的模型訓練過程中，輸入為不同姿態下的脊椎X 射線圖像。人體的脊塊部位特征基本類似，差異性較小，所以對于此類的訓練圖像，姿態信息即為核心特征信息，本文通過Encoder生成的隱空間及隱變量會最大限度地學習脊椎的姿態信息，并精確地回歸出姿態參數。

在自編碼器結構中，通常使用Decoder層來進行圖像的重建，從而訓練隱空間的生成過程。然而如果使用Decoder重建配準圖像，與CT 圖像中的脊椎信息沒有直接關系，缺少醫學的可解釋性，會對醫生的判斷造成干擾。本文中使用DRR 投影代替Decoder工作，DRR投影的計算原理與X射線圖像的生成原理類似，都是通過光射線穿過體素后的衰減情況來確定2D 圖像各像素值。相比之下，以CT 圖像為基礎投影生成的DRR 圖像更加滿足本文方向上圖像重建的質量要求，外觀與X射線圖像更加相似。

2.3 損失函數

對于神經網絡而言，性能除了取決于本身的結構之外，損失函數的設計也極其重要。為了更好地學習網絡參數，本文設計了一個基于“粗細”結合的損失函數。其中“粗配”包括預測的姿態參數與真實姿態r(i)的L1 范數損失LossL1。在訓練的過程中，通過姿態編碼器回歸出的姿態參數以長度為3的向量的形式存在，并在回歸完成前對其進行歸一化，使得姿態參數向量的3個旋轉姿態值在（-1，1），對應(-57.3°，57.3°)角度范圍，因此L1范數損失適合此類信息的計算。此外還包括X 射線圖像中提取的隱變量Zi與DRR 圖像中提取的隱變量使用余弦相似度計算的Losscosine。損失函數定義如下：

其中γ是一個正則化參數（在3.5節中會對其進行實驗分析取值）。實驗表明當初始姿態與真實姿態差異較大時，式（9）損失函數在收斂過程中更有效，但當DRR 圖像與X 射線圖像逐漸一致時，該損失函數對細微的差異敏感性變低。

在此基礎上，本文繼續引入一種用于“細配”的損失函數：基于梯度的歸一化互相關（Gradient-based Normalized Cross Correlation，GradNCC）［22］，基于梯度的度量最初通過微分進行轉換。水平和垂直的Sobel算子用于創建梯度圖像，并在圖像的兩個正交軸上表示透視強度的導數。然后計算DRR 的梯度圖像和X 光梯度圖像之間的歸一化互相關，該測量的最終值是這些歸一化相互關的平均值。梯度測量的優點是濾除圖像之間的低空間頻率差異，如軟組織引起的差異。在剛體的圖像中，邊緣含有豐富的特征信息，所以該損失集中于邊緣信息的相似性度量。將該損失函數可以表示為：

其中：IXI與IDRR分別對應X 射線圖像與DRR 圖像在區域(i，j)的強度值，與為重疊區域(i，j) ∈T內圖像的均值。

本文使用隨機梯度下降（Stochastic Gradient Descent，SGD）優化器進行優化，計算了Losscru最近10 次迭代的標準偏差（STandard Deviation，STD），設置了STD<3×10-3的停止準則，然后使用SGD 優化器切換到GradNCC 相似度，直至數據集訓練結束。

3 實驗與結果分析

3.1 實驗數據

實驗數據采用中國科學院計算技術研究所與北京積水潭醫院共同發布的CTSpine1K 脊椎數據集［23］。包括1 005 個CT 樣本（超過500 000 個CT 切片和11 000 個脊塊）。該數據集將所有DICOM 圖像重新格式化為NIfTI，以簡化數據處理。

本實驗中用于訓練的X 射線圖像通過對相應的CT 數據進行DRR 投影合成得到，這樣方便記錄對應真實的姿態參數。為了能夠充分考慮實際手術中醫生之所以需要配準的原因，通過在合成的X 射線圖像上進行高斯模糊，以使合成X 射線圖像的外觀與術中X 圖像一致，從而使其訓練的回歸器可以很好地推廣到真實X射線圖像上，整體流程如圖2所示。

圖2 合成X射線圖像過程Fig.2 Synthesis process of X-ray image

在脊椎手術過程中，醫生手術的部位主要是病人具體的某一兩塊脊椎。因此考慮手術的實際情況，本文使用每一張脊椎CT 數據的12 節胸椎部位實驗，通過使用ITK-SNAP 將12 節胸椎從第一節開始，每三節分割成一部分，如圖3 所示。這樣共將人體胸椎分割為四部分，可以根據手術要求確認需配準的對應部分。實驗中，本文采用100 位病人CT 圖像作為訓練集，對人體胸椎每三節分割并進行實驗，共得400 張部分脊椎CT。對其在（-20°，20°）旋轉姿態范圍內，遵循均勻分布合成204 800 張X 射線圖像作為訓練圖像，圖像大小為128×128。訓練階段，使用對應脊椎部位的合成X 射線圖像進行訓練，每個病人的脊椎CT 對應約2 000 張合成X 射線圖像，進行10 折交叉驗證測試。

圖3 脊椎CT圖像及分割示例圖像Fig.3 Spine CT images and segmentation example images

3.2 實驗環境與配置

實驗在一臺配備Intel Core i7-4790k CPU、16 GB RAM 和NVIDIA Tesla P100 GPU 的工作站上進行。神經網絡是使用開源深度學習框架PyTorch 實現。訓練階段中，使用SGD optimizer 對網絡進行訓練，循環學習率在1E-6 和1E-4 之間，每100 步，動量為0.99，Batch size 大小選擇為8，訓練了2 000 次迭代，直到收斂。

3.3 評價指標

為了評估最終的配準精度，本文采用平均目標配準誤差（mean Target Registration Error，mTRE）以及平均絕對誤差（Mean Absolute Error，MAE）。將mTRE 小于2 mm 且MAE 小于0.1 的結果視為配準成功，配準成功的百分比定義為成功率。為了準確地評估本文方法的配準實時性，同樣在測試的過程中記錄每次配準所需時間。

3.4 實驗結果分析

本文在每位病人的CT 測試集上進行200 次配準測試，共完成20 000次配準測試。隨機展示了兩組在不同部分脊椎上的配準結果，如圖4 所示。第一列與第三列為加上高斯模糊合成X 射線圖像，第二列與第四列為對應前列的通過模型預測姿態投影生成的配準圖像，其中每一行為不同分段的部分三節脊椎。由表1 可見，同樣展示了四部分脊椎在合成X 射線圖像真實旋轉的姿態參數與預測旋轉參數上的誤差。

圖4 脊椎配準結果可視化Fig.4 Visualization of spine registration results

表1 不同部位脊椎的平均旋轉角度誤差 單位：（°）Tab.1 Average rotation angle error of spine at different parts unit：（°）

這里選取了三種基于優化的方法（Opt-GC［3］、Opt-GO［3］、Opt-NGI［4］）以及三種基于學習的方法（Bresenham［11］、Fourier［8］和MLP［14］）作為對比方法與本文模型進行了比較。

Opt-GC［3］、Opt-GO［3］、Opt-NGI［4］：三種基于優化的方法主要通過迭代評估每次DRR圖像與X射線圖像的相似性優化投影參數，從而在多次優化后獲取精度最高的配準圖像。三種方法采用針對不同特征使用不同的圖像相似度計算方法，充分考慮醫學圖像特點，其最終配準精度滿足醫生對配準的要求。但迭代配準過程所需時間難以滿足醫生對術中實時的要求。

Bresenham［11］：該方法以加快DRR 投影為主要研究方向，提出基于 Bresenham 直線生成算法改進的模式強度與梯度相結合的混合配準算法，提高投影速度。最終采用改進的鮑威爾優化算法對參數進行迭代優化，完成配準過程。

Fourier［8］：該方法基于主方向傅里葉變換算子的模板匹配初始化方法，可避免接近真值的初值需求，并顯著擴大了捕獲范圍。在確保配準精度的前提下，大幅縮短配準時間。

MLP［14］：該方法以多層感知機為主要網絡結構，通過學習提取的手工特征預測脊椎姿態，避免了多次投影迭代過程。

本文在同樣的測試數據集以及同樣設備環境中進行200次實驗，并在配準精度、配準成功率、配準耗時上進行量化比較，進一步在mTRE 和MAE 兩個指標上，將本文方法和其他對比方法進行了配對t 檢驗，結果如表2 所示。在mTRE 和MAE 兩個指標上，本文方法顯著優于MLP 方法，與其他五種對比方法沒有統計意義上顯著差異（P>0.05），可以滿足外科醫生在實際手術中的配準精度指標。在配準成功率上，本文方法明顯高于其他對比方法，體現出較好的魯棒性。在配準耗時上，本文方法優于對比方法，基本可以滿足2 s 內顯示圖像的臨床需求。因此本文方法能在保證較高配準精度的同時有效提高配準速度，滿足臨床快速配準的要求。而三種基于學習的方法（Bresenham、Fourier 和MLP）雖然在配準速度上較傳統方法有所加快，但Bresenham 與Fourier 方法仍然無法達到術中實時要求；其中MLP 回歸特征的方法雖然時間性能滿足要求，但卻以犧牲配準精度為代價。

表2 本文方法與其他方法的實驗結果對比Tab.2 Comparison of experimental results of the proposed method and other methods

本文模型通過姿態編碼器直接預測脊椎姿態，相較于基于優化的方法以及Bresenham 方法，解決了耗費大量時間的迭代優化問題，整個配準僅進行一次DRR 計算。Fourier 方法雖相較于傳統方法配準時間有大幅減少，但同樣分級配準過程仍計算量較大，對于實時配準的要求仍無法滿足。針對MLP 直接通過提取圖像特征預測姿態信息，姿態無法作為核心信息被提取的問題，本文姿態編碼器通過提取X 射線圖像中含有核心姿態信息的隱空間來預測姿態信息，從而避免圖像中其他非核心信息的干擾，保證了配準精度。

3.5 正則化參數分析

圖5 展示了在不同正則化參數γ值下的模型訓練的mTRE 測試實驗結果。本文首先在0～0.5 每隔0.01 對γ進行賦值訓練，根據圖5 折線圖中實線可知，γ在0.32～0.33 取值，模型預測結果的mTRE 值較低。因此，為得到更精確的參數值，本文繼續在0.32～0.33 以0.001 的間隔繼續對γ進行賦值訓練，圖5 中虛線結果表明：當γ取0.323 時，對Losscru訓練效果更好。

圖5 不同γ取值下的mTRE結果Fig.5 mTRE results under different γ values

3.6 消融實驗

為了展示引入的三種損失函數組合的有效性，在同樣的整個測試集上對LossL1、Losscosine、Losscru、LossNCC、LossL1+LossNCC、Losscosine+LossNCC以及本文采用的LossNCC+Losscru分別進行實驗，在配準的角度誤差以及配準精度方面記錄實驗情況。

如表3 與圖6 所示，實驗結果表明，在“粗細”結合的損失函數的約束下，模型訓練效果較好。模型測試配準結果的平均目標配準誤差僅有1.1 mm，平均絕對誤差僅有0.05。在旋轉的3 個角度的誤差上也控制在了1.1°左右。同樣在此實驗中進行t 檢驗，在最終的配準精度指標方面，LossNCC+Losscru明顯優于其他組合損失（P<0.05）。

表3 不同損失訓練的相似度結果Tab.3 Similarity results of different loss training

圖6 不同損失的旋轉角度誤差Fig.6 Rotation angle errors of different losses

4 結語

2D/3D 醫學圖像配準是脊椎骨科手術導航中的關鍵技術。針對傳統的基于迭代優化的方法無法實時配準的問題，提出一個改進的姿態編碼器的配準模型，設計了針對姿態回歸的編碼器［24］。通過該編碼器對不同姿態的X 射線圖像生成隱空間表示，提取圖像姿態特征信息，繼而回歸脊椎部位的旋轉姿態參數，僅通過一次投影即可獲得配準圖像，避免了大量迭代姿態搜索過程。在真實的CTSpine1K 數據集上進行了實驗，結果表明，本文模型在配準精度、配準時間、配準成功率指標上取得了較好的結果，配準時間優于傳統方法，在Intel Core i7-4790k CPU、16 GB RAM、NVIDIA Tesla P100 GPU 硬件條件下可以滿足術中實時配準的要求。本文目前的實驗結果僅限于脊柱部位，其他非剛體配準［25］的效果有待進一步驗證，以及如何在更低成本的嵌入式硬件平臺上實現2D/3D 的實時配準是我們下一步的研究方向。