徐少康，張戰(zhàn)成，2*，姚浩男，鄒智偉，張寶成

（1.蘇州科技大學(xué) 電子與信息工程學(xué)院，江蘇 蘇州 215009；2.蘇州科技大學(xué) 蘇州市虛擬現(xiàn)實(shí)智能交互及應(yīng)用技術(shù)重點(diǎn)實(shí)驗(yàn)室，江蘇 蘇州 215009；3.中部戰(zhàn)區(qū)總醫(yī)院 骨科，武漢 430070）

0 引言

2D/3D 醫(yī)學(xué)圖像配準(zhǔn)是三維手術(shù)導(dǎo)航系統(tǒng)中的關(guān)鍵技術(shù)之一，在脊椎手術(shù)中，椎弓根螺釘植入是臨床脊柱外科手術(shù)中最常用的內(nèi)固定方法，術(shù)中精確植入螺釘?shù)耐瑫r(shí)，不損傷血管和神經(jīng)是一項(xiàng)挑戰(zhàn)性的工作［1］。醫(yī)生術(shù)中通過(guò)X 射線圖像（2D）對(duì)病人脊椎位置進(jìn)行分析，查看螺釘?shù)闹踩肭闆r。但2D 圖像的清晰度能夠反映出的信息較少，想要精確地分析出植入螺釘位置信息較為困難。而3D 圖像可以提供此類的多維信息，如果在手術(shù)過(guò)程中多次實(shí)施CT 成像會(huì)對(duì)醫(yī)生以及病人造成傷害，將術(shù)前3D 浮動(dòng)圖像對(duì)照術(shù)中2D 參考圖像進(jìn)行空間轉(zhuǎn)換達(dá)到坐標(biāo)系相統(tǒng)一，從3D 浮動(dòng)圖像中引導(dǎo)出對(duì)應(yīng)位置的清晰的2D 圖像，即可更精確地引導(dǎo)手術(shù)，提高植入椎弓根螺釘?shù)木_度與外科醫(yī)生的手術(shù)效率，同時(shí)減少醫(yī)生和病人的輻射傷害。

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

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

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

為了提高2D/3D醫(yī)學(xué)圖像配準(zhǔn)的速度，達(dá)到術(shù)中實(shí)時(shí)配準(zhǔn)的要求，受隱空間回歸啟發(fā)，將2D/3D 醫(yī)學(xué)圖像配準(zhǔn)問(wèn)題表達(dá)為基于自編碼的回歸網(wǎng)絡(luò)，并加入隱空間回歸約束，以捕獲X射線圖像中姿態(tài)信息，從而精確地回歸脊椎的旋轉(zhuǎn)的姿態(tài)參數(shù)。對(duì)比了傳統(tǒng)的優(yōu)化迭代配準(zhǔn)方法，本文方法避免了大量的DRR渲染時(shí)間，使得2D/3D圖像配準(zhǔn)的時(shí)間滿足實(shí)時(shí)的要求。

1 相關(guān)工作

1.1 DRR投影模型

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

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

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

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

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

1.2 三維轉(zhuǎn)換參數(shù)

一個(gè)剛體三維變換可以用一個(gè)含有6 個(gè)參數(shù)的矢量來(lái)參數(shù)化［10］。在本文方法中，通過(guò)3個(gè)平面內(nèi)和3個(gè)平面外的轉(zhuǎn)換參數(shù)來(lái)參數(shù)化這個(gè)三維變換。其中，平面內(nèi)的變換參數(shù)包含2個(gè)平移參數(shù)tx、ty和1個(gè)旋轉(zhuǎn)參數(shù)rz。平面內(nèi)的變換參數(shù)的結(jié)果看作為二維剛體變換。平面外的變換參數(shù)包括1 個(gè)面外平移參數(shù)tz和2 個(gè)面外旋轉(zhuǎn)參數(shù)rx和ry。平面外平移和旋轉(zhuǎn)的結(jié)果分別對(duì)應(yīng)著剛體圖像的縮放和形狀變化。

在脊椎手術(shù)中，為了能夠提示外科醫(yī)生在鉆探椎弓根螺釘導(dǎo)孔時(shí)可能引發(fā)的椎體皮質(zhì)骨折，通常會(huì)在植入椎弓根螺釘前向脊椎內(nèi)植入一枚探針用于探查脊椎管道情況，再選擇合適的椎弓根螺釘進(jìn)行植入。在進(jìn)行脊椎圖像配準(zhǔn)前，外科醫(yī)生可以直接通過(guò)探針植入的位置準(zhǔn)確地知道需要配準(zhǔn)的是第幾節(jié)脊椎位置，所以對(duì)于X射線圖像的固定平移拍攝，平移的三個(gè)參數(shù)是已知的，只需要計(jì)算出旋轉(zhuǎn)姿態(tài)參數(shù)rx、ry、rz即可完成配準(zhǔn)。

