Highly variable affine registration results for repeated runs

Hi all,

I’m working on a workflow that includes image registration, with allowance for rigid, affine or bspline transforms, and initialization using either moments, geometry or landmarks.

Some background:

My workflow was developing quite nicely in SimpleElastix (Selx). In the early days I managed to get from 0 to “something workable” in Selx but I struggled in Sitk. But I eventually ran into some roadblocks with Selx (e.g. bspline registration with landmarks for initial alignment), so I decided to make a move to Sitk.

The data:

DICOM MR images (T1 and T2).

Variability test:

Perform an affine registration using center-of-mass or geometrical center 10 times. For Sitk, also vary the sampling percentage (5%, 10%, 25%, 50%, 75%, 100%) for evaluating the metric. Each “test” consists of 10 repeated runs with a fixed initialization method and sampling percentage (for Sitk).

Gauge performance by performing image addition and subtraction of a single slice in fixed image and registered image (middle of the stack, same slice index for both images), or adding the (10) registered images (pixel-by-pixel, frame-by-frame) from any given test.

Result using Selx:

Registration takes 20-30 s per run using either COM or geometry. COM and geometry results look different, but results look identical visually when looking at any given initialization method (i.e. COM or geometry).

Following is a plot of the sum of all registered images using a geometrical centering for initialization:

Result using Sitk:

Registration takes 25-80 s per run depending on sampling percentage, using either COM or geometry. Results are highly variable when looking at any given initialization method with any given sampling percentage.

The final metric values are similar from run to run, for different sampling and different initializations, with a similar number of iterations to reach convergence. The registration always terminates due to convergence (as opposed to hitting a limit on the maximum number of iterations).

Following is a plot of the sum of all registered images using a geometrical centering for initialization, and with 5% sampling percentage:

and with 50% sampling percentage:

I estimate the maximum pixel shift to be approximately 3 pix (3 mm).

Sitk Registration method:

