VoxelMorph: A Learning Framework for Deformable Medical Image Registration(医学图像配准文章研读)

VoxelMorph: A Learning Framework for Deformable Medical Image Registration*(VoxelMorph：一种可变形医学图像配准的学习框架）
发表时间：2019年
发表期刊：IEEE Transactions on Medical Imaging

摘要

我们提出了 voxelmorph，一个基于快速学习的框架，用于可变形的成对医学图像配准。传统的配准方法为每一对图像优化一个目标函数，这对于大数据集或者丰富的变形模型都是耗时的。相对于这种方法和最近的基于学习建立的方法，我们将配准构造为一个函数，映射一对输入图像对到一个变形领域，对齐这些图像。我们通过卷积神经网络对函数进行参数化，并在一组图像上优化神经网络的参数。通过一对新的扫描，VoxelMorph通过直接计算函数来快速计算变形场。。在本文中，我们探讨了两种不同的训练策略。在第一个(无监督的)设置中，我们训练模型最大化基于图像强度的标准图像匹配目标函数。在第二种情况下，我们利用训练数据中可用的辅助分割。我们证明，无监督模型的准确性是可以与最先进的方法相媲美的，同时运行数量级更快。我们还表明，使用辅助数据训练的VoxelMorph提高了测试时间的配准精度，并评估了训练集大小对配准的影响。我们的方法有望加快医学图像分析和处理的速度，同时在基于学习的配准和应用中促进新的方向。我们的代码可以在 https://github.com/voxelmorph/voxelmorph 免费获取。

关键词

Registration、machine learning、convolutional neural networks

Introduction

配准是各种医学影像研究中的一项基本任务，几十年来一直是活跃的研究热点。在可变形配准中，在一对图像（例如3D磁共振（MR）脑扫描）之间建立密集的非线性对应。传统的配准方法在对配准映射实施约束的同时，通过将体素与相似外观对齐来解决每个图像对的优化问题。不幸的是，求解成对优化可能是计算密集型的，因此在实践中很慢。例如，在CPU上运行的最新算法可能需要数十分钟到数小时才能以高精度 [1] – [3] 配准一对扫描。最近的GPU实现将此运行时间缩短到仅几分钟，但每次配准都需要GPU [4] 。

我们提出了一种新的配准方法，从一个图像集合中学习参数化的配准函数。我们使用一个卷积神经网络(cnn)来实现这个函数，它接受两个n维的输入图像，并输出一个映射，该映射将一个图像中的所有体素映射到另一个图像。网络的参数，即卷积核权重，可以只使用感兴趣的数据集中的一个训练集来优化。该过程学习一种公共表示形式，该公共表示形式允许从同一分布中对一对新的图像进行对齐。本质上，我们在训练阶段用一个全局函数优化代替每个测试图像对求解的高昂的优化代价。一个新的测试扫描对的配准是通过简单地评估在给定的图像上的学习函数来实现的，从而导致快速的配准，甚至在 cpu 上。我们将我们的方法实现为一个通用框架，voxelmorph，可以在https://github.com/voxelmorph/voxelmorph上找到。

在基于学习的VoxelMorph框架中，我们可以自由地采用任何可微的目标函数，在本文中我们提出了两种可能的选择。 第一种方法，我们称之为无监督的，仅使用输入图像对和模型计算的配准场。 与传统的图像配准算法类似，这种损失函数量化了两幅图像强度之间的差异和变形的空间规律性。

在整个研究中，我们使用配准 3D MR 脑部扫描的示例。 然而，我们的方法广泛适用于医学成像领域内外的其他配准任务。 我们在一个包含 3,500 多次扫描的多研究数据集上评估我们的工作，其中包含来自不同年龄组的健康和患病大脑的图像。 我们的无监督模型实现了与最先进的配准相当的准确性，同时减少了几个数量级的时间。 使用 CPU 进行 VoxelMorph 配准只需不到一分钟，而在 GPU 上配准则不到一秒，相比之下，最先进的基线在 CPU 上需要数十分钟到两个多小时。

本文扩展了 2018 年计算机视觉和模式识别国际会议上提出的工作的初步版本 [6]。我们在这项工作的基础上扩展了分析，并引入了一个辅助学习模型，该模型可以在训练期间使用解剖分割来改进对分割图不可得的新测试图像对的配准。我们专注于使用两个损失函数和各种设置对 VoxelMorph 算法的行为进行全面分析，如下所示。我们在更多数据集以及基于图集和主题对主题的配准上测试了无监督方法。然后，我们探讨了在训练期间可以使用不同类型和数量的解剖区域分割作为辅助信息的情况，并评估在不可用分割的情况下对测试数据配准的影响。我们提出了一项实证分析，量化了训练集大小对准确性的影响，并展示了特定于实例的优化如何改善结果。最后，我们对超参数选择进行敏感性分析，并将我们的模型解释为摊销优化。

论文编排如下。第二节介绍医学影像配准，第三节介绍相关工作。第4节介绍了我们的方法。第五部分给出了磁共振成像数据的实验结果。我们在第6节中讨论对结果的见解和结论。

Background

在传统的图像配准公式中，一个 (moving或者source) 图像被扭曲以与第二个 (fixed或者target) 图像对齐。图1显示了取自3D MRI体积的二维冠状切片样本，并概述了几个解剖结构的边界。由于自然解剖大脑的变化和健康状态的差异，不同受试者之间存在显著的差异。可变形配准可以比较扫描之间的结构。这样的分析对于了解不同人群的变异性或疾病患者大脑解剖结构随时间的演变很有用。可变形配准策略通常涉及两个步骤:用于全局对齐的初始仿射变换，然后是具有更多自由度的慢得多的可变形变换。我们专注于后一步，其中我们计算所有体素的密集非线性对应关系。

Fig.1 仿射对准后MRI脑数据集中的冠状切片示例。每列是不同的扫描（对象），每行是不同的冠状切片。用不同颜色勾勒出一些解剖区域：亮/暗蓝的左/右白质，黄/红的左/右心室，紫/绿的左/右海马。扫描之间存在显著的结构差异，因此需要一个可变形的配准步骤来分析扫描间变化。

大多数现有的可变形配准算法基于能量函数迭代地优化变换 [7]。令f和m分别表示fixed image和moving image，并令 φ 是将f的坐标映射到m的坐标的配准场。优化问题可以写为:

其中，moφ 表示由φ进行变形的m，函数Lsim（）衡量两个输入之间的图像相似度，Lsmooth（）施加正则化，λ是正则化权衡参数。

Φ，Lsim和Lsmooth有几种常见的公式。通常，φ 的特征是位移矢量场u指定每个体素从f到m的矢量偏移: φ = Id+u，其中Id是恒等变换 [8]。，微分变换模型 φ通过速度矢量场的积分，保留拓扑并保持变换的可逆性 [9]。用于Lsim的常用指标包括强度均方误差、互信息 [10] 和互相关 [11]。当体积具有变化的强度分布和对比度时，后两者特别有用。Lsmooth强制执行空间平滑变形，通常建模为u的空间梯度的函数。

传统算法针对每个图像对优化公式（1）。当配准许多图像时，这很昂贵，例如作为人口分析的一部分。相反，我们假设一个场可以通过数据的参数化函数来计算。我们通过在图像对数据集上最小化（1）形式的预期能量来优化函数参数。本质上，我们用共享参数的全局优化代替了变形场的特定对优化，在其他领域被称为摊销[12]-[15]。一旦估计了全局函数，就可以通过在给定的图像对上评估函数来产生一个场。在本文中，我们使用基于位移的向量场表示，并重点介绍学习框架的各个方面及其优势。然而，我们最近证明了在类似 VoxelMorph 的框架中也可以使用基于速度的表示，该框架也包含在我们的代码库 [5] 中。

Related Work

Medical Image Registration (Non-Learning-Based)

在三维医学图像配准方面有广泛的工作 [8] 、 [9] 、 [11] 、 [16] - [21] 。一些研究在位移矢量场空间内进行优化。这些包括弹性类型模型 [8] 、 [22] 、 [23] 、统计参数映射 [24] 、b样条自由变形 [25] 、离散方法 [17] 、 [18] 和Demons [19] 、 [26] 。不同的拓扑变换在各种计算解剖学研究中都取得了显著的成功。常用的公式包括Large Diffeomorphic Distance Metric Mapping（LDDMM）[9] 、[21] 、[27] - [32] 、DARTEL[16] 、diffeomorphic demons[33] 和标准对称归一化（SyN）[11] 。所有这些基于非学习的方法优化了每个图像对的能量函数(上述方法都是基于非学习的方法，它们通过优化图像对的能量函数实现配准)，导致配准缓慢。最近基于GPU的算法建立在这些概念的基础上，以将算法运行时间减少到几分钟，但要求GPU可用于每次配准[4]， [34] 。

Medical Image Registration (Learning-Based)

最近有几篇论文提出神经网络学习医学图像配准的函数。其中大部分依赖于ground truth 变形场[35] - [39]，其要么通过模拟变形和变形图像获得，要么在扫描对上运行经典配准方法。一些人还使用图像相似性来帮助引导配准 [35] 。虽然有监督的方法是一个很有前途的方向，但是通过常规的配准工具得到的ground truth 变形场作为地面真值可能很难获得，并且可能限制所学习的变形类型。相反，VoxelMorph是无监督的，并且还能够在训练期间利用辅助信息，例如分割（如果这些信息可用）。

最近的两篇论文[40]，[41]首次提出了基于无监督学习的图像配准方法。两者都提出了一个由CNN和空间变换函数[42]组成的神经网络，可以将图像相互扭曲。然而，这两种初始方法仅在有限的图像子集上演示，例如3D子区域[41]或2D切片[40]，并且仅支持小的变换[40]。

最近的一种方法提出了一种分割驱动的损失函数，用于在同一受试者中配准不同的成像模式——T2w MRI 和 3D 超声[43]、[44]。作者证明，仅基于分割图的损失函数可以导致准确的主体内跨模态配准网络。与这项工作平行，在我们的一个实验中，我们展示了在主体到图谱配准训练期间使用分割图。我们分析了不同解剖标签可用性对整体配准质量的影响，并评估了分割和基于图像的损失组合在各种场景中的表现。我们发现基于分割的损失可能会有所帮助，例如，如果输入分割标签与我们评估的相同（与 [43] 和 [44] 一致）。我们还表明，基于图像的损失和平滑度损失仍然是必要的，特别是当我们评估训练期间未观察到的标签上的配准精度时，并鼓励变形规则。

2D Image Alignment

光流估计是二维图像的一个相关配准问题。光流算法返回一个描述二维图像之间小位移的密集位移矢量场。传统的光流方法通常使用变分方法 [45] - [47] 来解决类似于 (1) 的优化问题。更好地处理大位移或外观剧烈变化的扩展包括基于特征的匹配 [48] 、 [49] 和最近邻场[50] 。

近年来，已经提出了几种使用神经网络进行光流估计的基于学习的方法[51]-[56]。 这些算法将一对图像作为输入，并使用卷积神经网络来学习从数据中捕获光流概念的图像特征。 其中一些工作需要以ground truth流场 [52]、[53]、[55]、[56] 的形式进行监督，而我们建立在一些使用无监督目标 [51]、[54] 的基础上。 空间变换层使神经网络能够在不需要监督标签的情况下执行全局参数 2D 图像对齐 [42] 和密集空间变换 [54]、[57]、[58]。 密集估计的另一种方法是使用 CNN 来匹配图像块 [59]-[62]。 这些方法需要对patch进行详尽的匹配，导致运行时间变慢。

我们基于这些想法，扩展了空间变换以实现n维图像配准，并进一步展示了在训练期间利用图像分割如何在测试时提高配准精度。

Method

设f，m是n维空间域Ω⊂Rⁿ 上定义的两个图像体。在本文的其余部分中，我们将重点讨论n=3的情况，但我们的方法和实现与维度无关。为简单起见，我们假设f和m包含单通道灰度数据。我们还假设f和m作为预处理步骤是仿射对齐的，因此图像体之间不对齐的唯一来源是非线性的。许多软件包可用于快速仿射对齐。

我们使用卷积神经网络 (CNN) 对函数 gθ( f, m) = u 进行建模，其中 θ 是网络参数，即卷积层的卷积核。 f 和 m 之间的位移场 u 实际上存储在 n + 1 维图像中。也就是说，对于每个体素 p ∈Ω ， u( p) 是一个位移，使得 f( p) 和 m ◦ φ 对应于相似的解剖位置，其中映射 φ = Id + u 是使用恒等变换和u。

图 2 概述了我们的方法。该网络将 f 和 m 作为输入，并使用一组参数 θ 计算 φ。我们使用空间变换函数将 m 扭曲到 m ◦ φ，从而能够评估 m ◦φ 和 f 的相似性。在测试期间给定未见过的图像 f 和 m，我们通过评估 gθ( f, m) 获得配准场。

Fig.2 方法概述。我们学习函数gθ (?,?) 的参数θ，并将3D图像体m配准到第二个fixed volume f。在训练过程中，我们使用空间变换函数对m和φ进行变换。可选地，可以在训练期间利用诸如解剖分割sf、sm的辅助信息(蓝色框)。

我们使用 (单元素) 随机梯度下降来通过使用训练数据集最小化预期损失函数来找到最佳参数 θ。在这项工作中，我们提出了两个无监督损失函数。第一个捕获图像相似性和场平滑度，而第二个还利用解剖分割。我们在下一节中详细描述我们的CNN体系结构和两个损失函数。

VoxelMorph CNN Architecture

在本节中，我们将描述我们实验中使用的特定架构，但强调各种架构可能同样工作良好，并且确切的架构不是我们的重点。 gθ(·,·) 的参数化基于类似于 UNet [63]、[64] 的卷积神经网络架构，它由具有跳跃连接的编码器和解码器部分组成。

图 3 描绘了 VoxelMorph 中使用的网络，该网络采用通过将 m 和 f 连接到 2 通道 3D 图像中形成的单个输入。 在我们的实验中，输入的大小为 160 × 192 × 224 × 2，但框架不受特定大小的限制。 我们在编码器和解码器阶段应用 3D 卷积，内核大小为 3，步长为 2。每个卷积后面都有一个参数为 0.2 的 LeakyReLU 层。 卷积层捕获输入图像对的层次特征，用于估计 φ。

在编码器中，我们使用跨步卷积将每一层的空间维度减少一半。因此，编码器的连续层对输入进行较粗略的配准操作，类似于传统图像配准工作中使用的图像金字塔。

在解码阶段，我们在上采样，卷积和级联跳过连接之间进行切换，这些连接将在编码阶段学习的特征直接传播到生成配准的层。解码器的连续层在更精细的空间尺度上运行，从而实现精确的解剖对齐。最小层的卷积核的感受野应至少与f和m中的对应体素之间的最大预期位移一样大。在我们的体系结构中，最小层在输入图像 (1/16)³的图像提上施加卷积。

Spatial Transformation Function

所提出的方法部分通过最小化m° φ 和f之间的差异来学习最佳参数值。为了使用基于标准梯度的方法，我们基于空间变压器网络 [42] 构造了一个可微运算来计算m °φ。

对于每个体素 p，我们计算 m 中的（子像素）体素位置 p‘ = p + u§。因为图像值仅在整数位置定义，我们对八个相邻体素的值进行线性插值：

其中 Z(p’) 是 p ‘的邻近体素，并且 d 在Ω的维度上进行迭代。 因为我们可以计算梯度或次梯度（在0处当0的次梯度是0时用绝对值，不是很理解）， 我们可以在优化过程中反向传播误差。

Loss Functions

在本节中，我们提出了两个损失函数：一个无监督损失 Lus，它仅使用输入图像和生成的配准场来评估模型，以及一个辅助损失 La，它也在训练时利用解剖分割。

Unsupervised Loss Function

无监督损失 Lus(·,·,·) 由两部分组成：惩罚外观差异的 Lsim 和惩罚 φ 的局部空间变化的 Lsmooth：

其中 λ 是正则化参数。我们对 Lsim 的两个常用函数进行了试验。第一个是均方体素差异，适用于 f 和 m 具有相似的图像强度分布和局部对比度时：

第二个是f和m的局部互相关，这对扫描和数据集中发现的强度变化更为稳健 [11] 。设f^hat（p）和[m^hat °φ]（p）表示局部平均强度的图像：
pi在p附近的n³ 图像体上迭代，n=9在我们的实验中。f和m°?的局部互相关写为：

较高的CC表示更好的对齐，产生损失函数：Lsim(f，m，φ)=−CC(f，m◦φ)。

将Lsim最小化将鼓励m°φ近似f，但可能产生物理上不现实的非光滑φ。本文在空间梯度位移u上利用扩散正则函数促进光滑位移场φ:

使用相邻体素之间的差异来近似空间梯度。

Auxiliary Data Loss Function

在这里，我们描述了VoxelMorph如何利用训练期间可用的辅助信息，而不是在测试期间。解剖分割图有时在训练期间可用，可以由人类专家或自动算法进行注释。分割图将每个体素指定给解剖结构。如果一个配准域φ表示精确的解剖对应关系，则f和m◦φ中与同一解剖结构相对应的区域应重叠良好。

让S^kf, S^km分别表示f和m°φ中结构k的体素。们使用Dice分数量化结构k的体积重叠 ：

Dice 得分为 1 表示解剖结构完美匹配，得分为 0 表示没有重叠。 我们将所有结构 k ∈[1, K] 上的分割损失 Lseg 定义为:

单独的 Lseg 并不鼓励图像外观的平滑度和一致性，然而这对于良好的配准至关重要。 因此，我们将 Lseg 与 (4) 结合以获得目标：

其中，γ是一个正则化参数。

在我们的实验中，使用了仿射变换对齐图像，我们证明了损失 (10) 可以带来显著的改善。通常，根据任务的不同，也可以按照 [43] 中介绍的多尺度方式计算这种损失，这取决于初始对准的质量。

由于解剖标签是分类的，因此计算 Sm ◦ φ 的线性插值的简单实现是不合适的，并且 (8) 的直接实现可能不适用于自动微分框架。 我们将 Sf 和 Sm 设计为具有 K 个通道的图像体，其中每个通道是指定特定结构的空间域的二进制掩码。 我们通过使用线性插值对 sm 的每个通道进行空间变换来计算 sm ◦ φ。 然后，我们分别通过乘加 Sf 和 Sm ◦ φ 来计算 (8) 的分子和分母。

Amortized Optimization Interpretation

我们的方法将变形场 φ 上的特定对优化替换为函数 gθ(·,·) 的函数参数 θ 的全局优化。 这个过程有时被称为摊销优化[66]。 因为函数 gθ(·,·) 的任务是估计任意两个图像之间的配准，所以参数 θ 全局共享的事实作为自然正则化。 我们在V-C节（正则化分析）中证明了这一点。此外，该函数输出的变形的质量和泛化能力将取决于所训练的数据。实际上，所得到的变形可被解释为简单地近似或初始化最佳变形φ*，并且所得到的差异有时被称为摊销间隙 [15] 、 [66] 。如果需要，可以使用任何特定于实例的优化来改进这个初始变形场。 在我们的实验中，我们通过将产生的位移 u 作为模型参数来实现这一点，并使用梯度下降独立地微调每个特定扫描的变形。 本质上，这实现了传统配准的自动微分版本，使用 VoxelMorph 输出作为初始化。 但是，大多数情况下，我们发现初始变形（VoxelMorph 输出）已经与最先进的结果一样准确。 我们在 V-D 节中介绍的实验中探讨了这些方面。

Experiments

我们在脑MRI配准任务上演示了我们的方法。我们首先(V-B节)提出了一系列基于地图集的配准实验，在这些实验中，我们计算了地图集或参考体与我们数据集中的每个图像体之间的配准场。基于图集的配准是人口分析中的一种常见表述，其中主体间配准是一个核心问题。 地图集代表参考或平均体积，通常通过联合和重复对齐大脑 MR 体积数据集并将它们平均在一起 [67] 来构建。我们使用外部数据集[1], [68] 计算的图集 。每个输入体积对由图集（图像f）和来自数据集的体积（图像m）组成。

图 4 显示了所有示例使用相同固定图集的示例图像对。 在第二个实验（第 V-C 节）中，我们进行了超参数敏感性分析。 在第三个实验（第 V-D 部分）中，我们研究了训练集大小对配准的影响，并展示了特定于实例的优化。 在第四个实验（第 V-E 节）中，我们在包含手动分割的数据集上展示结果。 在下一个实验（第 V-F 部分）中，我们使用随机的训练对象对作为输入来训练 VoxelMorph，并测试未见过的测试对象对之间的配准。 最后（第 V-G 节），我们对使用辅助分割数据的配准进行了实证分析。 本文中所有描绘大脑的图都显示了 2D 切片，但所有配准都是在 3D 中完成的。

Fig.4 从输入对（第 1-2 列）中提取的示例 MR 冠状切片，并为 VoxelMorph使用不同的损失函数生成的 m ◦ φ。 我们覆盖了一些结构的边界：心室（蓝色/深绿色）、丘脑（红色/粉红色）和海马体（浅绿色/橙色）。 良好的配准将使 m◦ φ 中的结构看起来与 f 中的结构相似。 我们的模型能够处理结构形状的各种变化，包括第 2 行和第 3 行脑室的扩张/收缩，以及第 4 行海马体的拉伸。

Experimental Setup

Dataset

我们使用了来自八个公开可用数据集的3731个T1-weighted脑MRI扫描的大规模，多站点，多研究数据集: OASIS [69]，ABIDE [70]，ADHD200 [71]，MCIC [72]，PPMI [73]，HABS [74]，Harvard GSP [75]，和FreeSurfer Buckner40 [1]。每个数据集的采集细节、受试者年龄范围和健康状况都不同。 所有扫描都被重新采样到具有 1mm 各向同性体素的 256×256×256 网格。 我们执行标准预处理步骤，包括使用 FreeSurfer [1] 对每次扫描进行仿射空间归一化和大脑提取，并将生成的图像裁剪为 160 × 192 × 224。所有的磁共振成像都是用FreeSurfer进行解剖分割的，我们使用视觉检查进行质量控制，以捕捉分割结果和仿射对齐中的严重错误。我们包括所有受试者的所有解剖结构，其体积至少为100个体素，形成30个结构。我们使用生成的分割图来评估我们的配准，如下所述。我们将数据集分为3231、250和250个图像体，分别用于训练、验证和测试集，尽管我们强调在任何阶段都不使用任何受监督的信息。此外，Buckner40数据集仅用于测试，使用手动分割。

Evaluation Metrics

为这些数据获得密集的ground truth配准并没有明确定义，因为许多配准场可以产生类似的扭曲图像。我们首先使用解剖分割的体积重叠来评估我们的方法。如果配准场 φ 表示准确的对应关系，则 f 和 m ◦ φ 中对应于相同解剖结构的区域应该很好地重叠（示例见图 4）。我们使用Dice评分量化结构之间的体积重叠 (8) 。我们还评价了变形场的规律性。具体地，Jacobian矩阵
捕捉了p周围φ的局部性质。我们计算了所有|Jφ §| <= 0的非背景体素，这些体素的变形不是可微分的 [16] 。

Baseline Methods

我们使用对称归一化 (SyN)作为第一基线方法 [11]，这是一项比较研究 [2] 中表现最好的配准算法。 我们在公开可用的高级标准化工具 (ANTs) 软件包 [3] 中使用 SyN 实现，并使用互相关相似性度量。在我们处理医学图像的整个过程中，我们发现默认的ANTs平滑度参数对于将ANTs应用于我们的数据来说是次优的。我们使用跨多个数据集的宽参数扫描获得了改进的参数，并在这些实验中使用了这些参数。具体来说，我们使用0.25的SyN步长，高斯参数 (9,0.2)，在三个尺度上，每个尺度最多201次迭代。我们还使用NiftyReg软件包作为第二个基线方法。不幸的是，目前还没有一个GPU实现，而是我们构建了一个多线程CPU版本。我们搜索了各种参数设置以获得改进的参数，并使用 CC 损失函数、5 的网格间距和 500 次迭代。

VoxelMorph Implementation

我们使用带有tensorflow后端 [77] 的Keras [76] 实现了我们的方法。我们将2D线性插值空间变压器层扩展到n-d，这里使用n = 3。我们使用的是学习率为10⁻⁴的ADAM优化器 [78]。虽然我们的实现允许小批量随机梯度下降，但在我们的实验中，每个训练批次都由一对图像体组成。我们的实现包括150,000迭代的默认。我们的代码和模型参数可在https://github.com/voxelmorph/voxelmorph 在线获得。

Atlas-Based Registration

在这个实验中，我们训练VoxelMorph进行基于图集的配准。我们训练具有不同λ正则化参数的独立VoxelMorph网络。然后在验证集上选择优化Dice分数的网络，并在测试集上报告结果。

表 I 显示了针对仅全局仿射对齐、ANTs 和 NiftyReg 以及具有不同损失的 VoxelMorph 的基线为所有受试者和结构计算的平均 Dice 分数。 VoxelMorph 变体在 Dice方面的性能与 ANTs 和 NiftyReg 相当，并且明显优于仿射对齐。来自我们算法的扭曲图像的示例视觉结果显示在图 4 和图6中。VoxelMorph 能够处理各种结构的显著形状变化。

表1
第一次实验的 Affine Alignment、ANTs、NiftyReg 和 VoxelMorph 的平均Dice分数和运行时间结果。括号中是结构和受试者之间的标准偏差。计算所有结构和受试者的平均Dice分数。时间是在预处理后计算的。我们的网络在Dice得分上与ANTs和 NiftyReg 产生可比的结果，而操作数量级更快。我们还显示了每种方法中具有非正雅可比行列式的体素的数量和百分比，我们的大脑中有 520 万个体素。所有方法都显示不到 1% 的此类体素。

图 5 将每个结构的 Dice 分数显示为箱线图。 为了便于可视化，我们将来自两个半球的相同结构的 Dice 分数平均为一个分数，例如，将左右海马分数平均。 VoxelMorph 模型在所有结构上都实现了与 ANTs 和 NiftyReg 相当的 Dice 测量，在某些结构（如侧脑室）上表现稍好，而在其他结构（如海马）上表现更差。

Fig.5 第一个 (无监督) 实验的ANTs，NiftyReg和VoxelMorph结果的各种解剖结构的Dice得分的箱线图。对于这种可视化，我们将左右脑半球的Dice得分平均为一个得分。结构按平均ANTs Dice得分排序。

表I包括Jacobian行列式为非正值的体素计数。我们发现，所有方法都会导致此类体素的小岛变形，但在绝大多数体素 (99.4% - 99.9%) 处是异形的。图6和11 (在补充材料中) 示出了几个示例VoxelMorph变形场。VoxelMorph对微分变形没有明确的约束，但是在这种情况下，平滑度损失会导致通常平滑且行为良好的结果。ANTs和NiftyReg包括可以强制实施或强烈鼓励微分变形的实现，但是在我们的参数搜索过程中，这些负面影响运行时间或结果。在这项工作中，我们使用产生最佳Dice得分的配置运行了基线实现，结果也产生了良好的变形规律性。

Fig.6 在无监督实验（第V-B节）中，通过将moving image（第1列）配准到fixed image（第2列）来提取变形场φ的示例（第4-5列）。扭曲的图像体m◦φ如第3列所示。每个空间维度中的置换都映射到第4列中的每个RGB颜色通道。VoxelMoyph（MSE）产生的变形场在大脑中是平滑的，即使在记录明显不同的moving image和fixed image时也是如此。

Runtime:

表 1 显示了使用 Intel Xeon (E5-2680) CPU 和 NVIDIA TitanX GPU 的运行时结果。我们报告了仿射对齐预处理步骤之后计算的运行时间，所有提出的方法共享，即使在 CPU 上也只需要几分钟。 ANTs 需要两个或更多小时的 CPU 时间，而 NiftyReg对于给定设置需要大约 39 分钟。 ANTs 的运行时间差异很大，因为它的收敛性取决于对齐任务的难度。与 ANTs 相比，使用 VoxelMorph 配准两个图像在 CPU 上的速度平均快 150 倍，比 NiftyReg 快 40 倍。使用 GPU 时，VoxelMorph 会在一秒钟内计算出配准。据我们所知，没有针对 GPU 的公开可用的 ANTs 实现。 SyN 算法很可能会从 GPU 实现中受益，但 VoxelMorph 的主要优势来自不需要对每个测试对进行优化，这可以从 CPU 比较中看出。不幸的是，NiftyReg GPU 版本在所有可用存储库历史的当前源代码中都不可用。

Regularization Analysis

图 7 显示了平滑正则化参数 λ 的不同值下验证集的平均 Dice 分数。 结果在很大范围的 λ 值上平滑变化，说明我们的模型对 λ 的选择是稳健的。 有趣的是，即使设置 λ = 0（对配准不强制进行显式正则化），也会导致对仿射配准的显著改进。 这可能是因为最优网络参数 θ 需要很好地配准训练集中的所有对，从而为函数 gθ(·,·) 生成隐式数据集正则化。

Fig.7 使用不同正则化参数λ的VoxelMorph在验证集上的Dice得分。

Training Set Size and Instance-Specific Optimization

我们评估了训练集大小对准确性的影响，以及摊销和特定于实例的优化之间的关系。 由于 MSE 和 CC 执行基于atlas 的类似配准，因此在本节中我们使用 MSE。 我们在训练数据集中不同大小的子集上训练 VoxelMorph，并在以下方面报告 Dice 分数：(1) 训练子集，(2) 保留的测试集，以及 (3) 当每个变形进一步针对每个测试图像对单独优化时的测试集。 我们通过在每个测试对上使用100次梯度下降迭代从 VoxelMorph 获得的位移 u 进行微调来执行 (3) ，这在 GPU 上花费了 23.7 ± 0.4 秒，在单线程 CPU 上花费了 628.0 ± 4.2 秒。

图8显示了我们的结果。与较大的训练集大小相比，较小的训练集大小（10次扫描）导致训练和测试Dice分数略低。然而，当使用100次扫描或完整数据集进行训练时，Dice得分没有显著差异。无论训练集大小如何，进一步优化每个测试图像对上的VoxelMorph参数都能获得更好的测试Dice得分，与最新技术水平相当。

Manual Anatomical Delineations

由于大多数数据集都无法进行手动分割，FreeSurfer分割的可用性使得上述实验范围更广。在这个实验中，我们使用已经在V-B节中训练过的VoxelMorph模型，在包含39次扫描的（未见过的）Buckner40数据集上测试配准。该数据集包含专家手动描绘的先前实验中使用的相同解剖结构，我们在这里用于评估。我们还使用实例特定优化计算VoxelMorph，如第V-D节所述。如表II所示，Dice评分结果表明，使用互相关损失的VoxelMorph与使用相同损失函数的ANTs和NiftyReg表现类似，这与我们评估FreeSurfer分割的第一个实验一致。与之前的实验类似，具有实例特定优化的VoxelMorph进一步改进了结果。在该数据集上，使用具有MSE损失的VoxelMorph的结果获得的分数略低，但通过特定于实例的优化程序得到了改进，与ANTs和NiftyReg相当。

表2
手工标注实验结果。我们展示了Affine、ANTS、NiftyReg和VoxelMorph，其中“inst”表示额外的特定于实例的优化，如第V-D节所述。计算所有结构和受试者的平均Dice分数，括号中包含结构和受试者的标准差。

Subject-to-Subject Registration

在这个实验中，我们训练VoxelMorph以进行受试者到受试者的配准。由于每次配准都有更多的可变性，我们将每个网络层的特征数量增加一倍。我们还使用实例特定优化VoxelMorph，如第V-D节所述。表III给出了要配准的250个随机选择的测试对的平均测试Dice分数。与文献一致，我们发现与使用MSE损失相比，归一化互相关损失导致更稳健的结果。VoxelMorph（特征数量加倍）Dice得分与ANTs相当，略低于NiftyReg，而VoxelMorph实例特定优化的结果与两个基线相当。

表3
使用Affine、ANTs和VoxelMorph变体对受试者比对的结果，其中 “X2” 是指我们将特征数量增加一倍以解决任务固有变异性增加的模型，以及 “inst”表示额外的特定于实例的优化。

Registration with Auxiliary Data

在本节中，我们在使用损失函数 (10) 的训练期间使用分割图时评估 VoxelMorph。 因为 MSE 和 CC 在基于图集的配准中表现相似，所以在本节中，我们使用 λ = 0.02 的 MSE。 我们在两个实际场景中对我们的模型进行了评估：（1）当解剖结构标签的子集在训练期间可用时，以及（2）当粗分割标签在训练期间可用时。 我们使用与之前的实验相同的训练/验证/测试拆分。

Training With a Subset of Anatomical Labels

在许多实际环境中，可能无法获得所有结构的训练分割。因此，我们首先考虑分割仅可用于30个结构的子集的情况。我们指的是观察到的分割结构，其余的是未观察到的结构。我们考虑了三种情况：观察到一个、15个（一半）和30个（全部）结构分割。前两个实验基本上模拟了不同数量的部分观察到的分割。对于每个实验，我们在观察到的结构的不同子集上训练单独的模型，如下所示。对于单结构分割，我们手动选择了四个重要的结构进行四折实验（每折一个）：海马、大脑皮层、大脑白质和心室。在第二个实验中，我们随机选择了30个结构中的15个，每5个结构都有不同的选择。对于观察到的标签的每一折和每个子集，我们在训练时使用分割图，并在不使用分割图的测试对上显示结果。

图9 a - c显示了在(10) 中扫描辅助正则化权衡参数γ时观察到的标签和未观察到的标签的Dice分数。我们使用FreeSurfer注释训练模型，并在使用FreeSurfer注释的常规测试集（顶部）和使用手动注释的Buckner40测试集（底部）上显示结果。极值γ=0（或logγ=-∞）和γ=∞作为理论极值，γ=0对应于无监督VoxelMorph，γ=∞对应于仅用辅助标签训练的VoxelMorph，没有光滑性和图像匹配客观项。

一般来说，对于Dice得分，具有辅助数据的VoxelMorph显著优于(四种设置中最大p值<10⁻⁹)无监督的VoxelMorph(相当于γ=0或logγ=−∞)和可观察到的结构上的ANTs。观察到的标签上的Dice得分通常随着γ的增加而增加。

有趣的是，与非监督变量相比，VoxelMorph（利用辅助数据训练）在γ值范围内提高了未观察到的结构的Dice分数（见图9a-b)，即使在训练期间未明确观察到这些分割。当在评估阶段使用的所有结构都在训练期间观察时，我们发现在较高的γ值下具有良好的Dice结果（图9 c。).当γ较大时，未观测结构的配准精度开始下降，范围为logγ∈[-3，- 2] 。这可以解释为在这个范围内模型开始与观察到的结构过度匹配，也就是说，它继续提高观察到的结构的Dice分数，同时损害其他结构的配准精度（图9 c)。

