基于正交视角X线图像重建的3D/2D配准方法

弥 佳，周宇佳，冯前进

南方医科大学生物医学工程学院//广东省医学图像处理重点实验室//广东省医学成像与诊断技术工程实验室，广东 广州510515

术前3D CT图像和术中2D X线图像的精确配准是椎间盘射频消融、腰椎椎间融合和经皮椎弓根置钉等脊柱手术导航的关键技术［1-3］。由于脊柱手术中有限的切口难以清晰显示感兴趣的解剖结构，其往往需要借助不同角度的X线透视，以获取病人内部的解剖结构影像，进而确认置钉或导针与目标及周围组织的相对位置。然而，X线透视图像只能提供二维信息，不足以清晰地显示病变部位的复杂结构及其空间关系［4-6］，导致定位的准确性严重依赖外科医生的主观经验。

3D/2D图像配准有望提升术中引导的质量［7］。通过将3D术前CT扫描与2D术中X线图像统一在同一坐标系中［1］，可以获得感兴趣区域的3D空间信息，并更好地显示骨结构和软组织。

3D/2D图像配准通常被建模为两个连续的基本任务。首先，考虑到X射线图像是投影图像，其中的每个像素值是由对应X射线路径上物质的吸收决定的，因此CT图像中的体素与X线图像的像素之间不存在一一对应关系，在配准前需要建立两种图像之间的维度对应。第二个任务就是估计6自由度的位姿参数，即x、y、z轴的平移和旋转偏移量。

目前大多数3D/2D配准都是使用前向投影策略建立维度对应，并基于优化方法迭代求解。在每次迭代中，使用当前最佳的位姿估计对CT图像投影得到数字重建X光片（DRRs）［8］，并根据DRR和目标X线图像的相似性测度得到位姿更新［1］。基于图像灰度的相似性测度考虑了所有像素的对应关系，从而可以获得亚像素配准精度。近年来，CNN的发展使得基于学习的医学图像配准方法成为可能。早期有研究提出使用CNN直接回归给定的DRR和X线图像对之间的位姿参数［9］。考虑到3D位姿的参数空间，另有研究提出，使用CNN学习从图像对到黎曼流形上的姿态更新的映射。Gao等［10］引入了投影空间变换（ProST）模块，该模块将空间变换网络（STN）［11］推广到投影几何中，并使用测地梯度实现了端到端可微的DRR渲染和位姿更新。然而，从3D到2D的投影过程不可避免地会导致空间信息丢失，在较大或较复杂的变换下，需要较为精确的初始化来保证后续参数估计和优化能得到有意义的更新。此外，迭代过程中需要多次DRR渲染，导致此类方法相对较慢。

重建（或反投影）策略可以避免对3D CT图像投影带来的信息损失。有研究通过整合多视图2D图像来重建3D图像，以实现3D/3D图像之间的配准［3］。但是在2D投影图像数量极少的情况下，基于数学模型的重建质量较低，需要特殊设计的相似性度量来优化后续位姿参数估计过程。随着深度学习的发展，有研究提出了使用少量X线投影图像重建3D CT图像的方法［12,13］。

基于特征的方法通过最小化3D/2D图像中相应特征之间的距离来确定位姿关系，如解剖标志点［14］、曲线［15］或曲面［16］。通过直接建立标志点之间的对应关系，可以避免投影或重建。理想情况下，若已知3D/2D标志点对应关系，则存在位姿的闭合解析解。Bier等［17］提出了由CNN回归器组成的多阶段的序列网络，自动检测X射线图像中的解剖标志点。后续有研究提出了一种自动推断患者特异性标志点的方法用于3D/2D配准［18,19］，该方法结合了预训练的二维图像标志点检测网络和基于多视角反投影的微调。另有研究通过三角测量层，将3D信息引入到训练过程中［20］。具体来说，该方法首先从成对的X射线和DRR图像中跟踪多视图2D标志点以定位3D标志点，然后通过三角测量层配准两组3D标志点，实现了端到端学习。基于特征的3D/2D图像配准方法可以避免DRR的迭代生成，提高计算效率；稀疏特征可以降低计算复杂度，但精度将受限至像素级精度。并且，此类方法精度受限于如何定义标记点对应关系，人工勾画解剖标志费时；由于X射线的透射特性，自动检测相应的标志具有挑战性。

