基于神經(jīng)網(wǎng)絡(luò)的雙X射線影像2D-3D配準(zhǔn)算法

沈延延，馮漢升

1.安徽大學(xué)電氣工程與自動(dòng)化學(xué)院，安徽合肥230601；2.中國科學(xué)院等離子體物理研究所，安徽合肥230031

前言

現(xiàn)代醫(yī)學(xué)中，放療是惡性腫瘤的重要治療手段。但放射線在轟擊病灶的同時(shí)也會(huì)損害腫瘤組織周圍的正常細(xì)胞，因此需要在治療過程中使用圖像引導(dǎo)放療（Ⅰmage-Guided Radiotherapy, ⅠGRT）技術(shù)保證放療精確程度，從而提高放療對(duì)腫瘤細(xì)胞的殺滅效果，并降低副作用。2D-3D 醫(yī)學(xué)圖像配準(zhǔn)技術(shù)可以進(jìn)行患者擺位誤差的計(jì)算，是極其重要的一環(huán)。

傳統(tǒng)基于灰度的2D-3D 配準(zhǔn)算法通常采用迭代優(yōu)化的方式，通過搜索CT 的空間姿態(tài)并進(jìn)行投影生成數(shù)字重建放射影像（Digital Reconstruction Radiograph,DRR），對(duì)比DRR圖像與治療中獲取的X射線影像的相似度來判斷當(dāng)前CT 姿態(tài)與患者姿態(tài)間的差異［1］。由于DRR 圖像的生成非常耗時(shí)，導(dǎo)致基于灰度的2D-3D 配準(zhǔn)算法較慢。針對(duì)這一問題，Mu［2］采取塊投影法，利用已有的DRR 圖像幫助生成新的DRR 圖像，減少DRR 圖像生成的耗時(shí)。Tornai等［3］則將GPU 并行計(jì)算應(yīng)用到2D-3D 配準(zhǔn)中加快計(jì)算過程。而Pan 等［4］采取閾值分割來排除不感興趣區(qū)域，減小計(jì)算數(shù)據(jù)的大小。Lei 等［5］則利用一對(duì)正交X 射線影像上兩個(gè)單面板的粗配準(zhǔn)結(jié)合作為2D-3D 配準(zhǔn)初始化減小迭代次數(shù)。Ghafurian等［6］提出一種新的圖像梯度概率密度直方圖作為圖像的特征，避免搜索優(yōu)化6 個(gè)自由度的變換向量。晉青鵬等［7］分析解剖特征與DRR 圖像變化的關(guān)系進(jìn)行非剛性配準(zhǔn)。同時(shí)隨著深度學(xué)習(xí)的發(fā)展，Miao 等［8-10］最先提出并改進(jìn)使用卷積神經(jīng)網(wǎng)絡(luò)（Convolutional Neural Network,CNN）從圖像殘差中計(jì)算CT 空間剛體變換參數(shù)。Pinheiro 等［11］則是將CNN 回歸與傳統(tǒng)方法結(jié)合到一起，迭代地對(duì)新生成的DRR 圖像使用CNN 回歸。Pei 等［12］同樣使用CNN 并采取多尺度特征融合的方法綜合提取局部和總體特征進(jìn)行配準(zhǔn)。

本研究提出一種基于CNN的雙X射線影像配準(zhǔn)方法，利用兩個(gè)具有一定夾角的X射線影像在平面上的配準(zhǔn)，將空間剛體變換參數(shù)的一部分轉(zhuǎn)化為平面上的近似剛體變換，同時(shí)利用3個(gè)神經(jīng)網(wǎng)絡(luò)分組進(jìn)行回歸預(yù)測(cè)的方法簡(jiǎn)化2D-3D配準(zhǔn)，從而達(dá)到快速進(jìn)行2D-3D配準(zhǔn)的目的。

1 方法

1.1 幾何分解

