ENH: Add SetCombinationTransform(transform) to TransformixFilter

SuperElastix · Dec 7, 2022 · a2e140d · a2e140d
1 parent 2b81bc8
commit a2e140d
Show file tree

Hide file tree

Showing 8 changed files with 130 additions and 11 deletions.
diff --git a/Core/Kernel/elxElastixBase.h b/Core/Kernel/elxElastixBase.h
@@ -303,7 +303,7 @@ class ElastixBase
 
   /** Empty ApplyTransform()-function to be overridden. */
   virtual int
-  ApplyTransform() = 0;
+  ApplyTransform(bool doReadTransform) = 0;
 
   /** Function that is called at the very beginning of ElastixTemplate::Run().
    * It checks the command line input arguments.

diff --git a/Core/Kernel/elxElastixTemplate.h b/Core/Kernel/elxElastixTemplate.h
@@ -217,7 +217,7 @@ class ITK_TEMPLATE_EXPORT ElastixTemplate final : public ElastixBase
   Run() override;
 
   int
-  ApplyTransform() override;
+  ApplyTransform(bool doReadTransform) override;
 
   /** The Callback functions. */
   int

diff --git a/Core/Kernel/elxElastixTemplate.hxx b/Core/Kernel/elxElastixTemplate.hxx
@@ -255,7 +255,7 @@ ElastixTemplate<TFixedImage, TMovingImage>::Run()
 
 template <class TFixedImage, class TMovingImage>
 int
