ENH: Add tests to the module.

CMakeLists.txt

@@ -1,8 +1,6 @@
 cmake_minimum_required(VERSION 3.10.2)
 project(NDReg)
 
-set(NDReg_LIBRARIES NDReg)
-
 if(NOT ITK_SOURCE_DIR)
   find_package(ITK REQUIRED)
   list(APPEND CMAKE_MODULE_PATH ${ITK_CMAKE_DIR})

include/itkMetamorphosisImageRegistrationMethodv4.h

@@ -71,17 +71,15 @@ public ImageToImageFilter<TInputImage, TOutputImage>
   ~ConstantImageFilter(){}
 
   /** Does the real work. */
-  void ThreadedGenerateData(  const typename OutputImageType::RegionType& outputRegionForThread, ThreadIdType threadId)
+  void DynamicThreadedGenerateData(const typename OutputImageType::RegionType& outputRegionForThread) override
   {
-    ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
     OutputPixelType constant = this->GetConstant();
     OutputImageIteratorType  it(this->GetOutput(),outputRegionForThread);
 
     while(!it.IsAtEnd())
     {
       it.Set(constant);
       ++it;
-      progress.CompletedPixel();
     }
   }
 
@@ -225,7 +223,7 @@ public TimeVaryingVelocityFieldImageRegistrationMethodv4<TFixedImage, TMovingIma
   void IntegrateRate();
   FieldPointer GetMetricDerivative(FieldPointer field, bool useImageGradients);
   void UpdateControls();
-  void StartOptimization();
+  void StartOptimization() override;
   void GenerateData() override;
   void PrintSelf(std::ostream& os, Indent indent) const override;
 

include/itkMetamorphosisImageRegistrationMethodv4.hxx

@@ -18,7 +18,6 @@
 #ifndef itkMetamorphosisImageRegistrationMethodv4_hxx
 #define itkMetamorphosisImageRegistrationMethodv4_hxx
 #include "itkMetamorphosisImageRegistrationMethodv4.h"
-#include "vcl_legacy_aliases.h"
 
 namespace itk
 {
@@ -56,11 +55,11 @@ MetamorphosisImageRegistrationMethodv4()
 
   using FixedGradientPixelType = typename ImageMetricType::FixedImageGradientImageType::PixelType;
   m_FixedImageConstantGradientFilter = FixedImageConstantGradientFilterType::New();
-  m_FixedImageConstantGradientFilter->SetConstant(NumericTraits<FixedGradientPixelType>::One);
+  m_FixedImageConstantGradientFilter->SetConstant(NumericTraits<FixedGradientPixelType>::OneValue());
 
   using MovingGradientPixelType = typename ImageMetricType::MovingImageGradientImageType::PixelType;
   m_MovingImageConstantGradientFilter = MovingImageConstantGradientFilterType::New();
-  m_MovingImageConstantGradientFilter->SetConstant(NumericTraits<MovingGradientPixelType>::One);
+  m_MovingImageConstantGradientFilter->SetConstant(NumericTraits<MovingGradientPixelType>::OneValue());
 
   this->SetMetric(MeanSquaresImageToImageMetricv4<FixedImageType, MovingImageType>::New());
 }
