From 99f04d7df49ae97105e8c7360f26ab5016c1e56d Mon Sep 17 00:00:00 2001
From: Niels Dekker <N.Dekker@lumc.nl>
Date: Sun, 5 Sep 2021 10:50:45 +0200
Subject: [PATCH 1/2] ENH: Add static member function
 `BSplineKernelFunction::FastEvaluate(u)`

Equivalent to the existing `BSplineKernelFunction::Evaluate(u)` member
function. Allows writing faster code, because it is non-virtual, and
because it does not require the creation of a `BSplineKernelFunction`
object.
---
 .../Common/include/itkBSplineKernelFunction.h | 32 ++++++++++++-------
 1 file changed, 20 insertions(+), 12 deletions(-)

diff --git a/Modules/Core/Common/include/itkBSplineKernelFunction.h b/Modules/Core/Common/include/itkBSplineKernelFunction.h
index 73e2e3a971e..430ecee6e30 100644
--- a/Modules/Core/Common/include/itkBSplineKernelFunction.h
+++ b/Modules/Core/Common/include/itkBSplineKernelFunction.h
@@ -60,11 +60,19 @@ class ITK_TEMPLATE_EXPORT BSplineKernelFunction : public KernelFunctionBase<TRea
   /** Enum of for spline order. */
   static constexpr unsigned int SplineOrder = VSplineOrder;
 
+  /** Evaluate the function. Faster than the `Evaluate` member function, because it is static (while `Evaluate` is
+   * virtual). */
+  static TRealValueType
+  FastEvaluate(const TRealValueType u)
+  {
+    return Self::Evaluate(Dispatch<VSplineOrder>(), u);
+  }
+
   /** Evaluate the function. */
   TRealValueType
   Evaluate(const TRealValueType & u) const override
   {
-    return this->Evaluate(Dispatch<VSplineOrder>(), u);
+    return Self::FastEvaluate(u);
   }
 
 protected:
@@ -86,8 +94,8 @@ class ITK_TEMPLATE_EXPORT BSplineKernelFunction : public KernelFunctionBase<TRea
   {};
 
   /** Zeroth order spline. */
-  inline TRealValueType
-  Evaluate(const Dispatch<0> &, const TRealValueType & u) const
+  inline static TRealValueType
+  Evaluate(const Dispatch<0> &, const TRealValueType & u)
   {
     const TRealValueType absValue = itk::Math::abs(u);
     if (absValue < static_cast<TRealValueType>(0.5))
@@ -105,8 +113,8 @@ class ITK_TEMPLATE_EXPORT BSplineKernelFunction : public KernelFunctionBase<TRea
   }
 
   /** First order spline */
-  inline TRealValueType
-  Evaluate(const Dispatch<1> &, const TRealValueType & u) const
+  inline static TRealValueType
+  Evaluate(const Dispatch<1> &, const TRealValueType & u)
   {
     const TRealValueType absValue = itk::Math::abs(u);
     if (absValue < NumericTraits<TRealValueType>::OneValue())
@@ -120,8 +128,8 @@ class ITK_TEMPLATE_EXPORT BSplineKernelFunction : public KernelFunctionBase<TRea
   }
 
   /** Second order spline. */
-  inline TRealValueType
-  Evaluate(const Dispatch<2> &, const TRealValueType & u) const
+  inline static TRealValueType
+  Evaluate(const Dispatch<2> &, const TRealValueType & u)
   {
     const TRealValueType absValue = itk::Math::abs(u);
     if (absValue < static_cast<TRealValueType>(0.5))
@@ -144,8 +152,8 @@ class ITK_TEMPLATE_EXPORT BSplineKernelFunction : public KernelFunctionBase<TRea
   }
 
   /**  Third order spline. */
-  inline TRealValueType
-  Evaluate(const Dispatch<3> &, const TRealValueType & u) const
+  inline static TRealValueType
+  Evaluate(const Dispatch<3> &, const TRealValueType & u)
   {
     const TRealValueType absValue = itk::Math::abs(u);
     if (absValue < NumericTraits<TRealValueType>::OneValue())
@@ -169,10 +177,10 @@ class ITK_TEMPLATE_EXPORT BSplineKernelFunction : public KernelFunctionBase<TRea
   }
 
   /** Unimplemented spline order */
