基于姿态编码器的2D/3D脊椎医学图像实时配准方法

徐少康，张战成，2*，姚浩男，邹智伟，张宝成

（1.苏州科技大学 电子与信息工程学院，江苏 苏州 215009；2.苏州科技大学 苏州市虚拟现实智能交互及应用技术重点实验室，江苏 苏州 215009；3.中部战区总医院 骨科，武汉 430070）

0 引言

2D/3D 医学图像配准是三维手术导航系统中的关键技术之一，在脊椎手术中，椎弓根螺钉植入是临床脊柱外科手术中最常用的内固定方法，术中精确植入螺钉的同时，不损伤血管和神经是一项挑战性的工作［1］。医生术中通过X 射线图像（2D）对病人脊椎位置进行分析，查看螺钉的植入情况。但2D 图像的清晰度能够反映出的信息较少，想要精确地分析出植入螺钉位置信息较为困难。而3D 图像可以提供此类的多维信息，如果在手术过程中多次实施CT 成像会对医生以及病人造成伤害，将术前3D 浮动图像对照术中2D 参考图像进行空间转换达到坐标系相统一，从3D 浮动图像中引导出对应位置的清晰的2D 图像，即可更精确地引导手术，提高植入椎弓根螺钉的精确度与外科医生的手术效率，同时减少医生和病人的辐射伤害。

传统的基于迭代回归优化的2D/3D 配准方法［2-4］：从一个假设的姿态开始，通过对CT（Computed Tomography）进行模拟虚拟X 射线的衰减创建数字重建影像（Digitally Reconstructed Radiographs，DRR），评估DRR 图像与术中X 射线图像的相似性，并优化搜索下一个潜在的姿态；然后继续为这个新姿态生成DRR，直到相似性达到最大时，完成配准。然而，在该传统方法中涉及多次的DRR 成像以及大量的图像相似度计算，计算成本都很高。因此该类方法无法满足手术过程中对实时配准的需求。此外，在配准优化过程中，优化器容易陷入局部最优值，导致配准的成功率无法满足临床应用的精度要求。

近年来，基于机器学习的配准方法已经逐渐运用于医学图像配准任务当中［5-9］。Miao 等［10］提出基于卷积神经网络回归方法，分层学习姿势估计（Pose Estimation Hierarchical Learning，PEHL），以实现具有较大捕获范围的实时2D/3D 配准。为了能够提高2D/3D配准速度，Pan等［11］提出了一种基于阈值分割的方法来加速配准过程，通过使用基于Bresenham方法的加速DRR 算法，以减少生成数字重建射线照片的时间，使得光线投射的时间下降到0.937 s；但该方法仅是提高了DRR 的计算速度，并没有真正避免迭代计算的过程，在整个配准的时间中，仍然耗费了91.301 s。Wohlhart 等［12］提出训练卷积神经网络（Convolutional Neural Network，CNN），从距离图像中学习姿势区分描述符，并使用最近邻进行姿势估计。虽然使用该方法可以实现全局姿态估计，但精度相对较低，即角度误差小于5°的成功率低于60%。Liao等［13］通过用于跟踪和三角测量的兴趣点网络（Point-Of-Interest Network for Tracking and Triangulation，POINT2）来解决多视图2D/3D 刚性配准问题。在2D 和3D 中建立干预前和干预内数据之间的点对点对应关系，然后在匹配点之间进行形状对齐以估计CT的位姿，从而避免了通常昂贵且不可靠的迭代姿势搜索；但该方法需要来自至少两个视图的X 射线可用。Gouveia 等［14］从X 射线图像中提取了一个手工制作的特征，并训练了一个多层感知器（Multi Layer Perceptron，MLP）回归器来估计三维变换参数；然而该方法所能达到的准确度远低于使用基于优化的方法所能达到的准确度。

Sundermeyer 等［15］提出了一种基于RGB 的实时目标检测和6D姿态估计方法，提供了一个由隐空间中的样本定义的对象位姿的隐空间表示。数据的隐空间表示包含表示原始数据点所需的重要信息，隐空间表示用于将更复杂的原始数据形式（如图像、视频）转换为更“易于处理”和分析的简单表示［16-18］。Gu 等［19］提出将隐空间分解为两个次隐空间：模态特定的隐空间和姿势特定的隐空间，用于三维手部姿态估计，在基于单图像的手部姿态估计方面的性能明显优于目前最先进的方法。Jiang 等［20］提出了一种通用的自适应对抗隐空间框架，该框架主要由自适应隐空间发生器（Adaptive Latent Space Generator，ALSG）和基于约束的对抗自动编码器（Constrained Antiversarial AutoEncoder，CAAE）两部分组成。通过自适应映射器和对抗学习方法，建立ALSG，获得真实的潜在空间分布。这些研究表明隐空间可以较好地重构姿态信息。

