From b41eb0747905d314bab1ab209a01ad12d3ea8922 Mon Sep 17 00:00:00 2001
From: Niels Dekker <N.Dekker@lumc.nl>
Date: Thu, 21 Mar 2024 12:51:46 +0100
Subject: [PATCH] ENH: Add elastix parameter
 "UseSamplingForGroupwiseRegistration"

Allows elastix users to specify that the image sampling is used for group-wise registration:

    (UseSamplingForGroupwiseRegistration "true")

Aims to implement a generalization to the ReducedFullSampler, pull request https://github.com/SuperElastix/elastix/pull/52 proposed by Mathias Polfliet in 2018.
---
 Common/GTesting/itkImageFullSamplerGTest.cxx  | 38 +++++++++++++++++++
 Common/ImageSamplers/itkImageSamplerBase.h    |  9 +++++
 Common/ImageSamplers/itkImageSamplerBase.hxx  |  6 +++
 .../elxImageSamplerBase.hxx                   |  4 ++
 4 files changed, 57 insertions(+)

diff --git a/Common/GTesting/itkImageFullSamplerGTest.cxx b/Common/GTesting/itkImageFullSamplerGTest.cxx
index b423897ef..1d346ebd4 100644
--- a/Common/GTesting/itkImageFullSamplerGTest.cxx
+++ b/Common/GTesting/itkImageFullSamplerGTest.cxx
@@ -185,4 +185,42 @@ GTEST_TEST(ImageFullSampler, ExactlyEqualVersusSlightlyDifferentMaskImageDomain)
   EXPECT_FALSE(samplesOnExactlyEqualImageDomains.empty());
 
   EXPECT_EQ(samplesOnExactlyEqualImageDomains, samplesOnSlightlyDifferentImageDomains);
+}
+
+
+GTEST_TEST(ImageFullSampler, UseForGroupwiseRegistration)
+{
+  using PixelType = std::uint8_t;
+  static constexpr auto Dimension = 3U;
+  using ImageType = itk::Image<PixelType, Dimension>;
+  using ImageFullSamplerType = itk::ImageFullSampler<ImageType>;
+
+  std::mt19937 randomNumberEngine{};
+  const auto   imageDomain = CreateRandomImageDomain<Dimension>(randomNumberEngine);
+  const auto   image = CreateImageFilledWithSequenceOfNaturalNumbers<PixelType>(imageDomain);
+  elx::DefaultConstruct<ImageFullSamplerType> sampler{};
+
+  sampler.SetInput(image);
+  sampler.UseForGroupwiseRegistration();
+  sampler.Update();
+
+  const auto & output = DerefRawPointer(sampler.GetOutput());
+
+  const auto reducedRegion = [imageDomain] {
+    auto size = imageDomain.size;
+    size.back() = 1;
+    return itk::ImageRegion<Dimension>{ imageDomain.index, size };
+  }();
+
+  const itk::ImageRegionRange imageRegionRange(*image, reducedRegion);
+  const std::size_t           numberOfSamples{ output.size() };
+
+  ASSERT_EQ(numberOfSamples, imageRegionRange.size());
+
+  auto imageRegionIterator = imageRegionRange.cbegin();
+  for (std::size_t i{}; i < numberOfSamples; ++i)
+  {
+    EXPECT_EQ(output[i].m_ImageValue, *imageRegionIterator);
+    ++imageRegionIterator;
+  }
 }
\ No newline at end of file
diff --git a/Common/ImageSamplers/itkImageSamplerBase.h b/Common/ImageSamplers/itkImageSamplerBase.h
index 284a141b5..55cc343c7 100644
--- a/Common/ImageSamplers/itkImageSamplerBase.h
+++ b/Common/ImageSamplers/itkImageSamplerBase.h
@@ -209,6 +209,13 @@ class ITK_TEMPLATE_EXPORT ImageSamplerBase
   /** Allows disabling the use of multi-threading, by `SetUseMultiThread(false)`. */
   itkSetMacro(UseMultiThread, bool);
 
+  /** Tells that the sampler is used for groupwise registration. */
+  void
+  UseForGroupwiseRegistration()
+  {
+    m_UseForGroupwiseRegistration = true;
+  }
+
 protected:
   /** The constructor. */
   ImageSamplerBase();
@@ -264,6 +271,8 @@ class ITK_TEMPLATE_EXPORT ImageSamplerBase
 
   InputImageRegionType m_CroppedInputImageRegion{};
   InputImageRegionType m_DummyInputImageRegion{};
+
+  bool m_UseForGroupwiseRegistration{ false };
 };
 
 } // end namespace itk
diff --git a/Common/ImageSamplers/itkImageSamplerBase.hxx b/Common/ImageSamplers/itkImageSamplerBase.hxx
index ca87803a7..2ba2734f6 100644
--- a/Common/ImageSamplers/itkImageSamplerBase.hxx
+++ b/Common/ImageSamplers/itkImageSamplerBase.hxx
@@ -301,6 +301,12 @@ ImageSamplerBase<TInputImage>::CropInputImageRegion()
    * InputImageRegion and the BoundingBoxRegion.
    */
   m_CroppedInputImageRegion = m_InputImageRegion;
+
+  if (m_UseForGroupwiseRegistration)
+  {
+    m_CroppedInputImageRegion.GetModifiableSize().back() = 1;
+  }
+
   if (!m_Mask.IsNull())
   {
     /** Get a handle to the input image. */
diff --git a/Core/ComponentBaseClasses/elxImageSamplerBase.hxx b/Core/ComponentBaseClasses/elxImageSamplerBase.hxx
index 8a07fcabb..c217b5d8c 100644
--- a/Core/ComponentBaseClasses/elxImageSamplerBase.hxx
+++ b/Core/ComponentBaseClasses/elxImageSamplerBase.hxx
@@ -36,6 +36,10 @@ ImageSamplerBase<TElastix>::BeforeRegistrationBase()
   const Configuration & configuration = Deref(Superclass::GetConfiguration());
   ITKBaseType &         sampler = GetSelf();
   sampler.SetUseMultiThread(configuration.RetrieveParameterValue(true, "UseMultiThreadingForSamplers", 0, false));
+  if (configuration.RetrieveParameterValue(false, "UseSamplingForGroupwiseRegistration", 0, false))
+  {
+    sampler.UseForGroupwiseRegistration();
+  }
 }
 
 /**