WIP: COMP: Upgrade ITK from v5.2.0 to v5.3rc02

Upgraded to ITK version 5.3 (RC02), which was released on 11 November 2021.

Including upgrade to C++14, and various performance improvements:

InsightSoftwareConsortium/ITK@6a8569e Use the faster `TransformPhysicalPointToContinuousIndex` overload
InsightSoftwareConsortium/ITK@eb6ac88 Use the faster `TransformPhysicalPointToIndex` overload in filter
InsightSoftwareConsortium/ITK@0539a2c Remove protected `itk::Transform` data member `m_DirectionChange`
InsightSoftwareConsortium/ITK@eec9fe6 Remove unnecessary `IdentityTransform::m_ZeroJacobian`
InsightSoftwareConsortium/ITK@9961ccd Use FastEvaluate in MattesMutualInformationImageToImageMetric + v4
InsightSoftwareConsortium/ITK@9745409 Use FixedArray for table within BSplineInterpolationWeightFunction
InsightSoftwareConsortium/ITK@c23944b Remove BSplineInterpolationWeightFunction Kernel, use FastEvaluate
InsightSoftwareConsortium/ITK@bc7c5df Use FixedArray for BSplineBaseTransform ParameterIndexArrayType
InsightSoftwareConsortium/ITK@9bf745b Use FixedArray for BSplineInterpolationWeightFunction OutputType
InsightSoftwareConsortium/ITK@c64a58d Make `ResampleImageFilter::LinearThreadedGenerateData` faster

Release Notes: https://github.com/InsightSoftwareConsortium/ITK/releases/tag/v5.3rc02

Follow-up to pull request #475 commit 6803b26 "COMP: Upgrade ITK from v5.1.1 to v5.2.0" (merged on 31 May 2021).

.github/workflows/ElastixGitHubActions.yml

@@ -18,21 +18,21 @@ jobs:
           - os: ubuntu-18.04
             c-compiler: "gcc"
             cxx-compiler: "g++"
-            itk-git-tag: "v5.2.0"
+            itk-git-tag: "v5.3rc02"
             cmake-build-type: "Release"
             ANNLib: "libANNlib-5.0.so"
             ANNLib2: "libANNlib-5.0.so.1"
           - os: windows-2019
             c-compiler: "cl.exe"
             cxx-compiler: "cl.exe"
-            itk-git-tag: "v5.2.0"
+            itk-git-tag: "v5.3rc02"
             cmake-build-type: "Release"
             ANNLib: "ANNlib-5.0.dll"
             vcvars64: "C:/Program Files (x86)/Microsoft Visual Studio/2019/Enterprise/VC/Auxiliary/Build/vcvars64.bat"
           - os: macos-10.15
             c-compiler: "clang"  
             cxx-compiler: "clang++"
-            itk-git-tag: "v5.2.0"
+            itk-git-tag: "v5.3rc02"
             cmake-build-type: "Release"
             ANNLib: "libANNlib-5.0.1.dylib"
             ANNLib2: "libANNlib-5.0.dylib"

CMakeLists.txt

@@ -29,7 +29,7 @@ if(ELASTIX_USE_OPENCL)
 endif()
 
 # Find ITK.
