Registration Overview

The goal of registration is to estimate the transformation which maps points from one image to the corresponding points in another image. The transformation estimated via registration is said to map points from the fixed image coordinate system to the moving image coordinate system.

SimpleITK provides a configurable multi-resolution registration framework, implemented in the ImageRegistrationMethod class. In addition, a number of variations of the Demons registration algorithm are implemented independently from this class as they do not fit into the framework.

Actual Code

Code illustrating various aspects of the registration framework can be found in the set of examples which are part of the SimpleITK distribution and in the SimpleITK Jupyter notebook repository.

Initialization and Center of Rotation

The task of registration is formulated using non-linear optimization which requires an initial estimate. The two most common initialization approaches are (1) Use the identity transform (a.k.a. forgot to initialize). (2) Align the physical centers of the two images (see CenteredTransformInitializerFilter). If after initialization there is no overlap between the images, registration will fail. The closer the initialization transformation is to the actual transformation, the higher the probability of convergence to the correct solution.

If your registration involves the use of a global domain transform (described here), you should also set an appropriate center of rotation. In many cases you want the center of rotation to be the physical center of the fixed image (the CenteredTransformCenteredTransformInitializerFilter ensures this). This is of significant importance for registration convergence due to the non-linear nature of rotation. When the center of rotation is far from our physical region of interest (ROI), a small rotational angle results in a large displacement. Think of moving the pivot/fulcrum point of a lever. For the same rotation angle, the farther you are from the fulcrum the larger the displacement. For numerical stability we do not want our computations to be sensitive to very small variations in the rotation angle, thus the ideal center of rotation is the point which minimizes the distance to the farthest point in our ROI:

\[p_{center} = \underset{p_{rotation}} {\arg\min}\ dist(p_{rotation}, \{p_{roi}\})\]

Without additional knowledge we can only assume that the ROI is the whole fixed image. If your ROI is only in a sub region of the image, a more appropriate point would be the center of the oriented bounding box of that ROI.

ImageRegistrationMethod

To create a specific registration instance using the ImageRegistrationMethod you need to select several components which together define the registration instance:

  1. Transformation.

  2. Similarity metric.

  3. Optimizer.

  4. Interpolator.

Transform

The type of transformation defines the mapping between the two images. SimpleITK supports a variety of global and local transformations. The available transformations include:

The parameters modified by the registration framework are those returned by the transforms GetParameters() method. This requires special attention when the using a composite transform, as the specific parameters vary based on the content of your composite transformation.

Similarity Metric

The similarity metric reflects the relationship between the intensities of the images (identity, affine, stochastic…). The available metrics include:

In the ITKv4 and consequentially in SimpleITK all similarity metrics are minimized. For metrics whose optimum corresponds to a maximum, such as mutual information, the metric value is negated internally. The selection of a similarity metric is done using the ImageRegistrationMethod’s SetMetricAsX() methods.

Optimizer

The optimizer is selected using the SetOptimizerAsX() methods. When selecting the optimizer you will also need to configure it (e.g. set the number of iterations). The available optimizers include:

Interpolator

SimpleITK has a large number of interpolators. In most cases linear interpolation, the default setting, is sufficient. Unlike the similarity metric and optimizer, the interpolator is set using the SetInterpolator method which receives a parameter indicating the interpolator type.

Features of Interest

Transforms and image spaces

While the goal of registration, as defined above, refers to a single transformation and two images, the ITKv4 registration and the SimpleITK ImageRegistrationMethod provide additional flexibility in registration configuration.

From a coordinate system standpoint ITKv4 introduced the virtual image domain, making registration a symmetric process so that both images are treated similarly. As a consequence the ImageRegistrationMethod has methods for setting three transformations:

1. SetInitialTransform \(T_o\) - composed with the moving initial transform, maps points from the virtual image domain to the moving image domain, modified during optimization.

2. SetFixedInitialTransform \(T_f\) - maps points from the virtual image domain to the fixed image domain, never modified.