@@ -141,19 +140,19 @@ InitializeKernels(TimeVaryingImagePointer kernel, TimeVaryingImagePointer invers
   for(KIt.GoToBegin(), LIt.GoToBegin(); !KIt.IsAtEnd(); ++KIt,++LIt)
   {
     typename TimeVaryingImageType::IndexType  k = KIt.GetIndex();  // Get the frequency index
-    double  A, B;
+    double  A;
 
     // For every dimension accumulate the sum in A
     unsigned int i;
     for(i = 0, A = gamma; i < ImageDimension; i++)
     {
-      A += 2 * alpha * vcl_pow(size[i],2) * ( 1.0-cos(2*vnl_math::pi*k[i]/size[i]) );
-      //A += 2 * alpha * vcl_pow(size[i]/m_Scale,2) * ( 1.0-cos(2*vnl_math::pi*k[i]/size[i]) );
-      //A += 2 * alpha * vcl_pow(spacing[i],-2) * ( 1.0 - cos(2*vnl_math::pi*k[i]*spacing[i]) );
-      //A += alpha * vcl_pow(2*vnl_math::pi*k[i]*spacing[i],2);
+      A += 2 * alpha * std::pow(size[i],2) * ( 1.0-std::cos(2*vnl_math::pi*k[i]/size[i]) );
+      //A += 2 * alpha * std::pow(size[i]/m_Scale,2) * ( 1.0-std::cos(2*vnl_math::pi*k[i]/size[i]) );
+      //A += 2 * alpha * std::pow(spacing[i],-2) * ( 1.0 - std::cos(2*vnl_math::pi*k[i]*spacing[i]) );
+      //A += alpha * std::pow(2*vnl_math::pi*k[i]*spacing[i],2);
     }
 
-    KIt.Set(vcl_pow(A,-2)); // Kernel
+    KIt.Set(std::pow(A,-2)); // Kernel
     LIt.Set(A);             // "Inverse" kernel
   }
 }
@@ -177,7 +176,7 @@ Initialize()
   for(unsigned int i = 0; i < ImageDimension; i++)
   {
     velocityIndex[i] = fixedRegion.GetIndex()[i];
-    velocitySize[i] =  vcl_floor(fixedRegion.GetSize()[i]*m_Scale + 1.0001);
+    velocitySize[i] =  Math::floor(fixedRegion.GetSize()[i]*m_Scale + 1.0001);
     velocityOrigin[i] = fixedImage->GetOrigin()[i];
     velocitySpacing[i] = fixedImage->GetSpacing()[i] / m_Scale;
     for(unsigned int j = 0; j < ImageDimension; j++)
@@ -205,7 +204,7 @@ Initialize()
 
   using PadderType = FFTPadImageFilter<TimeVaryingFieldType>;
   typename PadderType::Pointer padder = PadderType::New();
-  padder->SetSizeGreatestPrimeFactor(vnl_math_min((unsigned int)5, greatestPrimeFactor));
+  padder->SetSizeGreatestPrimeFactor(std::min((unsigned int)5, greatestPrimeFactor));
   padder->SetInput(velocity);
   padder->Update();
   velocity = padder->GetOutput();
@@ -214,7 +213,7 @@ Initialize()
   velocitySize = velocity->GetLargestPossibleRegion().GetSize();
   velocityRegion.SetSize(velocitySize);
   velocity->SetRegions(velocityRegion);
-  velocity->FillBuffer(NumericTraits<VectorType>::Zero);
+  velocity->FillBuffer(NumericTraits<VectorType>::ZeroValue());
 
   // Initialize displacement, /phi_{10}
   this->m_OutputTransform->SetVelocityField(velocity);
@@ -260,13 +259,13 @@ Initialize()
   m_Rate->SetRegions(velocityRegion);
   m_Rate->CopyInformation(velocity);
   m_Rate->Allocate();
-  m_Rate->FillBuffer(NumericTraits<VirtualPixelType>::Zero);
+  m_Rate->FillBuffer(NumericTraits<VirtualPixelType>::ZeroValue());
 
   // Initialize bias, B = 0
   m_Bias->CopyInformation(m_VirtualImage);
   m_Bias->SetRegions(virtualRegion);
   m_Bias->Allocate();
-  m_Bias->FillBuffer(NumericTraits<VirtualPixelType>::Zero);
+  m_Bias->FillBuffer(NumericTraits<VirtualPixelType>::ZeroValue());
 
   // Initialize forward image I(1)
   using MovingCasterType = CastImageFilter<MovingImageType, VirtualImageType>;
@@ -357,8 +356,8 @@ CalculateNorm(TimeVaryingImagePointer image)
   calculator->Update();
 
   // sumOfSquares = (var(x) + mean(x)^2)*length(x)
-  double sumOfSquares = (calculator->GetVariance()+vcl_pow(calculator->GetMean(),2))*image->GetLargestPossibleRegion().GetNumberOfPixels();
-  return vcl_sqrt(sumOfSquares*m_VoxelVolume*m_TimeStep);
+  double sumOfSquares = (calculator->GetVariance()+std::pow(calculator->GetMean(),2))*image->GetLargestPossibleRegion().GetNumberOfPixels();
+  return std::sqrt(sumOfSquares*m_VoxelVolume*m_TimeStep);
 }
 
 template<typename TFixedImage, typename TMovingImage>
