ENH: Make random numbers multithreading registrations deterministic

Stopped using the global instance of MersenneTwisterRandomVariateGenerator. Stopped calling its `GetInstance()`.

Added a deterministic default-constructed `m_RandomVariateGenerator` to ElastixBase.

Added pointers to this generator to AdvancedImageToImageMetric, ImageRandomSamplerBase, and CMAEvolutionStrategyOptimizer. Added `SetRandomVariateGenerator` member functions to these classes (using a default-constructed generator when `SetRandomVariateGenerator` is not yet called).

Also removed the one `MersenneTwisterRandomVariateGenerator::New()` call from elastix, which was for a local generator in `ImageRandomSamplerBase::GenerateRandomNumberList()`. Instead, added `m_Seed` to ImageRandomSamplerBase, and used that seed for a default-constructed local generator.

Aims to make the results of running multiple registrations parallel (multi-threaded) within a single process deterministic.

Triggered by pull request InsightSoftwareConsortium/ITK#5287 "Deterministic multithreading usage of itkMersenneTwisterRandomVariateGenerator.cxx", Michal Meszaros, Mar 21, 2025.

Author: N-Dekker

Common/CostFunctions/itkAdvancedImageToImageMetric.h

@@ -35,6 +35,8 @@
 #include "itkAdvancedBSplineDeformableTransform.h"
 #include "itkAdvancedCombinationTransform.h"
 
+#include "elxDefaultConstruct.h"
+
 #include <cassert>
 #include <memory> // For unique_ptr.
 #include <typeinfo>
@@ -317,6 +319,12 @@ class ITK_TEMPLATE_EXPORT AdvancedImageToImageMetric : public ImageToImageMetric
   virtual void
   BeforeThreadedGetValueAndDerivative(const TransformParametersType & parameters) const;
 
+  void
+  SetRandomVariateGenerator(Statistics::MersenneTwisterRandomVariateGenerator & randomVariateGenerator)
+  {
+    m_RandomVariateGenerator = &randomVariateGenerator;
+  }
+
 protected:
   /** Constructor. */
   AdvancedImageToImageMetric();
@@ -598,6 +606,19 @@ class ITK_TEMPLATE_EXPORT AdvancedImageToImageMetric : public ImageToImageMetric
     itkExceptionMacro("Intentionally left unimplemented!");
   }
 
+  Statistics::MersenneTwisterRandomVariateGenerator &
+  GetRandomVariateGenerator()
+  {
+    return *m_RandomVariateGenerator;
+  }
+
+  // Note: Bypasses logical const-correctness
+  Statistics::MersenneTwisterRandomVariateGenerator &
+  GetMutableRandomVariateGenerator() const
+  {
+    return *m_RandomVariateGenerator;
+  }
+
   // Protected using-declaration, to avoid `-Woverloaded-virtual` warnings from GCC (GCC 11.4) or clang (macos-12).
   using Superclass::SetTransform;
 
@@ -627,6 +648,10 @@ class ITK_TEMPLATE_EXPORT AdvancedImageToImageMetric : public ImageToImageMetric
 
   MovingImageDerivativeScalesType m_MovingImageDerivativeScales{ MovingImageDerivativeScalesType::Filled(1.0) };
 
+  mutable elastix::DefaultConstruct<Statistics::MersenneTwisterRandomVariateGenerator>
+                                                      m_DefaultRandomVariateGenerator{};
+  Statistics::MersenneTwisterRandomVariateGenerator * m_RandomVariateGenerator{ &m_DefaultRandomVariateGenerator };
+
   // Private using-declarations, to avoid `-Woverloaded-virtual` warnings from GCC (GCC 11.4) or clang (macos-12).
   using Superclass::TransformPoint;
 

Common/GTesting/itkImageRandomCoordinateSamplerGTest.cxx

@@ -43,14 +43,14 @@ GTEST_TEST(ImageRandomCoordinateSampler, CheckImageValuesOfSamples)
   using ImageType = itk::Image<PixelType>;
   using SamplerType = itk::ImageRandomCoordinateSampler<ImageType>;
 
-  // Use a fixed seed, in order to have a reproducible sampler output.
-  DerefSmartPointer(MersenneTwisterRandomVariateGenerator::GetInstance()).SetSeed(1);
-
   const auto image =
     CreateImageFilledWithSequenceOfNaturalNumbers<PixelType>(ImageType::SizeType::Filled(minimumImageSizeValue));
 
+  elx::DefaultConstruct<MersenneTwisterRandomVariateGenerator> randomVariateGenerator{};
+  randomVariateGenerator.SetSeed(1);
   elx::DefaultConstruct<SamplerType> sampler{};
 
+  sampler.SetRandomVariateGenerator(randomVariateGenerator);
   const size_t numberOfSamples{ 3 };
   sampler.SetNumberOfSamples(numberOfSamples);
   sampler.SetInput(image);
