传统的配准方法[2-5]通过定义相似性度量函数和正则化约束，将配准问题转换为迭代优化问题，通过不断地优化使配准后的图像与固定图像之间的相似性最大，从而得到最优配准图像。然而，基于优化的传统方法通常需要在高维空间上进行迭代的梯度运算，导致运算成本昂贵。医学影像配准在应用中由于缺乏标签数据在训练与评估中面临较大问题，与此同时由于不同模态图像强度分布未知且复杂导致其极具挑战。
最近，基于深度学习来进行可变形图像配准的方法得到广泛发展，无监督学习[6]、弱监督学习[7]发展迅速，此外通过学习图像内在特征[8-11]学习配准参数也已经取得高效性能，这些方法采用卷积神经网络求解用于配准的网络模型，直接估计输入图像对之间的目标位移场。与传统的配准方法相比，基于深度学习的图像配准技术极大缩短了运行时间。这些方法在单模态配准问题上实现了快速配准和可比较的配准精度，如VoxelMorph方法[6]作为无监督单模态图像配准框架，在脑部数据集上取得了优越的性能。国内外****基于深度学习方法对多模态配准问题进行了大量研究，但由于多模态配准局限性，该领域仍有较大的发展潜力。对多模态图像配准来说，由于不同模态图像的强度分布未知且复杂，大多已有多模态图像配准方法严重依赖标签数据，例如真实的变形场或解剖标志，基于弱监督学习实现多模态配准的研究工作[7]，通过从解剖标签中包含的更高级别的对应信息中推断体素级别的转换，但获取真实解剖标签仍有较大困难。相关研究[8]探索了无监督的深度学习方法，基于图像转换方法将图像嵌入到低维的潜在空间中并在其中学习变形场，该方法在2D的多模态配准问题上取得了令人满意的性能。此外，部分研究[9]通过潜在空间特征编码的无监督域自适应技术打破不同模态之间的配准障碍，但是此方法在多模态实验上表现仍不够优异。
为应对以上挑战，本文提出一种完全无监督学习的深度多模态可变形图像配准方法，该方法在训练过程中不需要任何真实的变形场或对齐的多模态图像对。首先通过对不同模态的图像进行基于损失映射量的特征学习，然后基于最大后验概率对变形场进行估计。值得一提的是，本文借助空间转换器[12]和可微分的互信息损失函数最大化图像间的相似性，同时在配准字段上施加平滑约束，从而实现无监督学习。和以往方法相比，本文的主要贡献如下：
1) 提出了一种基于无监督学习的深度多模态可变形图像配准框架，该框架不需要任何对齐的图像对或解剖学标记数据。
2) 所提出的框架针对多模态配准问题首先分别学习不同模态下的图像信息，提取不同分辨率下的图像特征，并通过构建更具辨别力的特征损失映射量，然后对其进行处理，利用最大后验概率求解变形场。
3) 本文方法具有普适性，所提出的框架和损失函数可以转移到各种不同模态的医学图像应用中。
在脑部MRI T1、MRI T2以及CT-MRI上的3D多模态图像配准任务上，将所提出的方法与现有先进的多模态配准方法进行了比较。此外，还在最新的COVID-19的3D肺部CT的单模态配准任务上展示了本文方法的配准性能。大量的实验结果证明，相比于现有最好的配准方法，本文方法能够在更短的运行时间内得到具有竞争力的配准精度，并且在多模态与单模态医学图像配准中表现同样优异。
1 相关知识 不同图像配准允许图像变换具有不同的自由度，刚性变换和仿射变换自由度低，通常用作全局对齐的初始变换。可变形图像配准是一种非线性配准过程，旨在图像之间建立密集的非线性空间对应关系，允许更高的转换自由度。本文关注可变形的多模态配准问题。令来自不同模态的移动图像M与固定图像F作为输入，φ作为配准字段，也称作变形场。典型的可变形配准问题可以表示为
 | (1) |