模擬拍攝雙X射線影像時(shí)的機(jī)械結(jié)構(gòu)，建立如下空間坐標(biāo)系O 及兩個(gè)X 射線影像所在的平面坐標(biāo)系O1、O2。其中射線源與射線間繞Y 軸的夾角為θre，且該兩個(gè)射線都平行于平面X-O-Z。幾何示意如圖1所示。

圖1 空間坐標(biāo)系與平面坐標(biāo)系Fig.1 Spatial and plane coordinate systems

根據(jù)該幾何關(guān)系可將空間剛體變換的6 個(gè)參數(shù)中的5 個(gè)（包括沿X、Y、Z 軸的平移Tx、Ty、Tz，以及沿X、Y 軸的旋轉(zhuǎn)Rx、Ry）分解到兩個(gè)X 射線影像所在平面坐標(biāo)系中，反之也可由兩個(gè)平面坐標(biāo)系計(jì)算空間坐標(biāo)系的Tx、Ty、Tz、Rx、Ry。

首先，從空間參數(shù)分解到平面參數(shù)：

然后，從平面參數(shù)反演出除Rz以外的5 個(gè)空間參數(shù)Tx、Ty、Tz、Rx、Ry：

射線源到DRR平面的距離即源像距（SⅠD），射線源到等中心點(diǎn)的距離即源軸距（SAD），上述公式中放縮系數(shù)scale1與scale2的計(jì)算如下：

在該坐標(biāo)系下，模擬拍攝一對(duì)夾角為θre的X 射線影像，生成兩幅DRR圖像，由坐標(biāo)系分解可以得到在兩個(gè)平面上的X-DRR 圖像對(duì)配準(zhǔn)結(jié)果，計(jì)算出空間參數(shù)的Tx、Ty、Tz、Rx、Ry。

1.2 分步配準(zhǔn)

對(duì)應(yīng)射線下的X-DRR圖像對(duì)之間的差異主要由CT 姿態(tài)的空間變換造成，因此從X-DRR 圖像對(duì)的殘差到空間的姿態(tài)變換是存在著一個(gè)高度非線性的映射。而神經(jīng)網(wǎng)絡(luò)在理論上是可以無限逼近任意非線性函數(shù)的，因此本研究采用神經(jīng)網(wǎng)絡(luò)對(duì)空間剛體變換參數(shù)進(jìn)行回歸預(yù)測(cè)。為簡(jiǎn)化配準(zhǔn)問題的復(fù)雜程度，將配準(zhǔn)分解成兩個(gè)步驟進(jìn)行。由前文分析可知?jiǎng)傮w變換參數(shù)中的5個(gè)組成元素反映在平面變換中，根據(jù)參數(shù)類型的不同將其分為兩組分別進(jìn)行回歸，包含平面內(nèi)參數(shù)Pin={tx1,tx2,ty1,ty2,r1,r2} ，平面外參數(shù)Pout={R2} 。

由于參數(shù)Tx、Ty、Tz、Rx、Ry并非獨(dú)立作用在兩個(gè)平面上，其中平移參數(shù)在一個(gè)平面上表現(xiàn)為平移，在另一個(gè)平面上為放縮。例如Tx在平面坐標(biāo)系o1上表現(xiàn)為平移tx1，在平面坐標(biāo)系o2上表現(xiàn)為放縮系數(shù)scale2。旋轉(zhuǎn)參數(shù)在一個(gè)DRR 圖像上表現(xiàn)為旋轉(zhuǎn)，同時(shí)使得射線源到另一個(gè)DRR 圖像像素點(diǎn)的光路所經(jīng)過的CT體素點(diǎn)發(fā)生改變，進(jìn)而造成對(duì)應(yīng)DRR像素值改變以及圖像輪廓的形變（如Rx與r2以及DRR1）。因此X-DRR圖像對(duì)之間是一種非線性變換。當(dāng)初始位置到目標(biāo)位置的差值較小，即算法收斂范圍較小時(shí)可近似成剛體變換。Pin是兩個(gè)平面上X-DRR 圖像對(duì)之間的近似剛體變換參數(shù)。

