elastix 4.9 performance degradation

[mstaring]

elastix 4.9 has worse performance than 4.8.

Part of the problem is fixed in this commit:
eab944a

For the combo metric the problem persists, yielding registration times of 5 minutes in elastix 4.9 while being 1.5 minutes in elastix 4.8.

Recompiling with ITK 4.8.2 instead of ITK 4.13 showed no degradation, leading to the conclusion that a change in ITK inbetween those two versions caused the performance degradation.

[mstaring]

To be clear, this only happens with the MultiMetricMultiResolutionRegistration, not with the MultiResolutionRegistration

[mstaring]

@N-Dekker found the root cause:
InsightSoftwareConsortium/ITK#350

He is working on a fix.

[N-Dekker]

FYI, the fix that I proposed for ITK is at InsightSoftwareConsortium/ITK#351 ("PERF: Remove SystemInformation data from ResourceProbe...")

[mstaring]

Root cause fixed by @N-Dekker 👍 InsightSoftwareConsortium/ITK@747f3d1

Next step would be to upgrade elastix to ITK5.

[mstaring]

explanation is that the multi-metric created a itk::TimeProbe in every iteration, which suddenly took quite a bit of overhead (about 0.1 second in every iteration).

[mstaring]

This will be fixed with the new release with ITK5 support. Closing

[N-Dekker]

FYI, ITK4.13.2 also includes the ResourceProbe fix, that caused this large performance issue. Just released: https://discourse.itk.org/t/itk-4-13-2-has-been-released/1696

[rcorredorj]