为了提高2D/3D医学图像配准的速度，达到术中实时配准的要求，受隐空间回归启发，将2D/3D 医学图像配准问题表达为基于自编码的回归网络，并加入隐空间回归约束，以捕获X射线图像中姿态信息，从而精确地回归脊椎的旋转的姿态参数。对比了传统的优化迭代配准方法，本文方法避免了大量的DRR渲染时间，使得2D/3D图像配准的时间满足实时的要求。

1 相关工作

1.1 DRR投影模型

DRR 图像生成算法有许多种：抛雪球法（Splatting）、光线投射法（Ray-Casting）、光流场法（Light-Field）等。在2D/3D 医学图像配准中主要使用光线投射算法。Gao 等［21］提出的投影空间变换器模块（Projective Spatial Transformer，ProST）可将空间变换器推广到投影几何，从而实现可微分的体积渲染。

从三维立体场景到投影透射图像Im的映射可以建模为Im=T(θ)V，其中T(θ)是依赖于姿态参数θ的变换矩阵，V是3D CT体积图像。在基于强度的2D/3D配准中，寻求检索姿势参数θ，以便从3D CT扫描模拟的图像尽可能趋近于采集的图像。投影空间变换通过学习控制点G的变换来扩展标准投影几何，每条投影光线在CT 内有K个均匀间隔的控制点。给定θ∈SE(3)，通过仿射变换矩阵T(θ)获得一组变换后的控制点，投影平面长宽为M×N。T(θ)通常由三个平移参数tx、ty、tz和三个旋转参数rx、ry、rz来参数化，可以写成齐次坐标下的4×4矩阵：

其中R(θ)是控制CT 体积绕原点旋转的旋转矩阵。由于这些控制点位于体积内、位于体素之间，因此对控制点的值进行插值。

其中Gs∈RM×N×K，最后得到一个二维图像I∈RM×N沿每条射线积分，通过“堆叠”(m，n)处的Gs的k维度来实现。

上述过程利用了空间变换网格，将投影操作简化为一系列线性变换，使得2D/3D 配准的DRR 图像生成过程变得快速精确。

1.2 三维转换参数

一个刚体三维变换可以用一个含有6 个参数的矢量来参数化［10］。在本文方法中，通过3个平面内和3个平面外的转换参数来参数化这个三维变换。其中，平面内的变换参数包含2个平移参数tx、ty和1个旋转参数rz。平面内的变换参数的结果看作为二维刚体变换。平面外的变换参数包括1 个面外平移参数tz和2 个面外旋转参数rx和ry。平面外平移和旋转的结果分别对应着刚体图像的缩放和形状变化。

在脊椎手术中，为了能够提示外科医生在钻探椎弓根螺钉导孔时可能引发的椎体皮质骨折，通常会在植入椎弓根螺钉前向脊椎内植入一枚探针用于探查脊椎管道情况，再选择合适的椎弓根螺钉进行植入。在进行脊椎图像配准前，外科医生可以直接通过探针植入的位置准确地知道需要配准的是第几节脊椎位置，所以对于X射线图像的固定平移拍摄，平移的三个参数是已知的，只需要计算出旋转姿态参数rx、ry、rz即可完成配准。

2 脊椎姿态回归方法

2D/3D 配准的目标是搜索CT 图像的6 个未知自由度姿态，以该姿态对CT 体积模拟投射X 射线在平面上产生DRR与X 射线图像进行比较。因此2D/3D 配准的问题最终可以归结到如何从X射线图像中寻找最佳CT投影姿态的问题。

2.1 配准模型结构

训练模型结构如图1 所示，前半部分由Encoder 加全连接层组成，本文Encoder 采用卷积自编码器架构。该结构由4 个卷积组成，每层使用3×3的卷积核，步长为2，在编码过程中应用2D 卷积。在实验中，输入图像大小为128×128×1，但框架不受特定大小的限制。每个卷积后面都有一个参数为0.2 的LeakyReLU 层。模型的输入为脊椎X 射线图像，经过Encoder将X射线图像特征映射到隐空间并从中提取隐变量。然后通过接入128 → 64 → 32 → 16 → 3 全连接层，回归出对应的姿态参数向量。通过投影得到DRR 图像，即配准图像。在训练过程中，将配准图像输入前者的Encoder 模块，从配准图像中继续提取隐变量，用于辅助损失的计算（具体计算方式见2.3节）。