按難易順序先后回歸Pin和Pout，并且根據(jù)Pin回歸結(jié)果對(duì)DRR 圖像進(jìn)行對(duì)應(yīng)的剛體變換，進(jìn)而盡可能地消除已知參數(shù)對(duì)圖像的影響，使得在后續(xù)網(wǎng)絡(luò)中僅包含未知參數(shù)造成的圖像差異。其算法流程如圖2所示。

圖2 分步回歸Fig.2 Step-by-step regression

1.3 網(wǎng)絡(luò)結(jié)構(gòu)

在第1.2節(jié)中使用了3個(gè)獨(dú)立訓(xùn)練的CNN回歸網(wǎng)絡(luò)。其中網(wǎng)絡(luò)1、2的目的分別是回歸兩個(gè)平面上的參數(shù)，因此它們的網(wǎng)絡(luò)框架完全一致。輸入為X-DRR圖像對(duì)的殘差圖像，為了減小輸入圖像大小，將512×512的殘差圖像分割為互有重疊部分的4個(gè)258×258的區(qū)域，重疊在一起后輸入一個(gè)僅少量卷積層組成的擴(kuò)展網(wǎng)絡(luò)CNN-Pre。將擴(kuò)增的特征圖輸入由Dense Block［13］組成的特征提取網(wǎng)絡(luò)，并且每一個(gè)Dense Block的輸入都是由之前所有的Dense Block的輸出經(jīng)過不同的池化層將特征尺寸統(tǒng)一后連接再經(jīng)過1×1卷積層的特征融合操作得到。最后將提取到的所有特征進(jìn)行全局平均池化，再把池化結(jié)果輸入全連接層，輸出患者擺位誤差。網(wǎng)絡(luò)結(jié)構(gòu)如圖3所示，網(wǎng)絡(luò)1的CNN-Pre Block內(nèi)部結(jié)構(gòu)見表1。

網(wǎng)絡(luò)結(jié)構(gòu)經(jīng)過前兩個(gè)網(wǎng)絡(luò)的回歸，對(duì)DRR1、DRR2圖像獨(dú)立進(jìn)行了剛體變換來消除已知參數(shù)的影響。由于前一步的回歸精度問題，必然存在參數(shù)誤差，而且空間剛體變換參數(shù)并非獨(dú)立作用于兩個(gè)平面上，故而在回歸Pout前對(duì)兩個(gè)平面上圖像進(jìn)行獨(dú)立的剛體變換不足以完全消除已知參數(shù)的影響，即輸入Pout回歸網(wǎng)絡(luò)的兩個(gè)平面上的殘差圖像I1、I2中還包含Pin回歸誤差以及平面剛體變換帶來的干擾。因此需要一個(gè)擬合能力更強(qiáng)的回歸網(wǎng)絡(luò)，在網(wǎng)絡(luò)1、2的基礎(chǔ)上進(jìn)行改進(jìn)出網(wǎng)絡(luò)3。

圖3 網(wǎng)絡(luò)主要框架Fig.3 Main architecture of network

通過對(duì)網(wǎng)絡(luò)1、2 在其CNN-Pre 部分進(jìn)行結(jié)構(gòu)修改來提高該網(wǎng)絡(luò)的特征提取能力，修改后的CNN-Pre如圖4 所示。其中，Conv2D-1 使用深度可分離卷積［14］；Conv2D-2 使用1×1 的卷積核減少特征數(shù)量；Conv2D-4使用5×5空洞卷積增大感受野［15］；Conv2D-5使用1×1 卷積核進(jìn)行特征融合；其它卷積核為3×3。所有卷積前都添加Batch Normalizetion［16］層，激活函數(shù)為ReLu［17］，同時(shí)Conv2D-1、Conv2D-2、Conv2D-3、Conv2D-6使用2×2的平均池化。