3. SetMovingInitialTransform \(T_m\)- maps points from the virtual image domain to the moving image domain, never modified.

The transformation that maps points from the fixed to moving image domains is thus:

\[p_{moving}=T_o(T_m(T_f^{-1}(p_{fixed})))\]

Multi Resolution Framework

The ImageRegistrationMethod supports multi-resolution, pyramid, registration via two methods SetShrinkFactorsPerLevel and SetSmoothingSigmasPerLevel. The former receives the shrink factors to apply when moving from one level of the pyramid to the next and the later receives the sigmas to use for smoothing when moving from level to level. Sigmas can be specified either in voxel units or physical units (default) using SetSmoothingSigmasAreSpecifiedInPhysicalUnits.

Sampling

For many registration tasks one can use a fraction of the image voxels to estimate the similarity measure. Aggressive sampling can significantly reduce the registration runtime. The ImageRegistration method allows you to specify how/if to sample the voxels, SetMetricSamplingStrategy, and if using a sampling, what percentage, SetMetricSamplingPercentage.

The registration framework supports three sampling strategies:

  1. NONE - use all voxels, sampled points are the voxel centers.

  2. REGULAR - sample every n-th voxel while traversing the image in scan-line order, then within each voxel randomly perturb from center.

  3. RANDOM - sample image voxels with replacement using a uniform distribution, then within each voxel randomly perturb from center.

When using the REGULAR or RANDOM sampling strategies, running the same registration code multiple times will yield different results. To remove the randomness, set the pseudo-random number generator’s seed in the SetMetricSamplingPercentage method to a constant. The default seed value is wall clock time.

Note that using the RANDOM sampling strategy with a 100% sampling rate is not equivalent to using the sampling strategy of NONE. Given an image with N voxels, the former randomly selects N voxels with repetition and perturbs the points within each voxel, the latter uses the centers of all N voxels. Thus, for repeated random sampling with 100% rate, different samples are produced and likely none of them is of the centers of all N voxels.

Combining a mask with sampling is done using a rejection approach. First a sample is generated and then it is accepted or rejected if it is inside or outside the mask. This may cause problems when the mask region occupies a very small region in the original image. Because the sampling only discards data,the sample rate may be reduced from the requested one. For some similarity metrics (e.g. mutual information) this can result in an insufficient number of samples for metric value computation, leading to registraiton failure. Other metrics are more robust to small sample sizes (e.g. mean squares), but they all suffer from it. In such cases it is better to use a cropped version of the image for registration, possibly the mask’s bounding box, instead of the original image with a mask.

Scaling in Parameter Space

The ITKv4 framework introduced automated methods for estimating scaling factors for non-commensurate parameter units. These change the step size per parameter so that the effect of a unit of change has similar effects in physical space (think rotation of 1 radian and translation of 1 millimeter). The relevant methods are SetOptimizerScalesFromPhysicalShift, SetOptimizerScalesFromIndexShift and SetOptimizerScalesFromJacobian. In many cases this scaling is what determines if the the optimization converges to the correct optimum.

Observing Registration Progress

The ImageRegistrationMethod enables you to observe the registration process as it progresses. This is done using the Command-Observer pattern, associating callbacks with specific events. To associate a callback with a specific event use the AddCommand method.

Reproducibility

Generally speaking, repeated registrations of the same datasets will yield slightly different results. This is associated with the registration framework’s similarity metric computation implementation which utilizes both randomization and multi-threading.

The primary source of variability is the use of randomization. To eliminate randomization variability, set the seed parameter for the SetMetricSamplingPercentage to a fixed value (e.g. 42).

The secondary, and much more minor, source of variability has to do with the multithreading implementation, different order of operations in each registration run. To eliminate multithreading variability, set the number of threads to one via the SetGlobalDefaultNumberOfThreads method. From a practical standpoint, in the tradeoff between full reproducibility and computational efficiency with minor variability in results, computational efficiency is most often more important.