基于深度學習的2D/3D 圖像配準方法在脊柱微創(chuàng)導航手術中的應用

邢珍珍，顏立祥，張麟華

1.太原工業(yè)學院，山西 030008；2.電子科技大學

脊柱類疾病是比較常見的一種疾病。近年來，脊柱類疾病的發(fā)病率不斷攀升。據《2005—2017 年中國疾病負擔研究報告》[1]顯示，頸椎、腰椎類疾病是導致國人健康壽命損失的首要疾病之一。手術是治療脊柱類疾病的常見方式，隨著醫(yī)學影像技術的高速發(fā)展，影像導航手術逐漸成為脊椎疾病手術的研究熱點。醫(yī)生可以結合病人不同成像模式的影像獲取病變區(qū)域的信息[2]。術中可以通過X-ray 等2D 圖像實時跟蹤手術器械和人體組織的相對位置。3D 圖像一般包括計算機斷層掃描(CT)、錐形束CT(CBCT)、磁共振成像(MRI)等[3]，可以提供更加精確的信息，但其掃描操作難，只能在術前獲取。因此，將術前3D 影像數據和術中2D影像數據置于同一坐標系中，進行圖像配準，可以在術中提供圖像引導，提高手術的成功率[4-5]。傳統(tǒng)的2D/3D 配準方法主要包括基于灰度信息的配準和基于特征信息的配準[6-7]。其中基于灰度信息的配準算法，通過優(yōu)化圖像的灰度相似度來確定圖像間的變換參數，對于存在形變或灰度變化明顯的圖像，配準效果較差[8]?；谔卣餍畔⒌呐錅仕惴?，通過提取圖像中點、線等特征完成配準，過程簡單、速度較快[9]，但是容易丟失圖像的灰度、梯度等有價值信息。深度神經網絡近年來逐漸應用于醫(yī)學圖像配準領域[10]。其中，卷積神經網絡(convolutional neural networks,CNN)廣泛用于特征點檢測[11-12]、特征描述符[13]及圖像配準[14-15]中。本研究提出一種基于卷積神經網絡的2D/3D 脊柱CT的層級配準方法，為臨床脊柱微創(chuàng)手術中的3D 導航問題提供技術支持。

1 數據與方法

1.1 實驗數據 采用開源數據集LIDC-IDRI(The Lung Image Database Consortium)，美國國家癌癥研究所對1 010 例病人進行研究，制作了1 018 個研究樣本。每個樣本包含胸部醫(yī)學CT 圖像和對應的診斷結果標注。根據需求，本實驗下載了300 個CT 圖像，每個圖像包含數量不等的DCM 格式的序列圖像。實驗使用數據集中的3D CT 圖像生成的數字重建放射影像（DRR）圖像作為浮動圖像，通過對DRR 圖像進行隨機形變來仿真手術中的X-ray 影像，并將其作為參考圖像。CT 圖像分辨率為512×512×N，投影圖像的分辨率 為256×256×1。其 中，N表 示 該CT 圖 像 包 含DCM 序列的個數。實驗所用CT 圖像見圖1。

圖1 CT 圖像

1.2 研究方法

1.2.1 配準問題描述 采用一種基于卷積神經網絡的2D/3D 脊柱CT 的層級配準方法，圖像配準需要將術前3D 圖像通過投影算法生成模擬X-ray 的2D 圖像，與術中圖像實現維度的統(tǒng)一[16]。該過程通過DRR技術實現[17]。DRR 是醫(yī)療圖像配準的一個重要的前置步驟。在進行2D/3D 配準前，需要對術前CT 圖像進行空間變換，獲取不同角度下的DRR 圖像，通過計算DRR 圖像和術中X-ray 圖像的相似性測度值來構建目標函數，并對目標函數進行優(yōu)化，獲取最優(yōu)空間變換參數，使得待配準圖像的相似性測度值達到極值。配準公式如下。

其中：Tg表示配準后的空間變換，IR是參考圖像，IM是浮動圖像，T 表示空間變換，P表示DRR 投影，函數S是IM與IR間的相似性測度。IM與IR完全配準時，IR與IM生成的DRR 圖像在像素值和空間位置上均相同，此時滿足公式(2)。

1.2.2 配準流程 本研究的配準方法包含兩個部分。首先，利用深度學習網絡進行粗配準，其次利用參數優(yōu)化方法進行精配準。配準流程：首先，將3D CT 數據集產生形變生成DRR，通過浮動圖像和參考圖像來訓練深度學習網絡。利用回歸方法學習投影參數之間的關系，將待配準的術前CT 圖像生成的DRR 和術中2D的X-ray 圖像輸入網絡，結合學到的關系和幾何變換得到相關參數，完成粗配準。其次，將待配準的CT 圖像進行手動分割，完成精配準。計算分割后浮動圖像和參考圖像之間的相似性測度，通過Adam 參數優(yōu)化算法更新參數，完成單塊配準。循環(huán)單塊精配準過程，完成多塊椎骨的配準。流程圖見圖2。

