Unstable iteration for registration

Algorithms
1

Hi all,

I am using itk.RegularStepGradientDescentOptimizerv4 to perform a multiresolution translation registration between two images. I ran my code in the same computation environment for the same two images at different times, but I got a slightly different results (please see the attached figures) during the iterations. I checked the results and found the difference mostly starting from the middle level of a image pyramid (I used 3-level image pyramid), and then the difference would be increased and increased, something like “butterfly effect”. Actually, I have applied “ReturnBestParametersAndValueOn” for itk.RegularStepGradientDescentOptimizerv4. So I wonder is the phenomenon normal? If it is normal, how can I avoid this kind of unstable iteration as much as possible? Thanks in advance.

Regards,

Zhuangming Shen

Here is some parts of my code.

translation_registration = itk.ImageRegistrationMethodv4[itk.Image.F3, itk.Image.F3].New()

translation_interpolator_fixed_image = itk.LinearInterpolateImageFunction[itk.Image.F3, itk.D].New()
translation_interpolator_moving_image = itk.LinearInterpolateImageFunction[itk.Image.F3, itk.D].New()

translation_metric = itk.MattesMutualInformationImageToImageMetricv4[itk.Image.F3, itk.Image.F3].New()
translation_metric.SetFixedInterpolator(translation_interpolator_fixed_image)
translation_metric.SetMovingInterpolator(translation_interpolator_moving_image)
translation_metric.SetFixedImageMask(fixed_object)
translation_metric.SetNumberOfHistogramBins(32)

translation_transform = itk.TranslationTransform[itk.D, 3].New()

translation_optimizer_scales = itk.OptimizerParametersitk.D
translation_optimizer_scales.SetSize(translation_transform.GetNumberOfParameters())
translation_optimizer_scales.SetElement(0, 1.0) # x translation
translation_optimizer_scales.SetElement(1, 1.0) # y translation
translation_optimizer_scales.SetElement(2, 1.0) # z translation

translation_optimizer = itk.RegularStepGradientDescentOptimizerv4[itk.D].New()
translation_optimizer.SetScales(translation_optimizer_scales)
translation_optimizer.SetNumberOfIterations(100)
translation_optimizer.SetLearningRate(6.0)
translation_optimizer.SetMinimumStepLength(1e-1)
translation_optimizer.SetDoEstimateLearningRateAtEachIteration(True)
translation_optimizer.ReturnBestParametersAndValueOn()

translation_registration.SetFixedImage(caster_fixed_image.GetOutput())
translation_registration.SetMovingImage(caster_moving_image.GetOutput())
translation_registration.SetMetric(translation_metric)
translation_registration.SetOptimizer(translation_optimizer)
translation_registration.SetInitialTransform(translation_transform)
translation_registration.SetShrinkFactorsPerLevel([4, 2, 1])
translation_registration.SetMetricSamplingPercentagePerLevel([1,1,1])

def iterationUpdate():
current_parameters = translation_registration.GetOutput().Get().GetParameters()
print "[%d] optimizer value: %f current step length: %f current parameters: %f %f %f " %(translation_optimizer.GetCurrentIteration(), translation_optimizer.GetValue(), translation_optimizer.GetCurrentStepLength(), current_parameters.GetElement(0), current_parameters.GetElement(1), current_parameters.GetElement(2))

iteration_command = itk.PyCommand.New()
iteration_command.SetCommandCallable(iterationUpdate)
translation_optimizer.AddObserver(itk.IterationEvent(), iteration_command)

translation_registration.Update()

Fig1
Fig1.PNG718×403 75.3 KB

Fig2
Fig2.PNG721×386 70 KB

2

Have you tried translation_metric.ReinitializeSeed(12345);?

The documentation.

3

Hi dzenanz,

I added translation_registration.MetricSamplingReinitializeSeed(12345) to the code but still got non-deterministic results. It seems still existing random sampling in other places of the code. Any suggestion.

Regards,

Zhuangming Shen

Fig3
Fig3.PNG719×362 65.5 KB

Fig4
Fig4.PNG717×367 69.9 KB

4

What happens if you use only one thread? What happens if you use a different metric?

5

Hi blowekamp,

Thank you for your reply. But I’m not sure how can I configure to use one thread? I just added translation_registration.SetNumberOfThreads(1) to the code, but it still ran multithreads according to Linux “Top”.

The ITK 4.12.1 I used now only supports two metric for python-based image registration v4, i.e. MattesMutualInformationImageToImageMetricv4 and MeanSquaresImageToImageMetricv4. Because I need to align multimodal images, I have to use MattesMutualInformationImageToImageMetricv4.

Regards,

Zhuangming Shen

6

You tried option 3 from here. You can try options 1 or 2.

7

Hi dzenanz,

I tried option 1 and the code used only one thread. Now it shows deterministic results. So it seems the multi-threads cause the unstable registration. Can ITK guys fix the issue otherwise using one thread is too slow for practical medical image registration application?

Regards,

Zhuangming Shen

8

Having the metric algorithms be deterministic while being parallel would somewhat hurt their performance. I suppose that the decision was made to have the high performance for the parallel case, and determinism for single-threaded case. Determinism is usually required for debugging only.

@blowekamp @matt.mccormick @fbudin Is this the case, or is non-determinism with multiple threads but fixed seed a bug?

9

In the MattesMutualInformation metric, I am aware that the order of accumulation is non-deterministic with multiple threads. The magnitude of the difference seems plausible for this to be the source. I would try just setting the metric to use a single thread to verify the metric is the source of the non-determinism.

You may want to try the JointHistogramMutualInformation metric I don’t believe it has the same “optimization” as the mattes metric. It does not sound like it’s available in ITK Python, but it is available in SimpleITK along with all the other ITKv4 metrics and optimizers. There is not reason it can’t be wrapped in ITK python with the appropriate application of time and effort.

IMHO determinism or bounds of error is needed for reproducible research.

10

Hi dzenanz and blowekamp,

Thanks for your promptly response. I think Matt has wrapped JoinHistogramMutualInformation metric (http://review.source.kitware.com/#/c/22738/) for me. I will try it but I still expect you could solve the problem.

Regards,

Zhuangming Shen

11

Hi blowekamp,

You’re right. I tried JointHistogramMutualInformation metric, the results for different runs are consistent. Thanks for your suggestion.

Regards,

Zhuangming Shen

12

Hi Zhuangming,

To confirm, considerable effort and improvements were made to the v4 registration similarity matching metrics to give consistent results. Setting the random number generator seed and using the v4 metrics will help ensure consistency.

Regards,
Matt

13

Hi Matt,

I believe you have done a lot of work to make v4 registration more robust and efficient than v3 registration. Hope v4 registration can be better and better.

Regards,

Zhuangming Shen