图1 配准网络框架Fig.1 Framework of registration network

式（5）中：ICT为CT图像，IX为X射线图像，φ为Encoder+MLP姿态回归过程模型，Proj 为DRR 投影，Ireg为配准图像。编码器用F表示，该网络将原始数据X映射到隐空间F中(φ：X→F)。通过用激活函数σ传递的标准神经网络函数表示，其中Z是隐变量。

2.2 隐空间生成

在隐空间特征表示中，相似的样本图像之间的差别特征会被作为非核心信息被剔除，只有其中核心特征信息会被保留并且学习，所以在数据特征点被映射到隐空间后，相似特征点距离会更近。在本文的模型训练过程中，输入为不同姿态下的脊椎X 射线图像。人体的脊块部位特征基本类似，差异性较小，所以对于此类的训练图像，姿态信息即为核心特征信息，本文通过Encoder生成的隐空间及隐变量会最大限度地学习脊椎的姿态信息，并精确地回归出姿态参数。

在自编码器结构中，通常使用Decoder层来进行图像的重建，从而训练隐空间的生成过程。然而如果使用Decoder重建配准图像，与CT 图像中的脊椎信息没有直接关系，缺少医学的可解释性，会对医生的判断造成干扰。本文中使用DRR 投影代替Decoder工作，DRR投影的计算原理与X射线图像的生成原理类似，都是通过光射线穿过体素后的衰减情况来确定2D 图像各像素值。相比之下，以CT 图像为基础投影生成的DRR 图像更加满足本文方向上图像重建的质量要求，外观与X射线图像更加相似。

2.3 损失函数

对于神经网络而言，性能除了取决于本身的结构之外，损失函数的设计也极其重要。为了更好地学习网络参数，本文设计了一个基于“粗细”结合的损失函数。其中“粗配”包括预测的姿态参数与真实姿态r(i)的L1 范数损失LossL1。在训练的过程中，通过姿态编码器回归出的姿态参数以长度为3的向量的形式存在，并在回归完成前对其进行归一化，使得姿态参数向量的3个旋转姿态值在（-1，1），对应(-57.3°，57.3°)角度范围，因此L1范数损失适合此类信息的计算。此外还包括X 射线图像中提取的隐变量Zi与DRR 图像中提取的隐变量使用余弦相似度计算的Losscosine。损失函数定义如下：

其中γ是一个正则化参数（在3.5节中会对其进行实验分析取值）。实验表明当初始姿态与真实姿态差异较大时，式（9）损失函数在收敛过程中更有效，但当DRR 图像与X 射线图像逐渐一致时，该损失函数对细微的差异敏感性变低。

在此基础上，本文继续引入一种用于“细配”的损失函数：基于梯度的归一化互相关（Gradient-based Normalized Cross Correlation，GradNCC）［22］，基于梯度的度量最初通过微分进行转换。水平和垂直的Sobel算子用于创建梯度图像，并在图像的两个正交轴上表示透视强度的导数。然后计算DRR 的梯度图像和X 光梯度图像之间的归一化互相关，该测量的最终值是这些归一化相互关的平均值。梯度测量的优点是滤除图像之间的低空间频率差异，如软组织引起的差异。在刚体的图像中，边缘含有丰富的特征信息，所以该损失集中于边缘信息的相似性度量。将该损失函数可以表示为：

其中：IXI与IDRR分别对应X 射线图像与DRR 图像在区域(i，j)的强度值，与为重叠区域(i，j) ∈T内图像的均值。

本文使用随机梯度下降（Stochastic Gradient Descent，SGD）优化器进行优化，计算了Losscru最近10 次迭代的标准偏差（STandard Deviation，STD），设置了STD<3×10-3的停止准则，然后使用SGD 优化器切换到GradNCC 相似度，直至数据集训练结束。

3 实验与结果分析

3.1 实验数据

实验数据采用中国科学院计算技术研究所与北京积水潭医院共同发布的CTSpine1K 脊椎数据集［23］。包括1 005 个CT 样本（超过500 000 个CT 切片和11 000 个脊块）。该数据集将所有DICOM 图像重新格式化为NIfTI，以简化数据处理。

