ENH: Add support ComputeZYX parameter to EulerStackTransform (for 3D)

Supports specifying parameter `(ComputeZYX "true")`, when the transform is "EulerStackTransform" and the image stack is 4D (meaning that the images are 3D).

Tested that ComputeZYX = true yields a different result than ComputeZYX = false.

Triggered by a mail from Laurens Beljaards, subject "Elastix vraag", 10 June 2024.

Common/GTesting/elxTransformIOGTest.cxx

@@ -455,10 +455,16 @@ struct WithDimension
       { { "DeformationFieldFileName", { expectedDeformationFieldFileName } },
         { "DeformationFieldInterpolationOrder", { expectedZero } } });
     WithElastixTransform<EulerStackTransform>::Test_CreateTransformParameterMap_for_default_transform(
-      { { "CenterOfRotationPoint", ParameterValuesType(NDimension - 1, expectedZero) },
-        { "NumberOfSubTransforms", { expectedZero } },
-        { "StackOrigin", { expectedZero } },
-        { "StackSpacing", { expectedOne } } });
+      (NDimension == 4)
+        ? ParameterMapType{ { "CenterOfRotationPoint", ParameterValuesType(NDimension - 1, expectedZero) },
+                            { "ComputeZYX", { expectedFalse } },
+                            { "NumberOfSubTransforms", { expectedZero } },
+                            { "StackOrigin", { expectedZero } },
+                            { "StackSpacing", { expectedOne } } }
+        : ParameterMapType{ { "CenterOfRotationPoint", ParameterValuesType(NDimension - 1, expectedZero) },
+                            { "NumberOfSubTransforms", { expectedZero } },
+                            { "StackOrigin", { expectedZero } },
+                            { "StackSpacing", { expectedOne } } });
 
     WithElastixTransform<EulerTransformElastix>::Test_CreateTransformParameterMap_for_default_transform(
       (NDimension == 3)

Components/Transforms/EulerStackTransform/elxEulerStackTransform.hxx

@@ -124,6 +124,13 @@ EulerStackTransform<TElastix>::ReadFromFile()
 
     m_DummySubTransform->SetCenter(RDcenterOfRotationPoint);
 
+    if constexpr (ReducedSpaceDimension == 3)
+    {
+      // For 3D images, retrieve and set ComputeZYX.
+      m_DummySubTransform->SetComputeZYX(
+        configuration.RetrieveParameterValue(m_DummySubTransform->GetComputeZYX(), "ComputeZYX", 0, false));
+    }
+
     /** Set stack transform parameters. */
     m_StackTransform->SetNumberOfSubTransforms(m_NumberOfSubTransforms);
     m_StackTransform->SetStackOrigin(m_StackOrigin);
@@ -149,10 +156,19 @@ EulerStackTransform<TElastix>::CreateDerivedTransformParameterMap() const -> Par
 {
   const auto & itkTransform = *m_StackTransform;
 
-  return { { "CenterOfRotationPoint", Conversion::ToVectorOfStrings(m_DummySubTransform->GetCenter()) },
-           { "StackSpacing", { Conversion::ToString(itkTransform.GetStackSpacing()) } },
-           { "StackOrigin", { Conversion::ToString(itkTransform.GetStackOrigin()) } },
-           { "NumberOfSubTransforms", { Conversion::ToString(itkTransform.GetNumberOfSubTransforms()) } } };
+  ParameterMapType parameterMap{
+    { "CenterOfRotationPoint", Conversion::ToVectorOfStrings(m_DummySubTransform->GetCenter()) },
+    { "StackSpacing", { Conversion::ToString(itkTransform.GetStackSpacing()) } },
+    { "StackOrigin", { Conversion::ToString(itkTransform.GetStackOrigin()) } },
+    { "NumberOfSubTransforms", { Conversion::ToString(itkTransform.GetNumberOfSubTransforms()) } }
+  };
+
+  if constexpr (ReducedSpaceDimension == 3)
+  {
+    parameterMap["ComputeZYX"] = { Conversion::ToString(m_DummySubTransform->GetComputeZYX()) };
+  }
+
+  return parameterMap;
 
 } // end CreateDerivedTransformParameterMap()
 
@@ -253,6 +269,13 @@ EulerStackTransform<TElastix>::InitializeTransform()
   /** Set the translation to zero */
   m_DummySubTransform->SetTranslation(ReducedDimensionOutputVectorType());
 
+  if constexpr (ReducedSpaceDimension == 3)
+  {
+    // For 3D images, retrieve and set ComputeZYX.
+    m_DummySubTransform->SetComputeZYX(
+      configuration.RetrieveParameterValue(m_DummySubTransform->GetComputeZYX(), "ComputeZYX", 0, false));
+  }
+
   /** Set all subtransforms to a copy of the dummy Translation sub transform. */
   m_StackTransform->SetAllSubTransforms(*m_DummySubTransform);
 

Core/Main/GTesting/itkElastixRegistrationMethodGTest.cxx

@@ -34,6 +34,7 @@
 #include <itkCompositeTransform.h>
 #include <itkDisplacementFieldTransform.h>
 #include <itkEuler2DTransform.h>
+#include <itkEuler3DTransform.h>
 #include <itkImage.h>
 #include <itkIndexRange.h>
 #include <itkFileTools.h>
@@ -78,7 +79,9 @@ using elx::CoreMainGTestUtilities::ImageDomain;
 using elx::CoreMainGTestUtilities::TypeHolder;
 using elx::CoreMainGTestUtilities::minimumImageSizeValue;
 using elx::GTestUtilities::MakeMergedMap;
+
 using itk::Deref;
+using itk::Statistics::MersenneTwisterRandomVariateGenerator;
 
 
 template <typename TImage>
@@ -2375,3 +2378,61 @@ GTEST_TEST(itkElastixRegistrationMethod, SetAndGetInitialTransform)
   EXPECT_EQ(registration.GetInitialTransform(), nullptr);
   EXPECT_EQ(registration.GetExternalInitialTransform(), nullptr);
 }