Fig.9 在训练期间使用辅助数据时的测试扫描结果。顶部：常规测试集的FreeSurfer分割上的测试。底部：在Buckner40测试集的手动分割上测试相同的模型。我们测试具有不同数量的观察到的标签 (a - c)，以及具有更粗的分割图 (d) 。误差条表示受试者之间的标准偏差。每个图中所有标签的最左边的数据点对应于γ=0，表示不使用辅助数据（无监督）的VoxelMorph结果。通过将图像和平滑项设置为0来实现γ = ∞。我们展示了具有最佳参数的ANTs结果的Dice分数，不使用分割图进行比较。

Training With Coarse Labels

我们考虑只有粗略标签可用的场景，例如当所有白质被分割为一个结构时。 这种情况可以评估辅助数据如何在粗略描绘的结构内以更精细的尺度影响解剖配准。 为此，我们将 30 个结构合并为四大组：白质、灰质、脑脊液 (CSF) 和脑干，并评估原始结构配准的准确性。

图9d (顶部) 显示了不同γ下原始30个结构上的平均Dice得分。在 γ 为0.01的情况下，我们在FreeSurfer分割上获得的平均Dice得分为0.78 ± 0.03。与没有辅助信息的VoxelMorph相比，这大约是3个Dice点的改进 (p-value < 10^-10)。