因此，本文提出了一种基于重建的3D/2D配准方法，用于术中2D X线和术前3D CT图像的位姿参数估计（图1）。该方法对正交视角的X线图像进行重建，得到对应的3D表达，然后在流形空间中与3D术前CT图像进行配准，得到最终的6自由度的位姿参数。主要贡献如下：通过重建能避免对3D CT图像投影造成的信息损失，且在3D空间中完成配准能充分利用空间信息，实现更鲁棒的配准。将维度对应与位姿估计集成到统一的网络框架中，实现了端到端的训练。使用流形测地距离的损失函数，更好地适应位姿的参数空间，进一步提升模型的鲁棒性。

图1 本文方法概览Fig.1 Overview of the proposed method.

1 方法

1.1 实验数据

考虑到成对的X线和CT图像的收集和标注较为困难，本研究使用了大型公共CT数据集CTSpine1k［21］来进行仿真和训练。CTSpine1k包含了4个不同的开放获取的数据集，并对所有数据进行了椎骨的体素级标注。考虑到本研究集中在腰椎，因此，排除了HNSCC-3DCT-RT子数据集的头颈CT扫描和COVID-19子数据集的胸部CT扫描。其他低质量的图像也被排除在外。数据集已被划分为训练集、公共测试集和私有测试集。进行上述数据排除之后，训练集、公共测试集和私人测试集分别有546例、181例和169例。为简单起见，将CT扫描图像重采样至各向同性的1.25 mm×1.25 mm×1.25 mm分辨率。以图像中包含的所有腰椎的平均质心作为中心进行裁剪，得到128×128×128大小的3D图像。本研究首先利用DRR技术对原始CT图像进行投影，得到一对正交视角下的X线图像，分别是正位（AP）和侧位（Lat）视图。对CT原始图像进行随机刚性变换，以模拟术中与术前图像的位姿差异，并作为配准的金标准。为评估模型对于不同程度的位姿差异的配准表现，设置两个训练集，从不同的均匀分布中随机采样刚性变换的参数。训练集A）：旋转量∈(-10◦，10◦)，平移量∈(-20 mm,20 mm)；B）：旋转量∈(-45◦，45◦)，平移量∈(-50 mm,50 mm)。6个分量均独立采样。

1.2 模型

本研究的方法包含重建模块G和位姿估计模块R，前者完成2D图像到3D的重建，以实现待配准图像之间的维度对应，后者则完成待配准图像之间的相对位姿参数回归，得到最终的配准结果。考虑输入的X线图像xAP，xLat和术前CT图像ypre，为估计位姿变换Tgt，该方法可以描述为：

1.2.1 重建模块 重建模块（图2）基于X2CT-GAN 构建［13］，由两个分别用于不同视角的编码器-解码器网络和一个用于视角间融合的解码器组成。对于输入的两个正交视角X射线图像，由2D编码器分别编码，并转换为3D特征，在各自的3D解码器中解码。视角间融合解码器将两个正交视角的解码器网络融合在一起，重构出最终的3D图像。鉴于DenseNet在特征提取过程中的显著优势，这两个编码器都是使用DenseNet-121构建的［22］。DenseNet使用密集的连接重用特征，模型易于训练且参数使用效率高。每个解码器包含5个下采样块，每个下采样块首先使用一个步长为2的卷积层对特征进行下采样，然后通过一个密集连接块，再经过一个通道数减半的1×1卷积输出。解码器主要使用上采样块搭建，每个上采样块先使用一个卷积层压缩通道数，再经过2个卷积层之后使用步长为2的转置卷积扩大特征图的尺寸。每个卷积层后接一个归一化层和ReLU激活函数层。

图2 重建模块网络结构Fig.2 Network structure of the reconstruction module.