@@ -84,9 +84,11 @@ GTEST_TEST(ImageRandomCoordinateSampler, SetSeedMakesRandomizationDeterministic)
   for (const SamplerType::SeedIntegerType seed : { 0, 1 })
   {
     const auto generateSamples = [seed, image] {
+      elx::DefaultConstruct<MersenneTwisterRandomVariateGenerator> randomVariateGenerator{};
+      randomVariateGenerator.SetSeed(seed);
       elx::DefaultConstruct<SamplerType> sampler{};
 
-      DerefSmartPointer(MersenneTwisterRandomVariateGenerator::GetInstance()).SetSeed(seed);
+      sampler.SetRandomVariateGenerator(randomVariateGenerator);
       sampler.SetInput(image);
       sampler.Update();
       return std::move(Deref(sampler.GetOutput()).CastToSTLContainer());
@@ -114,9 +116,10 @@ GTEST_TEST(ImageRandomCoordinateSampler, HasSameOutputWhenUsingMultiThread)
     CreateImageFilledWithSequenceOfNaturalNumbers<PixelType>(ImageType::SizeType::Filled(minimumImageSizeValue));
 
   const auto generateSamples = [image](const bool useMultiThread) {
-    DerefSmartPointer(MersenneTwisterRandomVariateGenerator::GetInstance()).SetSeed(1);
+    elx::DefaultConstruct<MersenneTwisterRandomVariateGenerator> randomVariateGenerator{};
 
     elx::DefaultConstruct<SamplerType> sampler{};
+    sampler.SetRandomVariateGenerator(randomVariateGenerator);
     sampler.SetUseMultiThread(useMultiThread);
     sampler.SetInput(image);
     sampler.Update();

Common/GTesting/itkImageRandomSamplerSparseMaskGTest.cxx

@@ -48,9 +48,6 @@ GTEST_TEST(ImageRandomSamplerSparseMask, CheckImageValuesOfSamples)
   using ImageType = itk::Image<PixelType, Dimension>;
   using MaskSpatialObjectType = itk::ImageMaskSpatialObject<Dimension>;
 
-  // Use a fixed seed, in order to have a reproducible sampler output.
-  DerefSmartPointer(MersenneTwisterRandomVariateGenerator::GetInstance()).SetSeed(1);
-
   const auto imageSize = ImageType::SizeType::Filled(minimumImageSizeValue);
   const auto image = CreateImageFilledWithSequenceOfNaturalNumbers<PixelType>(imageSize);
 
@@ -61,9 +58,12 @@ GTEST_TEST(ImageRandomSamplerSparseMask, CheckImageValuesOfSamples)
   maskSpatialObject->SetImage(maskImage);
   maskSpatialObject->Update();
 
+  elx::DefaultConstruct<MersenneTwisterRandomVariateGenerator> randomVariateGenerator{};
+  randomVariateGenerator.SetSeed(1);
   elx::DefaultConstruct<itk::ImageRandomSamplerSparseMask<ImageType>> sampler{};
 
   const size_t numberOfSamples{ 3 };
+  sampler.SetRandomVariateGenerator(randomVariateGenerator);
   sampler.SetInput(image);
   sampler.SetMask(maskSpatialObject);
   sampler.SetNumberOfSamples(numberOfSamples);
@@ -104,9 +104,11 @@ GTEST_TEST(ImageRandomSamplerSparseMask, SetSeedMakesRandomizationDeterministic)
   for (const SamplerType::SeedIntegerType seed : { 0, 1 })
   {
     const auto generateSamples = [seed, image, maskSpatialObject] {
-      elx::DefaultConstruct<SamplerType> sampler{};
+      elx::DefaultConstruct<MersenneTwisterRandomVariateGenerator> randomVariateGenerator{};
+      elx::DefaultConstruct<SamplerType>                           sampler{};
 
-      DerefSmartPointer(MersenneTwisterRandomVariateGenerator::GetInstance()).SetSeed(seed);
+      randomVariateGenerator.SetSeed(seed);
+      sampler.SetRandomVariateGenerator(randomVariateGenerator);
       sampler.SetInput(image);
       sampler.SetMask(maskSpatialObject);
       sampler.Update();
@@ -144,8 +146,9 @@ GTEST_TEST(ImageRandomSamplerSparseMask, HasSameOutputWhenUsingMultiThread)
   maskSpatialObject->Update();
 
   const auto generateSamples = [image, maskSpatialObject](const bool useMultiThread) {
-    DerefSmartPointer(MersenneTwisterRandomVariateGenerator::GetInstance()).SetSeed(1);
-    elx::DefaultConstruct<SamplerType> sampler{};
+    elx::DefaultConstruct<MersenneTwisterRandomVariateGenerator> randomVariateGenerator{};
+    elx::DefaultConstruct<SamplerType>                           sampler{};
+    sampler.SetRandomVariateGenerator(randomVariateGenerator);
     sampler.SetUseMultiThread(useMultiThread);
     sampler.SetInput(image);
     sampler.SetMask(maskSpatialObject);

Common/ImageSamplers/itkImageRandomCoordinateSampler.h

@@ -131,8 +131,7 @@ class ITK_TEMPLATE_EXPORT ImageRandomCoordinateSampler : public ImageRandomSampl
     return interpolator;
   }();
 
-  RandomGeneratorPointer m_RandomGenerator{ RandomGeneratorType::GetInstance() };
-  InputImageSpacingType  m_SampleRegionSize{ itk::MakeFilled<InputImageSpacingType>(1.0) };
+  InputImageSpacingType m_SampleRegionSize{ itk::MakeFilled<InputImageSpacingType>(1.0) };
 
   /** Generate the two corners of a sampling region, given the two corners
    * of an image. If UseRandomSampleRegion=false, the smallesPoint and largestPoint

Common/ImageSamplers/itkImageRandomCoordinateSampler.hxx

@@ -207,10 +207,12 @@ ImageRandomCoordinateSampler<TInputImage>::GenerateRandomCoordinate(
   const InputImageContinuousIndexType & largestContIndex,
   InputImageContinuousIndexType &       randomContIndex)
 {
+  Statistics::MersenneTwisterRandomVariateGenerator & randomVariateGenerator = Superclass::GetRandomVariateGenerator();
+
   for (unsigned int i = 0; i < InputImageDimension; ++i)
   {
     randomContIndex[i] = static_cast<InputImagePointValueType>(
-      this->m_RandomGenerator->GetUniformVariate(smallestContIndex[i], largestContIndex[i]));
+      randomVariateGenerator.GetUniformVariate(smallestContIndex[i], largestContIndex[i]));
   }
 } // end GenerateRandomCoordinate()
 
@@ -268,7 +270,6 @@ ImageRandomCoordinateSampler<TInputImage>::PrintSelf(std::ostream & os, Indent i
   Superclass::PrintSelf(os, indent);
 
   os << indent << "Interpolator: " << this->m_Interpolator.GetPointer() << std::endl;
-  os << indent << "RandomGenerator: " << this->m_RandomGenerator.GetPointer() << std::endl;
 
 } // end PrintSelf()
 

Common/ImageSamplers/itkImageRandomSamplerBase.h

@@ -20,6 +20,7 @@
 
 #include "itkImageSamplerBase.h"
 #include <itkMersenneTwisterRandomVariateGenerator.h>
+#include "elxDefaultConstruct.h"
 #include <optional>
 
 namespace itk
@@ -94,6 +95,13 @@ class ITK_TEMPLATE_EXPORT ImageRandomSamplerBase : public ImageSamplerBase<TInpu
     return m_OptionalSeed;
   }
 
+  void
+  SetRandomVariateGenerator(Statistics::MersenneTwisterRandomVariateGenerator & randomVariateGenerator)
+  {
+    m_RandomVariateGenerator = &randomVariateGenerator;
+  }
+
+
   /** The input image dimension. */
   itkStaticConstMacro(InputImageDimension, unsigned int, Superclass::InputImageDimension);
 
@@ -108,6 +116,12 @@ class ITK_TEMPLATE_EXPORT ImageRandomSamplerBase : public ImageSamplerBase<TInpu
   void
   GenerateRandomNumberList();
 
+  Statistics::MersenneTwisterRandomVariateGenerator &
+  GetRandomVariateGenerator()
+  {
+    return *m_RandomVariateGenerator;
+  }
+
   /** PrintSelf. */
   void
   PrintSelf(std::ostream & os, Indent indent) const override;
@@ -117,6 +131,10 @@ class ITK_TEMPLATE_EXPORT ImageRandomSamplerBase : public ImageSamplerBase<TInpu
 
 private:
   std::optional<SeedIntegerType> m_OptionalSeed{};
+  SeedIntegerType                m_Seed{ 121212 + 1 };
+
+  elastix::DefaultConstruct<Statistics::MersenneTwisterRandomVariateGenerator> m_DefaultRandomVariateGenerator{};
+  Statistics::MersenneTwisterRandomVariateGenerator * m_RandomVariateGenerator{ &m_DefaultRandomVariateGenerator };
 };
 
 } // end namespace itk

Common/ImageSamplers/itkImageRandomSamplerBase.hxx

@@ -45,24 +45,19 @@ template <typename TInputImage>
 void
 ImageRandomSamplerBase<TInputImage>::GenerateRandomNumberList()
 {
-  /** Create a random number generator. Also used in the ImageRandomConstIteratorWithIndex. */
-  const auto localGenerator = Statistics::MersenneTwisterRandomVariateGenerator::New();
-
-  if (m_OptionalSeed)
-  {
-    localGenerator->SetSeed(*m_OptionalSeed);
-  }
+  elastix::DefaultConstruct<Statistics::MersenneTwisterRandomVariateGenerator> randomVariateGenerator{};
+  randomVariateGenerator.SetSeed(m_OptionalSeed.value_or(++m_Seed));
 
   /** Clear the random number list. */
   this->m_RandomNumberList.clear();
   this->m_RandomNumberList.reserve(this->m_NumberOfSamples);
 
   /** Fill the list with random numbers. */
   const auto numPixels = static_cast<double>(this->GetCroppedInputImageRegion().GetNumberOfPixels());
-  localGenerator->GetVariateWithOpenRange(numPixels - 0.5); // dummy jump
+  randomVariateGenerator.GetVariateWithOpenRange(numPixels - 0.5); // dummy jump
   for (unsigned long i = 0; i < this->m_NumberOfSamples; ++i)
   {
-    const double randomPosition = localGenerator->GetVariateWithOpenRange(numPixels - 0.5);
+    const double randomPosition = randomVariateGenerator.GetVariateWithOpenRange(numPixels - 0.5);
     this->m_RandomNumberList.push_back(randomPosition);
   }
 }

Common/ImageSamplers/itkImageRandomSamplerSparseMask.h

@@ -94,7 +94,6 @@ class ITK_TEMPLATE_EXPORT ImageRandomSamplerSparseMask : public ImageRandomSampl
   void
   GenerateData() override;
 
-  RandomGeneratorPointer     m_RandomGenerator{ RandomGeneratorType::GetInstance() };
   InternalFullSamplerPointer m_InternalFullSampler{ InternalFullSamplerType::New() };
 
 private:

Common/ImageSamplers/itkImageRandomSamplerSparseMask.hxx

@@ -75,6 +75,7 @@ ImageRandomSamplerSparseMask<TInputImage>::GenerateData()
   const ImageSampleContainerType & allValidSamples = Deref(this->m_InternalFullSampler->GetOutput());
   unsigned long                    numberOfValidSamples = allValidSamples.Size();
 
+  Statistics::MersenneTwisterRandomVariateGenerator & randomVariateGenerator = Superclass::GetRandomVariateGenerator();
 
   /** If desired we exercise a multi-threaded version. */
   if (Superclass::m_UseMultiThread)
@@ -84,7 +85,7 @@ ImageRandomSamplerSparseMask<TInputImage>::GenerateData()
 
     for (unsigned int i = 0; i < Superclass::m_NumberOfSamples; ++i)
     {
-      m_RandomIndices.push_back(m_RandomGenerator->GetIntegerVariate(numberOfValidSamples - 1));
+      m_RandomIndices.push_back(randomVariateGenerator.GetIntegerVariate(numberOfValidSamples - 1));
     }
 
     auto & samples = sampleContainer.CastToSTLContainer();
@@ -101,7 +102,7 @@ ImageRandomSamplerSparseMask<TInputImage>::GenerateData()
 
   for (unsigned int i = 0; i < Superclass::m_NumberOfSamples; ++i)
   {
-    unsigned long randomIndex = this->m_RandomGenerator->GetIntegerVariate(numberOfValidSamples - 1);
+    unsigned long randomIndex = randomVariateGenerator.GetIntegerVariate(numberOfValidSamples - 1);
     sampleVector.push_back(allValidSamples.ElementAt(randomIndex));
   }
 
@@ -157,7 +158,6 @@ ImageRandomSamplerSparseMask<TInputImage>::PrintSelf(std::ostream & os, Indent i
   Superclass::PrintSelf(os, indent);
 
   os << indent << "InternalFullSampler: " << this->m_InternalFullSampler.GetPointer() << std::endl;
-  os << indent << "RandomGenerator: " << this->m_RandomGenerator.GetPointer() << std::endl;
 
 } // end PrintSelf()
 

Common/ImageSamplers/itkMultiInputImageRandomCoordinateSampler.h

@@ -133,8 +133,7 @@ class ITK_TEMPLATE_EXPORT MultiInputImageRandomCoordinateSampler : public ImageR
     return interpolator;
   }();
 
-  RandomGeneratorPointer m_RandomGenerator{ RandomGeneratorType::GetInstance() };
-  InputImageSpacingType  m_SampleRegionSize{ itk::MakeFilled<InputImageSpacingType>(1.0) };
+  InputImageSpacingType m_SampleRegionSize{ itk::MakeFilled<InputImageSpacingType>(1.0) };
 
   /** Generate the two corners of a sampling region. */
   virtual void