From 8111c672508d89995550e63ba54ab871e324105b Mon Sep 17 00:00:00 2001 From: Niels Dekker Date: Fri, 14 Jun 2024 11:58:45 +0200 Subject: [PATCH] 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 | 14 +++-- .../elxEulerStackTransform.hxx | 31 ++++++++-- .../itkElastixRegistrationMethodGTest.cxx | 61 +++++++++++++++++++ 3 files changed, 98 insertions(+), 8 deletions(-) diff --git a/Common/GTesting/elxTransformIOGTest.cxx b/Common/GTesting/elxTransformIOGTest.cxx index b9fdce84c..93b68996c 100644 --- a/Common/GTesting/elxTransformIOGTest.cxx +++ b/Common/GTesting/elxTransformIOGTest.cxx @@ -455,10 +455,16 @@ struct WithDimension { { "DeformationFieldFileName", { expectedDeformationFieldFileName } }, { "DeformationFieldInterpolationOrder", { expectedZero } } }); WithElastixTransform::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::Test_CreateTransformParameterMap_for_default_transform( (NDimension == 3) diff --git a/Components/Transforms/EulerStackTransform/elxEulerStackTransform.hxx b/Components/Transforms/EulerStackTransform/elxEulerStackTransform.hxx index 027d8b206..88243c479 100644 --- a/Components/Transforms/EulerStackTransform/elxEulerStackTransform.hxx +++ b/Components/Transforms/EulerStackTransform/elxEulerStackTransform.hxx @@ -124,6 +124,13 @@ EulerStackTransform::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::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::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); diff --git a/Core/Main/GTesting/itkElastixRegistrationMethodGTest.cxx b/Core/Main/GTesting/itkElastixRegistrationMethodGTest.cxx index ec32cd92c..7db8706dc 100644 --- a/Core/Main/GTesting/itkElastixRegistrationMethodGTest.cxx +++ b/Core/Main/GTesting/itkElastixRegistrationMethodGTest.cxx @@ -34,6 +34,7 @@ #include #include #include +#include #include #include #include @@ -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 @@ -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; + + const auto numberOfImages = 4U; + const auto imageSizeValue = 6U; + const itk::Size imageStackSize{ imageSizeValue, imageSizeValue, imageSizeValue, numberOfImages }; + const auto imageStack = CreateImage(imageStackSize); + + const itk::ImageBufferRange imageBufferRange(*imageStack); + + std::mt19937 randomNumberEngine{}; + + std::generate(imageBufferRange.begin(), imageBufferRange.end(), [&randomNumberEngine] { + return std::uniform_real_distribution{ 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> 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)); +}