Regularity of Deformations

我们还通过计算Jacobian的行列式为非正的体素数来直观地评价变形场的规律性。表IV提供了所有γ值的定量规则性度量，表明VoxelMorph变形规则性作为γ的函数缓慢下降（在对数标度上显示），其中大约0.2%的体素在最低参数值处显示折叠，当γ=0.1时至多2.3%。在极值γ = ∞时，不鼓励平滑的模型变形表现出10-13%的折叠体素。因此，较低的γ值例如γ = 0.01提供了所有结构的高Dice分数的良好折衷，同时避免了高度不规则的变形场，并避免了如上所述的过拟合。图10示出了对于γ = 0.01和γ = ∞的变形场的示例，并且对于每个实验设置，我们在补充材料中包括更多的图。
表IV
使用 FreeSurfer、MSE 损失函数和平滑度参数 0.02 获得的辅助分割训练时变形场的规律性，使用计数和非正雅可比行列式体素数量的百分比来测量。

Fig.10. γ对扭曲图像和变形场的影响。我们展示了moving、fixed和扭曲的图像（第1-3列），其中叠加了在训练时刻观察到的结构。第4列和第5列显示了由此产生的变形场。在为观察到的结构提供更好的Dice分数的同时，通过γ=∞训练产生的变形场 远比使用γ=0.01时更不规则。同样，对于γ=∞，扭曲的图像在视觉上也不太连贯.