-find_package( ITK 5.2.0 REQUIRED COMPONENTS
+find_package( ITK 5.3 REQUIRED COMPONENTS
     ITKCommon
     ITKDisplacementField
     ITKDistanceMap

Common/Transforms/itkAdvancedBSplineDeformableTransform.hxx

@@ -306,9 +306,8 @@ AdvancedBSplineDeformableTransform<TScalarType, NDimensions, VSplineOrder>::Tran
 {
   /** Allocate memory on the stack: */
   const unsigned long                         numberOfWeights = WeightsFunctionType::NumberOfWeights;
-  typename WeightsType::ValueType             weightsArray[numberOfWeights];
   typename ParameterIndexArrayType::ValueType indicesArray[numberOfWeights];
-  WeightsType                                 weights(weightsArray, numberOfWeights, false);
+  WeightsType                                 weights;
   ParameterIndexArrayType                     indices(indicesArray, numberOfWeights, false);
 
   OutputPointType outputPoint;
@@ -394,9 +393,8 @@ AdvancedBSplineDeformableTransform<TScalarType, NDimensions, VSplineOrder>::GetJ
   /** Compute the number of affected B-spline parameters.
    * Allocate memory on the stack.
    */
-  const unsigned long             numberOfWeights = WeightsFunctionType::NumberOfWeights;
-  typename WeightsType::ValueType weightsArray[numberOfWeights];
-  WeightsType                     weights(weightsArray, numberOfWeights, false);
+  const unsigned long numberOfWeights = WeightsFunctionType::NumberOfWeights;
+  WeightsType         weights;
 
   /** Compute the weights. */
   IndexType supportIndex;
@@ -413,7 +411,7 @@ AdvancedBSplineDeformableTransform<TScalarType, NDimensions, VSplineOrder>::GetJ
   for (unsigned int d = 0; d < SpaceDimension; ++d)
   {
     unsigned long offset = d * SpaceDimension * numberOfWeights + d * numberOfWeights;
-    std::copy(weightsArray, weightsArray + numberOfWeights, jacobianPointer + offset);
+    std::copy_n(weights.data(), numberOfWeights, jacobianPointer + offset);
   }
 
   /** Compute the nonzero Jacobian indices.
@@ -463,9 +461,8 @@ AdvancedBSplineDeformableTransform<TScalarType, NDimensions, VSplineOrder>::Eval
   /** Compute the number of affected B-spline parameters.
    * Allocate memory on the stack.
    */
-  const unsigned long             numberOfWeights = WeightsFunctionType::NumberOfWeights;
-  typename WeightsType::ValueType weightsArray[numberOfWeights];
-  WeightsType                     weights(weightsArray, numberOfWeights, false);
+  const unsigned long numberOfWeights = WeightsFunctionType::NumberOfWeights;
+  WeightsType         weights;
 
   /** Compute the B-spline derivative weights. */
   IndexType supportIndex;
@@ -479,7 +476,7 @@ AdvancedBSplineDeformableTransform<TScalarType, NDimensions, VSplineOrder>::Eval
     const MovingImageGradientValueType mig = movingImageGradient[d];
     for (NumberOfParametersType i = 0; i < nnzjiPerDimension; ++i)
     {
-      imageJacobian[counter] = weightsArray[i] * mig;
+      imageJacobian[counter] = weights[i] * mig;
       ++counter;
     }
   }
@@ -523,13 +520,11 @@ AdvancedBSplineDeformableTransform<TScalarType, NDimensions, VSplineOrder>::GetS
 
   /** Compute the number of affected B-spline parameters. */
   /** Allocate memory on the stack: */
-  const unsigned long             numberOfWeights = WeightsFunctionType::NumberOfWeights;
-  typename WeightsType::ValueType weightsArray[numberOfWeights];
-  WeightsType                     weights(weightsArray, numberOfWeights, false);
+  const unsigned long numberOfWeights = WeightsFunctionType::NumberOfWeights;
+  WeightsType         weights;
 
   /** Array for CoefficientImage values */
-  typename WeightsType::ValueType coeffArray[numberOfWeights * SpaceDimension];
-  WeightsType                     coeffs(coeffArray, numberOfWeights * SpaceDimension, false);
+  WeightsType coeffs;
 
   IndexType supportIndex;
   this->m_DerivativeWeightsFunctions[0]->ComputeStartIndex(cindex, supportIndex);
@@ -633,12 +628,10 @@ AdvancedBSplineDeformableTransform<TScalarType, NDimensions, VSplineOrder>::GetS
   /** Helper variables. */
   /** Allocate memory on the stack: */
   const unsigned long numberOfWeights = WeightsFunctionType::NumberOfWeights;
-  WeightsValueType    weightsArray[numberOfWeights];
-  WeightsType         weights(weightsArray, numberOfWeights, false);
+  WeightsType         weights;
 
   /** Array for CoefficientImage values */
-  WeightsValueType coeffArray[numberOfWeights * SpaceDimension];
-  WeightsType      coeffs(coeffArray, numberOfWeights * SpaceDimension, false);
+  WeightsType coeffs;
 
   IndexType supportIndex;
   this->m_SODerivativeWeightsFunctions[0][0]->ComputeStartIndex(cindex, supportIndex);
@@ -761,9 +754,8 @@ AdvancedBSplineDeformableTransform<TScalarType, NDimensions, VSplineOrder>::GetJ
   /** Helper variables. */
 
   /** Allocate memory on the stack: */
-  const unsigned long             numberOfWeights = WeightsFunctionType::NumberOfWeights;
-  typename WeightsType::ValueType weightsArray[numberOfWeights];
-  WeightsType                     weights(weightsArray, numberOfWeights, false);
+  const unsigned long numberOfWeights = WeightsFunctionType::NumberOfWeights;
+  WeightsType         weights;
 
   IndexType supportIndex;
   this->m_DerivativeWeightsFunctions[0]->ComputeStartIndex(cindex, supportIndex);
@@ -784,7 +776,7 @@ AdvancedBSplineDeformableTransform<TScalarType, NDimensions, VSplineOrder>::GetJ
     this->m_DerivativeWeightsFunctions[i]->Evaluate(cindex, supportIndex, weights);
 
     /** Remember the weights. */
-    std::copy(weights.data_block(), weights.data_block() + numberOfWeights, weightVector + i * numberOfWeights);
+    std::copy_n(weights.data(), numberOfWeights, weightVector + i * numberOfWeights);
 
   } // end for i
 
@@ -870,12 +862,10 @@ AdvancedBSplineDeformableTransform<TScalarType, NDimensions, VSplineOrder>::GetJ
   /** Allocate weight on the stack. */
   typedef typename WeightsType::ValueType WeightsValueType;
   const unsigned long                     numberOfWeights = WeightsFunctionType::NumberOfWeights;
-  WeightsValueType                        weightsArray[numberOfWeights];
-  WeightsType                             weights(weightsArray, numberOfWeights, false);
+  WeightsType                             weights;
 
   /** Allocate coefficients on the stack. */
-  WeightsValueType coeffArray[numberOfWeights * SpaceDimension];
-  WeightsType      coeffs(coeffArray, numberOfWeights * SpaceDimension, false);
+  WeightsType coeffs;
 
   /** Copy values from coefficient image to linear coeffs array. */
   // takes considerable amount of time : 27% of this function. // with old region iterator, check with new
@@ -918,7 +908,7 @@ AdvancedBSplineDeformableTransform<TScalarType, NDimensions, VSplineOrder>::GetJ
      * weights at once for all dimensions */
 
     /** Remember the weights. */
-    std::copy(weights.data_block(), weights.data_block() + numberOfWeights, weightVector + i * numberOfWeights);
+    std::copy_n(weights.data(), numberOfWeights, weightVector + i * numberOfWeights);
 
     /** Reset coeffs iterator */
     itCoeffs = coeffs.begin();
@@ -1027,9 +1017,8 @@ AdvancedBSplineDeformableTransform<TScalarType, NDimensions, VSplineOrder>::GetJ
   /** Compute the number of affected B-spline parameters. */
 
   /** Allocate memory on the stack: */
-  const unsigned long             numberOfWeights = WeightsFunctionType::NumberOfWeights;
-  typename WeightsType::ValueType weightsArray[numberOfWeights];
-  WeightsType                     weights(weightsArray, numberOfWeights, false);
+  const unsigned long numberOfWeights = WeightsFunctionType::NumberOfWeights;
+  WeightsType         weights;
 
   IndexType supportIndex;
   this->m_SODerivativeWeightsFunctions[0][0]->ComputeStartIndex(cindex, supportIndex);
@@ -1157,12 +1146,10 @@ AdvancedBSplineDeformableTransform<TScalarType, NDimensions, VSplineOrder>::GetJ
   /** Allocate weight on the stack. */
   typedef typename WeightsType::ValueType WeightsValueType;
   const unsigned long                     numberOfWeights = WeightsFunctionType::NumberOfWeights;
-  WeightsValueType                        weightsArray[numberOfWeights];
-  WeightsType                             weights(weightsArray, numberOfWeights, false);
+  WeightsType                             weights;
 
   /** Allocate coefficients on the stack. */
-  WeightsValueType coeffArray[numberOfWeights * SpaceDimension];
-  WeightsType      coeffs(coeffArray, numberOfWeights * SpaceDimension, false);
+  WeightsType coeffs;
 
   /** Copy values from coefficient image to linear coeffs array. */
   // takes considerable amount of time : 27% of this function. // with old region iterator, check with new
@@ -1205,7 +1192,7 @@ AdvancedBSplineDeformableTransform<TScalarType, NDimensions, VSplineOrder>::GetJ
       this->m_SODerivativeWeightsFunctions[i][j]->Evaluate(cindex, supportIndex, weights);
 
       /** Remember the weights. */
-      std::copy(weights.data_block(), weights.data_block() + numberOfWeights, weightVector + count * numberOfWeights);
+      std::copy_n(weights.data(), numberOfWeights, weightVector + count * numberOfWeights);
       count++;
 
       /** Reset coeffs iterator */

Common/Transforms/itkBSplineInterpolationWeightFunctionBase.h

@@ -48,14 +48,17 @@ namespace itk
  */
 template <class TCoordRep = float, unsigned int VSpaceDimension = 2, unsigned int VSplineOrder = 3>
 class ITK_TEMPLATE_EXPORT BSplineInterpolationWeightFunctionBase
-  : public FunctionBase<ContinuousIndex<TCoordRep, VSpaceDimension>, Array<double>>
+  : public FunctionBase<ContinuousIndex<TCoordRep, VSpaceDimension>,
+                        FixedArray<double, Math::UnsignedPower(VSplineOrder + 1, VSpaceDimension)>>
 {
 public:
   /** Standard class typedefs. */
-  typedef BSplineInterpolationWeightFunctionBase                                   Self;
-  typedef FunctionBase<ContinuousIndex<TCoordRep, VSpaceDimension>, Array<double>> Superclass;
-  typedef SmartPointer<Self>                                                       Pointer;
-  typedef SmartPointer<const Self>                                                 ConstPointer;
+  typedef BSplineInterpolationWeightFunctionBase Self;
+  typedef FunctionBase<ContinuousIndex<TCoordRep, VSpaceDimension>,
+                       FixedArray<double, Math::UnsignedPower(VSplineOrder + 1, VSpaceDimension)>>
+                                   Superclass;
+  typedef SmartPointer<Self>       Pointer;
+  typedef SmartPointer<const Self> ConstPointer;
 
   /** Run-time type information (and related methods). */
   itkTypeMacro(BSplineInterpolationWeightFunctionBase, FunctionBase);
@@ -70,7 +73,7 @@ class ITK_TEMPLATE_EXPORT BSplineInterpolationWeightFunctionBase
   static constexpr unsigned long NumberOfWeights = Math::UnsignedPower(VSplineOrder + 1, VSpaceDimension);
 
   /** OutputType typedef support. */
-  typedef Array<double> WeightsType;
+  typedef FixedArray<double, NumberOfWeights> WeightsType;
 
   /** Index and size typedef support. */
   typedef Index<VSpaceDimension> IndexType;

Common/Transforms/itkCyclicBSplineDeformableTransform.hxx

@@ -319,9 +319,7 @@ CyclicBSplineDeformableTransform<TScalarType, NDimensions, VSplineOrder>::GetSpa
 
   /** Compute the number of affected B-spline parameters. */
   /** Allocate memory on the stack: */
-  const SizeValueType             numberOfWeights = WeightsFunctionType::NumberOfWeights;
-  typename WeightsType::ValueType weightsArray[numberOfWeights];
-  WeightsType                     weights(weightsArray, numberOfWeights, false);
+  WeightsType weights;
 
   IndexType supportIndex;
   this->m_DerivativeWeightsFunctions[0]->ComputeStartIndex(cindex, supportIndex);

Common/Transforms/itkRecursiveBSplineTransform.hxx

@@ -57,8 +57,7 @@ RecursiveBSplineTransform<TScalar, NDimensions, VSplineOrder>::TransformPoint(co
   OutputPointType outputPoint;
 
   /** Allocate weights on the stack: */
-  typename WeightsType::ValueType weightsArray1D[numberOfWeights];
-  WeightsType                     weights1D(weightsArray1D, numberOfWeights, false);
+  WeightsType weights1D;
 
   /** Check if the coefficient image has been set. */
   if (!this->m_CoefficientImages[0])
@@ -102,7 +101,7 @@ RecursiveBSplineTransform<TScalar, NDimensions, VSplineOrder>::TransformPoint(co
   /** Call the recursive TransformPoint function. */
   ScalarType displacement[SpaceDimension];
   RecursiveBSplineTransformImplementation<SpaceDimension, SpaceDimension, SplineOrder, TScalar>::TransformPoint(
-    displacement, mu, bsplineOffsetTable, weightsArray1D);
+    displacement, mu, bsplineOffsetTable, weights1D.data());
 
   // The output point is the start point + displacement.
   for (unsigned int j = 0; j < SpaceDimension; ++j)
@@ -156,10 +155,8 @@ RecursiveBSplineTransform<TScalar, NDimensions, VSplineOrder>::GetJacobian(
    * In contrast to the normal B-spline weights function, the recursive version
    * returns the individual weights instead of the multiplied ones.
    */
-  const unsigned int              numberOfWeights = RecursiveBSplineWeightFunctionType::NumberOfWeights;
-  typename WeightsType::ValueType weightsArray1D[numberOfWeights];
-  WeightsType                     weights1D(weightsArray1D, numberOfWeights, false);
-  IndexType                       supportIndex;
+  WeightsType weights1D;
+  IndexType   supportIndex;
   this->m_RecursiveBSplineWeightFunction->Evaluate(cindex, weights1D, supportIndex);
 
   /** Recursively compute the first numberOfIndices entries of the Jacobian.
@@ -168,7 +165,7 @@ RecursiveBSplineTransform<TScalar, NDimensions, VSplineOrder>::GetJacobian(
    */
   ParametersValueType * jacobianPointer = jacobian.data_block();
   RecursiveBSplineTransformImplementation<SpaceDimension, SpaceDimension, SplineOrder, TScalar>::GetJacobian(
-    jacobianPointer, weightsArray1D, 1.0);
+    jacobianPointer, weights1D.data(), 1.0);
 
   /** Compute the nonzero Jacobian indices.
    * Takes a significant portion of the computation time of this function.
@@ -217,10 +214,8 @@ RecursiveBSplineTransform<TScalar, NDimensions, VSplineOrder>::EvaluateJacobianW
    * In contrast to the normal B-spline weights function, the recursive version
    * returns the individual weights instead of the multiplied ones.
    */
-  const unsigned int              numberOfWeights = RecursiveBSplineWeightFunctionType::NumberOfWeights;
-  typename WeightsType::ValueType weightsArray1D[numberOfWeights];
-  WeightsType                     weights1D(weightsArray1D, numberOfWeights, false);
-  IndexType                       supportIndex;
+  WeightsType weights1D;
+  IndexType   supportIndex;
   this->m_RecursiveBSplineWeightFunction->Evaluate(cindex, weights1D, supportIndex);
 
   /** Recursively compute the inner product of the Jacobian and the moving image gradient.
@@ -234,7 +229,7 @@ RecursiveBSplineTransform<TScalar, NDimensions, VSplineOrder>::EvaluateJacobianW
   }
   ParametersValueType * imageJacobianPointer = imageJacobian.data_block();
   RecursiveBSplineTransformImplementation<SpaceDimension, SpaceDimension, SplineOrder, TScalar>::
-    EvaluateJacobianWithImageGradientProduct(imageJacobianPointer, migArray, weightsArray1D, 1.0);
+    EvaluateJacobianWithImageGradientProduct(imageJacobianPointer, migArray, weights1D.data(), 1.0);
 
   /** Setup support region needed for the nonZeroJacobianIndices. */
   RegionType supportRegion;
@@ -273,11 +268,8 @@ RecursiveBSplineTransform<TScalar, NDimensions, VSplineOrder>::GetSpatialJacobia
   }
 
   /** Create storage for the B-spline interpolation weights. */
-  const unsigned int              numberOfWeights = RecursiveBSplineWeightFunctionType::NumberOfWeights;
-  typename WeightsType::ValueType weightsArray1D[numberOfWeights];
-  WeightsType                     weights1D(weightsArray1D, numberOfWeights, false);
-  typename WeightsType::ValueType derivativeWeightsArray1D[numberOfWeights];
-  WeightsType                     derivativeWeights1D(derivativeWeightsArray1D, numberOfWeights, false);
+  WeightsType weights1D;
+  WeightsType derivativeWeights1D;
 
   double * weightsPointer = &(weights1D[0]);
   double * derivativeWeightsPointer = &(derivativeWeights1D[0]);
@@ -361,13 +353,9 @@ RecursiveBSplineTransform<TScalar, NDimensions, VSplineOrder>::GetSpatialHessian
   }
 
   /** Create storage for the B-spline interpolation weights. */
-  const unsigned int              numberOfWeights = RecursiveBSplineWeightFunctionType::NumberOfWeights;
-  typename WeightsType::ValueType weightsArray1D[numberOfWeights];
-  WeightsType                     weights1D(weightsArray1D, numberOfWeights, false);
-  typename WeightsType::ValueType derivativeWeightsArray1D[numberOfWeights];
-  WeightsType                     derivativeWeights1D(derivativeWeightsArray1D, numberOfWeights, false);
-  typename WeightsType::ValueType hessianWeightsArray1D[numberOfWeights];
-  WeightsType                     hessianWeights1D(hessianWeightsArray1D, numberOfWeights, false);
+  WeightsType weights1D;
+  WeightsType derivativeWeights1D;
+  WeightsType hessianWeights1D;
 
   double * weightsPointer = &(weights1D[0]);
   double * derivativeWeightsPointer = &(derivativeWeights1D[0]);
@@ -481,11 +469,8 @@ RecursiveBSplineTransform<TScalar, NDimensions, VSplineOrder>::GetJacobianOfSpat
   }
 
   /** Create storage for the B-spline interpolation weights. */