@@ -415,7 +414,7 @@ double
 MetamorphosisImageRegistrationMethodv4<TFixedImage, TMovingImage>::
 GetVelocityEnergy()
 {
-  return 0.5 * vcl_pow(CalculateNorm(ApplyKernel(m_InverseVelocityKernel,this->m_OutputTransform->GetVelocityField())),2); // 0.5 ||L_V V||^2
+  return 0.5 * std::pow(CalculateNorm(ApplyKernel(m_InverseVelocityKernel,this->m_OutputTransform->GetVelocityField())),2); // 0.5 ||L_V V||^2
 }
 
 template<typename TFixedImage, typename TMovingImage>
@@ -425,7 +424,7 @@ GetRateEnergy()
 {
   if(m_UseBias)
   {
-    return 0.5 * vcl_pow(m_Mu * CalculateNorm(ApplyKernel(m_InverseRateKernel,m_Rate)),2); // 0.5 \mu^2 ||L_R r||^2
+    return 0.5 * std::pow(m_Mu * CalculateNorm(ApplyKernel(m_InverseRateKernel,m_Rate)),2); // 0.5 \mu^2 ||L_R r||^2
   }
   else
   {
@@ -452,7 +451,7 @@ GetImageEnergy(VirtualImagePointer movingImage, MaskPointer movingMask)
   metric->SetVirtualDomainFromImage(m_VirtualImage);
   metric->Initialize();
 
-  return 0.5*vcl_pow(m_Sigma,-2) * metric->GetValue() * metric->GetNumberOfValidPoints() * m_VoxelVolume;         // 0.5 \sigma^{-2} ||I(1) - I_1||
+  return 0.5*std::pow(m_Sigma,-2) * metric->GetValue() * metric->GetNumberOfValidPoints() * m_VoxelVolume;         // 0.5 \sigma^{-2} ||I(1) - I_1||
 }
 
 template<typename TFixedImage, typename TMovingImage>
@@ -653,7 +652,7 @@ GetMetricDerivative(FieldPointer field, bool useImageGradients)
 
   typename FieldMultiplierType::Pointer multiplier1 = FieldMultiplierType::New();
   multiplier1->SetInput(metricDerivativeField);  // -dM(I(1) o \phi{t1}, I_1 o \phi{t1})
-  multiplier1->SetConstant(vcl_pow(m_Sigma,-2)); // 0.5 \sigma^{-2}
+  multiplier1->SetConstant(std::pow(m_Sigma,-2)); // 0.5 \sigma^{-2}
   multiplier1->Update();
 
   return multiplier1->GetOutput(); // p(t) \nabla I(t) = p(1, \phi{t1})  \nabla I(1, \phi{t1}) = -0.5 \sigma^{-2} -dM(I(1) o \phi{t1}, I_1 o \phi{t1})
@@ -679,7 +678,7 @@ UpdateControls()
     // Compute reverse mapping, \phi_{t1} by integrating velocity field, v(t).
     if(j == m_NumberOfTimeSteps-1)
     {
-      this->m_OutputTransform->GetModifiableDisplacementField()->FillBuffer(NumericTraits<VectorType>::Zero);
+      this->m_OutputTransform->GetModifiableDisplacementField()->FillBuffer(NumericTraits<VectorType>::ZeroValue());
     }
     else
     {
@@ -727,7 +726,7 @@ UpdateControls()
 
     typename TimeVaryingImageMultiplierType::Pointer multiplier1 = TimeVaryingImageMultiplierType::New();
     multiplier1->SetInput(ApplyKernel(m_RateKernel,rateJoiner->GetOutput())); // K_R[p]
-    multiplier1->SetConstant(-vcl_pow(m_Mu,-2));     // -\mu^-2
+    multiplier1->SetConstant(-std::pow(m_Mu,-2));     // -\mu^-2
 
     typename TimeVaryingImageAdderType::Pointer adder1 = TimeVaryingImageAdderType::New();
     adder1->SetInput1(m_Rate);                     // r

include/itkWrapExtrapolateImageFunction.h

@@ -109,7 +109,7 @@ class WrapExtrapolateImageFunction:
   }
 
 
-  void SetInputImage(const InputImageType*ptr)
+  void SetInputImage(const InputImageType*ptr) override
   {
     Superclass::SetInputImage(ptr);
     m_Interpolator->SetInputImage(this->GetInputImage());

test/CMakeLists.txt

@@ -0,0 +1,31 @@
+itk_module_test()
+
+set(NDRegTests
+  itkMetamorphosisImageRegistrationMethodv4Test.cxx
+  #itkTimeVaryingVelocityFieldSemiLagrangianIntegrationImageFilterTest.cxx
+  itkTimeVaryingVelocityFieldSemiLagrangianTransformTest.cxx
+  itkWrapExtrapolateImageFunctionTest.cxx
+  )
+
+CreateTestDriver(NDReg "${NDReg-Test_LIBRARIES}" "${NDRegTests}")
+
+itk_add_test(NAME itkMetamorphosisImageRegistrationMethodv4Test
+      COMMAND NDRegTestDriver
+    itkMetamorphosisImageRegistrationMethodv4Test ${ITK_TEST_OUTPUT_DIR}/itkMyFilterTestOutput.mha
+  )
+
+# Todo: Fix and re-enable
+#itk_add_test(NAME itkTimeVaryingVelocityFieldSemiLagrangianIntegrationImageFilterTest
+      #COMMAND NDRegTestDriver
+    #itkTimeVaryingVelocityFieldSemiLagrangianIntegrationImageFilterTest ${ITK_TEST_OUTPUT_DIR}/itkMyFilterTestOutput.mha
+  #)
+
+itk_add_test(NAME itkTimeVaryingVelocityFieldSemiLagrangianTransformTest
+      COMMAND NDRegTestDriver --without-threads
+    itkTimeVaryingVelocityFieldSemiLagrangianTransformTest ${ITK_TEST_OUTPUT_DIR}/itkNormalDistributionImageSourceTestOutput.mha
+  )
+
+itk_add_test(NAME itkWrapExtrapolateImageFunctionTest
+      COMMAND NDRegTestDriver
+    itkWrapExtrapolateImageFunctionTest ${ITK_TEST_OUTPUT_DIR}/itkMyFilterTestOutput.mha
+  )

test/itkMetamorphosisImageRegistrationMethodv4Test.cxx

@@ -0,0 +1,52 @@
+/*=========================================================================
+ *
+ *  Copyright Insight Software Consortium
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+
+#include "itkMetamorphosisImageRegistrationMethodv4.h"
+#include "itkImageFileWriter.h"
+#include "itkTestingMacros.h"
+
+
+int itkMetamorphosisImageRegistrationMethodv4Test( int argc, char * argv[] )
+{
+  if( argc < 2 )
+    {
+    std::cerr << "Missing parameters." << std::endl;
+    std::cerr << "Usage: " << argv[0];
+    std::cerr << " outputImage";
+    std::cerr << std::endl;
+    return EXIT_FAILURE;
+    }
+
+
+  const char * outputImageFileName  = argv[1];
+
+  const unsigned int Dimension = 2;
+  using PixelType = float;
+  using ImageType = itk::Image< PixelType, Dimension >;
+
+  using MetamorphosisImageRegistrationMethodv4Type = itk::MetamorphosisImageRegistrationMethodv4< ImageType, ImageType >;
+  MetamorphosisImageRegistrationMethodv4Type::Pointer metamorphosisImageRegistration =
+    MetamorphosisImageRegistrationMethodv4Type::New();
+
+  EXERCISE_BASIC_OBJECT_METHODS( metamorphosisImageRegistration, MetamorphosisImageRegistrationMethodv4,
+    TimeVaryingVelocityFieldImageRegistrationMethodv4 );
+
+
+  std::cout << "Test finished." << std::endl;
+  return EXIT_SUCCESS;
+}

test/itkTimeVaryingVelocityFieldSemiLagrangianIntegrationImageFilterTest.cxx

@@ -0,0 +1,55 @@
+/*=========================================================================
+ *
+ *  Copyright Insight Software Consortium
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+
+#include "itkTimeVaryingVelocityFieldSemiLagrangianIntegrationImageFilter.h"
+#include "itkImageFileWriter.h"
+#include "itkTestingMacros.h"
+
+
+int itkTimeVaryingVelocityFieldSemiLagrangianIntegrationImageFilterTest( int argc, char * argv[] )
+{
+  if( argc < 2 )
+    {
+    std::cerr << "Missing parameters." << std::endl;
+    std::cerr << "Usage: " << argv[0];
+    std::cerr << " outputImage";
+    std::cerr << std::endl;
+    return EXIT_FAILURE;
+    }
+
+  constexpr unsigned int Dimension = 3;
+
+  using PixelType = float;
+  using ImageType = itk::Image< PixelType, Dimension >;
+
+  using VectorType = itk::Vector< PixelType, Dimension >;
+  using TimeVaryingVelocityFieldType = itk::Image< VectorType, Dimension >;
+
+  using TimeVaryingVelocityFieldSemiLagrangianIntegrationImageFilterType =
+    itk::TimeVaryingVelocityFieldSemiLagrangianIntegrationImageFilter< TimeVaryingVelocityFieldType >;
+  TimeVaryingVelocityFieldSemiLagrangianIntegrationImageFilterType::Pointer
+    timeVaryingVelocityFieldSemiLagrangianIntegrationImageFilter =
+    TimeVaryingVelocityFieldSemiLagrangianIntegrationImageFilterType::New();
+
+  EXERCISE_BASIC_OBJECT_METHODS( timeVaryingVelocityFieldSemiLagrangianIntegrationImageFilter,
+    TimeVaryingVelocityFieldSemiLagrangianIntegrationImageFilter, TimeVaryingVelocityFieldIntegrationImageFilter );
+
+
+  std::cout << "Test finished." << std::endl;
+  return EXIT_SUCCESS;
+}

test/itkTimeVaryingVelocityFieldSemiLagrangianTransformTest.cxx

@@ -0,0 +1,50 @@
+/*=========================================================================
+ *
+ *  Copyright Insight Software Consortium
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+
+#include "itkTimeVaryingVelocityFieldSemiLagrangianTransform.h"
+#include "itkImageFileWriter.h"
+#include "itkTestingMacros.h"
+
+
+int itkTimeVaryingVelocityFieldSemiLagrangianTransformTest( int argc, char * argv[] )
+{
+  if( argc < 2 )
+    {
+    std::cerr << "Missing parameters." << std::endl;
+    std::cerr << "Usage: " << argv[0];
+    std::cerr << " outputImage";
+    std::cerr << std::endl;
+    return EXIT_FAILURE;
+    }
+
+  constexpr unsigned int Dimension = 2;
+
+  using ParametersValueType = double;
+
+  using TimeVaryingVelocityFieldSemiLagrangianTransformType =
+    itk::TimeVaryingVelocityFieldSemiLagrangianTransform< ParametersValueType, Dimension >;
+  TimeVaryingVelocityFieldSemiLagrangianTransformType::Pointer timeVaryingVelocityFieldSemiLagrangianTransform =
+    TimeVaryingVelocityFieldSemiLagrangianTransformType::New();
+
+  EXERCISE_BASIC_OBJECT_METHODS( timeVaryingVelocityFieldSemiLagrangianTransform,
+    TimeVaryingVelocityFieldSemiLagrangianTransform, TimeVaryingVelocityFieldTransform );
+
+
+  std::cout << "Test finished." << std::endl;
+  return EXIT_SUCCESS;
+}