Testing on Manual Segmentation Maps

我们还在上面使用的 Buckner40 数据集中的手动分割上测试了这些模型，得到了图 9（bottom）。 我们观察到与上述结论一致的行为，Dice 得分改进较小，这可能是由于 Buckner40 数据上获得的基线 Dice 得分较高。

Discussion and Conclusion

在Dice得分方面，具有无监督损失的VoxelMorph的性能与最新的ANTs和NiftyReg软件相当，同时在CPU上的计算时间从数小时减少到数分钟，在GPU上的计算时间减少不到一秒。VoxelMorph是可伸缩的，并且在训练过程中处理部分观察到的或粗略描绘的辅助信息，这可以导致Dice分数的提高，同时仍保留运行时的改进。

VoxelMorph执行摊销优化，学习全局函数参数，这些参数对于整个训练数据集都是最优的。如图8所示，数据集不需要很大: 只有100训练图像，VoxelMorph在训练和测试集上都获得了最先进的配准质量分数。特定于实例的优化进一步将VoxelMorph性能提高了一个Dice点。这是一个很小的增长，说明摊销优化可以获得几乎最佳的配准。

我们进行了一组完整的实验，证明为了合理选择γ，在训练期间解剖分割的可用性显著提高了VoxelMorph的测试配准性能（就Dice分数而言），同时提供平滑变形（例如，对于γ=0.01，小于0.5%的折叠体素）。性能增益取决于可用解剖分割的质量和数量。在训练过程中，如果有一个单一的标记解剖结构，那么在不影响其他解剖结构的情况下，测试对象配准该标记的准确性会提高。如果观察到一半或全部标签，甚至在训练时提供了一个粗略的分割，则在测试期间所有标签的配准精度都会提高。在本研究中，我们使用了一种辅助数据，而VoxelMorph可以利用其他辅助数据，例如不同的模式或解剖关键点。增加γ也会增加显示配准场折叠的体素数量。如[5]所述，通过对VoxelMorph使用不同的微分表示，可以缓解这种影响。

VoxelMorph是一种通用的学习模型，并不局限于特定的图像类型或解剖结构——它可能在其他医学图像配准应用中有用，如心脏Mr扫描或肺CT图像。通过适当的损失函数，如互信息，该模型还可以进行多模态配准。VoxelMorph有望显著加快医学图像的分析和处理管道，同时为基于学习的配准开辟了新的方向。