式中：
在大多数可变形图像配准设置中，仿射变换已考虑在内，图像之间的不匹配来源仅包括非线性的。在本文中，遵循这一假设，所有实验数据均已经过仿射变换。
2 基于深度学习的医学图像配准 2.1 方法流程 本文提出了一种完全无监督的深度多模态可变形图像配准方法，整体流程如图 1所示。首先，在特征提取的过程中分别对2种模态图像进行多尺度提取特征；其次，进行损失映射量计算，通过概率变形场生成模块对损失映射量的处理，可以由粗到细地得到变形场；最后，通过空间转换网络得到最终的配准图像，并使用评估方法对配准图像进行评估，利用损失函数完成无监督训练。
 |
| 图 1 由基于匹配量的特征学习和基于最大后验概率的变形场学习组成的无监督多模态配准框架 Fig. 1 Unsupervised multimodal registration framework composed of feature learning based on matching amount and deformation field learning based on maximum posterior probability |
| 图选项 |
2.2 网络设计
2.2.1 特征提取 本文提出的特征提取方法是基于金字塔模型进行的多尺度特征提取，借鉴了文献[13]中所提的设计思想。值得一提的是，为了对不同模态的图像提取到更多信息，本文利用不同的网络学习不同模态下的特征。本文中多模态实验在CT和MRI图像上进行。通过特征提取模块可以形成多层不同分辨率下的特征图像组成的特征金字塔，其中金字塔底部是输入图像，使用卷积进行下采样，在每次卷积之后进行激活，最后分别输出由小到大的金字塔形特征图像。在本文中共使用2个尺度的特征图像。通过提取特征金字塔的方式，可以得到关于原图像更多的特征信息，并且对不同低分辨率的图像操作使得整个网络更轻便快速。针对不同的模态图像使用不同特征提取器，这也是本文方法在实现多模态配准中的显著优点。
2.2.2 损失映射量计算 由于不同模态图像包含信息以及在成像方式上的巨大差别，因此本文在特征学习过程中引入更具辨识力的损失映射量，用于衡量不同模态特征的不匹配度，并在变形场生成中发挥作用，该模块借鉴光流领域中TV-L1能量模型[14]。具体来说，损失映射函数为L1范数，通过计算移动图像与固定图像在特征金字塔中相同层次的特征图像的L1范数来表示特征匹配中的误差，该模块定义如下：
 | (2) |
式中：i为处在特征金字塔的层次；fixed为固定图像特征；moving为移动图像特征。式(2)即定义为fixed图像特征和moving图像特征之间的数据残差项得到对应层次特征关联之间的匹配成本。
2.2.3 基于概率的变形场生成 变形场生成模块对提取到的多个层次的特征损失映射量进行处理，得到不同分辨率的变形场。由于不同模态图像保存的空间位置信息以及保存图像方式的不同，在配准过程中预测变形场有更大的挑战，本文提出将最大后验概率应用到多模态图像配准中以获得最佳变形场。
变形场生成网络构成由卷积、激活函数交替组成，可对多模态图像特征进行融合处理。随后对于不同层次的特征进行转换，通过求解最大后验概率得到所求当前尺度的变形场。和以往直接求得变形场相比，本文提出的变形场估计网络求解由变形场平均值和方差组成的近似后验概率参数，经过采样后输出当前层次的变形场。在预测得到变形场的过程中，设z为所求的变形场，对变形场概率分布进行建模：
 | (3) |
通过让
 | (4) |
所提出的方法旨在最大化后验概p(z|x; y)。具体来说，通过最大后验概率估计获得当前图像对(x, y)的最佳变形场φz。由于直接求后验概率p(z|x, y)并不容易，使用变分方法，并引入由ψ参数化的近似后验概率qψ(z|x; y)，在不断优化的过程中，将散度损失KL差异最小化：
 | (5) |
本文将近似后验概率qψ(z|x; y)建模为多元正态分布，通过ψ参数化神经网络，输出μz|x, y、
2.3 损失函数 损失函数衡量对模型输出的预测值与真实值之间相似性匹配的程度，以及变形场的光滑程度。通过不断优化损失函数从而指导模型的学习过程。
互信息是信息论中的重要概念，描述了2个系统之间的相关性，或互相包含信息的多少。在2幅图像中，通过熵以及联合熵的概念反映图像间信息相互包含的程度。互信息越大则表明2幅图像间的相关性越大。对于移动图像M、固定图像F来说，其互信息表示为
 | (6) |