网络中使用了三种特殊的连接来建立起2D空间到3D空间的联系。在单视角编码器-解码器网络中，连接A调整编码器输出的特征F∈ℝC×H×W的维度，得到伪三维特征，然后使用三维卷积层处理之后连接到解码器作为输入。为了将编码过程中不同分辨率的信息直接传递到解码器，分别在两个单视角的编码器中通过连接B将不同尺度的二维特征复制并卷积成伪三维特征，传递给对应尺度的解码层。最后，为了融合两个正交视角的信息，通过连接C对来自两个正交视角解码器的特征进行平均和融合。

1.2.2 位姿估计模块 位姿估计模块（图3）通过孪生网络完成特征提取，其两个支路分别用于术前CT图像和术中X线重建的伪CT图像。两条支路对应结构参数共享，作为一种隐式的约束以提取共有的特征进行位姿回归，同时也能减少网络的参数量避免过拟合。每个支路除了处理输入图像的第一个卷积层使用7×7的卷积核，其他卷积块都使用3×3的卷积核。每个卷积块包含一个步长为2的卷积层用于下采样，以及两个步长为1的卷积层用于增加特征通道数。每个卷积层后接一个归一化层和ReLU激活函数层。

图3 位姿估计模块网络结构Fig.3 Network structure of the pose estimation module.

经过5次下采样之后，对应于两个输入的特征被展平为1D特征，并进行拼接，作为全连接层的输入。全连接层建模高维特征之间的非线性关系，回归两个输入之间的相对位姿。

两个模块联合训练，将重建模块的输出作为位姿估计模块的其中一个输入，损失函数的梯度会回传至两个模块，同时对两个模块的参数进行优化。两个任务之间相互约束和监督，以达到共同提升的效果。

1.3 损失函数

为了保证重建图像与输入图像的结构一致性，并且考虑到CT图像包含了大量的细节信息，使用重建损失约束重建结果与金标准之间体素级别的一致性。以术中CT图像yintra作为金标准，重建损失基于均方误差定义：

进一步地，考虑重建结果与金标准在3个投影方向下的一致性。使用正投影得到冠状面、失状面及横断面三个方向的投影图像，比较金标准和重建图像的对应投影。投影损失基于L1范数定义：

其中Pi(·)表示在第i个方向上进行投影。

3D位姿包含3D空间中的旋转和平移，其集合为特殊欧式群SE(3)。但大多数基于深度学习的方法将其参数化为单独的旋转和平移，并使用L2范数在欧式空间中对其分别进行回归。考虑到SE(3)空间的固有结构，我们参考［23］将位姿参数化为SE(3)的李代数，即se(3)的元素，并计算两个位姿之间的测地距离作为回归损失函数。

网络权重通过反向传播损失函数相对于Tpred的梯度来更新。测地损失使用geomstats实现［24］。

总的损失函数定义为：

1.4 实验细节

本研究的框架是基于PyTorch［25］构建的。使用Adam求解器［26］对网络进行训练，初始学习率为2e-4，动量参数β1=0.5和β2=0.99。学习率在训练过程中线性衰减至0。训练在1 台Quadro RT X 8000 GPU 上进行。批大小为4。预测值和金标准值被转换为欧拉角，并使用平均绝对误差（MAE）评估旋转角度及平移距离的估计误差。所有的指标在测试集上进行平均。

1.5 对比方法

1.5.1 POINT2 POINT2 是一个用于2D 感兴趣点（POIs）跟踪和3D POIs配准的端到端框架［20］。该方法根据由术前CT生成的DRR中已知的POIs追踪术中X线图像上的POIs，基于多视角的追踪结果得到3D POIs，与术前CT上的POIs进行配准以估计术前、术中的位姿差异。实验中，使用存在位姿差异的一对CT图像生成两个视角的DRR对，每个DRR对中两张DRR图像分别作为术前和术中图像。从CT图像中随机选择20个3D POIs，然后投影到DRR图像上进行跟踪。在实验中，POINT2是从头开始训练的，因为它比该方法建议的两阶段方式收敛得更快。分别使用了U-Net和Dense-UNet（使用密集连接块构建编码器）作为POINT2的骨干网络，记为UNet-POINT2 和DenseUNet-POINT2。