-  const unsigned int              numberOfWeights = RecursiveBSplineWeightFunctionType::NumberOfWeights;
-  typename WeightsType::ValueType weightsArray1D[numberOfWeights];
-  WeightsType                     weights1D(weightsArray1D, numberOfWeights, false);
-  typename WeightsType::ValueType derivativeWeightsArray1D[numberOfWeights];
-  WeightsType                     derivativeWeights1D(derivativeWeightsArray1D, numberOfWeights, false);
+  WeightsType weights1D;
+  WeightsType derivativeWeights1D;
 
   double * weightsPointer = &(weights1D[0]);
   double * derivativeWeightsPointer = &(derivativeWeights1D[0]);
@@ -583,13 +568,9 @@ RecursiveBSplineTransform<TScalar, NDimensions, VSplineOrder>::GetJacobianOfSpat
   }
 
   /** Create storage for the B-spline interpolation weights. */
-  const unsigned int              numberOfWeights = RecursiveBSplineWeightFunctionType::NumberOfWeights;
-  typename WeightsType::ValueType weightsArray1D[numberOfWeights];
-  WeightsType                     weights1D(weightsArray1D, numberOfWeights, false);
-  typename WeightsType::ValueType derivativeWeightsArray1D[numberOfWeights];
-  WeightsType                     derivativeWeights1D(derivativeWeightsArray1D, numberOfWeights, false);
-  typename WeightsType::ValueType hessianWeightsArray1D[numberOfWeights];
-  WeightsType                     hessianWeights1D(hessianWeightsArray1D, numberOfWeights, false);
+  WeightsType weights1D;
+  WeightsType derivativeWeights1D;
+  WeightsType hessianWeights1D;
 
   double * weightsPointer = &(weights1D[0]);
   double * derivativeWeightsPointer = &(derivativeWeights1D[0]);

Testing/CI/Azure/ci.yml

@@ -1,5 +1,5 @@
 variables:
-  ITKv5_VERSION: v5.2.0
+  ITKv5_VERSION: v5.3rc02
   ITK_GIT_URL: https://github.com/InsightSoftwareConsortium/ITK
   ITK_SOURCE_DIR: $(Agent.BuildDirectory)/ITK-source
   ITK_BINARY_DIR: $(Agent.BuildDirectory)/ITK-build