式中：p(m, f)为M和F的联合概率分布函数；p(m)p(f)为M和F的边缘概率分布函数。
当2幅图像相似度越高或重合部分越大时，其相关性越大，互信息越大。对于多模态医学图像来说，不同模态图像在成像过程中存在显著差异，因此使用互信息评估两幅不同模态图像配准程度是有效且恰当的方式。通过互信息进行描述提取特征的匹配性[15-16]也是当前研究下的热点。
在本文的设计中，通过最大后验概率预测得到变形场，在预测得到变形场的过程中设z为所求的变形场，在此过程中涉及到KL散度损失，如式(5)所示，对其进行化简，将此部分损失定义如下：
 | (7) |
式中：Eq表示期望。最小化M和F图像之间的差异将鼓励M(φ)近似F，但可能导致不连续。所以在其空间梯度上使用扩散正则化器保证平滑：
 | (8) |
式中: Ω代表几何空间；Δ φ(p)为在点p处求梯度。
对于多模态图像配准，本文使用互信息损失、正则化项、KL散度设计损失函数，多模态配准使用的损失函数如下所示：
 | (9) |
式中：λ为正则化参数。
3 实验结果分析 本文使用Tensorflow框架设计整体网络。在多模态图像配准实验中，经历2次特征提取后，特征图像分辨率分别变为原图像的1/2和1/4。使用Adam优化算法，并将学习率设为0.000 1。由于实验环境内存大小的限制，在本文的所有实验中仅将批处理大小设置为1。
本文多模态医学图像配准实验在脑部MRI T1、MRI T2扫描数据、头部CT扫描数据。将所有配准实验中的数据分为3个部分，即训练数据、验证数据、测试数据，比例为0.75∶0.05∶0.2。对于头部MRI数据，首先将其进行预处理，进行仿射空间归一化和脑部提取，将所有的处理后数据大小设置为(160，160，160)。
此外，还在最新的肺部CT扫描数据进行实验，将COVID-19数据作为模板图像，验证本文方法在单模态图像配准中的可行性。对于肺部CT数据，将切片数据转换为三维数据，并将图像大小设置为(512，512，16)。
3.1 数据集 本文使用的实验数据包括以下几个数据集：
ADNI[17]：用于阿尔茨海默病早期诊断的方法，由阿尔茨海默病神经影像学倡议。
NSCLC[18-19]：美国癌症影像档案馆NSCLC-Radiomics馆藏的402幅CT扫描划定的左右胸部体积分割，胸部分割中包含肺部实质等数据。
CQ500[20]：收集了包含头部CT扫描以及来自各个中心临床报告的数据集。
COVID-19[21]：来自于COVID-19-CT-Seg-Benchmark，为新型冠状病毒肺炎分割数据。
3.2 评估指标
3.2.1 Dice Dice评分是一种集合相似度度量函数，通常用于计算2个样本的相似度。对于医学图像配准来说，如果移动图像通过和变形场拟合得到配准图像，那么配准图像与固定图像之间的医学解剖结构应该可以最大限度的重合，并可量化2幅图像之间的体积重叠，因此通过计算配准后图像与固定图像之间的Dice值来进行评估最终的配准效果：
 | (10) |
Dice得分为1时表示结构完全重叠，Dice得分为0表示无重叠，在评估中Dice得分越大越好。
3.2.2 平均幅度 为了量化变形场的平滑性，定性地显示了变形场的雅可比行列式，同时计算了雅可比行列式梯度的平均幅度(Grad Det-Jac)，即通过雅可比行列式梯度的平均幅度来定量表示变形场的平滑性，幅度越小，表示变形场变形越平滑。
3.2.3 均方误差 均方误差(Root Mean Square Error, RMSE)是一种通过对2幅图像中每个像素点对应关系的度量来计算图像相似度的方法。RMSE越小说明模型选择和拟合更好，数据预测也越成功。RMSE定义如下：
 | (11) |