2 脊椎姿態(tài)回歸方法

2D/3D 配準(zhǔn)的目標(biāo)是搜索CT 圖像的6 個(gè)未知自由度姿態(tài)，以該姿態(tài)對(duì)CT 體積模擬投射X 射線在平面上產(chǎn)生DRR與X 射線圖像進(jìn)行比較。因此2D/3D 配準(zhǔn)的問(wèn)題最終可以歸結(jié)到如何從X射線圖像中尋找最佳CT投影姿態(tài)的問(wèn)題。

2.1 配準(zhǔn)模型結(jié)構(gòu)

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

圖1 配準(zhǔn)網(wǎng)絡(luò)框架Fig.1 Framework of registration network

式（5）中：ICT為CT圖像，IX為X射線圖像，φ為Encoder+MLP姿態(tài)回歸過(guò)程模型，Proj 為DRR 投影，Ireg為配準(zhǔn)圖像。編碼器用F表示，該網(wǎng)絡(luò)將原始數(shù)據(jù)X映射到隱空間F中(φ：X→F)。通過(guò)用激活函數(shù)σ傳遞的標(biāo)準(zhǔn)神經(jīng)網(wǎng)絡(luò)函數(shù)表示，其中Z是隱變量。

2.2 隱空間生成

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

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

2.3 損失函數(shù)

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

其中γ是一個(gè)正則化參數(shù)（在3.5節(jié)中會(huì)對(duì)其進(jìn)行實(shí)驗(yàn)分析取值）。實(shí)驗(yàn)表明當(dāng)初始姿態(tài)與真實(shí)姿態(tài)差異較大時(shí)，式（9）損失函數(shù)在收斂過(guò)程中更有效，但當(dāng)DRR 圖像與X 射線圖像逐漸一致時(shí)，該損失函數(shù)對(duì)細(xì)微的差異敏感性變低。

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

其中：IXI與IDRR分別對(duì)應(yīng)X 射線圖像與DRR 圖像在區(qū)域(i，j)的強(qiáng)度值，與為重疊區(qū)域(i，j) ∈T內(nèi)圖像的均值。

本文使用隨機(jī)梯度下降（Stochastic Gradient Descent，SGD）優(yōu)化器進(jìn)行優(yōu)化，計(jì)算了Losscru最近10 次迭代的標(biāo)準(zhǔn)偏差（STandard Deviation，STD），設(shè)置了STD<3×10-3的停止準(zhǔn)則，然后使用SGD 優(yōu)化器切換到GradNCC 相似度，直至數(shù)據(jù)集訓(xùn)練結(jié)束。

3 實(shí)驗(yàn)與結(jié)果分析

3.1 實(shí)驗(yàn)數(shù)據(jù)

實(shí)驗(yàn)數(shù)據(jù)采用中國(guó)科學(xué)院計(jì)算技術(shù)研究所與北京積水潭醫(yī)院共同發(fā)布的CTSpine1K 脊椎數(shù)據(jù)集［23］。包括1 005 個(gè)CT 樣本（超過(guò)500 000 個(gè)CT 切片和11 000 個(gè)脊塊）。該數(shù)據(jù)集將所有DICOM 圖像重新格式化為NIfTI，以簡(jiǎn)化數(shù)據(jù)處理。

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

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

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

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

3.2 實(shí)驗(yàn)環(huán)境與配置

實(shí)驗(yàn)在一臺(tái)配備Intel Core i7-4790k CPU、16 GB RAM 和NVIDIA Tesla P100 GPU 的工作站上進(jìn)行。神經(jīng)網(wǎng)絡(luò)是使用開(kāi)源深度學(xué)習(xí)框架PyTorch 實(shí)現(xiàn)。訓(xùn)練階段中，使用SGD optimizer 對(duì)網(wǎng)絡(luò)進(jìn)行訓(xùn)練，循環(huán)學(xué)習(xí)率在1E-6 和1E-4 之間，每100 步，動(dòng)量為0.99，Batch size 大小選擇為8，訓(xùn)練了2 000 次迭代，直到收斂。

3.3 評(píng)價(jià)指標(biāo)

為了評(píng)估最終的配準(zhǔn)精度，本文采用平均目標(biāo)配準(zhǔn)誤差（mean Target Registration Error，mTRE）以及平均絕對(duì)誤差（Mean Absolute Error，MAE）。將mTRE 小于2 mm 且MAE 小于0.1 的結(jié)果視為配準(zhǔn)成功，配準(zhǔn)成功的百分比定義為成功率。為了準(zhǔn)確地評(píng)估本文方法的配準(zhǔn)實(shí)時(shí)性，同樣在測(cè)試的過(guò)程中記錄每次配準(zhǔn)所需時(shí)間。