Hi @N-Dekker and @mstaring,
Just to confirm. The "stable" binaries for elastix 4.9 here available (https://github.com/SuperElastix/elastix/releases/tag/4.9.0) were not compiled with the ITK version that solves this issue ? I'm asking because I was planning to upgrade from 4.8 to 4.9 and I notice differences not only in execution time, but also in the final registration results ... Do you consider that 4.9 a really stable version or better to consider 4.8 or 5.0 ?
Thanks!
RaC

[N-Dekker]

Thanks for your questions, Ricardo @rcorredorj

The "stable" binaries for elastix 4.9 here available were not compiled with the ITK version that solves this issue ?

Correct, these elastix 4.9 binaries still have the ITK performance issue. If you would want elastix 4.9 without this performance issue, you would have to build it yourself, using a fixed version of ITK4. (>= ITK v4.13.2)

Do you consider that 4.9 a really stable version or better to consider 4.8 or 5.0 ?

I don't think so. But is there a specific reason why you would be reluctant to use elastix 5.0?

[rcorredorj]

@N-Dekker thanks for your quick reply. Some of our applications running use ITK 4.11.0. It is going to be easier to push an upgrade from ITK 4.11.0 to 4.13.x. We actually tried to compile elastix 4.9 with ITK 4.11.0. That seems to compile properly, but I tried first to run a comparison between elastix 4.8 vs elastix 4.9 using the executables available in the releases in GitHub. However, I just noticed considerable differences on the output images using the same parameters files. Not sure now whether changing to ITK 5 or 4.13.2 will provide more consistent results with elastix 4.8 (Note: many reproducibility studies run already in our team with 4.8, this is why we were expecting to keep consistent results as much as possible)

[rcorredorj]

Actually ... I just run the same registration with the binaries in the release of elastix 5.0. I get the same results that I got with elastix 4.9, both being different compared to elastix 4.8. Then the question @N-Dekker is, were there fundamental changes in the registration algorithms and dependencies between 4.8 and later versions ?

I'm not doing anything really fancy in the registration. The parameters are practically from one of the example parameters files for an Affine registration and look as follow:

// C-style comments: //

// The internal pixel type, used for internal computations
// Leave to float in general.
// NB: this is not the type of the input images! The pixel
// type of the input images is automatically read from the
// images themselves.
// This setting can be changed to "short" to save some memory
// in case of very large 3D images.
(FixedInternalImagePixelType "float")
(MovingInternalImagePixelType "float")

// The dimensions of the fixed and moving image
// NB: This has to be specified by the user. The dimension of
// the images is currently NOT read from the images.
// Also note that some other settings may have to specified
// for each dimension separately.
(FixedImageDimension 3)
(MovingImageDimension 3)

// Specify whether you want to take into account the so-called
// direction cosines of the images. Recommended: true.
// In some cases, the direction cosines of the image are corrupt,
// due to image format conversions for example. In that case, you
// may want to set this option to "false".
(UseDirectionCosines "true")

// **************** Main Components **************************

// The following components should usually be left as they are:
(Registration "MultiResolutionRegistration")
(Interpolator "BSplineInterpolator")
(ResampleInterpolator "FinalBSplineInterpolator")
(Resampler "DefaultResampler")

// These may be changed to Fixed/MovingSmoothingImagePyramid.
// See the manual.
(FixedImagePyramid "FixedSmoothingImagePyramid")
(MovingImagePyramid "MovingSmoothingImagePyramid")

// The following components are most important:
// The optimizer AdaptiveStochasticGradientDescent (ASGD) works
// quite ok in general. The Transform and Metric are important
// and need to be chosen careful for each application. See manual.
(Optimizer "AdaptiveStochasticGradientDescent")
(Transform "AffineTransform")
//(Transform "EulerTransform")
//(Transform "SimilarityTransform")

(Metric "AdvancedMattesMutualInformation")

// ***************** Transformation **************************

// Scales the affine matrix elements compared to the translations, to make
// sure they are in the same range. In general, it's best to
// use automatic scales estimation:
(AutomaticScalesEstimation "true")

// Automatically guess an initial translation by aligning the
// geometric centers of the fixed and moving.
(AutomaticTransformInitialization "true")

// Whether transforms are combined by composition or by addition.
// In generally, Compose is the best option in most cases.
// It does not influence the results very much.
(HowToCombineTransforms "Compose")

// ******************* Similarity measure *********************

// Number of grey level bins in each resolution level,
// for the mutual information. 16 or 32 usually works fine.
// You could also employ a hierarchical strategy:
//(NumberOfHistogramBins 16 32 64)
(NumberOfHistogramBins 32)

// If you use a mask, this option is important.
// If the mask serves as region of interest, set it to false.
// If the mask indicates which pixels are valid, then set it to true.
// If you do not use a mask, the option doesn't matter.
(ErodeMask "false")

// ******************** Multiresolution **********************

// The number of resolutions. 1 Is only enough if the expected
// deformations are small. 3 or 4 mostly works fine. For large
// images and large deformations, 5 or 6 may even be useful.
(NumberOfResolutions 4)

// The downsampling/blurring factors for the image pyramids.
// By default, the images are downsampled by a factor of 2
// compared to the next resolution.
// So, in 2D, with 4 resolutions, the following schedule is used:
//(ImagePyramidSchedule 8 8 4 4 2 2 1 1 )
// And in 3D:
(ImagePyramidSchedule 8 8 8 4 4 4 2 2 2 1 1 1 )

// ******************* Optimizer ****************************

// Maximum number of iterations in each resolution level:
// 200-500 works usually fine for affine registration.
// For more robustness, you may increase this to 1000-2000.
(MaximumNumberOfIterations 1000)

// The step size of the optimizer, in mm. By default the voxel size is used.
// which usually works well. In case of unusual high-resolution images
// (eg histology) it is necessary to increase this value a bit, to the size
// of the "smallest visible structure" in the image:
//(MaximumStepLength 1.0)

// **************** Image sampling **********************

// Number of spatial samples used to compute the mutual
// information (and its derivative) in each iteration.
// With an AdaptiveStochasticGradientDescent optimizer,
// in combination with the two options below, around 2000
// samples may already suffice.
(NumberOfSpatialSamples 2048)

// Refresh these spatial samples in every iteration, and select
// them randomly. See the manual for information on other sampling
// strategies.
(NewSamplesEveryIteration "true")
(ImageSampler "Random")

// ************* Interpolation and Resampling ****************

// Order of B-Spline interpolation used during registration/optimisation.
// It may improve accuracy if you set this to 3. Never use 0.
// An order of 1 gives linear interpolation. This is in most
// applications a good choice.
(BSplineInterpolationOrder 1)

// Order of B-Spline interpolation used for applying the final
// deformation.
// 3 gives good accuracy; recommended in most cases.
// 1 gives worse accuracy (linear interpolation)
// 0 gives worst accuracy, but is appropriate for binary images
// (masks, segmentations); equivalent to nearest neighbor interpolation.
(FinalBSplineInterpolationOrder 3)

//Default pixel value for pixels that come from outside the picture:
(DefaultPixelValue 0)

// Choose whether to generate the deformed moving image.
// You can save some time by setting this to false, if you are
// only interested in the final (nonrigidly) deformed moving image
// for example.
(WriteResultImage "true")
// The pixel type and format of the resulting deformed moving image
(ResultImagePixelType "int")
(ResultImageFormat "nii")

[N-Dekker]

@rcorredorj

were there fundamental changes in the registration algorithms and dependencies between 4.8 and later versions ?

Thanks for the parameter file! Honestly I wasn't involved yet between elastix 4.8 and 4.9. But there appear some changes regarding BSpline at https://github.com/SuperElastix/elastix/releases/tag/4.9.0

Certain B-spline kernels used for interpolation had a small implementation mistake

HTH, Niels

[rcorredorj]

@N-Dekker Thank you so much for the fast replies.

@mstaring , I notice you have done those changes mentioned by Niels. 43687e5

May these changes in the B-spline kernels affect the B-Spline interpolators? (we are just using AffineTransform; not a RecursiveBSplineTransform that internally uses the BSplineDerivativeKernelFunction2 affected by changes in kernels)

[N-Dekker]

@rcorredorj I just discussed your question with my colleagues Marius (@mstaring) and Stefan (@stefanklein). We concluded that there were no intended fundamental changes in the registration algorithms, between elastix 4.8 an 5.0. Differences between the registration results of those versions are always possible, but they should only be caused by little code improvements and bug fixes. We do think that elastix 5.0 is stable. That's certainly the intention!

However, you might have a look at the log files from elastix 4.8 versus 4.9 or 5.0, to get more details about those difference. If you still suspect a fundamental change, and if you think there's something to be fixed, you may of course open a new issue at https://github.com/SuperElastix/elastix/issues/new

Thanks for your interesting questions, Ricardo!

[rcorredorj]

Thank you so much @N-Dekker .

I've been comparing the two versions elastix_04_8...4.9.0, but I don't have the deep knowledge of the code to know whether the core modules were affected by some bug fixes. I just checked the bug fix related to the bsplines, but didn't seem to me to touch the interpolators or anything that it's used in the simple paramaters that I'm using. It's a bit tricky because when computing the absolute difference between the result of 4.8 and 4.9 I get an image like this one (max abs. difference = 21)

The results to me seem quite affected, not just by a few numerical/floating-point issues. Between 4.9 and 5.0 I get the same result, though.

Considering this, you suggest anyways to use the 4.9 or newer, correct ?

[stefanklein]

HI, the numerical issues (or for example changes in random samples being picked in every iteration) could lead to small (usually) subvoxel differences in the final image registration result, because it's an iterative procedure. Then subtracting the two deformed images will lead to relatively large differences, more than you might expect. The difference images would be more clear if you compute the signed difference instead of the absolute difference. That makes it easier to see the residual misalignment between the two result images. In general we recommend to use the latest Elastix version. Stefan Op 6/22/2020 om 3:08 PM schreef Ricardo A. Corredor: Thank you so much @N-Dekker<https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2FN-Dekker&data=02%7C01%7Cs.klein%40erasmusmc.nl%7Ca5cb0830fe5c487f71c608d816ad5ed7%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C637284281158023380&sdata=LkQ8JmvwYz0sJ6wn%2BWNPzw0v5Uae02Hyhgpj7noB%2F4I%3D&reserved=0> . I've been comparing the two versions elastix_04_8...4.9.0<https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2FSuperElastix%2Felastix%2Fcompare%2Felastix_04_8...4.9.0&data=02%7C01%7Cs.klein%40erasmusmc.nl%7Ca5cb0830fe5c487f71c608d816ad5ed7%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C637284281158023380&sdata=usSw2vZ7Bi%2F7Qx7AhPAbwyFi7FNuAOaaSerC%2FBweFWI%3D&reserved=0>, but I don't have the deep knowledge of the code to know whether the core modules were affected by some bug fixes. I just checked the bug fix related to the bsplines, but didn't seem to me to touch the interpolators or anything that it's used in the simple paramaters that I'm using. It's a bit tricky because when computing the absolute difference between the result of 4.8 and 4.9 I get an image like this one (max abs. difference = 21) [image]<https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fuser-images.githubusercontent.com%2F2558598%2F85290731-c701c780-b499-11ea-90a3-da36b210f998.png&data=02%7C01%7Cs.klein%40erasmusmc.nl%7Ca5cb0830fe5c487f71c608d816ad5ed7%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C637284281158033376&sdata=U5WN1g1oMC8qnAxBjUL%2B%2FSp57v5v1c0JguC3tpIZAKw%3D&reserved=0> The results to me seem quite affected, not just by a few numerical/floating-point issues. Between 4.9 and 5.0 I get the same result, though. Considering this, you suggest anyways to use the 4.9 or newer, correct ? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub<https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2FSuperElastix%2Felastix%2Fissues%2F89%23issuecomment-647508120&data=02%7C01%7Cs.klein%40erasmusmc.nl%7Ca5cb0830fe5c487f71c608d816ad5ed7%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C637284281158033376&sdata=YigsypZA3OnJQfxa2bWy4PRDKkjt%2Fmo5yc6dQoIgrQg%3D&reserved=0>, or unsubscribe<https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FAAF2LNNSXVWSTTUOKMD7LJLRX5JVDANCNFSM4GFT6WTA&data=02%7C01%7Cs.klein%40erasmusmc.nl%7Ca5cb0830fe5c487f71c608d816ad5ed7%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C637284281158033376&sdata=S8JFU26Wya5HTtjWcXoBUc1OMHaByGGI0a9NefGh0B4%3D&reserved=0>.

…

-- Stefan Klein +31 10 7043442 http://www.bigr.nl/people/StefanKlein

[rcorredorj]

Hi,

Thanks @stefanklein. The substraction between the two images gives me a resulting image with a min and max of -21 and 21 respectively:

I'm just a bit surprised. I have 258949 voxels with differences < -2 (in intensity values) and 259618 voxels with differences > 2. Seems quite a lot for being numerical issues due to subvoxel operations (considering that are the same images and same parameters). The tricky part is that we are now trying to understand whether we can confirm that results with 4.9 are better than the ones with 4.8 (computing the exact metric is not that conclusive).

In any case, I understand from your comments that you suggest to upgrade to more recent versions. I guess one of the bug fixes between 4.8 and 4.9, had an impact on the metrics computed at each iteration inducing different results between the two versions. Maybe would be good to have such information somewhere in the documentation (in the FAQ? something related to compatibility across versions maybe?).

Thanks once again for your quick answers and support! All your replies have been very helpful. Thank you all