-  inline TRealValueType
-  Evaluate(const DispatchBase &, const TRealValueType &) const
+  inline static TRealValueType
+  Evaluate(const DispatchBase &, const TRealValueType &)
   {
-    itkExceptionMacro("Evaluate not implemented for spline order " << SplineOrder);
+    itkGenericExceptionMacro("Evaluate not implemented for spline order " << SplineOrder);
   }
 };
 } // end namespace itk

From c23944ba54205757c326e1520a4f43d83ac30c4e Mon Sep 17 00:00:00 2001
From: Niels Dekker <N.Dekker@lumc.nl>
Date: Sun, 5 Sep 2021 11:12:51 +0200
Subject: [PATCH 2/2] PERF: Remove BSplineInterpolationWeightFunction Kernel,
 use FastEvaluate

Removed the dynamically allocated `m_Kernel` data member from
`BSplineInterpolationWeightFunction` and replaced the virtual function
call `m_Kernel->Evaluate(x)` by a faster static member function call, to
the newly added `BSplineKernelFunction::FastEvaluate` member function.

May significantly improve the performance of
`BSplineInterpolationWeightFunction`.
---
 .../Common/include/itkBSplineInterpolationWeightFunction.h | 7 -------
 .../include/itkBSplineInterpolationWeightFunction.hxx      | 7 ++-----
 2 files changed, 2 insertions(+), 12 deletions(-)

diff --git a/Modules/Core/Common/include/itkBSplineInterpolationWeightFunction.h b/Modules/Core/Common/include/itkBSplineInterpolationWeightFunction.h
index 3c3a6c10f5d..7930a740917 100644
--- a/Modules/Core/Common/include/itkBSplineInterpolationWeightFunction.h
+++ b/Modules/Core/Common/include/itkBSplineInterpolationWeightFunction.h
@@ -20,7 +20,6 @@
 
 #include "itkFunctionBase.h"
 #include "itkContinuousIndex.h"
-#include "itkBSplineKernelFunction.h"
 #include "itkArray.h"
 #include "itkArray2D.h"
 #include "itkMath.h"
@@ -122,12 +121,6 @@ class ITK_TEMPLATE_EXPORT BSplineInterpolationWeightFunction
 
   /** Table mapping linear offset to indices. */
   TableType m_OffsetToIndexTable;
-
-  /** Interpolation kernel type. */
-  using KernelType = BSplineKernelFunction<Self::SplineOrder>;
-
-  /** Interpolation kernel. */
-  typename KernelType::Pointer m_Kernel;
 };
 } // end namespace itk
 
diff --git a/Modules/Core/Common/include/itkBSplineInterpolationWeightFunction.hxx b/Modules/Core/Common/include/itkBSplineInterpolationWeightFunction.hxx
index 062451ad6ba..17a6c1151be 100644
--- a/Modules/Core/Common/include/itkBSplineInterpolationWeightFunction.hxx
+++ b/Modules/Core/Common/include/itkBSplineInterpolationWeightFunction.hxx
@@ -19,6 +19,7 @@
 #define itkBSplineInterpolationWeightFunction_hxx
 
 #include "itkBSplineInterpolationWeightFunction.h"
+#include "itkBSplineKernelFunction.h"
 #include "itkImage.h"
 #include "itkMatrix.h"
 #include "itkMath.h"
@@ -44,10 +45,6 @@ BSplineInterpolationWeightFunction<TCoordRep, VSpaceDimension, VSplineOrder>::BS
     }
     ++counter;
   }
-
-
-  // Initialize the interpolation kernel
-  m_Kernel = KernelType::New();
 }
 
 /**
@@ -104,7 +101,7 @@ BSplineInterpolationWeightFunction<TCoordRep, VSpaceDimension, VSplineOrder>::Ev
 
     for (k = 0; k <= SplineOrder; ++k)
     {
-      weights1D[j][k] = m_Kernel->Evaluate(x);
+      weights1D[j][k] = BSplineKernelFunction<SplineOrder>::FastEvaluate(x);
       x -= 1.0;
     }
   }