1.5.2 X2CT-Reg X2CT-Reg 也是基于重建的3D/2D配准方法，但对于重建和位姿估计两个模块分开训练。前者的训练按照X2CT-GAN的训练策略进行［13］，增加判别器以实现对抗训练。为调整对抗损失相对于重建损失和投影损失对模型优化的影响，在损失函数中分别对这3项乘以权重0.1、10和10。在训练位姿估计模块时，冻结重建模块的权重，仅更新位姿估计模块的参数。此时的损失函数只有回归项。

1.6 基于3D椎骨分割任务的配准模型

在上述框架的基础上，本研究增加对图像的3D椎骨分割任务探究分割信息对于配准的影响。重建模块G得到双通道输出，分别是重建的CT图像和对应的椎骨分割图；使用另一个3D U-Net网络处理CT图像得到双通道输出，分别是术前CT图像的重建结果和对应的椎骨分割图。计算重建结果和对应椎骨分割图的哈达玛(Hadamard)积，作为位姿估计模块R的输入结果。网络仍然作为一个整体进行端对端训练，损失函数中加入用于监督分割任务的Dice损失。

其中，λ作为一个参数用于平衡Dice损失对整体训练的影响。Dice损失的定义为：

其中，对于预测的分割图pi∈Spred和真值分割图qi∈Sgt的N个体素求和。分别将X射线图像分割和CT图像分割的Dice损失记为

2 结果

2.1 总体性能

在两个训练集对应的测试集A)和测试集B)上分别对4 种方法进行测试。UNet-POINT2 旋转和平移的平均误差分别为3.973±2.731 和2.233±1.993。DenseUNet-POINT2 将旋转和平移误差分别降低了42.72%和34.82%。而基于重建的方法可以取得更好的性能。X2CT-Reg 旋转和平移的平均误差分别为0.128±0.112和0.186±0.160。本文提出的方法的旋转和平移误差分别为0.115±0.095和0.144±0.124。相较于X2CT-Reg提升了10.30%和22.53%（表1）。当位姿差增大时，POINT2的性能急剧下降，UNet-POINT2的旋转和平移误差分别为12.160±11.116和3.859±3.249，而DenseUNet-POINT2 的旋转和平移误差分别为10.411±11.214 和4.584±3.827。与测试集（A）相比，POINT2类的方法对位姿参数（尤其是旋转参数）的估计严重受限，即使使用了更为复杂的骨干网络DenseUNet，对其性能的提升也有限，旋转误差仅减少了14.39%，而平移误差增加了18.78%。基于重建的方法可以得到更好的结果，X2CT-Reg旋转和平移的误差分别为0.987±0.834和0.931±0.758。而本文提出的方法的误差分别为0.792±0.659和0.867±0.701。采用联合训练的策略，比X2CT-Reg 提高了19.70%和6.88%（表2）。

表2 测试集B)位姿估计误差量化结果比较Tab.2 Comparison of pose estimation error on test set B)

2.2 分割信息对位姿估计的影响

同样地，在测试集A)和测试集B)上对其性能进行评估（表3～6）。其中表3和表5分别展示了测试集A)和测试集B)上的位姿估计误差，表4和表6分别展示了对应的分割量化结果。

表3 引入分割信息时测试集A)位姿估计误差量化结果Tab.3 Quantification results of pose estimation error on test setA)when segmentation information is introduced

表4 测试集A)分割量化结果Tab.4 Quantification results of segmentation performance on test setA)

表5 引入分割信息时测试集B)位姿估计误差量化结果Tab.5 Quantification results of pose estimation error on test set B)when segmentation information introduced

表6 测试集B)分割量化结果Tab.6 Quantification results of segmentation performance on test set B)

