From 58e0a7bdd76f69a8615cb4a8856038c0892ea530 Mon Sep 17 00:00:00 2001 From: Niels Dekker Date: Fri, 28 Apr 2023 18:28:14 +0200 Subject: [PATCH] ENH: Convert the input images to the user-specified internal pixel type Both ElastixRegistrationMethod and TransformixFilter now convert their input images to the internal pixel type, specified by elastix/transformix parameters, "FixedInternalImagePixelType" and "MovingInternalImagePixelType" (which are "float" by default). The conversion is performed internally by `itk::CastImageFilter`. It is only performed when the pixel type of an input image is actually different from the specified internal pixel type. The converted images are cached in memory, to avoid repeating the very same (potentially expensive) image conversion multiple times. GoogleTest unit tests are included. Related to issue https://github.com/SuperElastix/elastix/issues/322 "TransformixFilter does not process non-float pixel type images" --- .../Main/GTesting/elxCoreMainGTestUtilities.h | 131 ++++++++++- .../itkElastixRegistrationMethodGTest.cxx | 212 ++++++++++++++++++ .../GTesting/itkTransformixFilterGTest.cxx | 185 ++++++++++++++- Core/Main/elxForEachSupportedImageType.h | 18 ++ Core/Main/elxLibUtilities.cxx | 19 +- Core/Main/elxLibUtilities.h | 59 +++++ Core/Main/itkElastixRegistrationMethod.h | 17 ++ Core/Main/itkElastixRegistrationMethod.hxx | 64 ++++-- Core/Main/itkTransformixFilter.hxx | 31 +-- 9 files changed, 688 insertions(+), 48 deletions(-) diff --git a/Core/Main/GTesting/elxCoreMainGTestUtilities.h b/Core/Main/GTesting/elxCoreMainGTestUtilities.h index 276d98b4d..7e6bc6196 100644 --- a/Core/Main/GTesting/elxCoreMainGTestUtilities.h +++ b/Core/Main/GTesting/elxCoreMainGTestUtilities.h @@ -37,6 +37,7 @@ #include // For begin and end. #include #include // For iota. +#include #include #include // For is_pointer, is_same, and integral_constant. #include @@ -50,6 +51,13 @@ namespace elastix namespace CoreMainGTestUtilities { +/// Eases passing a type as argument to a generic lambda. +template +struct TypeHolder +{ + using Type = TNested; +}; + /// Simple exception class, to be used by unit tests. class Exception : public std::exception { @@ -243,6 +251,14 @@ ParameterObject::Pointer inline CreateParameterObject( } +ParameterObject::Pointer inline CreateParameterObject(const ParameterObject::ParameterMapType & parameterMap) +{ + const auto parameterObject = ParameterObject::New(); + parameterObject->SetParameterMap(parameterMap); + return parameterObject; +} + + inline std::vector GetTransformParametersFromMaps(const std::vector & transformParameterMaps) { @@ -275,6 +291,11 @@ GetTransformParametersFromFilter(TFilter & filter) } +// ITK's RecursiveSeparableImageFilter "requires a minimum of four pixels along the dimension to be processed", at +// https://github.com/InsightSoftwareConsortium/ITK/blob/v5.3.0/Modules/Filtering/ImageFilterBase/include/itkRecursiveSeparableImageFilter.hxx#L226 +constexpr itk::SizeValueType minimumImageSizeValue{ 4 }; + + // The image domain. ITK calls it the "geometry" of an image. ("The geometry of an image is defined by its position, // orientation, spacing, and extent", according to https://itk.org/Doxygen52/html/classitk_1_1ImageBase.html#details). // The elastix manual (elastix-5.1.0-manual.pdf, January 16, 2023) simply calls it "the @@ -304,6 +325,14 @@ struct ImageDomain : size(size) {} + explicit ImageDomain(const ImageBaseType & image) + : direction(image.GetDirection()) + , index(image.GetLargestPossibleRegion().GetIndex()) + , size(image.GetLargestPossibleRegion().GetSize()) + , spacing(image.GetSpacing()) + , origin(image.GetOrigin()) + {} + // Constructor, allowing to explicitly specify all the settings of the domain. ImageDomain(const DirectionType & direction, const IndexType & index, @@ -332,16 +361,103 @@ struct ImageDomain AsParameterMap() const { return { + // Parameters in alphabetic order: { "Direction", elx::Conversion::ToVectorOfStrings(direction) }, { "Index", elx::Conversion::ToVectorOfStrings(index) }, + { "Origin", elx::Conversion::ToVectorOfStrings(origin) }, { "Size", elx::Conversion::ToVectorOfStrings(size) }, { "Spacing", elx::Conversion::ToVectorOfStrings(spacing) }, - { "Origin", elx::Conversion::ToVectorOfStrings(origin) }, }; } + + friend bool + operator==(const ImageDomain & lhs, const ImageDomain & rhs) + { + return lhs.direction == rhs.direction && lhs.index == rhs.index && lhs.size == rhs.size && + lhs.spacing == rhs.spacing && lhs.origin == rhs.origin; + } + + friend bool + operator!=(const ImageDomain & lhs, const ImageDomain & rhs) + { + return !(lhs == rhs); + } }; +template +int +GenerateRandomSign(TRandomNumberEngine & randomNumberEngine) +{ + return (randomNumberEngine() % 2 == 0) ? -1 : 1; +} + + +template +auto +CreateRandomImageDomain(std::mt19937 & randomNumberEngine) +{ + using ImageDomainType = ImageDomain; + + const auto createRandomDirection = [&randomNumberEngine] { + using DirectionType = typename ImageDomainType::DirectionType; + auto randomDirection = DirectionType::GetIdentity(); + + // For now, just a single random rotation + const auto randomRotation = std::uniform_real_distribution<>{ -M_PI, M_PI }(randomNumberEngine); + const auto cosRandomRotation = std::cos(randomRotation); + const auto sinRandomRotation = std::sin(randomRotation); + + randomDirection[0][0] = cosRandomRotation; + randomDirection[0][1] = sinRandomRotation; + randomDirection[1][0] = -sinRandomRotation; + randomDirection[1][1] = cosRandomRotation; + + return randomDirection; + }; + const auto createRandomIndex = [&randomNumberEngine] { + typename ImageDomainType::IndexType randomIndex{}; + std::generate(randomIndex.begin(), randomIndex.end(), [&randomNumberEngine] { + return std::uniform_int_distribution{ + std::numeric_limits::min() / 2, std::numeric_limits::max() / 2 + }(randomNumberEngine); + }); + return randomIndex; + }; + const auto createRandomSmallImageSize = [&randomNumberEngine] { + typename ImageDomainType::SizeType randomImageSize{}; + std::generate(randomImageSize.begin(), randomImageSize.end(), [&randomNumberEngine] { + return std::uniform_int_distribution{ minimumImageSizeValue, + 2 * minimumImageSizeValue }(randomNumberEngine); + }); + return randomImageSize; + }; + const auto createRandomSpacing = [&randomNumberEngine] { + typename ImageDomainType::SpacingType randomSpacing{}; + std::generate(randomSpacing.begin(), randomSpacing.end(), [&randomNumberEngine] { + // Originally tried the maximum interval from std::numeric_limits::min() to + // std::numeric_limits::max(), but that caused errors during inverse matrix computation. + return std::uniform_real_distribution{ 0.1, 10.0 }(randomNumberEngine); + }); + return randomSpacing; + }; + const auto createRandomPoint = [&randomNumberEngine] { + typename ImageDomainType::PointType randomPoint{}; + std::generate(randomPoint.begin(), randomPoint.end(), [&randomNumberEngine] { + return GenerateRandomSign(randomNumberEngine) * std::uniform_real_distribution{ + itk::SpacePrecisionType{}, std::numeric_limits::max() / 2.0 + }(randomNumberEngine); + }); + return randomPoint; + }; + + return ImageDomainType{ createRandomDirection(), + createRandomIndex(), + createRandomSmallImageSize(), + createRandomSpacing(), + createRandomPoint() }; +} + // Creates a test image, filled with zero. template auto @@ -357,11 +473,11 @@ CreateImage(const itk::Size & imageSize) // Creates a test image, filled with a sequence of natural numbers, 1, 2, 3, ..., N. template auto -CreateImageFilledWithSequenceOfNaturalNumbers(const itk::Size & imageSize) +CreateImageFilledWithSequenceOfNaturalNumbers(const ImageDomain & imageDomain) { using ImageType = itk::Image; const auto image = ImageType::New(); - image->SetRegions(imageSize); + imageDomain.ToImage(*image); image->Allocate(); const itk::ImageBufferRange imageBufferRange{ *image }; std::iota(imageBufferRange.begin(), imageBufferRange.end(), TPixel{ 1 }); @@ -369,6 +485,15 @@ CreateImageFilledWithSequenceOfNaturalNumbers(const itk::Size & } +// Creates a test image, filled with a sequence of natural numbers, 1, 2, 3, ..., N. +template +auto +CreateImageFilledWithSequenceOfNaturalNumbers(const itk::Size & imageSize) +{ + return CreateImageFilledWithSequenceOfNaturalNumbers(ImageDomain{ imageSize }); +} + + std::string GetDataDirectoryPath(); diff --git a/Core/Main/GTesting/itkElastixRegistrationMethodGTest.cxx b/Core/Main/GTesting/itkElastixRegistrationMethodGTest.cxx index bcd428c07..c79f4a841 100644 --- a/Core/Main/GTesting/itkElastixRegistrationMethodGTest.cxx +++ b/Core/Main/GTesting/itkElastixRegistrationMethodGTest.cxx @@ -23,8 +23,11 @@ #include +#include "GTesting/elxGTestUtilities.h" + #include "elxCoreMainGTestUtilities.h" #include "elxDefaultConstruct.h" +#include "elxForEachSupportedImageType.h" #include "elxTransformIO.h" // ITK header file: @@ -64,6 +67,7 @@ using elx::CoreMainGTestUtilities::CreateImage; using elx::CoreMainGTestUtilities::CreateImageFilledWithSequenceOfNaturalNumbers; using elx::CoreMainGTestUtilities::CreateParameterMap; using elx::CoreMainGTestUtilities::CreateParameterObject; +using elx::CoreMainGTestUtilities::CreateRandomImageDomain; using elx::CoreMainGTestUtilities::DerefRawPointer; using elx::CoreMainGTestUtilities::DerefSmartPointer; using elx::CoreMainGTestUtilities::FillImageRegion; @@ -73,6 +77,9 @@ using elx::CoreMainGTestUtilities::GetDataDirectoryPath; using elx::CoreMainGTestUtilities::GetNameOfTest; using elx::CoreMainGTestUtilities::GetTransformParametersFromFilter; using elx::CoreMainGTestUtilities::ImageDomain; +using elx::CoreMainGTestUtilities::TypeHolder; +using elx::CoreMainGTestUtilities::minimumImageSizeValue; +using elx::GTestUtilities::MakeMergedMap; template @@ -80,6 +87,7 @@ using ElastixRegistrationMethodType = itk::ElastixRegistrationMethod::Filled(minimumImageSizeValue); + const ImageDomain imageDomain(imageSize); + + elx::DefaultConstruct fixedImage{}; + imageDomain.ToImage(fixedImage); + fixedImage.Allocate(true); + + elx::DefaultConstruct movingImage{}; + imageDomain.ToImage(movingImage); + movingImage.Allocate(true); + + // Some "extreme" values to test if each of them is preserved during the transformation. + const std::array pixelValues{ PixelType{}, + PixelType{ 1 }, + std::numeric_limits::lowest(), + std::numeric_limits::min(), + PixelType{ std::numeric_limits::max() - 1 }, + std::numeric_limits::max() }; + std::copy(pixelValues.cbegin(), pixelValues.cend(), itk::ImageBufferRange(movingImage).begin()); + + // A dummy registration (that does not do any optimization). + elx::DefaultConstruct> registration{}; + + registration.SetParameterObject( + CreateParameterObject({ // Parameters in alphabetic order: + { "AutomaticParameterEstimation", { "false" } }, + { "FixedInternalImagePixelType", { ElxTypedef::FixedPixelTypeString } }, + { "MovingInternalImagePixelType", { ElxTypedef::MovingPixelTypeString } }, + { "ImageSampler", "Full" }, + { "MaximumNumberOfIterations", "0" }, + { "Metric", "AdvancedNormalizedCorrelation" }, + { "Optimizer", { "AdaptiveStochasticGradientDescent" } }, + { "ResampleInterpolator", { "FinalLinearInterpolator" } }, + { "Transform", "TranslationTransform" } })); + registration.SetFixedImage(&fixedImage); + registration.SetMovingImage(&movingImage); + registration.Update(); + + EXPECT_EQ(DerefRawPointer(registration.GetOutput()), movingImage); + }); +} + + +// Checks a zero-filled moving image with a random domain, having the same pixel type as any of the supported internal +// pixel types. +GTEST_TEST(itkElastixRegistrationMethod, CheckZeroFilledMovingImageWithRandomDomainHavingInternalPixelType) +{ + std::mt19937 randomNumberEngine{}; + + elx::ForEachSupportedImageType([&randomNumberEngine](const auto elxTypedef) { + using ElxTypedef = decltype(elxTypedef); + using ImageType = typename ElxTypedef::MovingImageType; + constexpr auto ImageDimension = ElxTypedef::MovingDimension; + + using PixelType = typename ImageType::PixelType; + + auto imageDomain = CreateRandomImageDomain(randomNumberEngine); + + // Reset index to avoid "FixedImageRegion does not overlap the fixed image buffered region" exceptions from + // itk::ImageToImageMetric::Initialize() + imageDomain.index = {}; + + // Create an image with values 1, 2, 3, ... N. We could have used arbitrary pixel values instead. + const auto fixedImage = CreateImageFilledWithSequenceOfNaturalNumbers(imageDomain); + + elx::DefaultConstruct movingImage{}; + imageDomain.ToImage(movingImage); + movingImage.Allocate(true); + + // A dummy registration (that does not do any optimization). + elx::DefaultConstruct> registration{}; + + registration.SetParameterObject( + CreateParameterObject({ // Parameters in alphabetic order: + { "AutomaticParameterEstimation", { "false" } }, + { "FixedInternalImagePixelType", { ElxTypedef::FixedPixelTypeString } }, + { "MovingInternalImagePixelType", { ElxTypedef::MovingPixelTypeString } }, + { "ImageSampler", "Full" }, + { "MaximumNumberOfIterations", "0" }, + { "Metric", "AdvancedNormalizedCorrelation" }, + { "Optimizer", { "AdaptiveStochasticGradientDescent" } }, + { "ResampleInterpolator", { "FinalLinearInterpolator" } }, + { "Transform", "TranslationTransform" } })); + registration.SetFixedImage(fixedImage); + registration.SetMovingImage(&movingImage); + registration.Update(); + + EXPECT_EQ(DerefRawPointer(registration.GetOutput()), movingImage); + }); +} + + +// Checks a minimum size moving image using any supported internal pixel type (which may be different from the input +// pixel type). +GTEST_TEST(itkElastixRegistrationMethod, CheckMinimumMovingImageUsingAnyInternalPixelType) +{ + const auto check = [](const auto inputPixelTypeHolder) { + elx::ForEachSupportedImageType([](const auto elxTypedef) { + using ElxTypedef = decltype(elxTypedef); + using InputPixelType = typename decltype(inputPixelTypeHolder)::Type; + using InputImageType = itk::Image; + + const ImageDomain imageDomain( + itk::Size::Filled(minimumImageSizeValue)); + + elx::DefaultConstruct fixedImage{}; + imageDomain.ToImage(fixedImage); + fixedImage.Allocate(true); + + const auto movingImage = CreateImageFilledWithSequenceOfNaturalNumbers(imageDomain); + + // A dummy registration (that does not do any optimization). + elx::DefaultConstruct> registration{}; + + registration.SetParameterObject( + CreateParameterObject({ // Parameters in alphabetic order: + { "AutomaticParameterEstimation", { "false" } }, + { "FixedInternalImagePixelType", { ElxTypedef::FixedPixelTypeString } }, + { "MovingInternalImagePixelType", { ElxTypedef::MovingPixelTypeString } }, + { "ImageSampler", "Full" }, + { "MaximumNumberOfIterations", "0" }, + { "Metric", "AdvancedNormalizedCorrelation" }, + { "Optimizer", { "AdaptiveStochasticGradientDescent" } }, + { "ResampleInterpolator", { "FinalLinearInterpolator" } }, + { "Transform", "TranslationTransform" } })); + registration.SetFixedImage(&fixedImage); + registration.SetMovingImage(movingImage); + registration.Update(); + + EXPECT_EQ(DerefRawPointer(registration.GetOutput()), DerefSmartPointer(movingImage)); + }); + }; + + check(TypeHolder{}); + check(TypeHolder{}); + check(TypeHolder{}); + check(TypeHolder{}); +} + + +// Checks a zero-filled moving image with a random domain, using any supported internal pixel type (which may be +// different from the input pixel type). +GTEST_TEST(itkElastixRegistrationMethod, CheckZeroFilledMovingImageWithRandomDomainUsingAnyInternalPixelType) +{ + std::mt19937 randomNumberEngine{}; + + const auto check = [&randomNumberEngine](const auto inputPixelTypeHolder) { + elx::ForEachSupportedImageType([&randomNumberEngine](const auto elxTypedef) { + using ElxTypedef = decltype(elxTypedef); + using InputPixelType = typename decltype(inputPixelTypeHolder)::Type; + using InputImageType = itk::Image; + + auto imageDomain = CreateRandomImageDomain(randomNumberEngine); + + // Reset index to avoid "FixedImageRegion does not overlap the fixed image buffered region" exceptions from + // itk::ImageToImageMetric::Initialize() + imageDomain.index = {}; + + // Create an image with values 1, 2, 3, ... N. We could have used arbitrary pixel values instead. + const auto fixedImage = CreateImageFilledWithSequenceOfNaturalNumbers(imageDomain); + + elx::DefaultConstruct movingImage{}; + imageDomain.ToImage(movingImage); + movingImage.Allocate(true); + + // A dummy registration (that does not do any optimization). + elx::DefaultConstruct> registration{}; + + registration.SetParameterObject( + CreateParameterObject({ // Parameters in alphabetic order: + { "AutomaticParameterEstimation", { "false" } }, + { "FixedInternalImagePixelType", { ElxTypedef::FixedPixelTypeString } }, + { "MovingInternalImagePixelType", { ElxTypedef::MovingPixelTypeString } }, + { "ImageSampler", "Full" }, + { "MaximumNumberOfIterations", "0" }, + { "Metric", "AdvancedNormalizedCorrelation" }, + { "Optimizer", { "AdaptiveStochasticGradientDescent" } }, + { "ResampleInterpolator", { "FinalLinearInterpolator" } }, + { "Transform", "TranslationTransform" } })); + registration.SetFixedImage(fixedImage); + registration.SetMovingImage(&movingImage); + registration.Update(); + + EXPECT_EQ(DerefRawPointer(registration.GetOutput()), movingImage); + }); + }; + + check(TypeHolder{}); + check(TypeHolder{}); + check(TypeHolder{}); + check(TypeHolder{}); +} diff --git a/Core/Main/GTesting/itkTransformixFilterGTest.cxx b/Core/Main/GTesting/itkTransformixFilterGTest.cxx index b22e6ba05..c54716915 100644 --- a/Core/Main/GTesting/itkTransformixFilterGTest.cxx +++ b/Core/Main/GTesting/itkTransformixFilterGTest.cxx @@ -27,6 +27,7 @@ #include "elxDefaultConstruct.h" #include "elxTransformIO.h" #include "GTesting/elxGTestUtilities.h" +#include "elxForEachSupportedImageType.h" #include @@ -66,13 +67,20 @@ using ParameterMapVectorType = elx::ParameterObject::ParameterMapVectorType; using elx::CoreMainGTestUtilities::CheckNew; using elx::CoreMainGTestUtilities::CreateImage; using elx::CoreMainGTestUtilities::CreateImageFilledWithSequenceOfNaturalNumbers; +using elx::CoreMainGTestUtilities::CreateParameterObject; +using elx::CoreMainGTestUtilities::CreateRandomImageDomain; using elx::CoreMainGTestUtilities::DerefRawPointer; using elx::CoreMainGTestUtilities::DerefSmartPointer; using elx::CoreMainGTestUtilities::FillImageRegion; using elx::CoreMainGTestUtilities::GetCurrentBinaryDirectoryPath; using elx::CoreMainGTestUtilities::GetDataDirectoryPath; using elx::CoreMainGTestUtilities::GetNameOfTest; +using elx::CoreMainGTestUtilities::GenerateRandomSign; +using elx::CoreMainGTestUtilities::ImageDomain; +using elx::CoreMainGTestUtilities::TypeHolder; +using elx::CoreMainGTestUtilities::minimumImageSizeValue; using elx::GTestUtilities::GeneratePseudoRandomParameters; +using elx::GTestUtilities::MakeMergedMap; template using DefaultConstructibleTransformixFilter = elx::DefaultConstruct>; @@ -90,15 +98,6 @@ ConvertToItkVector(const T & arg) } -auto -CreateParameterObject(const ParameterMapType & parameterMap) -{ - const auto parameterObject = CheckNew(); - parameterObject->SetParameterMap(parameterMap); - return parameterObject; -} - - template auto CreateDefaultDirectionParameterValues() @@ -417,6 +416,7 @@ Test_BSplineViaExternalTransformFile(const std::string & rootOutputDirectoryPath } } + } // namespace @@ -1292,3 +1292,170 @@ GTEST_TEST(itkTransformixFilter, ComputeSpatialJacobianDeterminantImage) EXPECT_EQ(matrix, expectedMatrix); } } + + +// Checks a minimum size moving image having the same pixel type as any of the supported internal pixel types. +GTEST_TEST(itkTransformixFilter, CheckMinimumMovingImageHavingInternalPixelType) +{ + elx::ForEachSupportedImageType([](const auto elxTypedef) { + using ElxTypedef = decltype(elxTypedef); + using ImageType = typename ElxTypedef::MovingImageType; + constexpr auto ImageDimension = ElxTypedef::MovingDimension; + + using PixelType = typename ImageType::PixelType; + using SizeType = itk::Size; + + const ImageDomain imageDomain( + itk::Size::Filled(minimumImageSizeValue)); + + elx::DefaultConstruct movingImage{}; + imageDomain.ToImage(movingImage); + movingImage.Allocate(true); + + // Some "extreme" values to test if each of them is preserved during the transformation. + const std::array pixelValues{ PixelType{}, + PixelType{ 1 }, + std::numeric_limits::lowest(), + std::numeric_limits::min(), + PixelType{ std::numeric_limits::max() - 1 }, + std::numeric_limits::max() }; + std::copy(pixelValues.cbegin(), pixelValues.cend(), itk::ImageBufferRange(movingImage).begin()); + + const ParameterMapType parameterMap = MakeMergedMap( + { // Parameters in alphabetic order: + { "FixedInternalImagePixelType", { ElxTypedef::FixedPixelTypeString } }, + { "MovingInternalImagePixelType", { ElxTypedef::MovingPixelTypeString } }, + { "ResampleInterpolator", { "FinalLinearInterpolator" } } }, + imageDomain.AsParameterMap()); + + const elx::DefaultConstruct> identityTransform{}; + elx::DefaultConstruct> transformixFilter{}; + + transformixFilter.SetMovingImage(&movingImage); + transformixFilter.SetTransform(&identityTransform); + transformixFilter.SetTransformParameterObject(CreateParameterObject(parameterMap)); + transformixFilter.Update(); + + EXPECT_EQ(DerefRawPointer(transformixFilter.GetOutput()), movingImage); + }); +} + + +// Checks a zero-filled moving image with a random domain, having the same pixel type as any of the supported internal +// pixel types. +GTEST_TEST(itkTransformixFilter, CheckZeroFilledMovingImageWithRandomDomainHavingInternalPixelType) +{ + std::mt19937 randomNumberEngine{}; + + elx::ForEachSupportedImageType([&randomNumberEngine](const auto elxTypedef) { + using ElxTypedef = decltype(elxTypedef); + using ImageType = typename ElxTypedef::MovingImageType; + constexpr auto ImageDimension = ElxTypedef::MovingDimension; + + using PixelType = typename ImageType::PixelType; + using SizeType = itk::Size; + + const auto imageDomain = CreateRandomImageDomain(randomNumberEngine); + + elx::DefaultConstruct movingImage{}; + imageDomain.ToImage(movingImage); + movingImage.Allocate(true); + + const ParameterMapType parameterMap = MakeMergedMap( + { // Parameters in alphabetic order: + { "FixedInternalImagePixelType", { ElxTypedef::FixedPixelTypeString } }, + { "MovingInternalImagePixelType", { ElxTypedef::MovingPixelTypeString } }, + { "ResampleInterpolator", { "FinalLinearInterpolator" } } }, + imageDomain.AsParameterMap()); + + const elx::DefaultConstruct> identityTransform{}; + elx::DefaultConstruct> transformixFilter{}; + + transformixFilter.SetMovingImage(&movingImage); + transformixFilter.SetTransform(&identityTransform); + transformixFilter.SetTransformParameterObject(CreateParameterObject(parameterMap)); + transformixFilter.Update(); + + EXPECT_EQ(DerefRawPointer(transformixFilter.GetOutput()), movingImage); + }); +} + + +// Checks a minimum size moving image using any supported internal pixel type (which may be different from the input +// pixel type). +GTEST_TEST(itkTransformixFilter, CheckMinimumMovingImageUsingAnyInternalPixelType) +{ + const auto check = [](const auto inputPixelTypeHolder) { + elx::ForEachSupportedImageType([](const auto elxTypedef) { + using ElxTypedef = decltype(elxTypedef); + using InputPixelType = typename decltype(inputPixelTypeHolder)::Type; + using InputImageType = itk::Image; + + const ImageDomain imageDomain( + itk::Size::Filled(minimumImageSizeValue)); + const auto movingImage = CreateImageFilledWithSequenceOfNaturalNumbers(imageDomain); + + const ParameterMapType parameterMap = MakeMergedMap( + { // Parameters in alphabetic order: + { "FixedInternalImagePixelType", { ElxTypedef::FixedPixelTypeString } }, + { "MovingInternalImagePixelType", { ElxTypedef::MovingPixelTypeString } }, + { "ResampleInterpolator", { "FinalLinearInterpolator" } } }, + imageDomain.AsParameterMap()); + + const elx::DefaultConstruct> identityTransform{}; + elx::DefaultConstruct> transformixFilter{}; + + transformixFilter.SetMovingImage(movingImage); + transformixFilter.SetTransform(&identityTransform); + transformixFilter.SetTransformParameterObject(CreateParameterObject(parameterMap)); + transformixFilter.Update(); + + EXPECT_EQ(DerefRawPointer(transformixFilter.GetOutput()), DerefSmartPointer(movingImage)); + }); + }; + + check(TypeHolder{}); + check(TypeHolder{}); + check(TypeHolder{}); +} + +// Checks a zero-filled moving image with a random domain, using any supported internal pixel type (which may be +// different from the input pixel type). +GTEST_TEST(itkTransformixFilter, CheckZeroFilledMovingImageWithRandomDomainUsingAnyInternalPixelType) +{ + std::mt19937 randomNumberEngine{}; + + const auto check = [&randomNumberEngine](const auto inputPixelTypeHolder) { + elx::ForEachSupportedImageType([&randomNumberEngine](const auto elxTypedef) { + using ElxTypedef = decltype(elxTypedef); + using InputPixelType = typename decltype(inputPixelTypeHolder)::Type; + using InputImageType = itk::Image; + + const auto imageDomain = CreateRandomImageDomain(randomNumberEngine); + elx::DefaultConstruct movingImage{}; + imageDomain.ToImage(movingImage); + movingImage.Allocate(true); + + const ParameterMapType parameterMap = MakeMergedMap( + { // Parameters in alphabetic order: + { "FixedInternalImagePixelType", { ElxTypedef::FixedPixelTypeString } }, + { "MovingInternalImagePixelType", { ElxTypedef::MovingPixelTypeString } }, + { "ResampleInterpolator", { "FinalLinearInterpolator" } } }, + imageDomain.AsParameterMap()); + + const elx::DefaultConstruct> identityTransform{}; + elx::DefaultConstruct> transformixFilter{}; + + transformixFilter.SetMovingImage(&movingImage); + transformixFilter.SetTransform(&identityTransform); + transformixFilter.SetTransformParameterObject(CreateParameterObject(parameterMap)); + transformixFilter.Update(); + + EXPECT_EQ(DerefRawPointer(transformixFilter.GetOutput()), movingImage); + }); + }; + + check(TypeHolder{}); + check(TypeHolder{}); + check(TypeHolder{}); +} diff --git a/Core/Main/elxForEachSupportedImageType.h b/Core/Main/elxForEachSupportedImageType.h index 60c4c8b14..5312dae8a 100644 --- a/Core/Main/elxForEachSupportedImageType.h +++ b/Core/Main/elxForEachSupportedImageType.h @@ -46,6 +46,24 @@ ForEachSupportedImageType(const TFunction & func) ForEachSupportedImageType(func, std::make_index_sequence()); } +template +bool +ForEachSupportedImageTypeUntilTrue(const TFunction & func, const std::index_sequence &) +{ + // Expand the sequence by a fold expression of the form (func(n) || func(n - 1) || ... || func(1)). + return (func(elx::ElastixTypedef{}) || ...); +} + + +/** Runs a function `func(ElastixTypedef{})`, for each supported image type from "elxSupportedImageTypes.h", + * until the function returns true. */ +template +bool +ForEachSupportedImageTypeUntilTrue(const TFunction & func) +{ + return ForEachSupportedImageTypeUntilTrue(func, std::make_index_sequence()); +} + } // namespace elastix diff --git a/Core/Main/elxLibUtilities.cxx b/Core/Main/elxLibUtilities.cxx index 7498b9aed..b59ece42f 100644 --- a/Core/Main/elxLibUtilities.cxx +++ b/Core/Main/elxLibUtilities.cxx @@ -28,7 +28,6 @@ LibUtilities::SetParameterValueAndWarnOnOverride(ParameterMapType & parameterMa const std::string & parameterName, const std::string & parameterValue) { - if (const auto found = parameterMap.find(parameterName); found == parameterMap.end()) { parameterMap[parameterName] = { parameterValue }; @@ -44,4 +43,22 @@ LibUtilities::SetParameterValueAndWarnOnOverride(ParameterMapType & parameterMa } } + +std::string +LibUtilities::RetrievePixelTypeParameterValue(const ParameterMapType & parameterMap, const std::string & parameterName) +{ + if (const auto found = parameterMap.find(parameterName); found != parameterMap.end()) + { + if (const auto & second = found->second; !second.empty()) + { + if (const auto & front = second.front(); !front.empty()) + { + return front; + } + } + } + return "float"; +} + + } // namespace elastix diff --git a/Core/Main/elxLibUtilities.h b/Core/Main/elxLibUtilities.h index 4fad1d7b1..81cf51269 100644 --- a/Core/Main/elxLibUtilities.h +++ b/Core/Main/elxLibUtilities.h @@ -18,8 +18,15 @@ #ifndef elxLibUtilities_h #define elxLibUtilities_h +#include "elxForEachSupportedImageType.h" + +#include +#include +#include + #include #include +#include // For is_same_v. #include @@ -35,6 +42,58 @@ SetParameterValueAndWarnOnOverride(ParameterMapType & parameterMap, const std::string & parameterName, const std::string & parameterValue); + +/** Retrieves the PixelType string value of the specified parameter. Returns "float" by default. */ +std::string +RetrievePixelTypeParameterValue(const ParameterMapType & parameterMap, const std::string & parameterName); + +template +itk::SmartPointer +CastToInternalPixelType(itk::SmartPointer inputImage, const std::string & internalPixelTypeString) +{ + if (inputImage == nullptr) + { + itkGenericExceptionMacro("The specified input image should not be null!"); + } + + itk::SmartPointer outputImage; + + elx::ForEachSupportedImageTypeUntilTrue([inputImage, &outputImage, &internalPixelTypeString](const auto elxTypedef) { + using ElxTypedef = decltype(elxTypedef); + + if constexpr (TInputImage::ImageDimension == ElxTypedef::MovingDimension) + { + using InternalImageType = typename ElxTypedef::MovingImageType; + + if (internalPixelTypeString == ElxTypedef::MovingPixelTypeString) + { + if constexpr (std::is_same_v) + { + outputImage = inputImage; + } + else + { + const auto castFilter = itk::CastImageFilter::New(); + castFilter->SetInput(inputImage); + castFilter->Update(); + outputImage = castFilter->GetOutput(); + } + return true; + } + } + return false; + }); + + if (outputImage == nullptr) + { + itkGenericExceptionMacro("Failed to cast to the specified internal pixel type \"" + << internalPixelTypeString + << "\". It may need to be added to the CMake variable ELASTIX_IMAGE_" + << TInputImage::ImageDimension << "D_PIXELTYPES."); + } + return outputImage; +} + } // namespace elastix::LibUtilities #endif diff --git a/Core/Main/itkElastixRegistrationMethod.h b/Core/Main/itkElastixRegistrationMethod.h index 7afd7b69f..d4d42ca3d 100644 --- a/Core/Main/itkElastixRegistrationMethod.h +++ b/Core/Main/itkElastixRegistrationMethod.h @@ -307,6 +307,23 @@ class ITK_TEMPLATE_EXPORT ElastixRegistrationMethod : public itk::ImageSource + std::vector + GetInputImages(const char * const inputTypeString) + { + std::vector images; + for (const auto & inputName : this->GetInputNames()) + { + if (this->IsInputOfType(inputTypeString, inputName)) + { + images.push_back(itkDynamicCastInDebugMode(this->ProcessObject::GetInput(inputName))); + } + } + return images; + } + /** Private using-declaration, just to avoid GCC compilation warnings: '...' was hidden [-Woverloaded-virtual] */ using Superclass::SetInput; diff --git a/Core/Main/itkElastixRegistrationMethod.hxx b/Core/Main/itkElastixRegistrationMethod.hxx index 9e2cdda93..5fea5a60b 100644 --- a/Core/Main/itkElastixRegistrationMethod.hxx +++ b/Core/Main/itkElastixRegistrationMethod.hxx @@ -77,6 +77,8 @@ template void ElastixRegistrationMethod::GenerateData() { + using elx::LibUtilities::CastToInternalPixelType; + using elx::LibUtilities::RetrievePixelTypeParameterValue; using elx::LibUtilities::SetParameterValueAndWarnOnOverride; // Force compiler to instantiate the image dimensions, otherwise we may get @@ -86,8 +88,6 @@ ElastixRegistrationMethod::GenerateData() const unsigned int fixedImageDimension = FixedImageDimension; const unsigned int movingImageDimension = MovingImageDimension; - DataObjectContainerPointer fixedImageContainer = DataObjectContainerType::New(); - DataObjectContainerPointer movingImageContainer = DataObjectContainerType::New(); DataObjectContainerPointer fixedMaskContainer = nullptr; DataObjectContainerPointer movingMaskContainer = nullptr; DataObjectContainerPointer resultImageContainer = nullptr; @@ -98,18 +98,6 @@ ElastixRegistrationMethod::GenerateData() // Split inputs into separate containers for (const auto & inputName : this->GetInputNames()) { - if (this->IsInputOfType("FixedImage", inputName)) - { - fixedImageContainer->push_back(this->ProcessObject::GetInput(inputName)); - continue; - } - - if (this->IsInputOfType("MovingImage", inputName)) - { - movingImageContainer->push_back(this->ProcessObject::GetInput(inputName)); - continue; - } - if (this->IsInputOfType("FixedMask", inputName)) { if (fixedMaskContainer.IsNull()) @@ -200,18 +188,46 @@ ElastixRegistrationMethod::GenerateData() const auto fixedImageDimensionString = std::to_string(fixedImageDimension); const auto fixedImagePixelTypeString = elx::PixelTypeToString(); const auto movingImageDimensionString = std::to_string(movingImageDimension); - const auto movingImagePixelTypeString = elx::PixelTypeToString(); + + const std::vector fixedInputImages = GetInputImages("FixedImage"); + const std::vector movingInputImages = GetInputImages("MovingImage"); + + // Cache the containers of internal images, to avoid casting the input images to the same internal pixel type more + // than once, in the `for` loop below here. + std::map> fixedInternalImageContainers; + std::map> movingInternalImageContainers; // Run the (possibly multiple) registration(s) for (unsigned int i = 0; i < parameterMapVector.size(); ++i) { auto & parameterMap = parameterMapVector[i]; - // Set image dimension and pixel type from input images (overrides user settings) + // Lambda to create a FixedImageContainer or a MovingImageContainer. + const auto createImageContainer = + [](const auto & inputImages, const auto internalPixelTypeString, auto & imageContainers) { + if (const auto found = imageContainers.find(internalPixelTypeString); found == imageContainers.end()) + { + const auto imageContainer = DataObjectContainerType::New(); + + for (const auto inputImage : inputImages) + { + const auto internalImage = CastToInternalPixelType(inputImage, internalPixelTypeString); + imageContainer->push_back(internalImage); + } + imageContainers[internalPixelTypeString] = imageContainer; + return imageContainer; + } + else + { + // There was an image container already for the specified internal pixel type. So just use that one! + return found->second; + } + }; + + // Set image dimension from input images (overrides user settings) SetParameterValueAndWarnOnOverride(parameterMap, "FixedImageDimension", fixedImageDimensionString); - SetParameterValueAndWarnOnOverride(parameterMap, "FixedInternalImagePixelType", fixedImagePixelTypeString); SetParameterValueAndWarnOnOverride(parameterMap, "MovingImageDimension", movingImageDimensionString); - SetParameterValueAndWarnOnOverride(parameterMap, "MovingInternalImagePixelType", movingImagePixelTypeString); + SetParameterValueAndWarnOnOverride(parameterMap, "ResultImagePixelType", fixedImagePixelTypeString); // Initial transform parameter files are handled via arguments and enclosing loop, not @@ -228,8 +244,14 @@ ElastixRegistrationMethod::GenerateData() // Set stuff we get from a previous registration elastixMain->SetInitialTransform(transform); - elastixMain->SetFixedImageContainer(fixedImageContainer); - elastixMain->SetMovingImageContainer(movingImageContainer); + elastixMain->SetFixedImageContainer( + createImageContainer(fixedInputImages, + RetrievePixelTypeParameterValue(parameterMap, "FixedInternalImagePixelType"), + fixedInternalImageContainers)); + elastixMain->SetMovingImageContainer( + createImageContainer(movingInputImages, + RetrievePixelTypeParameterValue(parameterMap, "MovingInternalImagePixelType"), + movingInternalImageContainers)); elastixMain->SetFixedMaskContainer(fixedMaskContainer); elastixMain->SetMovingMaskContainer(movingMaskContainer); elastixMain->SetResultImageContainer(resultImageContainer); @@ -256,8 +278,6 @@ ElastixRegistrationMethod::GenerateData() // Get stuff in order to put it in the next registration transform = elastixMain->GetFinalTransform(); - fixedImageContainer = elastixMain->GetFixedImageContainer(); - movingImageContainer = elastixMain->GetMovingImageContainer(); fixedMaskContainer = elastixMain->GetFixedMaskContainer(); movingMaskContainer = elastixMain->GetMovingMaskContainer(); resultImageContainer = elastixMain->GetResultImageContainer(); diff --git a/Core/Main/itkTransformixFilter.hxx b/Core/Main/itkTransformixFilter.hxx index 3fdc7971e..90bf78171 100644 --- a/Core/Main/itkTransformixFilter.hxx +++ b/Core/Main/itkTransformixFilter.hxx @@ -64,6 +64,8 @@ template void TransformixFilter::GenerateData() { + using elx::LibUtilities::CastToInternalPixelType; + using elx::LibUtilities::RetrievePixelTypeParameterValue; using elx::LibUtilities::SetParameterValueAndWarnOnOverride; // Force compiler to instantiate the image dimension, otherwise we may get @@ -153,15 +155,6 @@ TransformixFilter::GenerateData() const auto transformixMain = elx::TransformixMain::New(); m_TransformixMain = transformixMain; - // Setup transformix for warping input image if given - DataObjectContainerPointer inputImageContainer = nullptr; - if (!this->IsEmpty(this->GetMovingImage())) - { - inputImageContainer = DataObjectContainerType::New(); - inputImageContainer->InsertElement(0, const_cast(this->GetMovingImage())); - transformixMain->SetInputImageContainer(inputImageContainer); - } - // Get ParameterMap ParameterObjectPointer transformParameterObject = this->GetTransformParameterObject(); ParameterMapVectorType transformParameterMapVector = transformParameterObject->GetParameterMaps(); @@ -172,6 +165,22 @@ TransformixFilter::GenerateData() itkExceptionMacro("Empty parameter map in parameter object."); } + // Setup transformix for warping input image if given + DataObjectContainerPointer inputImageContainer = nullptr; + + if (const auto movingImage = itkDynamicCastInDebugMode(this->ProcessObject::GetInput("MovingImage")); + !Self::IsEmpty(movingImage)) + { + // Note that the internal pixel type is only retrieved from the very first transform parameter map. + const auto internalPixelTypeString = + RetrievePixelTypeParameterValue(transformParameterMapVector.front(), "MovingInternalImagePixelType"); + const auto internalImage = CastToInternalPixelType(movingImage, internalPixelTypeString); + + inputImageContainer = DataObjectContainerType::New(); + inputImageContainer->push_back(internalImage); + transformixMain->SetInputImageContainer(inputImageContainer); + } + if (m_Transform) { // Adjust the local transformParameterMap according to this m_Transform. @@ -251,11 +260,7 @@ TransformixFilter::GenerateData() auto & transformParameterMap = transformParameterMapVector[i]; SetParameterValueAndWarnOnOverride(transformParameterMap, "FixedImageDimension", movingImageDimensionString); - SetParameterValueAndWarnOnOverride( - transformParameterMap, "FixedInternalImagePixelType", movingImagePixelTypeString); SetParameterValueAndWarnOnOverride(transformParameterMap, "MovingImageDimension", movingImageDimensionString); - SetParameterValueAndWarnOnOverride( - transformParameterMap, "MovingInternalImagePixelType", movingImagePixelTypeString); SetParameterValueAndWarnOnOverride(transformParameterMap, "ResultImagePixelType", movingImagePixelTypeString); }