5.1 选择配准组件
执行配准时,必须仔细选择几个组件,如第2章所述。组件必须在参数文件中指定。 例如:
(Transform "BSplineTransform")
(Metric "AdvancedMattesMutualInformation")
. . .
在表5.1中给出了需要指定的组件的列表,并提出了一些建议。 “配准”组件在第2章中没有提及。配准组件用于连接所有其他组件,并实现了配准的多分辨率方面。 所以,可以说实际上实现了图2.2的块方案。 此外,第二章中没有明确提及“重新采样器”组件。它仅用于在配准后生成变形的运动图像。 目前,elastix中只有一个Resampler可用:DefaultResampler。 因此,该组件不再进一步讨论。
5.2 所有参数概述
elastix的所有可用组件的列表可以在以下位置找到:http://elastix.isi.uu.nl/doxygen/modules.html
可以在elastix网站上找到可以为每个elastix组件指定的所有参数的列表:
http://elastix.isi.uu.nl/doxygen/parameter.html
在该站点,您可以找到如何指定参数以及默认值。 我们试图提出合理的默认值,虽然默认值在所有情况下一定不会奏效。 维基上可以找到成功的参数文件集合:http://elastix.bigr.nl/wiki/index.php/Parameter_file_database
这可能会让您开始使用您的特定应用程序。
5.3 重要参数
按照与2.3节相同的顺序,我们将讨论每个组件的重要参数,并解释表5.1中提出的建议。
5.3.1 Registration
只需使用MultiResolutionRegistration方法,因为多分辨率是一个好主意。 如果您仍然认为您不需要所有这些多分辨率,您可以随时将NumberOfResolutions设置为1。您不必设置其他任何内容。 第5.3.7节更详细地讨论了决议的数量。
5.3.2 Metric
AdvancedMattesMutualInformation通常适用于单模和多模式图像。 它支持度量值和导数的快速计算,如果变换是通过利用其紧凑支持的B样条曲线。 您需要设置直方图箱的数量,这是计算联合直方图所需要的数量。 一个很好的价值取决于你的输入图像的动态范围,但以我们的经验32通常是可以的:
(NumberOfHistogramBins 32)
5.3.3 Sampler(取样)
RandomCoordinate采样器与StandardGradientDescent和AdaptiveStochasticGradientDescent优化器结合使用,这是推荐的优化例程。 这些优化方法可以与每次迭代随机选择的少量样本一起使用,参见第2.7节,这大大降低了配准时间。 将NumberOfSpatialSamples设置为3000。不要低于2000。与在体素网格上绘制样本的采样器(如随机采样器)相比,RandomCoordinate采样器避免了所谓的网格效应[Th'evenaz and Unser,2008]。
第5.3.6节讨论的随机采样器的一个重要选项是:
(NewSamplesEveryIteration "true")
其在每次迭代中强制选择新样本。
RandomCoordinate采样器的一个有趣的选项是UseRandomSampleRegion参数,与SampleRegionSize参数组合使用。 如果UseRandomSampleRegion设置为“false”
(默认值),采样器从整个图像域中抽取样本。 当设置为“true”时,采样器随机选择一个体素,然后在该体素周围的正方形邻域中选择剩余的样本。 邻域的大小由SampleRegionSize(物理坐标)决定。 3D图像的示例:
(ImageSampler "RandomCoordinate")
(NewSamplesEveryIteration "true")
(UseRandomSampleRegion "true")
(SampleRegionSize 50.0 50.0 50.0)
(NumberOfSpatialSamples 2000)
在每次迭代中,随机选择50^3 mm的正方形区域。 在该地区,根据均匀分布选择了2000个样本。 有效地获得了一种局部相似性度量,有时会提供更好的配准结果。 有关这种方法的更多信息见Klein el al. [2008]。 对于样本区域大小,尝试的合理值是≈图像大小的1/3。
5.3.4 Interpolator(插补)
在配准期间,使用LinearInterpolator。 在我们目前的实现中,它比第一阶B样条内插器快得多,即使它们在理论上是一样的。
我们建议使用更高质量的三阶B样条插值来生成产生的变形运动图像:
(ResampleInterpolator "FinalBSplineInterpolator")
(FinalBSplineInterpolationOrder 3)
5.3.5 Transform(转变)
这个选择取决于手头的应用。对于您希望没有非刚性变形的同一患者的图像,您可以考虑刚性变换,即选择EulerTransform。如果要补偿比例差异,请考虑仿射变换:AffineTransform。这两个转换需要旋转中心,可以由用户设置。默认情况下固定图像的几何中心,这是推荐的。需要设置的另一个参数是Scales。量表为变换参数μ的每个元素定义了在优化期间使用的缩放值。缩放用于将μ的元素置于相同的范围内(对应于旋转的参数通常具有比对应于平移的参数小得多的范围)。我们建议让elastix自动计算:(AutomaticScalesEstimation “true”)
。在做非刚性变换之前,始终要先进行刚性或仿射变换,以获得良好的初始对齐。
对于非刚性配准问题,elastix具有BSplineTransform。 B样条非刚体变换由控制点的统一网格定义。 该网格由网格节点之间的间距定义。 间距定义网格是多么密集,或者您可以模拟的转换的位置。 对于每个分辨率级别,您可以定义不同的网格间距。 这就是我们所说的多网格。 一般来说,我们建议从粗糙的B样条网格开始,即更全面的变换。 这样,较大的结构首先匹配,原因与为什么要开始一个刚性或仿射变换相同的原因。 在以后的分辨率中,您可以逐步细化转化; 这个想法是,您随后匹配较小的结构,达到最终精度。 最终网格间距用以下指定:
(FinalGridSpacingInPhysicalUnits 10.0 10.0 10.0)
数量与图像中的尺寸一样多。 大多数医疗图像的间距是以毫米为单位。 也可以以体素单位指定网格:
(FinalGridSpacingInVoxels 16.0 16.0 16.0)
如果最终的B样条网格间距选择为高,则不能匹配小结构。 另一方面,如果网格间距非常低,则可以匹配小结构,但是您可能允许转换具有太多的自由度。 这可能导致不规则的变换,特别是对于图像的均匀部分,因为在可以引导配准的区域中没有边缘(或其他信息)。 惩罚或正则化术语,见公式(2.2),可以帮助避免这些问题。 由于这取决于所需的精度,因此很难推荐最终网格间距的值。 但是我们可以尝试:如果您对较大的结构感兴趣,可以将其设置为32个体素,您可以下调到16或8个体素以便匹配较小的结构,甚至可以达到4。最后的选择可能需要一些正规化术语,除非你仔细地逐渐改进网格间距。
要指定多格式调度,请使用GridSpacingSchedule命令:
(NumberOfResolutions 4)
(FinalGridSpacingInVoxels 8.0 8.0)
(GridSpacingSchedule 6.0 6.0 4.0 4.0 2.5 2.5 1.0 1.0)
GridSpacingSchedule定义所有分辨率级别的乘法因子。 结合最终网格间距,确定所有分辨率级别的网格间距。在二维图像的情况下,上述计划指定6×8 = 48像素分辨率的0级网格,通过32和20像素,在最后的8级像素分辨率。 GridSpacingSchedule的默认值使用了2的幂次方案: (GridSpacingSchedule 8.0 8.0 4.0 4.0 2.0 1.0 1.0 1.0)
(用于2D图像)。
作为旁注:在(2.1)中最小化的参数的数量由μ的大小确定,即在通过控制点网格间隔的B样条变形变换的情况下。 如果您将间距加倍,则3D图像的参数数量将增加8倍。 对于2563图像和16个体素的网格间距,这将导致大约(256/16)3×3≈12.000个参数; 对于8个像素的网格间距,这几乎是100.000个参数。 参数的数量可以直接关系到内存消耗和配准时间,具体取决于具体的实现。
在大多数文献中,立体(3阶)B样条用于图像配准。 当然,其他spline orders也是可能的。 您可以尝试使用BSplineTransformSplineOrder选项。 支持1,2和3 次。 较低的次将减少计算时间,但可能导致较不平滑的变形。 对于1次,严格来说,变形场不再是可微分的。
作为B样条变换的替代方案,elastix包括SplineKernelTransform,它实现了薄板样条类型的变换。 有关此转换的更多信息,请参见第2.6节和第6.4节。
最后,在每个参数文件中包含以下行是明智的:
(HowToCombineTransforms "Compose")
到elastix版本4.2,默认情况下,如果此行被省略,使用“Add”用于向后兼容。 从elastix版本4.3,默认值已更改为“Compose”,这在大多数应用程序中更好。 参见2.6节,方程(2.18)和(2.19)进行说明。
5.3.6 Optimiser(优化)
StandardGradientDescent方法,参见方程(2.21)和(2.23)提供了执行快速注册的可能性,参见Klein et al[2007]。 关键思想是您使用每次迭代中新选择的体素(样本)的随机子集来计算成本函数导数。 使用参数NumberOfSpatialSamples指定要使用的样本数,请参见第5.3.3节。 通常,2000-5000就够了。 告诉优化器在每次迭代中选择新的样品是非常重要的:
(NewSamplesEveryIteration "true")
StandardGradientDescent方法的缺点是您需要调整增益因子ak的参数,请参见第5.3节。 方程(2.22)需要一个步长ak的选择,它在elastix方程(2.23)中定义。 图5.2给出了一些例子。
参数α和A定义了函数的衰减斜率。 对于参数α,我们建议值为0.6。 对于A,使用的顺序为50:
(SP alpha 0.6)
(SP A 50.0)
这将参数a(在elastix中称为SP_a)作为调整的最重要参数。 这是一个重要的参数,它可能意味着成功的一个好的选择和失败,如果没有! 如果a设置得太高,则迭代求解算法(2.22)变得不稳定,您可能会使图像变形超出识别。 如果一个设置太低,你将永远不会使其达到最佳状态,或者可能会停在非常小的附近局部最佳状态。 图5.2说明了这一点。
a一个很好的选择取决于用于配准的成本函数:为MSD提供一个很好的结果的a与为MI提供良好结果的a不一样。 最后,还取决于您在固定图像和运动图像之间所期望的变形量。 所以再次建议很难给予。 一般来说,,我们建议您考虑数量级,如果a = 10太小,尝试a = 100而不是a = 11。对于互信息,归一化相关系数和归一化的相互信息,您可以从a = 1000左右开始。对于均方差的度量,可以尝试小于1的值。 a选择的方式太大,可能会遇到“运动图像缓冲区外的地图太多”的错误信息。 此错误也可能有其他原因。 elastix网站的常见问题提供了有关此错误消息的更多信息。
对于每个分辨率,您可以指定不同的SP_a值,但是对于每个分辨率,可能会更容易地以相同的值开始。
(SP_a 1000.0 1000.0 1000.0)
或等效地:
(SP_a 1000.0)
与优化器相关的最后一个重要选项是最大迭代次数:
(MaximumNumberOfIterations 500)
在StandardGradientDescent的情况下,不仅最大值也是最小值,因为此优化器没有实现其他停止条件。 一般来说,迭代越多,配准结果越好。 但是,当然,更多的迭代需要更多的时间。 500是一个好的开始。 如果计算时间不是问题,请使用2000。 如果您匆忙,您可以尝试下200次迭代。 对于小的2D图像和刚性配准,甚至更少的迭代就足够了。 使用更多迭代的一个好处是,更广泛的SP_a给出了好的结果。 调节SP然后变得更容易。
你厌倦了调整a了吗?
开始使用AdaptiveStochasticGradientDescent优化器。 此优化器与StandardGradientDescent非常相似,但会自动估计SP_a的正确初始值。更多细节见Klein et al [2009]。 实际上,该优化器可以在许多应用程序中使用其默认设置。 用户必须指定迭代次数:
(Optimizer "AdaptiveStochasticGradientDescent")
(MaximumNumberOfIterations 500)
有一些可选的额外的参数,如SigmoidInitialTime和MaximumStepLength,这在文献[Klein et al,2009]中有解释。 除了自动计算步长之外,该优化器还实现了一种自适应步长机制,这通常使优化有点更强大。 这个优化器的唯一缺点是,在自由度很大(μ大的情况下)估计a是相对耗时的时间,但是我们正在解决这个问题[Qiao et al,2014]。
5.3.7 Image pyramids(图像金字塔)
FixedImagePyramid和MovingImagePyramid具有相同的选项。 下面关于FixedImagePyramid的描述同样适用于MovingImagePyramid。
使用FixedSmoothingImagePyramid,因为它不会丢弃有价值的信息,而且由于您没有使用FullSampler,所以下采样不会随时保存。 尽管对于大图像和许多分辨率级别,它可能消耗相当多的内存。 必须设置两个参数以定义多分辨率策略:分辨率数(NumberOfResolutions)和每个分辨率中使用的特定下采样计划(FixedImagePyramidSchedule)。 如果仅设置NumberOfResolutions,则将使用默认计划,将固定图像在每个维度中平滑2倍,从最后一个分辨率的σ= 0.5开始。 那个计划通常很好。 如果您有高度各向异性的数据,您可能需要在最大间距的方向上模糊较少。
一般来说3 resolutions是一个很好的起点。 如果固定和运动的图像最初是远离的,您可以将分辨率级别的数量增加到例如5或6。这样,图像更加模糊,并且更注意配准大型主导结构。
金字塔计划定义了每个方向x,y,z和每个分辨率级别的模糊量(以及使用FixedRecursiveImagePyramid的情况下的下采样)。 可以指定如下:
(NumberOfResolutions 4)
(FixedImagePyramidSchedule 8 8 4 4 2 2 1 1)
在这个例子中,使用了2D图像的4个分辨率。 在分辨率0时,图像在每个方向上的σ= 8/2体素模糊(σ是金字塔计划值的一半)。 在1级使用σ= 4/2,最后在最后一级4级,原始图像用于配准。 指定的固定和移动图像金字塔具有相同的时间表可以用一个命令:(ImagePyramidSchedule 4 4 2 2 2 1 1 1 1)
对于具有3个分辨率级别的3D图像,其中在z方向上执行较少的平滑。
5.4 Masks(面具,掩盖)
有时候,您有兴趣只对齐图像的一部分。 关注这一部分的一个可能性是裁剪图像。 然而,裁剪仅将感兴趣区域(ROI)限制为正方形(2D)或立方体(3D)。 如果您需要不规则形状的ROI,您可以使用掩码。 掩码是填满0和1的一个二进制图像。 如果使用掩码,则只对掩码内部的图像部分进行配准,即掩码为1。
你可以/应该使用掩码
- 当您的图像包含没有真正意义的人造边缘时。配准可能会对齐这些人造边缘,从而忽略有意义的边缘。 超声图像中的圆锥光束边缘是这样的人造边缘的示例。
- 当图像包含您的ROI附近的结构,可能会影响您的ROI内的配准。 这是例如匹配肺数据的情况。 通常情况下,您对肺部感兴趣,而不是肋骨是否对齐。 然而,肋骨是例如在CT中可以对相似性度量具有强烈影响的结构,特别是如果您使用MSD度量。 在这种情况下,肋骨框架可以以与肋骨框架在肺边界附近的血管结构为代价完全对准。 在这种情况下,如果您使用扩张的肺分割作为掩码,它将帮助您。
掩码可用于固定和运动图像。 固定图像掩码足以将注册聚焦在ROI上,因为从固定图像中抽取样本。 当您的运动图像包含ROI附近的废话灰色值时,您只需要为运动图像使用掩码。
如果您使用掩码来防止人为缘故障,您还需要设置参数:
(ErodeMask "true")
如果没有,那么当执行多分辨率时,由于平滑步骤,来自人造边缘的信息将流入您的ROI。 如果ROI周围的边缘有意义,例如 在肺的例子中,你应该将它设置为false,因为这个边缘将有助于引导配准。
在绘制样品时elastix弹出的常见异常是:“在合理的时间内找不到足够的图像样本。 也许掩码太小了。“这样做的可能原因是你的固定图像掩码太小了。 有关详细信息,请参阅常见问题。
5.5 故障排除
5.5.1 常见错误
常见问题解答中收集了一些常见的混淆和疑问来源,可以在这里找到
http://elastix.isi.uu.nl/FAQ.php
5.5.2 不好的初始对准
当两个图像之间的初始对齐非常糟糕时,您无法启动非刚性配准。 有时使它正确可能是麻烦的。 什么因素可以帮助你做到正确?
- 从具有低自由度的转换开始,即平移,刚性,相似性或仿射变换。 有时,图像真的很远,并且没有重叠开始(注意:物理空间中图像的位置由原点和体素间距决定;参见第2.2节)。 然后,将一个解决方案添加到您的参数文件中:
(AutomaticTransformInitialization "true")
该参数有助于自动估计上述转换的初始对准。 支持三种方法:将固定和运动图像的中心对齐的默认方法,对齐重心的方法以及简单地对齐图像来源的方法。 可以通过将以下一行添加到参数文件中来选择一种方法:
(AutomaticTransformInitializationMethod "GeometricalCenter")
(AutomaticTransformInitializationMethod "CenterOfGravity")
(AutomaticTransformInitializationMethod "Origins")
请注意,“Origins”目前仅适用于仿射变换。
- 您需要一个良好的多分辨率策略,即相当多的分辨率级别。 这样一来,大部分的平滑化正在发生,模糊所有的细节,从而将配准集中在主要的结构上。
- 使用更多的迭代。
- 采取更大的步骤。 将参数a设置得高,以使得变形的转换组件在每次迭代中需要几个体素的步长,最多可达10个。也许它会工作,在这种情况下,您将获得对齐的速度很快,但是步长a仍然是 很大,所以你马上再次离开对准。 如果发生这种情况,您应该让序列ak = a /(A + k)^α相对较快地衰减。 这可以通过设置a和A稍低一点来实现,参见图5.1。
- 如果您需要找到大的旋转,您可能需要为旋转采取更大的步骤,但不需要为变形。 这可以通过修改比例参数来实现,参见第5.3.5节:
(Scales 10000.0)
您可以将其设置为较低以采取较大的旋转步长。 如果设置为1.0,则为旋转方向转动大的步数(但旋转在径向中定义)。 如果你把它设置得很高(> 10^6),你根本不会旋转。 你可能不需要低于1000.0。 请注意,AutomaticScalesEstimation选项通常工作正常,因此不需要指定Scales。
5.5.3 内存消耗
临床图像的典型尺寸随着时间的推移而增加。 因此,内存效率将会更加成为一个问题。 小图像elastix消耗≈100 MB的内存,对于较大的图像对(256^3)和一些常见的组件,消耗可以是大约1-1.5 GB。 使用非常大的图像(400^3及以上),内存消耗可能会升高到旧版32位Windows系统的2 GB以上的限制。 或者您正在使用笔记本电脑。 或者你正在另一个消耗消费的主程序中使用elastix。 大图像怎么办?
- 给自己买一个拥有大内存的全新电脑。 现在计算机和操作系统都将是64位,所以你可以实际解决这么多的内存。
- elastix中的图像在内部默认情况下表示为一堆浮动类型的体素。 您可以将其修改为简短图像:
(FixedInternalImagePixelType "short")
(MovingInternalImagePixelType "short")
这样,您可以节省用于存储固定和移动图像的内存量以及它们的多分辨率金字塔。 这将以牺牲精确度为代价,但可能不是有害的。 这个选项对于elastix和transformix都是有用的。
- 更改配准过程中使用的插值器。 在使用B样条内插的情况下,请注意,它将内部的系数图像存储在双重类型中。 您还可以指定浮动版本:
(Interpolator "BSplineInterpolatorFloat")
这可以节省你的另一个内存大小的短图像。 此选项仅适用于elastix。 为了节省更多的内存,请使用LinearInterpolator。 - 更改重新采样图像时使用的插值器:
(ResampleInterpolator "FinalBSplineInterpolatorFloat")
这个选项对于elastix和transformix都是有用的。 然而,对于elastix,它将只会在配准结束时为您节省一些内存。 - 在配准期间使用下采样图像。 这可能不会太大地影响配准准确性。 配准后,您可以使用transformix将所得到的变换应用于原始的全尺寸运动图像。 有关详细信息,请参阅常见问题。