本实验中用于训练的X 射线图像通过对相应的CT 数据进行DRR 投影合成得到，这样方便记录对应真实的姿态参数。为了能够充分考虑实际手术中医生之所以需要配准的原因，通过在合成的X 射线图像上进行高斯模糊，以使合成X 射线图像的外观与术中X 图像一致，从而使其训练的回归器可以很好地推广到真实X射线图像上，整体流程如图2所示。

图2 合成X射线图像过程Fig.2 Synthesis process of X-ray image

在脊椎手术过程中，医生手术的部位主要是病人具体的某一两块脊椎。因此考虑手术的实际情况，本文使用每一张脊椎CT 数据的12 节胸椎部位实验，通过使用ITK-SNAP 将12 节胸椎从第一节开始，每三节分割成一部分，如图3 所示。这样共将人体胸椎分割为四部分，可以根据手术要求确认需配准的对应部分。实验中，本文采用100 位病人CT 图像作为训练集，对人体胸椎每三节分割并进行实验，共得400 张部分脊椎CT。对其在（-20°，20°）旋转姿态范围内，遵循均匀分布合成204 800 张X 射线图像作为训练图像，图像大小为128×128。训练阶段，使用对应脊椎部位的合成X 射线图像进行训练，每个病人的脊椎CT 对应约2 000 张合成X 射线图像，进行10 折交叉验证测试。

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

3.2 实验环境与配置

实验在一台配备Intel Core i7-4790k CPU、16 GB RAM 和NVIDIA Tesla P100 GPU 的工作站上进行。神经网络是使用开源深度学习框架PyTorch 实现。训练阶段中，使用SGD optimizer 对网络进行训练，循环学习率在1E-6 和1E-4 之间，每100 步，动量为0.99，Batch size 大小选择为8，训练了2 000 次迭代，直到收敛。

3.3 评价指标

为了评估最终的配准精度，本文采用平均目标配准误差（mean Target Registration Error，mTRE）以及平均绝对误差（Mean Absolute Error，MAE）。将mTRE 小于2 mm 且MAE 小于0.1 的结果视为配准成功，配准成功的百分比定义为成功率。为了准确地评估本文方法的配准实时性，同样在测试的过程中记录每次配准所需时间。

3.4 实验结果分析

本文在每位病人的CT 测试集上进行200 次配准测试，共完成20 000次配准测试。随机展示了两组在不同部分脊椎上的配准结果，如图4 所示。第一列与第三列为加上高斯模糊合成X 射线图像，第二列与第四列为对应前列的通过模型预测姿态投影生成的配准图像，其中每一行为不同分段的部分三节脊椎。由表1 可见，同样展示了四部分脊椎在合成X 射线图像真实旋转的姿态参数与预测旋转参数上的误差。

图4 脊椎配准结果可视化Fig.4 Visualization of spine registration results

表1 不同部位脊椎的平均旋转角度误差 单位：（°）Tab.1 Average rotation angle error of spine at different parts unit：（°）

这里选取了三种基于优化的方法（Opt-GC［3］、Opt-GO［3］、Opt-NGI［4］）以及三种基于学习的方法（Bresenham［11］、Fourier［8］和MLP［14］）作为对比方法与本文模型进行了比较。

Opt-GC［3］、Opt-GO［3］、Opt-NGI［4］：三种基于优化的方法主要通过迭代评估每次DRR图像与X射线图像的相似性优化投影参数，从而在多次优化后获取精度最高的配准图像。三种方法采用针对不同特征使用不同的图像相似度计算方法，充分考虑医学图像特点，其最终配准精度满足医生对配准的要求。但迭代配准过程所需时间难以满足医生对术中实时的要求。

Bresenham［11］：该方法以加快DRR 投影为主要研究方向，提出基于 Bresenham 直线生成算法改进的模式强度与梯度相结合的混合配准算法，提高投影速度。最终采用改进的鲍威尔优化算法对参数进行迭代优化，完成配准过程。

Fourier［8］：该方法基于主方向傅里叶变换算子的模板匹配初始化方法，可避免接近真值的初值需求，并显著扩大了捕获范围。在确保配准精度的前提下，大幅缩短配准时间。

MLP［14］：该方法以多层感知机为主要网络结构，通过学习提取的手工特征预测脊椎姿态，避免了多次投影迭代过程。

