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); }