From 510586750ad3514fd6f52c26b1925a53f503fe01 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C5=BEenan=20Zuki=C4=87?= Date: Fri, 9 Apr 2021 11:51:19 -0400 Subject: [PATCH] ENH: Add Transform::ApplyToImageMetadata method This is equivalent to Slicer's "Harden Transform" operation. --- Modules/Core/Transform/include/itkTransform.h | 22 +++++ .../Core/Transform/include/itkTransform.hxx | 46 +++++++++ .../itkResampleInPlaceImageFilterTest.cxx | 99 +++++++++++-------- 3 files changed, 127 insertions(+), 40 deletions(-) diff --git a/Modules/Core/Transform/include/itkTransform.h b/Modules/Core/Transform/include/itkTransform.h index 15832f4f586..cf19ca42b43 100644 --- a/Modules/Core/Transform/include/itkTransform.h +++ b/Modules/Core/Transform/include/itkTransform.h @@ -18,6 +18,7 @@ #ifndef itkTransform_h #define itkTransform_h +#include // For std::enable_if #include "itkTransformBase.h" #include "itkVector.h" #include "itkSymmetricSecondRankTensor.h" @@ -542,6 +543,27 @@ class ITK_TEMPLATE_EXPORT Transform : public TransformBaseTemplate + typename std::enable_if::type + ApplyToImageMetadata(TImage * image) const; + template + typename std::enable_if::type + ApplyToImageMetadata(SmartPointer image) const + { + this->ApplyToImageMetadata(image.GetPointer()); // Delegate to the raw pointer signature + } + protected: /** * Clone the current transform. diff --git a/Modules/Core/Transform/include/itkTransform.hxx b/Modules/Core/Transform/include/itkTransform.hxx index 1f02fddad80..d574cf948b2 100644 --- a/Modules/Core/Transform/include/itkTransform.hxx +++ b/Modules/Core/Transform/include/itkTransform.hxx @@ -21,6 +21,7 @@ #include "itkTransform.h" #include "itkCrossHelper.h" #include "vnl/algo/vnl_svd_fixed.h" +#include "itkMetaProgrammingLibrary.h" namespace itk { @@ -470,6 +471,51 @@ Transform::TransformS return outputTensor; } +template +template +typename std::enable_if::type +Transform::ApplyToImageMetadata(TImage * image) const +{ + using ImageType = TImage; + + if (!this->IsLinear()) + { + itkWarningMacro(<< "ApplyToImageMetadata was invoked with non-linear transform of type: " << this->GetNameOfClass() + << ". This might produce unexpected results."); + } + + typename Self::Pointer inverse = this->GetInverseTransform(); + + // transform origin + typename ImageType::PointType origin = image->GetOrigin(); + origin = inverse->TransformPoint(origin); + image->SetOrigin(origin); + + typename ImageType::SpacingType spacing = image->GetSpacing(); + typename ImageType::DirectionType direction = image->GetDirection(); + // transform direction cosines and compute new spacing + for (unsigned i = 0; i < ImageType::ImageDimension; i++) + { + Vector dirVector; + for (unsigned k = 0; k < ImageType::ImageDimension; k++) + { + dirVector[k] = direction[k][i]; + } + + dirVector *= spacing[i]; + dirVector = inverse->TransformVector(dirVector); + spacing[i] = dirVector.Normalize(); + + for (unsigned k = 0; k < ImageType::ImageDimension; k++) + { + direction[k][i] = dirVector[k]; + } + } + image->SetDirection(direction); + image->SetSpacing(spacing); +} + #if !defined(ITK_LEGACY_REMOVE) template diff --git a/Modules/Core/Transform/test/itkResampleInPlaceImageFilterTest.cxx b/Modules/Core/Transform/test/itkResampleInPlaceImageFilterTest.cxx index 3713a10a2c4..244ede35101 100644 --- a/Modules/Core/Transform/test/itkResampleInPlaceImageFilterTest.cxx +++ b/Modules/Core/Transform/test/itkResampleInPlaceImageFilterTest.cxx @@ -34,6 +34,49 @@ Validate(double input, double desired, double tolerance) return std::abs(input - desired) > tolerance * std::abs(desired); } +// Returns true if images are different +template +bool +imagesDifferent(ImageType * baselineImage, ImageType * outputImage) +{ + double tol = 1.e-3; // tolerance + + typename ImageType::PointType origin = outputImage->GetOrigin(); + typename ImageType::DirectionType direction = outputImage->GetDirection(); + typename ImageType::SpacingType spacing = outputImage->GetSpacing(); + + typename ImageType::PointType origin_d = baselineImage->GetOrigin(); + typename ImageType::DirectionType direction_d = baselineImage->GetDirection(); + typename ImageType::SpacingType spacing_d = baselineImage->GetSpacing(); + + // Image info validation + bool result = false; // no difference by default + for (unsigned int i = 0; i < ImageType::ImageDimension; ++i) + { + result = (result || Validate(origin[i], origin_d[i], tol)); + result = (result || Validate(spacing[i], spacing_d[i], tol)); + for (unsigned int j = 0; j < ImageType::ImageDimension; ++j) + { + result = (result || Validate(direction(i, j), direction_d(i, j), tol)); + } + } + + // Voxel contents validation + using ImageConstIterator = itk::ImageRegionConstIterator; + ImageConstIterator it1(outputImage, outputImage->GetRequestedRegion()); + ImageConstIterator it2(baselineImage, baselineImage->GetRequestedRegion()); + it1.GoToBegin(); + it2.GoToBegin(); + while (!it1.IsAtEnd()) + { + result = (result || Validate(it1.Get(), it2.Get(), tol)); + ++it1; + ++it2; + } + + return result; +} + int itkResampleInPlaceImageFilterTest(int argc, char * argv[]) { @@ -45,19 +88,12 @@ itkResampleInPlaceImageFilterTest(int argc, char * argv[]) return EXIT_FAILURE; } - double tol = 1.e-3; // tolerance - constexpr unsigned int Dimension = 3; using PixelType = short; using ImageType = itk::Image; using ImagePointer = ImageType::Pointer; - using ImagePointType = ImageType::PointType; - using ImageDirectionType = ImageType::DirectionType; - using ImageSpacingType = ImageType::SpacingType; - using ImageConstIterator = itk::ImageRegionConstIterator; using TransformType = itk::VersorRigid3DTransform; - using FilterType = itk::ResampleInPlaceImageFilter; // Read in input test image @@ -92,50 +128,33 @@ itkResampleInPlaceImageFilterTest(int argc, char * argv[]) // Set up the resample filter FilterType::Pointer filter = FilterType::New(); ITK_EXERCISE_BASIC_OBJECT_METHODS(filter, ResampleInPlaceImageFilter, ImageToImageFilter); + filter->SetInputImage(inputImage); ITK_TEST_SET_GET_VALUE(inputImage, filter->GetInputImage()); filter->SetRigidTransform(transform); ITK_TEST_SET_GET_VALUE(transform, filter->GetRigidTransform()); ITK_TRY_EXPECT_NO_EXCEPTION(filter->Update()); ImagePointer outputImage = filter->GetOutput(); - - // Get image info - ImagePointType origin = outputImage->GetOrigin(); - ImageDirectionType direction = outputImage->GetDirection(); - ImageSpacingType spacing = outputImage->GetSpacing(); + ITK_TRY_EXPECT_NO_EXCEPTION(itk::WriteImage(outputImage, argv[3])); // Read in baseline image ImagePointer baselineImage = nullptr; ITK_TRY_EXPECT_NO_EXCEPTION(baselineImage = itk::ReadImage(argv[2])); - ImagePointType origin_d = baselineImage->GetOrigin(); - ImageDirectionType direction_d = baselineImage->GetDirection(); - ImageSpacingType spacing_d = baselineImage->GetSpacing(); - // Image info validation - bool result = false; // test result default = no failure - for (unsigned int i = 0; i < Dimension; ++i) - { - result = (result || Validate(origin[i], origin_d[i], tol)); - result = (result || Validate(spacing[i], spacing_d[i], tol)); - for (unsigned int j = 0; j < Dimension; ++j) - { - result = (result || Validate(direction(i, j), direction_d(i, j), tol)); - } - } - - // Voxel contents validation - ImageConstIterator it1(outputImage, outputImage->GetRequestedRegion()); - ImageConstIterator it2(baselineImage, baselineImage->GetRequestedRegion()); - it1.GoToBegin(); - it2.GoToBegin(); - while (!it1.IsAtEnd()) - { - result = (result || Validate(it1.Get(), it2.Get(), tol)); - ++it1; - ++it2; - } - - ITK_TRY_EXPECT_NO_EXCEPTION(itk::WriteImage(outputImage, argv[3])); + // Now do comparisons + bool result = imagesDifferent(baselineImage, outputImage); + transform->ApplyToImageMetadata(inputImage); + result = (result || imagesDifferent(baselineImage, inputImage)); + + // Make sure we can invoke ApplyToImageMetadata via const/raw pointer + TransformType * rawPointerTransform = transform.GetPointer(); + rawPointerTransform->ApplyToImageMetadata(inputImage); + rawPointerTransform->ApplyToImageMetadata(inputImage.GetPointer()); + const TransformType * constRawPointerTransform = transform.GetPointer(); + constRawPointerTransform->ApplyToImageMetadata(inputImage); + TransformType::ConstPointer constPointerTransform = transform.GetPointer(); + constPointerTransform->ApplyToImageMetadata(inputImage); + constPointerTransform->ApplyToImageMetadata(inputImage.GetPointer()); return result; }