圖2 層級配準算法流程圖

1.2.3 粗配準 本研究對象為病人脊柱，屬于一種剛體器官，故通過卷積神經網絡對醫(yī)學圖像進行剛體配準[18]。卷積神經網絡可以通過卷積運算，學習到豐富的特征信息，進而對測試集進行分類、回歸等預測任務[19-20]。在剛體變換中，投影空間涉及3 個平移及3 個旋轉參數[21]。變換參數用向量T=(tx,ty,tz,rx,ry,rz)表示，其中，tx、ty、tz分別表示在X 軸、Y 軸、Z 軸上的平移參數，r、ry、rz分別表示繞X 軸、Y 軸、Z 軸的旋轉參數。通過卷積神經網絡學習剛體變換參數，實現粗配準。

1.2.3.1 特征提取 通過CNN 的卷積層和池化層提取參考圖像和浮動圖像的差值圖特征，通過全連接層來依次回歸形變參數，做回歸預測。可以通過目標區(qū)域的中心、寬度、高度和方向4 個參數(q,w,h,φ)對其進行定位，計算公式如下。

其 中：D是X 光 源 到DRR 平 面 的 距 離，w0是 目 標區(qū)域的大小。圖像配準的實質是預測變化參數，參數之間的變化用P 表示，公式如下，其中f(·)表示回歸算子，X(·)表示特征提取算子。

1.2.3.2 粗配準網絡結構 搭建CNN 深度學習網絡模型，如圖3 所示，共有8 層網絡結構。輸入層是大小為256×256 的浮動圖像和參考圖像，卷積層的卷積核為20 個，大小為5×5，不填充，步長為1。池化層采用大小為2×2 的最大池化，步長為2。采用ReLU 激活函數，對提取的特征進行非線性轉換。全連接層1 有1 024 個激活函數單元，輸出結點個數為1 024。全連接層2 輸出結點個數為128 個。最后1 層為輸出層，輸出結點個數為N個。

圖3 粗配準網絡結構示意圖

1.2.3.3 參數回歸 若同時回歸參數向量T的6 個形變參數，容易使得訓練陷入局部最優(yōu)[22]，且回歸精度較差。本研究采用按順序回歸的方法，首先回歸(tx,ty)，將圖3 網絡中最后全連接層的激活函數單元設置為2，進行回歸預測，更新差值圖的(tx,ty)，作為下一步網絡的輸入。用同樣的方法依次回歸(rz)，(tz)，(rx,ry)。最后進行組合，得到形變參數(tx,ty,tz,rx,ry,rz)。

采用均方誤差函數為損失函數，相對平移誤差，旋轉誤差更為明顯，所以，在損失函數中添加對旋轉參數的權重，損失函數公式如下。

每張浮動圖像與對應的參考圖像組合，構成一個訓練樣本。將訓練樣本依次輸入網絡模型，得到浮動圖像與參考圖像的差值圖像，如圖4 所示。通過學習差值圖的特征信息，輸出待配準圖像的形變參數。

圖4 參考圖像、浮動圖像和差值圖像

1.2.4 精配準 精配準主要是對粗配準得到的形變參數進行優(yōu)化，對椎骨影像進行分塊配準。采用手動分割的方式，通過矩形框標記，分割的每幅子圖僅包含1 塊椎骨。如圖5 所示方框區(qū)域所示。

圖5 圖像分割示意圖

將單塊椎骨圖生成DRR 圖像，作為浮動圖像，計算浮動圖像與參考圖像之間的Dice Loss 值，該值若小于預設閾值，則完成單塊椎骨圖的精配準，若不滿足條件，將Dice Loss 設置為Adam 參數優(yōu)化算法的目標函數，重復上述步驟，搜尋Dice Loss 值最小時的最優(yōu)精配準參數，完成單塊椎骨圖的配準。按同樣步驟完成所有單塊椎骨圖的精配準，從而實現脊柱CT 層級配準。

1.2.5 配準實驗

1.2.5.1 實驗一：單個對象的多形變配準 首先以數據集LIDC-IDRI 的第2 張CT 為實驗對象，模擬生成10種形變。采用Elastix的配準算法作為對照組，Elastix廣泛應用于醫(yī)學圖像的配準[23]。對于配準參數（tx，ty，tz，rx，ry，rz），令旋轉參數為10°，平移參數為20 mm。每次1 個參數形變時，選取3 組實驗，每次兩個參數進行形變時，選取3 組來簡化實驗。每次3 個參數形變，同樣選取3 組來簡化實驗。最后選取每個參數都有變化的1 組。共10 組。數據組形變參數（rx，ry，rz，tx，ty，tz）分別為(10,0,0,0,0,0)、（0,10,0,0,0,0）、（0，0，0，0，0，20）、（10，0，0，0，20，0）、（0,0,0,20,0,20）、（0,10,10,0,0,0）、（10，10，10,0,0,0）、（10,0,10,0,20,0）、（0,0,0,20,20,20）、（10,10,10,20,20,20）。按上述參數產生形變得到的DRR圖像為浮動圖像，（rx，ry，rz，tx，ty，tz）=（0,0,0,0,0,0）得到的DRR 圖像為參考圖像。