def SitkNonDefReg(FixIm, MovIm, Transform='affine', InitMethod='moments',
                  SamplingPercentage=5, NumIters=500):
    import SimpleITK as sitk
    from ipywidgets import interact, fixed
    def command_iteration(method):
        print(f"{method.GetOptimizerIteration():3} = {method.GetMetricValue():10.5f}")
        print("\t#: ", len(method.GetOptimizerPosition()))
    def command_multi_iteration(method):
        print("--------- Resolution Changing ---------")
    if InitMethod in ['geometry', 'geometricalcenter']:
        CTIF = sitk.CenteredTransformInitializerFilter.GEOMETRY
    elif InitMethod in ['moments', 'centerofgravity']:
        CTIF = sitk.CenteredTransformInitializerFilter.MOMENTS
    SamplingPercentage = SamplingPercentage/100
    LearningRate = 1.0
    #LearningRate = 5.0
    if Transform == 'rigid':
        SitkTx = sitk.Euler3DTransform()
    elif Transform == 'affine':
        SitkTx = sitk.AffineTransform(FixIm.GetDimension())
    """ Initial alignment. """   
    InitialTx = sitk.CenteredTransformInitializer(FixIm, MovIm,
                                                  SitkTx, CTIF)
    """ Registration. """
    RegMethod = sitk.ImageRegistrationMethod()
    RegMethod.SetInitialTransform(InitialTx, inPlace=False)

    """ Similarity metric settings. """
    """ Setup for the multi-resolution framework. """
    ShrinkFactors = [4,2,1]
    SmoothingSigmas = [2,1,1]
    """ Interpolator settings. """
    """ Optimizer settings. """
    # Connect all of the observers for plotting during registration:
    RegMethod.AddCommand(sitk.sitkStartEvent, StartPlot)
    RegMethod.AddCommand(sitk.sitkEndEvent, EndPlot)
    RegMethod.AddCommand(sitk.sitkMultiResolutionIterationEvent, UpdateMultiresIterations) 
    RegMethod.AddCommand(sitk.sitkIterationEvent, lambda: PlotValues(RegMethod))
    #FinalTx = RegMethod.Execute(sitk.Cast(FixIm, sitk.sitkFloat32), 
    #                            sitk.Cast(MovIm, sitk.sitkFloat32))
    FinalTx = RegMethod.Execute(FixIm, MovIm)
    """ Post registration analysis. """
    print(f'Final metric value: {RegMethod.GetMetricValue()}')
    print("Optimizer stopping condition:",
    #print(f"Final Iteration: {RegMethod.GetOptimizerIteration()}")
    RegIm = sitk.Resample(MovIm, FixIm, FinalTx, sitk.sitkLinear, 0.0, 
    return RegIm

I would very greatful if someone could check over the above method and suggest any ideas on where I’m going wrong, what could be causing the large variability and how to overcome it.

Hello @ctorti,

You are not doing anything wrong. What you are seeing are the effects of random sampling.
With the default settings every registration run uses a unique set of samples (wall clock is used as random seed). To remove this randomness you need to set the random seed in the SetMetricSamplingPercentage function to a constant. Any constant will do, but commonly you’ll see:

This should make all repeated runs behave the same.

This does not address the problem of “which of ten random runs should be taken as the correct solution” (taking the minimum amongst them is a common approach). You may find the registration initialization and panorama notebooks of interest. They both deal with registration initialization.


You can tune the number of samples, metric parameters, etc. and plot the metric value as a surface while perturbing two parameters around the optimum. You have found the good settings if the metric surface is smooth, with a single peak. If you see that the surface is noisy or the optimum is not a single peak but a ridge or there are multiple peaks then it indicates that that the metric is not good enough. If the metric surface plots look good then you can go on to choosing an optimizer and tuning its parameters (to make it move quickly, while preserving convergence, and make it quickly realize when it reached the optimum).

I can confirm your experience with SimpleElastix - its simple default rigid and rigid+bspline parameter sets work very well for all reasonable input images. In contrast, BRAINSTools (which is a thin wrapper over ITK’s registration framework) always requires hours of tuning for each registration problem. As a result, in recent years I completely switched to Elastix. I’ve tried ANTS, too, and it seem to have similarly good performance as Elastix. I’m sure that ITK’s built-in registration framework should be able to achieve similar performance as Elastix and ANTS, but somebody would need to take the time to thoroughly test and optimize it and come up with parameter sets that work well on all reasonable image inputs. Maybe you will be this person!

1 Like

Hi @zivy many thanks for your clarification and suggestions. I will re-visit those notebooks.

I was aware of the random nature of the sampling but what struck me as surprising was that increasing the percentage seemed to make no difference in terms of reducing the variability. Even if I use 100% of samples I still get large variations:

Is this to be expected?

I’m guessing that in this case the random seed changes the order of the samples fed into the metric, but all samples are utilized nonetheless.

Hi @lassoan thanks for your insight and suggestion to plot the metric surface. I also admire your optimism with regards to my abilities. :laughing: Many times I’ve gone in circular arguments when trying to decide whether to progress with SimpleElastix or SimpleITK.

Apologies for the length of my response but this is an interesting topic with a fair bit of complexity. So if you will forgive me…

Despite Kasper Marstal’s valiant efforts I think the people power behind the development of Sitk far exceeds that of Selx. And with numerous notebooks for examples and forums like this it’s possible to get support from the community. This doesn’t exist for Selx. But moreover I ran into some major roadblocks with Selx.

It might be helpful to take a moment to explain the bigger picture of what I’m trying to do: An add-on feature for the OHIF-Viewer that will enable a user to “copy” a ROI (contour(s) or segmentation(s)) from one DICOM series to another. The bulk of the code will likely remain in Python with some conversion from/to JavaScript at inputs/outputs to the Python code.

If the “source” and “target” series have different FrameOfReferenceUIDs registration will be performed. The images could be MR, CT, PET, etc. and might require rigid, affine or bspline transforms, and might require landmarks (e.g. MR to full-body CT) for initialization.

Once the user has patiently waited for a registration to be performed we don’t want them to wait again if unnecessary. So following registration the parameters will be stored as a DICOM Registration Object. If in the future the user wishes to perform a copy of a different ROI but for the same pair of series the parameters will be read from the DRO, and the “source” ROI will be transformed to the “target” domain skipping the computationally expensive registration step.

Here lies the main issues that resulted in the switch from Selx to Sitk:

  1. While affine/bspline registrations in Selx using a COM/geometric initialization took ~ 1 min, using landmarks (“CorrespondingPointsEuclideanDistanceMetric”) for initialization took 7-8 min. By comparison Sitk does the affine in about 1 min, and the bspline in 3-4 min. I tried tuning Selx parameters to try to reduce computational time but to no avail.
  2. Trying to get the deformation field following a bspline registration in Selx killed my kernel (known issue), which means I can’t get access to the parameters to store in a DRO
  3. If I perform the registration in Sitk v1.2 (which is what Selx is based on) I cannot SetCoefficientImage for BSpline Transforms (SetCoefficientImage was introduced in v2.0), which means I can’t use an existing DRO to transform/resample the “source” ROI

Perhaps if the import call to Selx was import SimpleElastix rather than import SimpleITK it might have been possible for me register using Selx and transform/resample with Sitk v2, but that would still have required me to get the deformation field. As far as I’m aware it isn’t possible/practical to import two different versions of a package at the same time.

So I’m in a difficult situation where I would benefit greatly from the ease and robustness of Selx (which is really needed if you’re developing a GUI tool to be used for a broad range of data by non-programmers) but as it stands my workflow is fragmented in Selx and complete (not completed) in Sitk.

You seem to be a very experienced developer with expertise on both sides of the Selx/Sitk fence. I wonder if you might have any advice on how you might approach this dilemma.

Let me just quickly address one point.

  • For all sampling strategies, there is a random perturbation so that the sample points are not at regular gird to reduce sample aliasing artifacts. So if you want reproducibility then set the seed parameter in the following function:
ImageRegistrationMethod::SetMetricSamplingPercentage (double percentage, unsigned int seed=sitkWallClock)
1 Like

Thanks @blowekamp.

I was able to get reproducible results by setting the seed parameter to a constant as suggested by Ziv and yourself. I completely take Ziv’s point that my new reproducible result is no “better” than any one of the variable results I previously had.

It would be interesting to know why setting the sampling to 100% didn’t eliminate the variability. Am I wrong in thinking that run #1 with 100% of samples should provide the same result as run #2 also with 100%?

And by this logic and in general, the variability should decrease as the sampling percentage increases?

Or is there subtlety in how the metric utilizes all 100% of samples, i.e. using a sub-set of the 100% samples at different iterations, but overall using all 100% (seems unlikely since knowing how many iterations will be performed is not known apriori, and hence, how to “bin” the 100% samples, or how much overlapping is required to avoid leaving out any sample prior to termination)?

Each same is perturbed by a random amount to avoid aliasing:

1 Like

Unfortunately, improving and maintaining classic image registration is not an academically rewarding topic (hard to say anything new and relevant) and it is hard to get funding. Everybody is onto deep learning now, so rather than spending time with supporting and further optimizing these very good classic methods, people now try to use deep learning for registration, essentially starting everything from scratch. Good for getting grants and easier and more interesting for researchers, but bad for users. Developers such as those that are here to help you are very busy with managing building, packaging, wrapping infrastructure, improving documentation, examples, etc. and we don’t have time to do thorough research, such as bring up ITK’s registration methods to the robustness of Elastix.

To use ITK, SimpleITK and automatic image registration tools, such as BRAINSTools, Elastix, Plastimatch, and ANTS, we build all these tools along with 3D Slicer and its virtual Python environment. Therefore, you can use all these libraries without any version conflicts, either using GUI, Python, C++, Jupyter notebook, or REST API. We take care of all the build and packaging issues on Windows, Mac, and Linux. Users just need to deal with finding good registration settings.

Taking all these into account, if you are not confident that you can make ITK’s registration more robust then I would recommend to experiment with Elastix (you know it works), ANTS (seems to have comparable performance to Elastix and still supported and developed), and maybe Plastimatch (it has nice landmark-guided registration options and it still has some developer support). Then, whichever you find that works best, build that from source (or use 3D Slicer’s build), to avoid version incompatibilities.

1 Like

Hi @ctorti,

There is a subtle but significant difference between using a uniform random sampling with 100% sampling rate and a sampling strategy of NONE.

Given an image with N voxels, the former randomly selects N voxels with repetition, the latter takes all N voxels. Thus, for a random sampling of 100% you still get different samples and likely none of them is of all N voxels.


You cannot expect to transfer registration from one image to another just by registration. It may work for special cases (such as intra-operative registration or registration between daily fractions) but when tissue is removed, structures appear (tumor, hematoma, gas etc. develops, implant is added, …), structures have sliding contact (prostate, lung, …), or move independently (such as heart valve leaflets) then you will not be able to propagate the segmentation from the previous segmentation. No matter how much input the user gives to guide the registration, it will be just impossible to achieve good enough segmentation with it.

Therefore, if you want to guarantee that your users can get acceptable segmentation results in all cases then you need to provide semi-automatic segmentation tools. However, if you have such segmentation tools then registration guidance tools would become redundant. Asking the user to provide inputs to guide the registration first, and then ask for inputs to fix the registration would be unnecessary and probably quite confusing and and time-consuming. Therefore, I would recommend to focus your attention on the interactive segmentation method, and don’t spend too much time with trying to get a very good registration (and I would not ask inputs from the user to just guide the registration, but would only ask inputs to directly guide the segmentation).

1 Like

Hi @lassoan and thanks for your valuable insight. I completely appreciate that ITK’s developers must have their hands full, which is all the more reason why newbies like me are so grateful to have their support! :grin:

Looking into what I can do in 3D Slicer and the Python interface has been on my list of things to do for some time. I’d really love to see a demo of sorts and to spend some time playing around with it. Finding time to do that is another issue in itself… :roll_eyes:

@zivy thanks for your explanation. I’d completely forgotten about the None strategy. It would be interesting to do a bit of experimenting with those different sampling strategies to see how things behave.

@lassoan you fill me with dispair. :grimacing:

By the way, I like your use of “propagating”. Maybe that’s a nicer way of putting it as opposed to “copying”.

I should be clear that I am not working on segmentation tools, but rather ROI propagating tools. Other members of my team are working on other (OHIF-Viewer) add-ons that includes the use of AI-assisted segmentation.

On a positive note, I’ve done a fair bit (not extensive) testing of my ROI-propagating workflow on (mostly MR to MR, and mostly of the brain) images, some of which include tumours. So far I’ve been reasonably pleased with the results.