表1 網(wǎng)絡(luò)1的CNN-Pre Block內(nèi)部結(jié)構(gòu)Tab.1 Internal structure of CNN-Pre Block of network 1

圖4 網(wǎng)絡(luò)3中CNN-Pre結(jié)構(gòu)圖Fig.4 Structure of CNN-Pre of network 3

對(duì)于網(wǎng)絡(luò)3 的輸入圖像需要進(jìn)行預(yù)處理。首先去除兩個(gè)殘差圖像四周邊緣的32個(gè)像素寬度丟棄平面變換時(shí)在圖像邊緣引入的誤差，然后再按照相互重疊32 個(gè)像素寬度的方式將剩余圖像分割為4 個(gè)256×256 的區(qū)塊，最后將兩個(gè)平面上殘差圖像分割結(jié)果重疊為8@256×256，作為網(wǎng)絡(luò)的輸入圖像。

2 實(shí)驗(yàn)

2.1 數(shù)據(jù)獲取

在本實(shí)驗(yàn)中，采用DRR 圖像模擬X 射線影像以獲取準(zhǔn)確參數(shù)用于訓(xùn)練。針對(duì)CT 尺寸為512×512×283，物理尺寸為（400×400×283）mm3的頭顱CT 生成約26萬組訓(xùn)練數(shù)據(jù)［18］。每一組數(shù)據(jù)包含在初始參數(shù)下生成的2幅DRR圖像，在目標(biāo)參數(shù)下生成的2幅模擬X 射線影像，所有圖像大小都是512×512。設(shè)置兩次拍攝時(shí)，繞Z 軸旋轉(zhuǎn)的夾角為90°。取95%的圖像數(shù)據(jù)作為訓(xùn)練集，剩余部分作為測(cè)試集。

2.2 實(shí)驗(yàn)結(jié)果

在上述訓(xùn)練數(shù)據(jù)下，在顯卡為Tesla P100、CPU為Ⅰntel E5-2680 的服務(wù)器上使用深度學(xué)習(xí)框架Pytorch 0.4 編寫訓(xùn)練。訓(xùn)練輪數(shù)Epoch 為32，每一輪中不重復(fù)地隨機(jī)抽取的mini Batch大小為64，直到所有組數(shù)據(jù)都被使用過。訓(xùn)練中采取神經(jīng)網(wǎng)絡(luò)常用的Adam 優(yōu)化方法［19］。損失函數(shù)采用L2 范數(shù)即剛體變換參數(shù)的均方誤差。網(wǎng)絡(luò)1 訓(xùn)練結(jié)果如圖5 所示，其中學(xué)習(xí)率LR按式（4）進(jìn)行調(diào)整。

為對(duì)比該配準(zhǔn)方法的配準(zhǔn)效果，所采用的對(duì)比算法為CUDA 加速迭代式的雙X 射線影像配準(zhǔn)（CUDA）。在第2.1 節(jié)所述頭顱數(shù)據(jù)上1 萬多幅測(cè)試圖像中隨機(jī)挑選50組，配準(zhǔn)結(jié)果如表2所示。

其中采用的配準(zhǔn)結(jié)果衡量標(biāo)準(zhǔn)為mTRE［20］與配準(zhǔn)參數(shù)平均誤差ME。

mTRE計(jì)算方法如下：

取CT 中包含N個(gè)點(diǎn)的點(diǎn)集Q，計(jì)算其在配準(zhǔn)結(jié)果的變換Treg與目標(biāo)真實(shí)變換Tgd下的對(duì)應(yīng)點(diǎn)距離誤差均值，同時(shí)取CT 數(shù)據(jù)的物理尺寸長(zhǎng)度1% 即2.83 mm 作為配準(zhǔn)成功的標(biāo)準(zhǔn)。ME作為配準(zhǔn)精度衡量標(biāo)準(zhǔn)，ME(P)計(jì)算方式如下：