將差值圖歸一化后輸入CNN 進行訓練。劃分100 對樣本為訓練集，10 對樣本為測試集。權重參數使用隨機梯度下降（SGD）方法進行學習和優(yōu)化，設置權重衰減系數為d=0.000 1，動量m=0.9。 設置batch_size=10，epoch_num=20，iteration=10。 其 中學習率策略如式(7)所示。當損失函數值小于0.01 或達到迭代次數時，停止訓練模型并對其進行保存，用于參數的預測。

其中：ω表示需要更新的參數，t表示迭代步數，α是學習率，St是梯度累積變量，S初始值為0，dx是損失函數對于ω 的梯度，β是衰減系數。設置α=0.001，ε=10-6，β=0.9，完成圖像的精配準。將分割區(qū)域圖像的像素矩陣值與原圖融合，得到最終的配準結果圖。

1.2.5.2 實驗二：多個對象的形變配準 上述實驗對數據集中第2 個病例的CT 影像模擬了10 種形變。為了證明DLIR 的魯棒性，本研究使用數據集中沒有參與過訓練和測試的第201 張到第220 張CT 圖像進行研究。設置這20 張CT 的形變參數(tx,ty,tz,rx,ry,rz)為(20,0,20,0,10,0)，同時采用Elastix 的配準算法作為對照組。

2 結果

2.1 評價標準 本研究用3 個配準精度指標對結果進行評價。以N×M 大小的兩幅圖像I1,I2為例，闡述配準評價指標。

2.1.1 平均絕對誤差（MAE） MAE 是直接衡量配準所得的參數和真值參數之間的誤差大小，能夠避免誤差相互抵消。MAE 的取值范圍為[0,+∞)，其值越小，說明模型預測越準確。公式如下：

其中：Treg(i)表示空間變換的第i個參數，Ttruth(i)表示空間變換真值的第i個參數，d表示維度。

2.1.2 歸一化互相關（NCC） NCC 用于歸一化待匹配目標原始像素間的相關度。NCC取值范圍為[-1,1]，若為-1 表示兩個樣本完全不相關，若為1 表示二者有很高的相關性。計算公式如下，其中，D(I)是圖像I的方差，σ(I1,I2)是圖像I1,I2的協(xié)方差。

聯合熵H(I1,I2)計算如下，其中，p ( i1,i2 )為兩圖聯合概率，Ni表示灰度值是i 的像素點個數，Ni1,i2為兩圖灰度值分別為i1,i2的像素對個數。

2.2 實驗結果 采用配準結果圖對實驗做定性分析，通過MAE、NCC 和NMI 評價標準對實驗做定量分析。為了方便起見，本研究所提深度學習圖像配準(deep learning image registration)算法用“DLIR”來表示，基于Elastix 的配準(elastix image registration)算法用“ExIR”來表示。在實驗一的10 組實驗中，從6 個配準參數形變1 個、2 個、3 個和6 個的4 組實驗中分別選取一組配準實驗結果圖，顯示配準實驗的浮動圖像、參考圖像、基于DLIR 和ExIR 的配準結果。見圖6。

圖6 實驗配準結果

從配準效果圖可以看出配準后的結果圖基本一致，下面對實驗結果進行定量分析，計算配準后的圖像與參考圖像之間的MAE、NCC 和NMI，數據見表1。

表1 配準實驗數據結果

實驗二對數據集第201 張～第220 張共20 張CT影像進行配準。平均實驗參數結果見表2。

表2 20 張CT 的平均實驗參數

3 討論

本研究提出了一種基于深度學習和參數優(yōu)化的2D/3D 層級配準方法，通過對某一病人的CT 圖像模擬生成了10 種形變，將基于Elastix 的配準方法作為基準實驗組，10 組數據中，MAE 指標下，DLIR 有10 組優(yōu)于ExIR；NCC 指標下，DLIR 有7 組優(yōu)于ExIR；NMI 指標下，DLIR 有10 組優(yōu)于ExIR。根據這10 組實驗結果，DLIR 組MAE、NCC 和NMI 的樣本平均值分別為0.039 8，0.834 7 和0.629 9，ExIR 組的上述3 項指標分別為0.135 1，0.643 2 和0.184 6，3 項指標上，DLIR 均優(yōu)于ExIR。另外，使用LIDI-IDRI 數據集中的20 張CT 的配準實驗結果顯示，表明提出的方法在MAE 指標上降低了0.097 8，在NCC 指標上提高了0.185 6，在NMI 指標上提高了0.445 6，表明本研究提出的方法3 項指標均優(yōu)于ExIR 方法。本研究所提方法既考慮到了脊柱剛體變換的全局信息，又考慮到了兩塊椎骨之間的局部信息，有效提高了配準精度。說明將深度學習應用于脊柱微創(chuàng)導航手術的CT 圖像配準中具有可行性，是一項非常具有前景的研究。