第二章 图像配准

本章介绍了基于elastix的基本配准概念。 更高级的配准主题将在第6章中讨论。
　　图像配准是医学影像领域的重要工具。 在许多临床情况下，为了分析患者的情况，制作了患者的几张图像。 这些图像采用例如X射线扫描仪，磁共振成像（MRI）扫描仪，计算机断层摄影（CT）扫描仪和超声波扫描仪来获取，其提供关于受试者解剖学的知识。 单一或多方式患者数据的组合通常会产生额外的临床信息，在单独的图像中不明显。 为此，必须找到图像之间的空间关系。 图像配准是从一个图像中的体素到另一个图像中的体素之间找到空间一对一映射的任务，见图2.1。 关于这个问题更多信息可参考Maintz and Viergever [1998], Lester and Arridge[1999], Hill et al. [2001], Hajnal et al. [2001], Zitov´a and Flusser [2003], Modersitzki [2004].

图2.1
　　以下部分介绍了配准过程的数学公式，并概述了一般配准方法的组成部分。 之后，在2.3-2.8节中，对每个组件进行了更详细的讨论。 对于每个组件，elastix使用的名称以打字机样式给出。 在2.9节中，讨论了评估配准结果的方法。

2.1 配准框架

配准过程涉及两个图像。变形浮动图像Im(x)以适应于固定图像If(x)。浮动图像和固定图像是d维的，各自被定义于自己的空间域：

。配准就是找到移位u(x)是的Im(x+u(x))在空间上与If(x)对齐。换句话说，就是找到转换T(x)=x+u(x)，使得Im(T(x))空间上与If(x)对齐。配准被定义为从固定图像到运动图像的映射，如：
对齐的质量由距离或相似性度量S定义，例如平方差之和（SSD），相关比或互信息（MI）度量。 由于这个问题对于非刚性变换T是不正确的，常常引入限制T的正则化或惩罚项P。
　　通常，配准问题被转换为优化问题，其中代价函数C被最小化为w.r.t. T：
其中γ重量与规律性相似（where γ weighs similarity against regularity.）。
为了解决上述最小化问题，基本上有两种方法：参数和非参数。 参考Fischer and Modersitzki [2004]概述了有关非参数的方法，但在本手册中没有讨论。elastix软件是基于参数方法。在参数化方法中，通过引入变换的参数化（模型）来限制可能的转换次数。 原来的优化问题就变成：
下标μ表示变换已被参数化。矢量μ包含“变换参数”的值。 例如，当将变换建模为2D刚体变换时，参数向量μ包含一个旋转角度以及x和y方向上的平移。 方程（2.3）也可以写为：
从这个方程式可以看出，原来的问题（2.1）已经被简化了。 不是优化“功能空间T”，我们现在优化μ。 其他变换模型的例子在2.6节中给出。
图2.2 图2.2显示了块方案中参数化配准算法的一般组成部分。 该计划是Ib'a~nez[2005]等人提出的计划的一个稍微扩展的版本。 方程（2.1） - （2.4）可以识别几个组成部分; 有些将在后面介绍。 首先，我们有图像。 需要定义图像的概念。 这在2.2节中完成。 然后我们有代价函数C或“度量”，它定义了对齐的质量。 如前所述，代价函数由相似性度量S和正则化项P组成。本章没有讨论正则化项P，而是在第六章中讨论。2.3节讨论相似性度量S。相似性度量的定义引入了采样器组件，它在2.4节中进行了说明。 变换模型Tμ的一些例子在2.6节中给出。 实际解决（2.4）问题的优化过程在第2.7节中进行了说明。 在优化期间，在非体素位置评估IM（Tμ（x））值，需要强度插值。 第2.5节介绍了插值器的选择。 方程式（2.1） - （2.4）中不清楚的另一件事是使用多分辨率策略来加速配准，并使其更加强大，参见第2.8节。

2.2 图像

由于图像配准是关于图像的，所以我们必须小心图像的意思。 我们采用 Insight Toolkit [Ib´a˜nez et al., 2005, p. 40]中的图像概念:

关于图像的附加信息被认为是强制性的。 特别地，与像素之间的物理间距和空间中的图像相对于某个世界坐标系的位置相关联的信息是非常重要的。 图像原点和间距是许多应用程序的基础。 例如，配准在物理坐标中执行。 不正确的间距和起始点将导致这些过程中不一致的结果。 没有空间信息的医学图像不应用于医学诊断，图像分析，特征提取，辅助放射治疗或图像引导手术。 换句话说，缺乏空间信息的医学图像不仅没有用，而且也是危险的。
　图2.3阐明了与itk :: Image相关联的主要几何概念。 在该图中，圆圈用于表示像素的中心。 假设像素的值作为位于像素中心的狄拉克三角函数存在。 在像素中心之间测量像素间隔，并且可以沿着每个维度不同。 图像原点与图像中第一个像素的坐标相关联。 像素被认为是围绕保持数据值的像素中心的矩形区域。 这可以被视为图像网格的Voronoi区域，如图的右侧所示。 图像值的线性插值在Delaunay区域内执行，该区域的边角是像素中心。

图2.3 因此，您应该注意，您使用的图像格式，能够存储相关信息（如MHD，DICOM）。一些图像格式，如BMP，不存储的原点和间距。这可能会导致严重的问题！
　　至于elastix版本4.2，在elastix中，图像方向（方向余弦）尚未完全支持。 从elastix 4.3，完全支持图像方向，但由于向后兼容性原因可以禁用图像方向。

2.3 指标

在文献中可以找到几种相似性度量的选择。 下面介绍一些常见的选择。 在括号中，给出了elastix中度量的名称：
　　均方差（MSD）：（AdvancedMeanSquares）MSD定义为：

ΩF为固定图像IF的域，|ΩF| 体素的数量。 给定变换T，可以通过循环固定图像中的体素，通过插值IF（xi），计算IM（Tμ（xi））并将平方差加到和来容易完成运算。

　　归一化相关系数（NCC）：（AdvancedNormalizedCrerelation）NCC定义为：
　　相互信息（MI）：（AdvancedMattesMutualInformation）MI定义为：
　　规范化互信息（NMI）：（NormalizedMualualInformation）NMI定义为NMI =（H（IF）+ H（IM））/ H（IF，IM），H表示熵。 这个表达式可以与MI的定义MI = H（IF）+ H（IM）-H（IF，IM）进行比较。 再次，由2.8定义的联合概率（使用B样条Parzen窗口），NMI可写为：

　　Kappa统计（KS）：（AdvancedKappaStatistic）KS定义为：

　　MSD测量是仅适用于具有相同强度分布的两个图像的度量，即对于来自相同模态的图像。 NCC不太严格，它假定固定和运动图像的强度值之间呈线性关系，因此可以更频繁地使用。 MI度量更为普遍：只要假设固定和运动图像强度的概率分布之间的关系。 对于MI，众所周知它不仅适用于单模态，而且适用于多模态图像对。 这种测量通常是图像配准的好选择。 NMI测量就像MI一样，适用于单模和多模态配准。 Studholme等[1999] 在某些情况下，似乎表明比MI更好的表现。 KS测量专门用于配准二进制图像（分段）。 它衡量分段的“重叠”。

2.4 图像采样器

在等式（2.5） - （2.8）中，我们观察到固定图像上的循环：