3.4 實(shí)驗(yàn)結(jié)果分析

本文在每位病人的CT 測(cè)試集上進(jìn)行200 次配準(zhǔn)測(cè)試，共完成20 000次配準(zhǔn)測(cè)試。隨機(jī)展示了兩組在不同部分脊椎上的配準(zhǔn)結(jié)果，如圖4 所示。第一列與第三列為加上高斯模糊合成X 射線圖像，第二列與第四列為對(duì)應(yīng)前列的通過(guò)模型預(yù)測(cè)姿態(tài)投影生成的配準(zhǔn)圖像，其中每一行為不同分段的部分三節(jié)脊椎。由表1 可見(jiàn)，同樣展示了四部分脊椎在合成X 射線圖像真實(shí)旋轉(zhuǎn)的姿態(tài)參數(shù)與預(yù)測(cè)旋轉(zhuǎn)參數(shù)上的誤差。

圖4 脊椎配準(zhǔn)結(jié)果可視化Fig.4 Visualization of spine registration results

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

這里選取了三種基于優(yōu)化的方法（Opt-GC［3］、Opt-GO［3］、Opt-NGI［4］）以及三種基于學(xué)習(xí)的方法（Bresenham［11］、Fourier［8］和MLP［14］）作為對(duì)比方法與本文模型進(jìn)行了比較。

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

Bresenham［11］：該方法以加快DRR 投影為主要研究方向，提出基于 Bresenham 直線生成算法改進(jìn)的模式強(qiáng)度與梯度相結(jié)合的混合配準(zhǔn)算法，提高投影速度。最終采用改進(jìn)的鮑威爾優(yōu)化算法對(duì)參數(shù)進(jìn)行迭代優(yōu)化，完成配準(zhǔn)過(guò)程。

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

MLP［14］：該方法以多層感知機(jī)為主要網(wǎng)絡(luò)結(jié)構(gòu)，通過(guò)學(xué)習(xí)提取的手工特征預(yù)測(cè)脊椎姿態(tài)，避免了多次投影迭代過(guò)程。

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

表2 本文方法與其他方法的實(shí)驗(yàn)結(jié)果對(duì)比Tab.2 Comparison of experimental results of the proposed method and other methods

本文模型通過(guò)姿態(tài)編碼器直接預(yù)測(cè)脊椎姿態(tài)，相較于基于優(yōu)化的方法以及Bresenham 方法，解決了耗費(fèi)大量時(shí)間的迭代優(yōu)化問(wèn)題，整個(gè)配準(zhǔn)僅進(jìn)行一次DRR 計(jì)算。Fourier 方法雖相較于傳統(tǒng)方法配準(zhǔn)時(shí)間有大幅減少，但同樣分級(jí)配準(zhǔn)過(guò)程仍計(jì)算量較大，對(duì)于實(shí)時(shí)配準(zhǔn)的要求仍無(wú)法滿足。針對(duì)MLP 直接通過(guò)提取圖像特征預(yù)測(cè)姿態(tài)信息，姿態(tài)無(wú)法作為核心信息被提取的問(wèn)題，本文姿態(tài)編碼器通過(guò)提取X 射線圖像中含有核心姿態(tài)信息的隱空間來(lái)預(yù)測(cè)姿態(tài)信息，從而避免圖像中其他非核心信息的干擾，保證了配準(zhǔn)精度。

3.5 正則化參數(shù)分析

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

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

3.6 消融實(shí)驗(yàn)

為了展示引入的三種損失函數(shù)組合的有效性，在同樣的整個(gè)測(cè)試集上對(duì)LossL1、Losscosine、Losscru、LossNCC、LossL1+LossNCC、Losscosine+LossNCC以及本文采用的LossNCC+Losscru分別進(jìn)行實(shí)驗(yàn)，在配準(zhǔn)的角度誤差以及配準(zhǔn)精度方面記錄實(shí)驗(yàn)情況。

如表3 與圖6 所示，實(shí)驗(yàn)結(jié)果表明，在“粗細(xì)”結(jié)合的損失函數(shù)的約束下，模型訓(xùn)練效果較好。模型測(cè)試配準(zhǔn)結(jié)果的平均目標(biāo)配準(zhǔn)誤差僅有1.1 mm，平均絕對(duì)誤差僅有0.05。在旋轉(zhuǎn)的3 個(gè)角度的誤差上也控制在了1.1°左右。同樣在此實(shí)驗(yàn)中進(jìn)行t 檢驗(yàn)，在最終的配準(zhǔn)精度指標(biāo)方面，LossNCC+Losscru明顯優(yōu)于其他組合損失（P<0.05）。

表3 不同損失訓(xùn)練的相似度結(jié)果Tab.3 Similarity results of different loss training

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

4 結(jié)語(yǔ)

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