-ElastixTemplate<TFixedImage, TMovingImage>::ApplyTransform()
+ElastixTemplate<TFixedImage, TMovingImage>::ApplyTransform(const bool doReadTransform)
 {
   /** Timer. */
   itk::TimeProbe timer;
@@ -307,7 +307,11 @@ ElastixTemplate<TFixedImage, TMovingImage>::ApplyTransform()
   auto & elxTransformBase = *(this->GetElxTransformBase());
 
   elxResamplerBase.ReadFromFile();
-  elxTransformBase.ReadFromFile();
+
+  if (doReadTransform)
+  {
+    elxTransformBase.ReadFromFile();
+  }
 
   /** Tell the user. */
   timer.Stop();

diff --git a/Core/Kernel/elxTransformixMain.cxx b/Core/Kernel/elxTransformixMain.cxx
@@ -45,6 +45,19 @@ namespace elastix
 
 int
 TransformixMain::Run()
+{
+  return RunWithTransform(nullptr);
+}
+
+/**
+ * **************************** RunWithTransform *****************************
+ *
+ * Assuming EnterCommandLineParameters has already been invoked.
+ * or that m_Configuration is initialized in another way.
+ */
+
+int
+TransformixMain::RunWithTransform(itk::TransformBase * const transform)
 {
   /** Set process properties. */
   this->SetProcessPriority();
@@ -113,7 +126,16 @@ TransformixMain::Run()
 
   elastixBase.SetResamplerContainer(this->CreateComponents("Resampler", "DefaultResampler", errorCode));
 
-  elastixBase.SetTransformContainer(this->CreateComponents("Transform", "", errorCode));
+  if (transform)
+  {
+    const auto transformContainer = elx::ElastixBase::ObjectContainerType::New();
+    transformContainer->push_back(transform);
+    elastixBase.SetTransformContainer(transformContainer);
+  }
+  else
+  {
+    elastixBase.SetTransformContainer(this->CreateComponents("Transform", "", errorCode));
+  }
 
   /** Check if all components could be created. */
   if (errorCode != 0)
@@ -135,7 +157,7 @@ TransformixMain::Run()
   /** ApplyTransform! */
   try
   {
-    errorCode = elastixBase.ApplyTransform();
+    errorCode = elastixBase.ApplyTransform(transform == nullptr);
   }
   catch (const itk::ExceptionObject & excp)
   {
@@ -185,10 +207,12 @@ TransformixMain::Run(const ArgumentMapType & argmap, const ParameterMapType & in
  */
 
 int
-TransformixMain::Run(const ArgumentMapType & argmap, const std::vector<ParameterMapType> & inputMaps)
+TransformixMain::Run(const ArgumentMapType &               argmap,
+                     const std::vector<ParameterMapType> & inputMaps,
+                     itk::TransformBase * const            transform)
 {
   this->EnterCommandLineArguments(argmap, inputMaps);
-  return this->Run();
+  return this->RunWithTransform(transform);
 } // end Run()
 
 

diff --git a/Core/Kernel/elxTransformixMain.h b/Core/Kernel/elxTransformixMain.h
@@ -19,6 +19,7 @@
 #define elxTransformixMain_h
 
 #include "elxElastixMain.h"
+#include <itkTransformBase.h>
 
 namespace elastix
 {
@@ -90,8 +91,8 @@ class TransformixMain : public ElastixMain
   Run(const ArgumentMapType & argmap, const ParameterMapType & inputMap) override;
 
   /** Run version for using transformix as library. */
-  virtual int
-  Run(const ArgumentMapType & argmap, const std::vector<ParameterMapType> & inputMaps);
+  int
+  Run(const ArgumentMapType & argmap, const std::vector<ParameterMapType> & inputMaps, itk::TransformBase * = nullptr);
 
   /** Get and Set input- and outputImage. */
   virtual void
@@ -106,6 +107,11 @@ class TransformixMain : public ElastixMain
    */
   int
   InitDBIndex() override;
+
+private:
+  /** Does Run with an optionally specified Transform. */
+  int
+  RunWithTransform(itk::TransformBase *);
 };
 
 } // end namespace elastix

diff --git a/Core/Main/GTesting/itkTransformixFilterGTest.cxx b/Core/Main/GTesting/itkTransformixFilterGTest.cxx
@@ -64,6 +64,7 @@ using ParameterMapVectorType = elx::ParameterObject::ParameterMapVectorType;
 
 // Using-declarations:
 using elx::CoreMainGTestUtilities::CheckNew;
+using elx::CoreMainGTestUtilities::CreateImage;
 using elx::CoreMainGTestUtilities::CreateImageFilledWithSequenceOfNaturalNumbers;
 using elx::CoreMainGTestUtilities::Deref;
 using elx::CoreMainGTestUtilities::DerefSmartPointer;
@@ -1022,6 +1023,84 @@ GTEST_TEST(itkTransformixFilter, SetTranslationTransform)
 }
 
 
+GTEST_TEST(itkTransformixFilter, SetCombinationTransform)
+{
+  constexpr auto ImageDimension = 2U;
+  using PixelType = float;
+  using ImageType = itk::Image<PixelType, ImageDimension>;
+  const itk::Size<ImageDimension> imageSize{ { 5, 6 } };
+
+  const auto numberOfPixels =
+    std::accumulate(imageSize.cbegin(), imageSize.cend(), std::size_t{ 1 }, std::multiplies<>{});
+
+  const auto fixedImage = CreateImage<PixelType>(imageSize);
+  const auto movingImage = CreateImage<PixelType>(imageSize);
+
+  std::mt19937 randomNumberEngine;
+
+  std::generate_n(fixedImage->GetBufferPointer(), numberOfPixels, [&randomNumberEngine] {
+    return std::uniform_real_distribution<>{ 1, 32 }(randomNumberEngine);
+  });
+  std::generate_n(movingImage->GetBufferPointer(), numberOfPixels, [&randomNumberEngine] {
+    return std::uniform_real_distribution<>{ 32, 64 }(randomNumberEngine);
+  });
+
+  EXPECT_NE(*movingImage, *fixedImage);
+
+  for (const bool useInitialTransform : { false, true })
+  {
+    const std::string initialTransformParameterFileName =
+      useInitialTransform ? (GetDataDirectoryPath() + "/Translation(1,-2)/TransformParameters.txt") : "";
+
+    for (const char * const transformName : { "AffineTransform",
+                                              "BSplineTransform",
+                                              "EulerTransform",
+                                              "RecursiveBSplineTransform",
+                                              "SimilarityTransform",
+                                              "TranslationTransform" })
+    {
+      elx::DefaultConstruct<itk::ElastixRegistrationMethod<ImageType, ImageType>> registration;
+      registration.SetFixedImage(fixedImage);
+      registration.SetMovingImage(movingImage);
+      registration.SetInitialTransformParameterFileName(initialTransformParameterFileName);
+      registration.SetParameterObject(CreateParameterObject({ // Parameters in alphabetic order:
+                                                              { "AutomaticTransformInitialization", { "false" } },
+                                                              { "ImageSampler", { "Full" } },
+                                                              { "MaximumNumberOfIterations", { "2" } },
+                                                              { "Metric", { "AdvancedNormalizedCorrelation" } },
+                                                              { "Optimizer", { "AdaptiveStochasticGradientDescent" } },
+                                                              { "ResampleInterpolator", { "FinalLinearInterpolator" } },
+                                                              { "Transform", { transformName } } }));
+      registration.Update();
+
+      const ImageType & registrationOutputImage = Deref(registration.GetOutput());
+
+      EXPECT_NE(registrationOutputImage, *fixedImage);
+      EXPECT_NE(registrationOutputImage, *movingImage);
+
+      const auto combinationTransform = registration.GetCombinationTransform();
+
+      EXPECT_NE(combinationTransform, nullptr);
+
+      elx::DefaultConstruct<itk::TransformixFilter<ImageType>> transformixFilter{};
+      transformixFilter.SetMovingImage(movingImage);
+      transformixFilter.SetCombinationTransform(combinationTransform);
+      transformixFilter.SetTransformParameterObject(
+        CreateParameterObject({ // Parameters in alphabetic order:
+                                { "Direction", CreateDefaultDirectionParameterValues<ImageDimension>() },
+                                { "Index", ParameterValuesType(ImageDimension, "0") },
+                                { "Origin", ParameterValuesType(ImageDimension, "0") },
+                                { "ResampleInterpolator", { "FinalLinearInterpolator" } },
+                                { "Size", ConvertToParameterValues(imageSize) },
+                                { "Spacing", ParameterValuesType(ImageDimension, "1") } }));
+      transformixFilter.Update();
+
+      ExpectEqualImages(Deref(transformixFilter.GetOutput()), registrationOutputImage);
+    }
+  }
+}
+
+
 // Tests that Update() throws an exception when the transform parameter object has zero parameter maps.
 GTEST_TEST(itkTransformixFilter, UpdateThrowsExceptionOnZeroParameterMaps)
 {

diff --git a/Core/Main/itkTransformixFilter.h b/Core/Main/itkTransformixFilter.h
@@ -103,6 +103,8 @@ class ITK_TEMPLATE_EXPORT TransformixFilter : public ImageSource<TMovingImage>
 
   using MeshType = Mesh<OutputImagePixelType, MovingImageDimension>;
 
+  using TransformType = Transform<double, MovingImageDimension, MovingImageDimension>;
+
   /** Typedefs for images of determinants of spatial Jacobian matrices, and images of spatial Jacobian matrices */
   using SpatialJacobianDeterminantImageType = itk::Image<float, MovingImageDimension>;
   using SpatialJacobianMatrixImageType =
@@ -225,6 +227,8 @@ class ITK_TEMPLATE_EXPORT TransformixFilter : public ImageSource<TMovingImage>
    * transform object, with additional information from the specified transform parameter object. */
   itkSetConstObjectMacro(Transform, TransformBase);
 
+  itkSetObjectMacro(CombinationTransform, TransformType);
+
   /** Computes the spatial Jacobian determinant for each pixel, and returns an image of the computed values.
   \note Before calling this member function, Update() must be called. */
   SmartPointer<SpatialJacobianDeterminantImageType>
@@ -288,6 +292,8 @@ class ITK_TEMPLATE_EXPORT TransformixFilter : public ImageSource<TMovingImage>
   typename MeshType::Pointer      m_OutputMesh{ nullptr };
 
   TransformBase::ConstPointer m_Transform;
+
+  SmartPointer<TransformType> m_CombinationTransform;
 };
 
 } // namespace itk

diff --git a/Core/Main/itkTransformixFilter.hxx b/Core/Main/itkTransformixFilter.hxx
@@ -265,7 +265,7 @@ TransformixFilter<TMovingImage>::GenerateData()
   unsigned int isError = 0;
   try
   {
-    isError = transformixMain->Run(argumentMap, transformParameterMapVector);
+    isError = transformixMain->Run(argumentMap, transformParameterMapVector, m_CombinationTransform);
 
     if (m_InputMesh)
     {