式中：F(p)和M(p)为固定图像与移动图像在p点处的值。
3.3 比较方法 医学影像研究界已经开发了基于深度学习的方法，并在包括图像配准在内的许多应用中实现了最先进的技术。选取性能表现优异的深度学习方法VoxelMorph[6]作为比较模型。同时，选取2个ANTs[22]传统方法中的对称归一化方法SyN[23]、Elastix[24]作为本文的对比方法。
3.4 实验结果 本文共进行3个方面的实验，接下来将详细介绍每种实验配置以及结果：
1) 同一主体MRI T1和MRI T2的多模态配准。
2) 将CT脑部图像作为atlas模板，使用任意3D MRI T2图像匹配到模板图像。
3) 将COVID-19图像作为模板，使用任意CT肺部图像匹配到模板图像。
3.4.1 MRI T1-MRI T2配准 在磁共振成像中，MRI T1观察解剖结构较好。很多病灶的MRI T2信号要强于周围的正常组织，故MRI T2显示组织病变较好。将MRI T1与MRI T2图像相配准可以很好地在较好的解剖结构中观察病变组织。将同一主体的MRI T1、MRI T2图像分别设置为固定图像与移动图像，表 1展示了分别将MRI T1、MRI T2作为固定图像时的RMSE配准结果，affined为不经过配准，直接使用评价指标度量的结果(粗体代表最优结果)，图 2为可视化效果。
表 1 MRI T1-MRI T2 RMSE配准结果 Table 1 MRI T1-MRI T2 RMSE registration results
| 方法 | 配准方向 | |
| MRI T2->MRI T1 | MRI T1->MRI T2 | |
| affined | 17.549(3.116) | 17.549(3.116) |
| VoxelMorph | 16.844(3.276) | 17.545(2.982) |
| Elastix | 17.731(3.140) | 17.717(3.111) |
| ANTs (SyN) | 16.932(3.141) | 17.206(3.078) |
| 本文 | 16.725(3.303) | 16.915(3.223) |
表选项
 |