如表3所示，增加分割任务(λ=0.1)时，术前图像和术中图像分割的DSC值分别为0.896±0.037和0.965±0.035。位姿估计的性能有所提升，旋转估计误差为0.103±0.080，平移估计误差为0.143±0.116，分别降低了10.43%和0.69%。增大λ为1.0，分割性能的提升并不显著，术前图像和术中图像的DSC值分别增大了0.78%和0.10%，但配准性能下降。对于测试集B)，当λ=0.1时分割结果较差，术前图像和术中图像分割的DSC值分别为0.899±0.027和0.887±0.050，此时的旋转估计误差和平移估计误差分别增大了18.18%和3.11%。增大λ为1.0，在术前图像和术中图像的DSC值分别增大了0.56%和7.78%的情况下，配准的旋转误差增大了6.82%，平移误差降低了9.11%。

3 讨论

术前3D CT图像和2D术中X线图像的精确配准是各种图像引导脊柱手术的关键技术之一，用于统一术前图像与术中图像的坐标系，为术中有限的成像条件下的2D影像提供补充的3D信息，以提升术中引导的质量，改善临床结果。

大多数3D/2D配准都是使用前向投影策略建立维度对应，并评估DRR和X线图像之间的相似性［9,10,27］。对CT图像投影造成了CT图像3D空间信息的损失，并且2D投影空间中的空间信息模糊，3D位姿估计困难。本文使用正交视角的X线图像重建3D图像，最终在3D空间中完成图像配准。以数据驱动的方式学习从2D投影图像到3D立体图像的空间信息补全，降低空间信息模糊对于位姿估计的影响，使得可以在单次前向推理中得到精度满足要求的配准结果。

基于特征的3D/2D图像配准方法［17-20］使用稀疏特征降低计算复杂度，但精度也将受限于特征检测。在2.1节所示的实验结果中，位姿差异增大时图像外观的变化显著，会使得POINT2对POIs的检测性能很大程度上降低，从而导致位姿估计误差增大。而基于重建的方法由于充分利用了3D空间信息，具有较强的鲁棒性。同时，相对于POIs，基于灰度的配准策略为姿态估计提供了更多特征，以产生亚像素级的精度。

此外，本方法将维度对应与位姿估计集成到统一的网络框架中，实现了端到端的训练。不同于先训练生成网络再进行配准的策略［28］，联合训练统一了两个模块的优化目标，使得重建过程也受到位姿估计的监督，以位姿估计为导向。同时网络也学习处理生成图像和真实图像外观之间存在的不一致，而不用额外人工设计相似性测度［3］，实现更高效的学习。

然而本文的算法中需要对3D CT图像进行逐像素的重建，其中包含了丰富的纹理、结构等细节信息，导致重建3D CT图像困难，而这些信息对于后续刚性配准来说可能是冗余的。分割信息的引入可以强调学习过程中对于刚性结构的关注，有可能提高后续位姿参数估计性能。通过2.2节的实验结果我们发现，分割信息的引入对位姿估计可以有一定的提升，但提升依赖于分割的性能，以及分割和配准任务在整体优化目标中的平衡。如表3和表4所示，对于测试集A)，分割任务性能较好时，可以对配准的性能有所提升。由于测试集A)数据的位姿采样自更小的范围，其分割难度较低，当对应于分割任务的损失函数权重增大，对分割性能的提升并不显著，反而会导致位姿估计相关的损失对于整体训练的影响减小，使得配准性能下降。但对于表5和表6所示的测试集B)，通过增大分割损失的权重可以有效地提升分割性能，继而优化位姿估计。因此，如何更有效地利用分割信息以及如何平衡多任务的学习，仍然需要进一步的研究。

本研究中使用的数据来自公共数据集，且采用标准的正侧位视角，其鲁棒性仍然需要通过真实的配对X射线和CT扫描图像进一步评估。

综上所述，本文提出了一种精确、鲁棒且端到端的3D/2D配准方法，速度满足实时性要求，有望成为提高脊柱手术导航精度的工具之一。该方法基于图像3D信息的重建，建立了3D/2D图像之间的维度对应，并联合位姿估计构成了一个统一的框架。通过联合训练，利用多任务之间的互相监督达到互相提升的目的。使用了流形测地距离的损失函数，更好地适应了位姿的参数空间，提升了模型的鲁棒性。