P與P′是包含N個(gè)平移配準(zhǔn)參數(shù)的向量，分別為配準(zhǔn)參數(shù)與真實(shí)參數(shù)，其中平移PT=[Tx,Ty,Tz]，旋轉(zhuǎn)PR=[Rx,Ry,Rz]，分別計(jì)算平移與旋轉(zhuǎn)的ME。

2.3 分析

實(shí)驗(yàn)結(jié)果表明，本文方法在ME 度量下，其平移參數(shù)回歸誤差為0.047 mm，相對(duì)迭代算法下降了81.20%；旋轉(zhuǎn)參數(shù)回歸誤差為0.057，相對(duì)CUDA 算法下降了68.51%，有顯著提升。在mTRE 標(biāo)準(zhǔn)下，本文方法配準(zhǔn)誤差為0.462 mm，僅為CUDA 算法的32.32%。因此本文方法的配準(zhǔn)精度在多個(gè)指標(biāo)下相對(duì)傳統(tǒng)算法都有較大提升。

圖5 網(wǎng)絡(luò)1訓(xùn)練結(jié)果Fig.5 Training result of network 1

表2 實(shí)驗(yàn)結(jié)果Tab.2 Experience results

X-DRR 圖像間的相似性度量是關(guān)于空間變換參數(shù)的高度非線性函數(shù)。因此在傳統(tǒng)迭代優(yōu)化過程中非常容易陷入局部最優(yōu)值，造成配準(zhǔn)誤差較大。而由于采用了CNN 回歸的方式，本文方法在配準(zhǔn)過程中沒有傳統(tǒng)迭代算法的優(yōu)化過程，依靠神經(jīng)網(wǎng)絡(luò)強(qiáng)大的擬合能力，從訓(xùn)練數(shù)據(jù)中尋找X-DRR 殘差圖像到空間參數(shù)的映射關(guān)系；又由于空間參數(shù)分步回歸，進(jìn)一步簡(jiǎn)化了這種映射，使得神經(jīng)網(wǎng)絡(luò)能回歸出更加精確的最優(yōu)解。

在50 組隨機(jī)數(shù)據(jù)中平均耗時(shí)僅為0.04 s，即達(dá)到25FPS 的處理速度，遠(yuǎn)快于CUDA 方法，達(dá)到實(shí)時(shí)配準(zhǔn)的標(biāo)準(zhǔn)。對(duì)于訓(xùn)練好的神經(jīng)網(wǎng)絡(luò)，在配準(zhǔn)中僅需要對(duì)于殘差圖像做一次單向的運(yùn)算，無需迭代生成DRR 圖像和計(jì)算圖像相似度，從而規(guī)避了傳統(tǒng)迭代優(yōu)化算法中最耗時(shí)的部分，具有非?？斓呐錅?zhǔn)能力。

3 結(jié)論

在ⅠGRT 系統(tǒng)的工作中，2D-3D 配準(zhǔn)速度的快慢直接影響患者的治療體驗(yàn)。本研究提出的利用雙X射線進(jìn)行空間坐標(biāo)分解，簡(jiǎn)化空間坐標(biāo)系變化與X-DRR殘差圖像映射關(guān)系；利用含有Dense Block 結(jié)構(gòu)的CNN 的強(qiáng)大非線性擬合能力，分步回歸出空間剛體變換參數(shù)。相對(duì)基于灰度的迭代優(yōu)化配準(zhǔn)算法，在配準(zhǔn)精度和速度都有明顯提升，尤其是其配準(zhǔn)速度已經(jīng)達(dá)到實(shí)時(shí)配準(zhǔn)的標(biāo)準(zhǔn)，在臨床治療中具有一定實(shí)用潛力。