| 图 2 MRI T2->MRI T1配准结果可视化 Fig. 2 Visualization of MRI T2->MRI T1 registration results |
| 图选项 |
3.4.2 CT-MRI配准 在多模态CT-MRI配准实验中，选择CT头部数据作为atlas图像，MRI图像作为移动图像。由于没有医学图像分割标签，使用RMSE作为评价指标。此实验在测试时，使MRI T2数据向模板图像进行配准。
表 2显示了将CT脑部数据作为模板图像进行实验的结果(粗体代表最优结果)，包括均值和标准差，进一步证明本文提出方法在多模态实验表现优异。
表 2 CT-MRI配准结果 Table 2 CT-MRI registration results
| 方法 | RMSE |
| affined | 29.932(0.819) |
| VoxelMorph | 28.196(4.374) |
| Elastix | 27.252(0.961) |
| ANTs (SyN) | 29.401(1.021) |
| 本文 | 25.886(0.374) |
表选项
3.4.3 CT-CT配准 为验证本文框架的普适性，在COVID-19数据上进行了单模态配准的实验。为了便于观察是否感染，将COVID-19作为atlas，未患病或疑似病人图像配准到atlas以便帮助医生确认是否感染病毒。训练阶段使用移动图像向atlas进行配准，表 3展示了CT配准结果，NCC代表图像的互相关程度(粗体代表最优结果)。
表 3 CT-CT配准结果 Table 3 CT-CT registration results
| 方法 | Dice | RMSE | NCC | Grad Det-Jac/10-3 |
| affined | 0.464(0.091) | 32.714(1.694) | 0.180(0.008) | |
| VoxelMorph | 0.636(0.109) | 21.086(2.313) | 0.224(0.016) | 4.497(1.058) |
| Elastix | 0.522(0.123) | 33.874(2.191) | 0.112(0.021) | |
| ANTs (SyN) | 0.611(0.233) | 31.382(4.169) | 0.133(0.063) | |
| 本文-L1_cost | 0.654(0.077) | 30.135(3.336) | 0.285(0.014) | 0.866(0.004) |
| 本文-Prob | 0.678(0.094) | 28.617(3.232) | 0.283(0.012) | 0.812(0.005) |
| 本文 | 0.708(0.086) | 27.461(3.162) | 0.287(0.014) | 0.812(0.004) |
表选项
实验结果表明，本文的框架在多个评价指标上取得领先性能。尤其是在Dice指标上，本文方法取得十分优异的性能。原因可能是因为在胸部扫描数据上，部分像素点在变形场中需要有较大的空间位置变换，本文方法利用多个分辨率图像，可以保留更多的空间信息，从而优于比较方法。其中本文-Prob方法直接进行变形场预测而不经过概率进行变形场学习，本文-L1_cost方法为去掉损失映射量模块的方法，通过与2组消融实验比较验证本文所提模块的有效性。
表 4为MRI T2配准到MRI T1图像的配准时间比较(粗体代表最优结果)。结果显示，本文方法在达到较高精度的同时比传统方法的时间更低，这是由于本文提出在多个低分辨率特征上进行操作，使得整个网络更加轻便，运行速度较其他深度学习方法更有优势。
表 4 MRI T2->MRI T1配准运行时间 Table 4 MRI T2->MRI T1 registration running time
| 方法 | 时间/s |
| VoxelMorph | 0.406 |
| Elastix | 6.595 |
| SyN | 20.038 |
| 本文 | 0.281 |
表选项
4 结论 1) 本文提出了一种基于深度学习的深度无监督多模态图像配准框架，不需要任何的标签数据。该框架由基于损失映射量的特征学习和基于最大后验概率的变形场学习组成。
2) 相较于监督方法，本文方法借助空间转换函数不需要Ground Truth(GT)，同时所提出的框架和损失函数可以转移到各种不同模态的医学图像应用中，具有普适性。
3) 大量的实验结果证明，相比于其他配准方法，本文方法能够在更短的时间内取得具有竞争力的配准精度。
参考文献
| [1] | MAINTZ J B A, VIERGEVER M A. A survey of medical image registration[J]. Medical Image Analysis, 1998, 2(1): 1-36. DOI:10.1016/S1361-8415(98)80001-7 |
| [2] | ASHBURNER J. A fast diffeomorphic image registration algorithm[J]. NeuroImage, 2007, 38(1): 95-113. DOI:10.1016/j.neuroimage.2007.07.007 |
| [3] | GUIZAR-SICAIROS M, THURMAN S T, FIENUP J R. Efficient subpixel image registration algorithms[J]. Optics Letters, 2008, 33(2): 156-158. DOI:10.1364/OL.33.000156 |
| [4] | MAES F, VANDERMEULEN D, SUETENS P. Medical image registration using mutual information[J]. Proceedings of the IEEE, 2003, 10(91): 1699-1722. |
| [5] | D'AGOSTINO E, MAES F, VANDERMEULEN D, et al. A viscous fluid model for multimodal non-rigid image registration using mutual information[J]. Medical Image Analysis, 2003, 7(4): 565-575. DOI:10.1016/S1361-8415(03)00039-2 |
| [6] | BALAKRISHNAN G, ZHAO A, SABUNCU M R, et al. An unsupervised learning model for deformable medical image registration[C]//IEEE Conference on Computer Vision and Pattern Recognition. Piscataway: IEEE Press, 2018: 9252-9260. |
| [7] | HU Y, MODAT M, GIBSON E, et al. Weakly-supervised convolutional neural networks for multimodal image registration[J]. Medical Image Analysis, 2018, 49: 1-13. DOI:10.1016/j.media.2018.07.002 |
| [8] | QIN C, SHI B, LIAO R, et al. Unsupervised deformable registration for multi-modal images via disentangled representations[C]//International Conference on Information Processing in Medical Imaging. Berlin: Springer, 2019: 249-261. |
| [9] | MAHAPATRA D, GE Z. Training data independent image registration using generative adversarial networks and domain adaptation[EB/OL]. [2020-07-21]. https://doi.org/10./olb/j.patcog.2019.107109. |
| [10] | ZHANG Y X, LIU R S, LI Z, et al. Coupling principled refinement with bi-directional deep estimation for robust deformable 3D medical image registration[C]//Proceedings of the 2020 IEEE 17th International Symposium on Biomedical Imaging. Piscataway: IEEE Press, 2020: 86-90. |
| [11] | LIU R S, LI Z, ZHANG Y X, et al. Bi-level probabilistic feature learning for deformable image registration[C]//International Joint Conference on Artificial Intelligence, 2020: 723-730. |
| [12] | JADERBERG M, SIMONYAN K, KAVUKCUOGLU K, et al. Spatial transformer networks[C]//Conference and Workshop on Neural Information Processing Systems. New York: Curran Associates, 2015: 2017-2025. |
| [13] | SUN D, YANG X, LIU M, et al. PWC-NET: CNNs for optical flow using pyramid, warping, and cost volume[C]//IEEE Conference on Computer Vision and Pattern Recognition. Piscataway: IEEE Press, 2018: 8934-8943. |
| [14] | ZACH C, POCK T, BISCHOF H. A duality based approach for real TV-L1 optical flow[J]. Pattern Recognition, 2007, 4713: 214-223. |
| [15] | 徐峻岭, 周毓明, 陈林, 等. 基于互信息的无监督特征选择[J]. 计算机研究与发展, 2012, 49(2): 158-168. XU J L, ZHOU Y M, CHEN L, et al. An unsupervised feature selection approach based on mutual information[J]. Journal of Computer Research and Development, 2012, 49(2): 158-168. (in Chinese) |
| [16] | MAHAPATRA D, ANTONY B, SEDAI S. Deformable medical image registration using generative adversarial networks[C]//Proceedings of the 2018 IEEE 15th International Symposium on Biomedical Imaging. Piscataway: IEEE Press, 2018: 1449-1453. |
| [17] | MUELLER S G, WEINER M W, THAL L J, et al. Ways toward an early diagnosis in Alzheimer's disease: The Alzheimer's disease neuroimaging initiative(ADNI)[J]. Alzheimer's & Dementia, 2005, 1(1): 55-66. |
| [18] | KISER K J, AHMED S, STIEB S M, et al. Data from the thoracic volume and pleural effusion segmentations in diseased lungs for benchmarking chest ct processing pipelines[EB/OL]. [2020-08-22]. https://doi.org/10.7937/tcia.2020.6c7y-gq39. |
| [19] | AERTS H J W L, WEE L, VELAZQUEZ E R, et al. Data from NSCLC-radiomics[EB/OL]. [2020-08-22]. https://doi.org/10.7937/K9/TCIA.2015.PF0M9REI. |
| [20] | CHILAMKURTHY S, GHOSH R, TANAMALA S, et al. Deep learning algorithms for detection of critical findings in head CT scans: A retrospective study[J]. Lancet, 2018, 392(10162): 2388-2396. DOI:10.1016/S0140-6736(18)31645-3 |
| [21] | MA J, GE C, WANG Y, et al. COVID-19 CT lung and infection segmentation dataset(Version1.0)[EB/OL]. [2020-08-22]. http://doi.org/10.5281/zenodo.3757476. |
| [22] | AVANTS B B, TUSTISON N, GANG S. Advanced normalization tools (ANTs)[J]. Insight, 2009, 2(365): 1-35. |
| [23] | AVANTS B B, EPSTEIN C L, GROSSMAN M, et al. Symmetric diffeomorphic image registration with cross-correlation: Evaluating automated labeling of elderly and neurodegenerative brain[J]. Medical Image Analysis, 2008, 12(1): 26-41. DOI:10.1016/j.media.2007.06.004 |
| [24] | KLEIN S, STARING M, MURPHY K, et al. Elastix: A toolbox for intensity based medical image registration[J]. IEEE Transactions on Medical Imaging, 2010, 29(1): 196-205. DOI:10.1109/TMI.2009.2035616 |
