From 0e0fd1364438060034a6d4bc6a500982642000f4 Mon Sep 17 00:00:00 2001
From: Niels Dekker <N.Dekker@lumc.nl>
Date: Wed, 20 Jan 2021 17:09:20 +0100
Subject: [PATCH] ENH: Add ElastixFilter unit tests UpdateSerially and
 UpdateInParallel

Originally based on a code snippet by Konstantinos Ntatsis at pull request #389 ("ENH: Make ElastixMain database creation + loading components thread safe"):
https://github.com/SuperElastix/elastix/pull/389#issuecomment-759744106
---
 Core/Main/GTesting/ElastixFilterGTest.cxx | 123 ++++++++++++++++++++++
 1 file changed, 123 insertions(+)

diff --git a/Core/Main/GTesting/ElastixFilterGTest.cxx b/Core/Main/GTesting/ElastixFilterGTest.cxx
index 5a988c969..8afc46018 100644
--- a/Core/Main/GTesting/ElastixFilterGTest.cxx
+++ b/Core/Main/GTesting/ElastixFilterGTest.cxx
@@ -31,6 +31,117 @@
 #include <utility> // For pair
 
 
+namespace
+{
+template <typename TImage>
+std::vector<std::string>
+GetTransformParameters(const elx::ElastixFilter<TImage, TImage> & filter)
+{
+  const auto transformParameterObject = filter.GetTransformParameterObject();
+
+  if (transformParameterObject != nullptr)
+  {
+    const auto & transformParameterMaps = transformParameterObject->GetParameterMap();
+
+    if (transformParameterMaps.size() == 1)
+    {
+      const auto & transformParameterMap = transformParameterMaps.front();
+      const auto   found = transformParameterMap.find("TransformParameters");
+
+      if (found != transformParameterMap.cend())
+      {
+        return found->second;
+      }
+    }
+  }
+  return {};
+}
+
+
+void
+UpdateFiltersInForLoop(const bool inParallel)
+{
+  constexpr auto ImageDimension = 2;
+  using ImageType = itk::Image<float, ImageDimension>;
+  using SizeType = itk::Size<ImageDimension>;
+
+  const auto                             imageSizeValue = 8;
+  const itk::ImageRegion<ImageDimension> imageRegion{ itk::Index<ImageDimension>::Filled(imageSizeValue / 2),
+                                                      SizeType::Filled(2) };
+
+  const auto               numberOfImages = 10;
+  ImageType::Pointer       images[numberOfImages];
+  std::vector<std::string> transformParameters[numberOfImages - 1];
+
+  const auto parameterObject = elx::ParameterObject::New();
+  parameterObject->SetParameterMap(
+    elx::ParameterObject::ParameterMapType{ // Parameters in alphabetic order:
+                                            { "ImageSampler", { "Full" } },
+                                            { "MaximumNumberOfIterations", { "2" } },
+                                            { "Metric", { "AdvancedNormalizedCorrelation" } },
+                                            { "Optimizer", { "AdaptiveStochasticGradientDescent" } },
+                                            { "Transform", { "TranslationTransform" } },
+                                            { "WriteFinalTransformParameters", { "false" } } });
+  for (auto & image : images)
+  {
+    image = ImageType::New();
+    image->SetRegions(SizeType::Filled(imageSizeValue));
+    image->Allocate(true);
+
+    const itk::Experimental::ImageRegionRange<ImageType> imageRegionRange{ *image, imageRegion };
+    std::fill(std::begin(imageRegionRange), std::end(imageRegionRange), 1.0f);
+  }
+
+  itk::SmartPointer<elx::ElastixFilter<ImageType, ImageType>> filters[numberOfImages - 1];
+  std::generate_n(filters, numberOfImages - 1, [] { return elx::ElastixFilter<ImageType, ImageType>::New(); });
+
+  for (auto i = 0; i < (numberOfImages - 1); ++i)
+  {
+    const auto & filter = filters[i];
+    filter->LogToConsoleOff();
+    filter->LogToFileOff();
+    filter->SetFixedImage(images[i]);
+    filter->SetMovingImage(images[i + std::size_t{ 1 }]);
+    filter->SetParameterObject(parameterObject);
+  }
+
+  if (inParallel)
+  {
+    // Note: The OpenMP implementation of GCC (including GCC 5.5 and GCC 10.2)
+    // does not seem to support value-initialization of a for-loop index by
+    // empty pair of braces, as in `for (int i{}; i < n; ++i)`.
+#pragma omp parallel for
+    for (auto i = 0; i < (numberOfImages - 1); ++i)
+    {
+      const auto & filter = filters[i];
+      filter->Update();
+      transformParameters[i] = GetTransformParameters(*filter);
+    }
+  }
+  else
+  {
+    for (auto i = 0; i < (numberOfImages - 1); ++i)
+    {
+      const auto & filter = filters[i];
+      filter->Update();
+      transformParameters[i] = GetTransformParameters(*filter);
+    }
+  }
+
+  // Test the resulting TransformParameters after the parallel for loop.
+  for (std::size_t i{}; i < (numberOfImages - 1); ++i)
+  {
+    EXPECT_EQ(transformParameters[i].size(), ImageDimension);
+
+    for (const auto & transformParameter : transformParameters[i])
+    {
+      EXPECT_EQ(transformParameter, "0");
+    }
+  }
+}
+
+} // namespace
+
 // Tests registering two small (5x6) binary images, which are translated with respect to each other.
 GTEST_TEST(ElastixFilter, Translation)
 {
@@ -108,3 +219,15 @@ GTEST_TEST(ElastixFilter, Translation)
     EXPECT_EQ(std::round(std::stod(transformParameters[i])), translationOffset[i]);
   }
 }
+
+
+GTEST_TEST(ElastixFilter, UpdateSerially)
+{
+  UpdateFiltersInForLoop(false);
+}
+
+
+GTEST_TEST(ElastixFilter, UpdateInParallel)
+{
+  UpdateFiltersInForLoop(true);
+}