From f691accd0264dd95018de0cf3e00df785a04416f Mon Sep 17 00:00:00 2001 From: Niels Dekker Date: Fri, 31 Dec 2021 10:24:57 +0100 Subject: [PATCH] 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: https://github.com/InsightSoftwareConsortium/ITK/commit/6a8569ef11cf4f788eeb659b532156479309d7a7 Use the faster `TransformPhysicalPointToContinuousIndex` overload https://github.com/InsightSoftwareConsortium/ITK/commit/eb6ac88c8209ff6cfb3f280bb25e8388f236e211 Use the faster `TransformPhysicalPointToIndex` overload in filter https://github.com/InsightSoftwareConsortium/ITK/commit/0539a2c4ddd2b189d1e48eaf5294ce5556efe732 Remove protected `itk::Transform` data member `m_DirectionChange` https://github.com/InsightSoftwareConsortium/ITK/commit/eec9fe67056bfb5c693a4915b86eac9271e1dc83 Remove unnecessary `IdentityTransform::m_ZeroJacobian` https://github.com/InsightSoftwareConsortium/ITK/commit/9961ccd672cb67ab2128390a49748104eb4459b6 Use FastEvaluate in MattesMutualInformationImageToImageMetric + v4 https://github.com/InsightSoftwareConsortium/ITK/commit/974540982b2926ee6dc4f3933c07962ff74ac9a6 Use FixedArray for table within BSplineInterpolationWeightFunction https://github.com/InsightSoftwareConsortium/ITK/commit/c23944ba54205757c326e1520a4f43d83ac30c4e Remove BSplineInterpolationWeightFunction Kernel, use FastEvaluate https://github.com/InsightSoftwareConsortium/ITK/commit/bc7c5dff2d6a29815bce50cbe7fd8bb5476a4a2b Use FixedArray for BSplineBaseTransform ParameterIndexArrayType https://github.com/InsightSoftwareConsortium/ITK/commit/9bf745bda51044c90208c7b0776e87d7c86adb11 Use FixedArray for BSplineInterpolationWeightFunction OutputType https://github.com/InsightSoftwareConsortium/ITK/commit/c64a58d42a4543eaf1ba20bd1f4c95c7d9658c00 Make `ResampleImageFilter::LinearThreadedGenerateData` faster Release Notes: https://github.com/InsightSoftwareConsortium/ITK/releases/tag/v5.3rc02 Follow-up to pull request https://github.com/SuperElastix/elastix/pull/475 commit 6803b26a73aea270cf7b75e0febc4d53ea21053e "COMP: Upgrade ITK from v5.1.1 to v5.2.0" (merged on 31 May 2021). --- .github/workflows/ElastixGitHubActions.yml | 6 +- CMakeLists.txt | 2 +- .../itkAdvancedBSplineDeformableTransform.hxx | 59 ++++++++----------- ...tkBSplineInterpolationWeightFunctionBase.h | 15 +++-- .../itkCyclicBSplineDeformableTransform.hxx | 4 +- .../itkRecursiveBSplineTransform.hxx | 55 ++++++----------- Testing/CI/Azure/ci.yml | 2 +- ...tkBSplineTransformPointPerformanceTest.cxx | 3 +- 8 files changed, 57 insertions(+), 89 deletions(-) diff --git a/.github/workflows/ElastixGitHubActions.yml b/.github/workflows/ElastixGitHubActions.yml index 0d3436b70..730d9de06 100644 --- a/.github/workflows/ElastixGitHubActions.yml +++ b/.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" diff --git a/CMakeLists.txt b/CMakeLists.txt index 6ffc80c6a..aa769705a 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/Common/Transforms/itkAdvancedBSplineDeformableTransform.hxx b/Common/Transforms/itkAdvancedBSplineDeformableTransform.hxx index c1de92d89..cff0aa7b4 100644 --- a/Common/Transforms/itkAdvancedBSplineDeformableTransform.hxx +++ b/Common/Transforms/itkAdvancedBSplineDeformableTransform.hxx @@ -306,9 +306,8 @@ AdvancedBSplineDeformableTransform::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::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::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::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::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::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::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::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::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::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::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::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::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::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 */ diff --git a/Common/Transforms/itkBSplineInterpolationWeightFunctionBase.h b/Common/Transforms/itkBSplineInterpolationWeightFunctionBase.h index 985b4e874..f0f419cbb 100644 --- a/Common/Transforms/itkBSplineInterpolationWeightFunctionBase.h +++ b/Common/Transforms/itkBSplineInterpolationWeightFunctionBase.h @@ -48,14 +48,17 @@ namespace itk */ template class ITK_TEMPLATE_EXPORT BSplineInterpolationWeightFunctionBase - : public FunctionBase, Array> + : public FunctionBase, + FixedArray> { public: /** Standard class typedefs. */ - typedef BSplineInterpolationWeightFunctionBase Self; - typedef FunctionBase, Array> Superclass; - typedef SmartPointer Pointer; - typedef SmartPointer ConstPointer; + typedef BSplineInterpolationWeightFunctionBase Self; + typedef FunctionBase, + FixedArray> + Superclass; + typedef SmartPointer Pointer; + typedef SmartPointer 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 WeightsType; + typedef FixedArray WeightsType; /** Index and size typedef support. */ typedef Index IndexType; diff --git a/Common/Transforms/itkCyclicBSplineDeformableTransform.hxx b/Common/Transforms/itkCyclicBSplineDeformableTransform.hxx index 858fc6af1..ec9c11e32 100644 --- a/Common/Transforms/itkCyclicBSplineDeformableTransform.hxx +++ b/Common/Transforms/itkCyclicBSplineDeformableTransform.hxx @@ -319,9 +319,7 @@ CyclicBSplineDeformableTransform::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); diff --git a/Common/Transforms/itkRecursiveBSplineTransform.hxx b/Common/Transforms/itkRecursiveBSplineTransform.hxx index 9098047ee..f78c7790b 100644 --- a/Common/Transforms/itkRecursiveBSplineTransform.hxx +++ b/Common/Transforms/itkRecursiveBSplineTransform.hxx @@ -57,8 +57,7 @@ RecursiveBSplineTransform::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::TransformPoint(co /** Call the recursive TransformPoint function. */ ScalarType displacement[SpaceDimension]; RecursiveBSplineTransformImplementation::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::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::GetJacobian( */ ParametersValueType * jacobianPointer = jacobian.data_block(); RecursiveBSplineTransformImplementation::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::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::EvaluateJacobianW } ParametersValueType * imageJacobianPointer = imageJacobian.data_block(); RecursiveBSplineTransformImplementation:: - 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::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::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::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::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]); diff --git a/Testing/CI/Azure/ci.yml b/Testing/CI/Azure/ci.yml index 61b0d8631..576b14e2d 100644 --- a/Testing/CI/Azure/ci.yml +++ b/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 diff --git a/Testing/itkBSplineTransformPointPerformanceTest.cxx b/Testing/itkBSplineTransformPointPerformanceTest.cxx index f7f9740cf..aff4617bb 100644 --- a/Testing/itkBSplineTransformPointPerformanceTest.cxx +++ b/Testing/itkBSplineTransformPointPerformanceTest.cxx @@ -61,9 +61,8 @@ class BSplineTransform_TEST : public AdvancedBSplineDeformableTransform