基于深度学习的2D/3D 图像配准方法在脊柱微创导航手术中的应用

邢珍珍，颜立祥，张麟华

1.太原工业学院，山西 030008；2.电子科技大学

脊柱类疾病是比较常见的一种疾病。近年来，脊柱类疾病的发病率不断攀升。据《2005—2017 年中国疾病负担研究报告》[1]显示，颈椎、腰椎类疾病是导致国人健康寿命损失的首要疾病之一。手术是治疗脊柱类疾病的常见方式，随着医学影像技术的高速发展，影像导航手术逐渐成为脊椎疾病手术的研究热点。医生可以结合病人不同成像模式的影像获取病变区域的信息[2]。术中可以通过X-ray 等2D 图像实时跟踪手术器械和人体组织的相对位置。3D 图像一般包括计算机断层扫描(CT)、锥形束CT(CBCT)、磁共振成像(MRI)等[3]，可以提供更加精确的信息，但其扫描操作难，只能在术前获取。因此，将术前3D 影像数据和术中2D影像数据置于同一坐标系中，进行图像配准，可以在术中提供图像引导，提高手术的成功率[4-5]。传统的2D/3D 配准方法主要包括基于灰度信息的配准和基于特征信息的配准[6-7]。其中基于灰度信息的配准算法，通过优化图像的灰度相似度来确定图像间的变换参数，对于存在形变或灰度变化明显的图像，配准效果较差[8]。基于特征信息的配准算法，通过提取图像中点、线等特征完成配准，过程简单、速度较快[9]，但是容易丢失图像的灰度、梯度等有价值信息。深度神经网络近年来逐渐应用于医学图像配准领域[10]。其中，卷积神经网络(convolutional neural networks,CNN)广泛用于特征点检测[11-12]、特征描述符[13]及图像配准[14-15]中。本研究提出一种基于卷积神经网络的2D/3D 脊柱CT的层级配准方法，为临床脊柱微创手术中的3D 导航问题提供技术支持。

1 数据与方法

1.1 实验数据 采用开源数据集LIDC-IDRI(The Lung Image Database Consortium)，美国国家癌症研究所对1 010 例病人进行研究，制作了1 018 个研究样本。每个样本包含胸部医学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 图像，与术中图像实现维度的统一[16]。该过程通过DRR技术实现[17]。DRR 是医疗图像配准的一个重要的前置步骤。在进行2D/3D 配准前，需要对术前CT 图像进行空间变换，获取不同角度下的DRR 图像，通过计算DRR 图像和术中X-ray 图像的相似性测度值来构建目标函数，并对目标函数进行优化，获取最优空间变换参数，使得待配准图像的相似性测度值达到极值。配准公式如下。

其中：Tg表示配准后的空间变换，IR是参考图像，IM是浮动图像，T 表示空间变换，P表示DRR 投影，函数S是IM与IR间的相似性测度。IM与IR完全配准时，IR与IM生成的DRR 图像在像素值和空间位置上均相同，此时满足公式(2)。

1.2.2 配准流程 本研究的配准方法包含两个部分。首先，利用深度学习网络进行粗配准，其次利用参数优化方法进行精配准。配准流程：首先，将3D CT 数据集产生形变生成DRR，通过浮动图像和参考图像来训练深度学习网络。利用回归方法学习投影参数之间的关系，将待配准的术前CT 图像生成的DRR 和术中2D的X-ray 图像输入网络，结合学到的关系和几何变换得到相关参数，完成粗配准。其次，将待配准的CT 图像进行手动分割，完成精配准。计算分割后浮动图像和参考图像之间的相似性测度，通过Adam 参数优化算法更新参数，完成单块配准。循环单块精配准过程，完成多块椎骨的配准。流程图见图2。

图2 层级配准算法流程图

1.2.3 粗配准 本研究对象为病人脊柱，属于一种刚体器官，故通过卷积神经网络对医学图像进行刚体配准[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 的卷积层和池化层提取参考图像和浮动图像的差值图特征，通过全连接层来依次回归形变参数，做回归预测。可以通过目标区域的中心、宽度、高度和方向4 个参数(q,w,h,φ)对其进行定位，计算公式如下。

其 中：D是X 光 源 到DRR 平 面 的 距 离，w0是 目 标区域的大小。图像配准的实质是预测变化参数，参数之间的变化用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 个形变参数，容易使得训练陷入局部最优[22]，且回归精度较差。本研究采用按顺序回归的方法，首先回归(tx,ty)，将图3 网络中最后全连接层的激活函数单元设置为2，进行回归预测，更新差值图的(tx,ty)，作为下一步网络的输入。用同样的方法依次回归(rz)，(tz)，(rx,ry)。最后进行组合，得到形变参数(tx,ty,tz,rx,ry,rz)。

采用均方误差函数为损失函数，相对平移误差，旋转误差更为明显，所以，在损失函数中添加对旋转参数的权重，损失函数公式如下。

每张浮动图像与对应的参考图像组合，构成一个训练样本。将训练样本依次输入网络模型，得到浮动图像与参考图像的差值图像，如图4 所示。通过学习差值图的特征信息，输出待配准图像的形变参数。

图4 参考图像、浮动图像和差值图像

1.2.4 精配准 精配准主要是对粗配准得到的形变参数进行优化，对椎骨影像进行分块配准。采用手动分割的方式，通过矩形框标记，分割的每幅子图仅包含1 块椎骨。如图5 所示方框区域所示。

图5 图像分割示意图

将单块椎骨图生成DRR 图像，作为浮动图像，计算浮动图像与参考图像之间的Dice Loss 值，该值若小于预设阈值，则完成单块椎骨图的精配准，若不满足条件，将Dice Loss 设置为Adam 参数优化算法的目标函数，重复上述步骤，搜寻Dice Loss 值最小时的最优精配准参数，完成单块椎骨图的配准。按同样步骤完成所有单块椎骨图的精配准，从而实现脊柱CT 层级配准。

1.2.5 配准实验

1.2.5.1 实验一：单个对象的多形变配准 首先以数据集LIDC-IDRI 的第2 张CT 为实验对象，模拟生成10种形变。采用Elastix的配准算法作为对照组，Elastix广泛应用于医学图像的配准[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）方法进行学习和优化，设置权重衰减系数为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，完成图像的精配准。将分割区域图像的像素矩阵值与原图融合，得到最终的配准结果图。

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的协方差。

联合熵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 讨论

本研究提出了一种基于深度学习和参数优化的2D/3D 层级配准方法，通过对某一病人的CT 图像模拟生成了10 种形变，将基于Elastix 的配准方法作为基准实验组，10 组数据中，MAE 指标下，DLIR 有10 组优于ExIR；NCC 指标下，DLIR 有7 组优于ExIR；NMI 指标下，DLIR 有10 组优于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 均优于ExIR。另外，使用LIDI-IDRI 数据集中的20 张CT 的配准实验结果显示，表明提出的方法在MAE 指标上降低了0.097 8，在NCC 指标上提高了0.185 6，在NMI 指标上提高了0.445 6，表明本研究提出的方法3 项指标均优于ExIR 方法。本研究所提方法既考虑到了脊柱刚体变换的全局信息，又考虑到了两块椎骨之间的局部信息，有效提高了配准精度。说明将深度学习应用于脊柱微创导航手术的CT 图像配准中具有可行性，是一项非常具有前景的研究。