ENH: Add Transform::ApplyToImageMetadata method

This is equivalent to Slicer's "Harden Transform" operation.

Author: dzenanz

Modules/Core/Transform/include/itkTransform.h

@@ -18,6 +18,7 @@
 #ifndef itkTransform_h
 #define itkTransform_h
 
+#include <type_traits> // For std::enable_if
 #include "itkTransformBase.h"
 #include "itkVector.h"
 #include "itkSymmetricSecondRankTensor.h"
@@ -542,6 +543,27 @@ class ITK_TEMPLATE_EXPORT Transform : public TransformBaseTemplate<TParametersVa
   itkLegacyMacro(virtual void ComputeInverseJacobianWithRespectToPosition(const InputPointType & x,
                                                                           JacobianType &         jacobian) const);
 
+  /** Apply this transform to an image without resampling.
+   *
+   * Updates image metadata (origin, spacing, direction cosines matrix) in place.
+   *
+   * Only available when input and output space are of the same dimension.
+   * Only works properly for linear transforms.
+   *
+   * The image parameter may be either a SmartPointer or a raw pointer.
+   * */
+  template <typename TImage>
+  typename std::enable_if<TImage::ImageDimension == NInputDimensions && TImage::ImageDimension == NOutputDimensions,
+                          void>::type
+  ApplyToImageMetadata(TImage * image) const;
+  template <typename TImage>
+  typename std::enable_if<TImage::ImageDimension == NInputDimensions && TImage::ImageDimension == NOutputDimensions,
+                          void>::type
+  ApplyToImageMetadata(SmartPointer<TImage> image) const
+  {
+    this->ApplyToImageMetadata(image.GetPointer()); // Delegate to the raw pointer signature
+  }
+
 protected:
   /**
    * Clone the current transform.

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<TParametersValueType, NInputDimensions, NOutputDimensions>::TransformS
   return outputTensor;
 }
 
+template <typename TParametersValueType, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+template <typename TImage>
+typename std::enable_if<TImage::ImageDimension == NInputDimensions && TImage::ImageDimension == NOutputDimensions,
+                        void>::type
+Transform<TParametersValueType, NInputDimensions, NOutputDimensions>::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<typename Self::ParametersValueType, ImageType::ImageDimension> 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 <typename TParametersValueType, unsigned int NInputDimensions, unsigned int NOutputDimensions>

Modules/Core/Transform/test/itkResampleInPlaceImageFilterTest.cxx

@@ -34,6 +34,49 @@ Validate(double input, double desired, double tolerance)
   return std::abs<double>(input - desired) > tolerance * std::abs<double>(desired);
 }
 
+// Returns true if images are different
+template <typename ImageType>
+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<ImageType>;
+  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<PixelType, Dimension>;
   using ImagePointer = ImageType::Pointer;
-  using ImagePointType = ImageType::PointType;
-  using ImageDirectionType = ImageType::DirectionType;
-  using ImageSpacingType = ImageType::SpacingType;
-  using ImageConstIterator = itk::ImageRegionConstIterator<ImageType>;
   using TransformType = itk::VersorRigid3DTransform<double>;
-
   using FilterType = itk::ResampleInPlaceImageFilter<ImageType, ImageType>;
 
   // 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<ImageType>(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<ImageType>(baselineImage, outputImage);
+  transform->ApplyToImageMetadata(inputImage);
+  result = (result || imagesDifferent<ImageType>(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;
 }