From 691c35869c01d63cb88f716dc9312cc485ead83b Mon Sep 17 00:00:00 2001 From: Niels Dekker Date: Wed, 24 May 2023 13:37:35 +0200 Subject: [PATCH] BUG: Add InitialTransformParameterObject maps to registration result Added the maps from the InitialTransformParameterObject of an ElastixRegistrationMethod to the output TransformParameterObject (at the begin, before the maps added during registration). Ensures that the initial transform is included, when the result from `registration.GetTransformParameterObject()` is passed to `transformix.SetTransformParameterObject`. Discussed with Marius (mstaring). --- .../itkElastixRegistrationMethodGTest.cxx | 96 ++++++++++++++++++- Core/Main/itkElastixRegistrationMethod.hxx | 6 +- 2 files changed, 96 insertions(+), 6 deletions(-) diff --git a/Core/Main/GTesting/itkElastixRegistrationMethodGTest.cxx b/Core/Main/GTesting/itkElastixRegistrationMethodGTest.cxx index 25e85a6a9..911bdb740 100644 --- a/Core/Main/GTesting/itkElastixRegistrationMethodGTest.cxx +++ b/Core/Main/GTesting/itkElastixRegistrationMethodGTest.cxx @@ -347,6 +347,85 @@ Expect_equal_output_SetInitialTransformParameterObject_and_Transformix_SetTransf } +// Expects transformix to yield the same output image as an elastix registration, when transformix has the +// "TransformParameterObject" taken from the elastix registration result. +template +void +Expect_equal_output_Transformix_SetTransformParameterObject_GetTransformParameterObject( + const ParameterMapVectorType & initialTransformParameterMaps, + const ImageDomain & fixedImageDomain, + const ImageDomain & movingImageDomain) +{ + using PixelType = float; + using ImageType = itk::Image; + + std::mt19937 randomNumberEngine{}; + + const auto fillImageBufferRandomly = [&randomNumberEngine](ImageType & image) { + const itk::ImageBufferRange imageBufferRange(image); + + std::generate(imageBufferRange.begin(), imageBufferRange.end(), [&randomNumberEngine] { + return std::uniform_real_distribution{ PixelType{ 1 }, PixelType{ 2 } }(randomNumberEngine); + }); + }; + + itk::Size movingImageSize; + std::iota(movingImageSize.begin(), movingImageSize.end(), 5U); + + elx::DefaultConstruct fixedImage{}; + fixedImageDomain.ToImage(fixedImage); + fixedImage.Allocate(true); + fillImageBufferRandomly(fixedImage); + + elx::DefaultConstruct movingImage{}; + movingImageDomain.ToImage(movingImage); + movingImage.Allocate(true); + fillImageBufferRandomly(movingImage); + + elx::DefaultConstruct registrationParameterObject{}; + + const ParameterMapType registrationParameterMap = CreateParameterMap({ + // Non-default parameters in alphabetic order: + NonDefaultRegistrationParameter({ "ImageSampler", { "Full" } }), // required + NonDefaultRegistrationParameter({ "MaximumNumberOfIterations", { "2" } }), + NonDefaultRegistrationParameter({ "Metric", { "AdvancedNormalizedCorrelation" } }), // default "" + NonDefaultRegistrationParameter({ "NumberOfResolutions", { "1" } }), + NonDefaultRegistrationParameter({ "Optimizer", { "AdaptiveStochasticGradientDescent" } }), // default "" + // RequiredRatioOfValidSamples as in the example in the elxMetricBase.h documentation. The FAQ even suggests 0.05: + // https://github.com/SuperElastix/elastix/wiki/FAQ/702a35cf0f5e0cf797b531fcbe3297ff9a9f3a18#i-am-getting-the-error-message-too-many-samples-map-outside-moving-image-buffer-what-does-that-mean + NonDefaultRegistrationParameter({ "RequiredRatioOfValidSamples", { "0.1" } }), + NonDefaultRegistrationParameter({ "Transform", { "BSplineTransform" } }), // default "" + }); + + registrationParameterObject.SetParameterMap(registrationParameterMap); + + elx::DefaultConstruct initialTransformParameterObject{}; + initialTransformParameterObject.SetParameterMaps(initialTransformParameterMaps); + + elx::DefaultConstruct> registration{}; + registration.SetParameterObject(®istrationParameterObject); + registration.SetInitialTransformParameterObject(&initialTransformParameterObject); + registration.SetFixedImage(&fixedImage); + registration.SetMovingImage(&movingImage); + registration.Update(); + + elx::DefaultConstruct> transformix{}; + transformix.SetTransformParameterObject(registration.GetTransformParameterObject()); + transformix.SetMovingImage(&movingImage); + transformix.Update(); + + const auto & transformixOutput = DerefRawPointer(transformix.GetOutput()); + + // Sanity checks, checking that our test is non-trivial. + EXPECT_NE(movingImage, fixedImage); + EXPECT_NE(transformixOutput, fixedImage); + EXPECT_NE(transformixOutput, movingImage); + + const auto & actualRegistrationOutput = DerefRawPointer(registration.GetOutput()); + EXPECT_EQ(actualRegistrationOutput, transformixOutput); +} + + template void Test_WriteBSplineTransformToItkFileFormat(const std::string & rootOutputDirectoryPath) @@ -951,6 +1030,7 @@ GTEST_TEST(itkElastixRegistrationMethod, SetInitialTransformParameterObject) { "Transform", { "TranslationTransform" } }, { "TransformParameters", { "0", "-2" } } } } }) { + const auto numberOfInitialTransformParameterMaps = initialTransformParameterMaps.size(); initialTransformParameterObject.SetParameterMaps(initialTransformParameterMaps); // Do the test for a few possible translations. @@ -966,10 +1046,16 @@ GTEST_TEST(itkElastixRegistrationMethod, SetInitialTransformParameterObject) const auto & transformParameterMaps = DerefRawPointer(registration.GetTransformParameterObject()).GetParameterMaps(); - ASSERT_EQ(transformParameterMaps.size(), numberOfRegistrationParameterMaps); + ASSERT_EQ(transformParameterMaps.size(), + numberOfInitialTransformParameterMaps + numberOfRegistrationParameterMaps); + + for (std::size_t i{}; i < numberOfInitialTransformParameterMaps; ++i) + { + EXPECT_EQ(transformParameterMaps[i], initialTransformParameterMaps[i]); + } // All registration parameter maps, except for the first one, should just have a zero-translation. - for (unsigned int i{ 1 }; i < numberOfRegistrationParameterMaps; ++i) + for (auto i = numberOfInitialTransformParameterMaps + 1; i < numberOfRegistrationParameterMaps; ++i) { const auto transformParameters = ConvertStringsToVectorOfDouble(transformParameterMaps[i].at("TransformParameters")); @@ -977,8 +1063,8 @@ GTEST_TEST(itkElastixRegistrationMethod, SetInitialTransformParameterObject) } // Together the initial translation and the first registration should have the actual image translation. - const auto transformParameters = - ConvertStringsToVectorOfDouble(transformParameterMaps.front().at("TransformParameters")); + const auto transformParameters = ConvertStringsToVectorOfDouble( + transformParameterMaps[numberOfInitialTransformParameterMaps].at("TransformParameters")); EXPECT_EQ(initialTranslation + ConvertToOffset(transformParameters), actualTranslation); } } @@ -1116,6 +1202,8 @@ GTEST_TEST(itkElastixRegistrationMethod, SetInitialTransformParameterObjectVersu { translationTransformParameterMap, bsplineTransformParameterMap }, createRandomImageDomain(), createRandomImageDomain()); + Expect_equal_output_Transformix_SetTransformParameterObject_GetTransformParameterObject( + { translationTransformParameterMap }, createRandomImageDomain(), createRandomImageDomain()); } } diff --git a/Core/Main/itkElastixRegistrationMethod.hxx b/Core/Main/itkElastixRegistrationMethod.hxx index 5fea5a60b..d7613f2d5 100644 --- a/Core/Main/itkElastixRegistrationMethod.hxx +++ b/Core/Main/itkElastixRegistrationMethod.hxx @@ -92,8 +92,10 @@ ElastixRegistrationMethod::GenerateData() DataObjectContainerPointer movingMaskContainer = nullptr; DataObjectContainerPointer resultImageContainer = nullptr; ElastixMainObjectPointer transform = nullptr; - ParameterMapVectorType transformParameterMapVector; - FlatDirectionCosinesType fixedImageOriginalDirection; + ParameterMapVectorType transformParameterMapVector = m_InitialTransformParameterObject + ? m_InitialTransformParameterObject->GetParameterMaps() + : ParameterMapVectorType{}; + FlatDirectionCosinesType fixedImageOriginalDirection; // Split inputs into separate containers for (const auto & inputName : this->GetInputNames())