。 直到现在，我们假设循环遍历固定图像的所有体素。 一般来说，这不是必需的。 一个子集可能就足够了[Th'evenaz and Unser，2000，Klein et al, 2007]。 可以以不同的方式选择子集：随机，网格等。采样器组件表示此过程。
　　通常使用以下采样器：

• Full：（Full）全采样器只需选择固定图像的所有体元坐标xi。
• Grid：（Grid）网格采样器定义固定图像上的常规网格，并选择网格上的坐标xi。 实际上，网格采样器因此下降了固定图像（不在平滑之前）。 网格的大小（或等效地，下采样因子，其是原始固定图像大小除以网格大小）是用户输入。
• Random：（Random）随机取样器从坐标为xi的固定图像中随机选择用户指定数量的体素。 每个体素都有相同的机会被选中。 样品不一定只选一次。
• Random Coordinate：（RandomCoordinate）随机坐标采样器与随机取样器相似。 它还随机选择用户指定的坐标数xi。 然而，随机坐标采样器不限于体元位置。 也可以选择体素之间的坐标。 当然这些位置的灰度值IF（xi）必须通过插值获得。
  　　虽然乍一看，完整的采样器似乎是最明显的选择，实际上并不总是被使用，因为它在大图像中的计算成本。 随机采样器与随机优化方法相结合特别有用[Klein et al，2007]。 另见第2.7节。 使用随机坐标采样器使成本函数C更为平滑的μ函数，这使得优化问题（2.4）更容易解决。 这已经在Th'evenaz and Unser [2008]中指出。

2.5 插补细分器

如前所述，优化的价值我在（Tµ（x））在非体素的位置进行评估，其强度需要插值。几种插值方法存在，质量和速度不同。如图2.4所示。

• Nearest neighbour: (NearestNeighborInterpolator) 这是最简单的技术，质量低，需要的资源很少。 返回距离最近的体素的强度。
• Linear: (LinearInterpolator) 返回值是周围体素的加权平均值，每个体素的距离取为体重。
• N-th order B-spline: (BSplineInterpolator or BSplineInterpolatorFloat for a memory efficient version) 订单越高，质量越好，还需要更多的计算时间。 实际上，最近邻（N = 0）和线性插值（N = 1）也属于这一类。 有关详细信息，请参阅Unser [1999]。
  　　在配准期间，一阶B样条插值（即线性插值）通常给出令人满意的结果。 质量和速度之间是一个很好的平衡。 为了产生最终结果，即配准的变形结果，通常需要一个更高阶的内插，我们建议N = 3。最后的结果是由所谓的ResampleInterpolator在elastix中产生的。 可以使用上述任一项，但您需要使用Final添加名称，例如：FinalLinearInterpolator。

2.6 变换

关于变换的频繁混乱是它的方向。在elastix变换中，变换被定义为从固定图像域到运动图像域的坐标映射：

。混乱通常源于短语：“运动图像变形以适合固定图像”。虽然人们可以说像这样的图像配准，但是这样的短语并不意味着反映数学基础：一个使运动图像变形，但变换仍然是从固定图像到运动图像。当尝试计算变形运动图像（配准结果）IM（Tμ（x））（此过程通常称为重采样）时，其原因变得清晰。如果从移动图像到固定图像来定义变换，则不会将固定图像域中的所有体素映射到（例如，在缩放的情况下），并且在变形的运动图像中将出现孔。通过定义的转换，重采样非常简单：循环固定图像域ΩF中的所有体素x，计算其映射位置y =Tμ（x），在y处插入运动图像，并将其填充到输出图像中的x位置。
　　用于Tμ的变换模型决定了您可以处理的固定图像和运动图像之间的变形类型。 为了增加灵活性，这些是平移，刚性，相似性，仿射，非刚性B样条和非刚性薄板样条转换。

• Translation: (TranslationTransform) 平移定义为：
• Rigid: (EulerTransform，欧拉转换) 刚性变换定义为：

  矩阵R是旋转矩阵（即正交和正确的），c是旋转中心，并且t再次变换。 图像被视为刚体，可以平移和旋转，但不能缩放/拉伸。 旋转矩阵由欧拉角（2D中的一个，3D中的三个）参数化。 参数矢量μ由欧拉角（rad）和平移矢量（向量）组成。 在2D中，给出长度为3的矢量（向量）：μ=（θz，tx，ty）T，其中θz表示围绕垂直于图像的轴的旋转。 在3D中，给出长度为6的向量：μ=（θx，θy，θz，tx，ty，tz）T。 旋转中心不是μ的一部分; 它是一个固定的设置，通常是图像的中心。

• Similarity: (SimilarityTransform) 相似性变换定义为：

  s为标量，R为旋转矩阵。 这意味着图像被视为一个对象，它可以各向同性地转换，旋转和缩放。 旋转矩阵在2D中以角度参数化，并且通过3D中的所谓的“versor”参数化（也可以使用欧拉角度）。 参数向量μ由角/度，平移向量和各向同性缩放因子组成。 在2D中，给出长度为4的矢量：μ=（s，θz，tx，ty）T。 在3D中，这给出了长度为7的向量：μ=（q1，q2，q3，tx，ty，tz，s）T，其中q1，q2和q3是主语的元素。 当你需要这种转换时，很少有这种情况。

• Affine: (AffineTransform) 仿射变换定义为： 其中矩阵A没有限制。 这意味着图像可以被平移，旋转，缩放和剪切。 参数向量μ由矩阵元素aij和平移向量形成。 在2D中，给出长度为6的矢量：μ=（a11，a12，a21，a22，tx，ty）T。 在3D中，这给出了长度为12的向量。
  我们还实现了仿射变换的另一种风格，具有相同的含义，但使用另一个参数化。 代替由矩阵元素+平移形成的μ，它由d旋转，d剪切因子，d尺度和d平移形成。 定义如下： R，G和S分别为旋转，剪切和缩放矩阵。 它可以使用AffineDTITransform进行选择，因为它首先使用DTI成像，尽管它可以在其他地方使用。 它目前仅在3D中实现，并且还具有12个参数。
• B-splines: (BSplineTransform) 对于非刚性变换的类别，B样条[Rueckert et
  al，1999]通常用作参数化：
  其中xk是控制点，β3（x）是立方三维B样条多项式[Unser，1999]，pk是B样条系数向量（松散地说是控制点位移），σ是B样条控制点间距， Nx在x的B样条的紧凑支持下的所有控制点的集合。 控制点xk定义在常规网格上，覆盖在固定图像上。 在这方面，我们谈论“放在固定图像上的控制点网格”，以及关于“移动的控制点”。 注意，Tμ（xk） 不等于 xk + pk，一个常见的误解。 因此，调用pk控制点位移实际上有些误导。 另请注意，控制点网格与网格图像采样器使用的网格完全无关，请参见第2.4节。
  控制点网格由控制点之间的空间量定义，σ=（σ1，...，σd）（d为图像尺寸），每个方向可以不同。 B样条具有局部支持（| Nx |很小），这意味着一个点的变换只能由几个周围的控制点计算。这对于建模局部变换和快速计算都是有益的。参数μ由B样条系数pk形成。控制点P =（P1，...，Pd）的数量通过M =（P1×...×Pd）×d来确定参数M的数量。 Pi依次由图像尺寸s和B样条网格间距确定，即Pi≈si /σi（我们使用≈，因为一些附加控制点放置在图像外部）。对于3D图像，M≈10000个参数不是异常情况，M可以容易地增长到105 -106。参数向量（2D图像）由以下组成：μ=（p1x，p2x，...，pP1，p1y，p2y，...，pP2）T。
• Thin-plate splines: (SplineKernelTransform) 薄板样条函数是非刚性变换的另一个众所周知的表示。 薄板样条是Davis等人[1997]，布鲁克斯和阿贝尔[2007]更通用的基于类的变换类的实例 。 该变换基于固定和运动图像中的一组K个对应的地标： 。 该变换表示为仿射分量和非刚性分量的和：
  其中G（r）是基函数，ck是与每个地标对应的系数。系数ck和A和t的元素从地标位移dk = 计算。基函数G（r）的具体选择决定了“物理行为”。最常使用的G（r）选择导致薄板样条函数，但另一个有用的替代方法是弹性体样条Davis et al [1997]。样条内核转换通常比B样条更低效（因为它们缺乏B样条的紧凑支持属性），但是它们允许在放置控制点 时具有更大的灵活性。移动的地标形成参数向量μ。需要两个地标集来定义转换。注意，为了执行配准，用户仅给出固定的界标位置;移动的地标被初始化为与身份变换相对应的固定地标，并且随后被优化。参数向量（对于2D图像）组成如下： 。注意与B样条变换相比的μ的排序差异。
  
  有关不同变换的说明，请参见图2.5。 选择适合您需求的变换：如果您希望基本问题包含局部变形，请选择非刚性转换，如果仅需要补偿姿态差异，请选择刚性转换。 要初始化非刚性配准问题，首先执行刚性或仿射。 初始刚性或仿射配准 的结果与以下两种方式之一的非刚性变换 组合： 后一种方法通常是优选的，因为它使得几个后配准分析任务更加直接。

2.7 优化器

为了求解优化问题（2.4），即获得最优变换参数向量μ，通常采用迭代优化策略：

其中dk是迭代k处的“搜索方向”，ak是沿搜索方向控制步长的标量增益因子。 优化过程如图2.6所示。 Klein et al [2007]给出了文献提供的各种优化例程的概述。 例子有准牛顿（QN），非线性共轭梯度（NCG），梯度下降（GD）和罗宾斯 - 蒙罗（RM）。 梯度下降和罗宾斯·蒙罗将在下面讨论。 有关我们参考的其他优化方法的详细信息请参阅 [Klein et al., 2007, Nocedal and Wright, 1999]。

• Gradient descent (GD): (StandardGradientDescent or RegularStepGradientDescent)梯度下降优化方法将搜索方向作为成本函数的负梯度：

  
  
  

  其中g（μk）=∂C/∂μ在当前位置μk评估。 增益因子ak存在若干选择。
  它可以例如由线搜索或通过使用k的预定义函数来确定。

• Robbins-Monro (RM): (StandardGradientDescent or FiniteDifferenceGradientDescent) RM优化方法通过近似g~k代替成本函数g（μk）的导数的计算。
  

  
  近似值可能更快地计算，但可能会降低GD方案的收敛性质，因为每次迭代都会产生近似误差g（μk） - g~k。 Klein et al. [2007]表明，仅使用固定图像中的一个小的随机子集（≈2000）可以显着加速配准，而不会影响配准的准确性。 2.4节中描述的随机或随机取样器是随机抽取体素的采样器的示例。 重要的是每个迭代k选择固定图像体素的新子集，使得近似误差具有零平均值。 RM方法通常与ak组合作为k的预定衰减函数：
  
  其中a> 0，A≥1，0≤α≤1是用户定义的常数。 根据我们的经验，合理的选择是α≈0.6，A约为用户定义的最大迭代次数的10％。 整体增益a的选择取决于μ和g的预期范围，因此是问题特定的。 在我们的经验中，配准结果对这些参数的小扰动不是很敏感。 5.3.6节给出了更多的建议。
  　　请注意，GD和RM其实非常相似。 使用完整采样器运行RM（见第2.4节），而不是随机采样器，相当于执行GD。 我们建议在GD上使用RM，因为它的速度要快得多，而且不会影响精度。 在这种情况下，参数a是要为应用程序调整的参数。 StandardGradientDescent的更高级版本是AdaptiveStochasticGradientDescent，它需要较少的参数设置，并且趋向于更加强大Klein et al. [2009]。
  　　elastix中的其他优化器有：FullSearch，ConjugateGradient，ConjugateGradientFRPR，QuasiNewtonLBFGS，RSGDEachParameterApart，SimultaneousPerturbation，CMAEvolutionStrategy。

2.8 多分辨率

对于多分辨率策略的良好概述，请参阅Lester and Arridge [1999]。 区分两种分层方法：减少数据复杂度，降低转换复杂度。

2.8.1 数据复杂度

通常使用具有较低复杂度的图像来开始配准过程，例如平滑的图像并且可能被下采样的图像。 这增加了配准成功的机会。 一系列具有增加的平滑度的图像称为刻度空间。 如果图像不仅平滑，而且进行了下采样，则数据不仅复杂度较低，而且数据量实际上减少了。 在这种情况下，我们谈论一个“金字塔”。 然而，令人困惑的是，我们也使用金字塔这个字来指代一个尺度空间。 文献中发现了几个尺度空间或金字塔，其中包括高斯和拉普拉斯金字塔，形态尺度空间，样条和小波金字塔。 高斯金字塔是最常见的金字塔。 在elastix我们有：

• Gaussian pyramid: (FixedRecursiveImagePyramid and MovingRecursiveImagePyramid)应用平滑和下采样
• Gaussian scale space: (FixedSmoothingImagePyramid and MovingSmoothingImagePyramid) 适用平滑和无下采样
• Shrinking pyramid: (FixedShrinkingImagePyramid and MovingShrinkingImagePyramid) 不适用于平滑，但仅适用于抽样。

图2.7显示了具有和不具有下采样的高斯金字塔。 与完全采样器（参见第2.4节）相结合，使用下采样的金字塔将在第一个分辨率级别中节省大量时间，因为图像包含少量的体素。 与随机取样器或RandomCoordinate组合，下采样步骤不是必需的，因为随机采样器无论如何都选择用户定义的样本数量，与图像大小无关。

2.8.2 转型的复杂性

第二个多分辨率策略是以较少的转换模式的自由度开始配准。 变换的自由度等于参数矢量μ的长度（元素数）。

第2.6节已经提到了一个例子：在非刚性（B样条）配准之前使用刚性变换。 我们甚至可以使用三级策略：首先是刚性的，然后仿射，然后是非刚性的B样条。

另一个例子是增加转型模型中的自由度。 通过B样条变换，通常很好的做法是开始使用粗略控制点网格进行配准，只能对粗糙变形进行建模。 在随后的分辨率中，B样条网格逐渐细化，从而引入了匹配较小结构的能力。 见5.3.5节。

2.9 评估配准

您如何验证您的配准是否成功？ 这是一个困难的问题。 一般来说，你不知道每个体素应该映射到哪个体素。 这里有一些提示：

• 变形的运动图像IM（Tμ（x））应该看起来与固定图像IF（x）相似。 因此，在观看者中并排比较图像。 您也可以使用棋盘视图或可拖动的十字架将两个图像叠加在一起。 除了看起来相似，还要检查变形的运动图像是否与运动图像具有相同的纹理。 变形图像中突然模糊的区域可能表示该区域的变形太大。
• 对于单模态图像数据，您可以检查差异图像。 完美的配准将导致差异图像没有任何边缘，只是噪音。
• 计算配准后分割解剖结构的重叠。 重叠越好，配准越好。 请注意，这需要您（手动）分段数据中的结构。 为了测量重叠，通常使用骰子相似系数（DSC）：
  

  其中X和Y表示二进制标签图像，和| ·| 表示等于1的体素数。较高的DSC表示较好的对应关系。 值1表示完全重叠，值0表示完全不重叠。 坦尼姆系数（TC）也经常使用。 它与DSC = 2TC /（TC + 1）相关。 另见Crum et al [2006]。 重要的是要认识到，分段结构的表面体积比影响了您通常获得的重叠值[Rohlfing et al，2004]。 DSC = 0.8的值对于复杂血管结构的重叠将是非常好的。 对于大的球形物体，重叠度<0.9通常不是很好。 什么是好的，当然取决于你的应用程序。

• 在您知道对应的点之间计算配准后的距离。 您可以通过在固定和运动图像中手动点击相应点来获得相应的点。 Murphy等人采用半自动化方法，更少的时间选择。 Murphy et al [2008]，旨在找到肺部相应的点。 理想情况下，登记已经发现与地面事实相同的信件。
• 通过查看Tμ（x）的Jacobian的行列式来检查变形场。小于1的值表示局部压缩，大于1的值表示局部膨胀，1表示体积保持。 测量是量化的：值为1.1意味着体积增加10％。 如果这个值大大偏离1，你可能会担心（但如果这是您期望的应用程序，可能不会）。 如果是负面的，你在你的转型中有“折叠”，你肯定应该担心。
• 检查收敛，通过计算每次迭代的精确度量值（而不是近似值，当您进行随机抽样时），并绘制它。 例如，对于MSD度量，度量值越低，配准越好。
• 不要使用图像相似性来评估您的配准。 Torsten Rohlfing可以解释为什么Rohlfing [2012]。

2.10 可视化配准

elastix是一个命令行程序，不做可视化。 它需要输入固定和运动图像，并且在配准结束时生成输出（结果）图像。 但是，通常情况下，您需要视觉检查最终结果。 为此，您可以使用外部查看器。 这样一个观察者并没有elastix包，而是一个独立的应用程序，具有可视化的专用功能。 我们在表2.1中列出了一些可视化工具。 所有这些都是免费提供的，有时甚至是开源的。 列表并不详尽。