本文在同样的测试数据集以及同样设备环境中进行200次实验，并在配准精度、配准成功率、配准耗时上进行量化比较，进一步在mTRE 和MAE 两个指标上，将本文方法和其他对比方法进行了配对t 检验，结果如表2 所示。在mTRE 和MAE 两个指标上，本文方法显著优于MLP 方法，与其他五种对比方法没有统计意义上显著差异（P>0.05），可以满足外科医生在实际手术中的配准精度指标。在配准成功率上，本文方法明显高于其他对比方法，体现出较好的鲁棒性。在配准耗时上，本文方法优于对比方法，基本可以满足2 s 内显示图像的临床需求。因此本文方法能在保证较高配准精度的同时有效提高配准速度，满足临床快速配准的要求。而三种基于学习的方法（Bresenham、Fourier 和MLP）虽然在配准速度上较传统方法有所加快，但Bresenham 与Fourier 方法仍然无法达到术中实时要求；其中MLP 回归特征的方法虽然时间性能满足要求，但却以牺牲配准精度为代价。

表2 本文方法与其他方法的实验结果对比Tab.2 Comparison of experimental results of the proposed method and other methods

本文模型通过姿态编码器直接预测脊椎姿态，相较于基于优化的方法以及Bresenham 方法，解决了耗费大量时间的迭代优化问题，整个配准仅进行一次DRR 计算。Fourier 方法虽相较于传统方法配准时间有大幅减少，但同样分级配准过程仍计算量较大，对于实时配准的要求仍无法满足。针对MLP 直接通过提取图像特征预测姿态信息，姿态无法作为核心信息被提取的问题，本文姿态编码器通过提取X 射线图像中含有核心姿态信息的隐空间来预测姿态信息，从而避免图像中其他非核心信息的干扰，保证了配准精度。

3.5 正则化参数分析

图5 展示了在不同正则化参数γ值下的模型训练的mTRE 测试实验结果。本文首先在0～0.5 每隔0.01 对γ进行赋值训练，根据图5 折线图中实线可知，γ在0.32～0.33 取值，模型预测结果的mTRE 值较低。因此，为得到更精确的参数值，本文继续在0.32～0.33 以0.001 的间隔继续对γ进行赋值训练，图5 中虚线结果表明：当γ取0.323 时，对Losscru训练效果更好。

图5 不同γ取值下的mTRE结果Fig.5 mTRE results under different γ values

3.6 消融实验

为了展示引入的三种损失函数组合的有效性，在同样的整个测试集上对LossL1、Losscosine、Losscru、LossNCC、LossL1+LossNCC、Losscosine+LossNCC以及本文采用的LossNCC+Losscru分别进行实验，在配准的角度误差以及配准精度方面记录实验情况。

如表3 与图6 所示，实验结果表明，在“粗细”结合的损失函数的约束下，模型训练效果较好。模型测试配准结果的平均目标配准误差仅有1.1 mm，平均绝对误差仅有0.05。在旋转的3 个角度的误差上也控制在了1.1°左右。同样在此实验中进行t 检验，在最终的配准精度指标方面，LossNCC+Losscru明显优于其他组合损失（P<0.05）。

表3 不同损失训练的相似度结果Tab.3 Similarity results of different loss training

图6 不同损失的旋转角度误差Fig.6 Rotation angle errors of different losses

4 结语

2D/3D 医学图像配准是脊椎骨科手术导航中的关键技术。针对传统的基于迭代优化的方法无法实时配准的问题，提出一个改进的姿态编码器的配准模型，设计了针对姿态回归的编码器［24］。通过该编码器对不同姿态的X 射线图像生成隐空间表示，提取图像姿态特征信息，继而回归脊椎部位的旋转姿态参数，仅通过一次投影即可获得配准图像，避免了大量迭代姿态搜索过程。在真实的CTSpine1K 数据集上进行了实验，结果表明，本文模型在配准精度、配准时间、配准成功率指标上取得了较好的结果，配准时间优于传统方法，在Intel Core i7-4790k CPU、16 GB RAM、NVIDIA Tesla P100 GPU 硬件条件下可以满足术中实时配准的要求。本文目前的实验结果仅限于脊柱部位，其他非刚体配准［25］的效果有待进一步验证，以及如何在更低成本的嵌入式硬件平台上实现2D/3D 的实时配准是我们下一步的研究方向。