邢珍珍，顏立祥，張麟華

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

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

1 數(shù)據(jù)與方法

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

圖1 CT 圖像

1.2 研究方法

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

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

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

圖2 層級(jí)配準(zhǔn)算法流程圖

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

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

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

1.2.3.2 粗配準(zhǔn)網(wǎng)絡(luò)結(jié)構(gòu) 搭建CNN 深度學(xué)習(xí)網(wǎng)絡(luò)模型，如圖3 所示，共有8 層網(wǎng)絡(luò)結(jié)構(gòu)。輸入層是大小為256×256 的浮動(dòng)圖像和參考圖像，卷積層的卷積核為20 個(gè)，大小為5×5，不填充，步長為1。池化層采用大小為2×2 的最大池化，步長為2。采用ReLU 激活函數(shù)，對(duì)提取的特征進(jìn)行非線性轉(zhuǎn)換。全連接層1 有1 024 個(gè)激活函數(shù)單元，輸出結(jié)點(diǎn)個(gè)數(shù)為1 024。全連接層2 輸出結(jié)點(diǎn)個(gè)數(shù)為128 個(gè)。最后1 層為輸出層，輸出結(jié)點(diǎn)個(gè)數(shù)為N個(gè)。

圖3 粗配準(zhǔn)網(wǎng)絡(luò)結(jié)構(gòu)示意圖

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

采用均方誤差函數(shù)為損失函數(shù)，相對(duì)平移誤差，旋轉(zhuǎn)誤差更為明顯，所以，在損失函數(shù)中添加對(duì)旋轉(zhuǎn)參數(shù)的權(quán)重，損失函數(shù)公式如下。

每張浮動(dòng)圖像與對(duì)應(yīng)的參考圖像組合，構(gòu)成一個(gè)訓(xùn)練樣本。將訓(xùn)練樣本依次輸入網(wǎng)絡(luò)模型，得到浮動(dòng)圖像與參考圖像的差值圖像，如圖4 所示。通過學(xué)習(xí)差值圖的特征信息，輸出待配準(zhǔn)圖像的形變參數(shù)。

圖4 參考圖像、浮動(dòng)圖像和差值圖像

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

圖5 圖像分割示意圖

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

1.2.5 配準(zhǔn)實(shí)驗(yàn)

1.2.5.1 實(shí)驗(yàn)一：單個(gè)對(duì)象的多形變配準(zhǔn) 首先以數(shù)據(jù)集LIDC-IDRI 的第2 張CT 為實(shí)驗(yàn)對(duì)象，模擬生成10種形變。采用Elastix的配準(zhǔn)算法作為對(duì)照組，Elastix廣泛應(yīng)用于醫(yī)學(xué)圖像的配準(zhǔn)[23]。對(duì)于配準(zhǔn)參數(shù)（tx，ty，tz，rx，ry，rz），令旋轉(zhuǎn)參數(shù)為10°，平移參數(shù)為20 mm。每次1 個(gè)參數(shù)形變時(shí)，選取3 組實(shí)驗(yàn)，每次兩個(gè)參數(shù)進(jìn)行形變時(shí)，選取3 組來簡化實(shí)驗(yàn)。每次3 個(gè)參數(shù)形變，同樣選取3 組來簡化實(shí)驗(yàn)。最后選取每個(gè)參數(shù)都有變化的1 組。共10 組。數(shù)據(jù)組形變參數(shù)（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）。按上述參數(shù)產(chǎn)生形變得到的DRR圖像為浮動(dòng)圖像，（rx，ry，rz，tx，ty，tz）=（0,0,0,0,0,0）得到的DRR 圖像為參考圖像。

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

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

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

2 結(jié)果

2.1 評(píng)價(jià)標(biāo)準(zhǔn) 本研究用3 個(gè)配準(zhǔn)精度指標(biāo)對(duì)結(jié)果進(jìn)行評(píng)價(jià)。以N×M 大小的兩幅圖像I1,I2為例，闡述配準(zhǔn)評(píng)價(jià)指標(biāo)。

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

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

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

聯(lián)合熵H(I1,I2)計(jì)算如下，其中，p ( i1,i2 )為兩圖聯(lián)合概率，Ni表示灰度值是i 的像素點(diǎn)個(gè)數(shù)，Ni1,i2為兩圖灰度值分別為i1,i2的像素對(duì)個(gè)數(shù)。

2.2 實(shí)驗(yàn)結(jié)果 采用配準(zhǔn)結(jié)果圖對(duì)實(shí)驗(yàn)做定性分析，通過MAE、NCC 和NMI 評(píng)價(jià)標(biāo)準(zhǔn)對(duì)實(shí)驗(yàn)做定量分析。為了方便起見，本研究所提深度學(xué)習(xí)圖像配準(zhǔn)(deep learning image registration)算法用“DLIR”來表示，基于Elastix 的配準(zhǔn)(elastix image registration)算法用“ExIR”來表示。在實(shí)驗(yàn)一的10 組實(shí)驗(yàn)中，從6 個(gè)配準(zhǔn)參數(shù)形變1 個(gè)、2 個(gè)、3 個(gè)和6 個(gè)的4 組實(shí)驗(yàn)中分別選取一組配準(zhǔn)實(shí)驗(yàn)結(jié)果圖，顯示配準(zhǔn)實(shí)驗(yàn)的浮動(dòng)圖像、參考圖像、基于DLIR 和ExIR 的配準(zhǔn)結(jié)果。見圖6。

圖6 實(shí)驗(yàn)配準(zhǔn)結(jié)果

從配準(zhǔn)效果圖可以看出配準(zhǔn)后的結(jié)果圖基本一致，下面對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行定量分析，計(jì)算配準(zhǔn)后的圖像與參考圖像之間的MAE、NCC 和NMI，數(shù)據(jù)見表1。

表1 配準(zhǔn)實(shí)驗(yàn)數(shù)據(jù)結(jié)果

實(shí)驗(yàn)二對(duì)數(shù)據(jù)集第201 張～第220 張共20 張CT影像進(jìn)行配準(zhǔn)。平均實(shí)驗(yàn)參數(shù)結(jié)果見表2。

表2 20 張CT 的平均實(shí)驗(yàn)參數(shù)

3 討論

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