+
+
+// Tests that ComputeZYX = true yields a different result than ComputeZYX = false.
+GTEST_TEST(itkElastixRegistrationMethod, EulerStackTransformComputeZYX)
+{
+  using PixelType = float;
+  static constexpr auto ImageDimension = 4U;
+  using ImageType = itk::Image<PixelType, ImageDimension>;
+
+  const auto                      numberOfImages = 4U;
+  const auto                      imageSizeValue = 6U;
+  const itk::Size<ImageDimension> imageStackSize{ imageSizeValue, imageSizeValue, imageSizeValue, numberOfImages };
+  const auto                      imageStack = CreateImage<PixelType>(imageStackSize);
+
+  const itk::ImageBufferRange<ImageType> imageBufferRange(*imageStack);
+
+  std::mt19937 randomNumberEngine{};
+
+  std::generate(imageBufferRange.begin(), imageBufferRange.end(), [&randomNumberEngine] {
+    return std::uniform_real_distribution<PixelType>{ PixelType{ 0 }, PixelType{ 2 } }(randomNumberEngine);
+  });
+
+  const auto doRegistration = [imageStack](const bool computeZYX) {
+    // Use a fixed seed, in order to have a reproducible sampler output.
+    DerefSmartPointer(MersenneTwisterRandomVariateGenerator::GetInstance()).SetSeed(1);
+    MersenneTwisterRandomVariateGenerator::ResetNextSeed();
+
+    elx::DefaultConstruct<ElastixRegistrationMethodType<ImageType>> registration{};
+    registration.SetFixedImage(imageStack);
+    registration.SetMovingImage(imageStack);
+    registration.SetParameterObject(CreateParameterObject({ // Parameters in alphabetic order:
+                                                            { "AutomaticTransformInitialization", "false" },
+                                                            { "ComputeZYX", elx::Conversion::BoolToString(computeZYX) },
+                                                            { "ImageSampler", "RandomCoordinate" },
+                                                            { "Interpolator", "ReducedDimensionBSplineInterpolator" },
+                                                            { "MaximumNumberOfIterations", "3" },
+                                                            { "Metric", "AdvancedNormalizedCorrelation" },
+                                                            { "Optimizer", "AdaptiveStochasticGradientDescent" },
+                                                            { "Transform", "EulerStackTransform" } }));
+    registration.Update();
+    return GetTransformParametersFromFilter(registration);
+  };
+
+  const auto transformParameters = doRegistration(true);
+
+  const auto expectedNumberofParametersPerImage = itk::Euler3DTransform<>::New()->GetParameters().size();
+  EXPECT_EQ(transformParameters.size(), numberOfImages * expectedNumberofParametersPerImage);
+
+  // Sanity check: Expect at least one non-zero parameter value, otherwise the test is probably trivial.
+  EXPECT_TRUE(std::any_of(
+    transformParameters.cbegin(), transformParameters.cend(), [](const double parameter) { return parameter != 0.0; }));
+
+  // Sanity check: Expect the same result twice, when running with computeZYX = true twice.
+  EXPECT_EQ(transformParameters, doRegistration(true));
+
+  // Expect that the result is different that, when running with computeZYX = false.
+  EXPECT_NE(transformParameters, doRegistration(false));
+}