PERF: Re-use jacobian rather than instantiation

A temporary jacobian value was always created for and resized at every iteration of the program call. This was increadibly in-efficient during registrations and was the #1 time consuming portion of the code for many registration examples. The new function: ComputeJacobianWithRespectToParametersCachedTemporaries was created to allow for improved performance. Change-Id: Ic405b9d26c5fdf3660941372e9bb9814c2792143
BRAINSia · May 10, 2014 · 0779906 · 0779906
1 parent 010c924
commit 0779906
Show file tree

Hide file tree

Showing 14 changed files with 73 additions and 31 deletions.
diff --git a/.../Compatibility/V3Compatibility/include/itkIterativeInverseDeformationFieldImageFilter.hxx b/.../Compatibility/V3Compatibility/include/itkIterativeInverseDeformationFieldImageFilter.hxx
@@ -51,7 +51,7 @@ void IterativeInverseDeformationFieldImageFilter< TInputImage, TOutputImage >
     {
     itkExceptionMacro("\n Input is missing.");
     }
-  if ( !TInputImage::ImageDimension == TOutputImage::ImageDimension )
+  if ( ! ( TInputImage::ImageDimension == TOutputImage::ImageDimension ) )
     {
     itkExceptionMacro("\n Image Dimensions must be the same.");
     }

diff --git a/Modules/Core/Transform/include/itkCompositeTransform.h b/Modules/Core/Transform/include/itkCompositeTransform.h
@@ -358,13 +358,13 @@ class CompositeTransform :
    * Compute the Jacobian with respect to the parameters for the compositie
    * transform using Jacobian rule. See comments in the implementation.
    */
-  virtual void ComputeJacobianWithRespectToParameters(const InputPointType  & p, JacobianType & j) const;
+  virtual void ComputeJacobianWithRespectToParameters(const InputPointType  & p, JacobianType & j) const ITK_OVERRIDE;
 
   /**
    * Expanded interface to Compute the Jacobian with respect to the parameters for the compositie
    * transform using Jacobian rule. This version takes in temporary variables to avoid excessive constructions.
    */
-  virtual void ComputeJacobianWithRespectToParametersCachedTemporaries( const InputPointType & p, JacobianType & outJacobian, JacobianType & jacobianWithRespectToPosition ) const;
+  virtual void ComputeJacobianWithRespectToParametersCachedTemporaries( const InputPointType & p, JacobianType & outJacobian, JacobianType & jacobianWithRespectToPosition ) const ITK_OVERRIDE;
 
 protected:
   CompositeTransform();

diff --git a/Modules/Core/Transform/include/itkCompositeTransform.hxx b/Modules/Core/Transform/include/itkCompositeTransform.hxx
@@ -671,7 +671,7 @@ CompositeTransform<TScalar, NDimensions>
    * N cols = total number of parameters in the selected sub transforms. */
   outJacobian.SetSize( NDimensions, this->GetNumberOfLocalParameters() );
   JacobianType jacobianWithRespectToPosition(NDimensions, NDimensions);
-  this->ComputeJacobianWithRespectToParametersCachedTemporaries( const InputPointType & p, JacobianType & outJacobian, JacobianType & jacobianWithRespectToPosition )
+  this->ComputeJacobianWithRespectToParametersCachedTemporaries( p, outJacobian, jacobianWithRespectToPosition );
 }
 
 
@@ -682,13 +682,9 @@ void
 CompositeTransform<TScalar, NDimensions>
 ::ComputeJacobianWithRespectToParametersCachedTemporaries( const InputPointType & p, JacobianType & outJacobian, JacobianType & jacobianWithRespectToPosition ) const
 {
-  /* Returns a concatenated MxN array, holding the Jacobian of each sub
-   * transform that is selected for optimization. The order is the same
-   * as that in which they're applied, i.e. reverse order.
-   * M rows = dimensionality of the transforms
-   * N cols = total number of parameters in the selected sub transforms. */
-  outJacobian.SetSize( NDimensions, this->GetNumberOfLocalParameters() );
-  JacobianType jacobianWithRespectToPosition(NDimensions, NDimensions);
+  //NOTE: This must have been done outside of outJacobian.SetSize( NDimensions, this->GetNumberOfLocalParameters() );
+  //NOTE: assert( outJacobian.GetSize == ( NDimensions, this->GetNumberOfLocalParameters() ) )
+  //NOTE: assert( jacobianWithRespectToPosition.GetSize == (NDimensions, NDimensions) )
 
   NumberOfParametersType offset = NumericTraits< NumberOfParametersType >::Zero;
 

diff --git a/Modules/Core/Transform/include/itkTransform.h b/Modules/Core/Transform/include/itkTransform.h
@@ -496,12 +496,20 @@ class Transform : public TransformBaseTemplate< TScalar >
    *  will most likely occur during multi-threading.
    *  To avoid repeatitive memory allocation, pass in 'jacobian' with its size
    *  already set. */
-  virtual void ComputeJacobianWithRespectToParameters(const InputPointType  & itkNotUsed(p), JacobianType & itkNotUsed(jacobian) ) const
+  virtual void ComputeJacobianWithRespectToParameters(const InputPointType  & itkNotUsed(p), JacobianType & itkNotUsed(jacobian) ) const = 0;
+#if 0
   {
     itkExceptionMacro(
       "ComputeJacobianWithRespectToParamters( InputPointType, JacobianType"
       " is unimplemented for " << this->GetNameOfClass() );
   }
+#endif
+
+  virtual void ComputeJacobianWithRespectToParametersCachedTemporaries(const InputPointType  & p, JacobianType & jacobian, JacobianType & itkNotUsed(jacobianWithRespectToPosition) ) const
+  {
+    //NOTE: default implementation is not optimized, and just falls back to original methods.
+    this->ComputeJacobianWithRespectToParameters(p, jacobian);
+  }
 
 
   /** This provides the ability to get a local jacobian value

diff --git a/Modules/Core/Transform/test/itkMultiTransformTest.cxx b/Modules/Core/Transform/test/itkMultiTransformTest.cxx
@@ -126,12 +126,20 @@ class MultiTransformTestTransform :
   typedef typename Superclass::TransformTypePointer         TransformTypePointer;
   /** InverseTransform type. */
   typedef typename Superclass::InverseTransformBasePointer  InverseTransformBasePointer;
+  typedef typename Superclass::InputPointType             InputPointType;
+  typedef typename Superclass::JacobianType               JacobianType;
 
-  typename Superclass::OutputPointType TransformPoint( const typename Superclass::InputPointType & point ) const
+  typename Superclass::OutputPointType TransformPoint( const InputPointType & point ) const
   {
     return point;
   }
 
+  virtual void ComputeJacobianWithRespectToParameters(const InputPointType  & itkNotUsed(p), JacobianType & itkNotUsed(jacobian) ) const ITK_OVERRIDE
+    {
+    itkExceptionMacro(
+      "ComputeJacobianWithRespectToParamters( InputPointType, JacobianType"
+      " is unimplemented for " << this->GetNameOfClass() );
+    }
 protected:
   MultiTransformTestTransform(){};
   virtual ~MultiTransformTestTransform(){};

diff --git a/...clude/itkANTSNeighborhoodCorrelationImageToImageMetricv4GetValueAndDerivativeThreader.hxx b/...clude/itkANTSNeighborhoodCorrelationImageToImageMetricv4GetValueAndDerivativeThreader.hxx
@@ -503,10 +503,14 @@ ANTSNeighborhoodCorrelationImageToImageMetricv4GetValueAndDerivativeThreader< TD
       }
 
     /* Use a pre-allocated jacobian object for efficiency */
-    JacobianType & jacobian = this->m_GetValueAndDerivativePerThreadVariables[threadId].MovingTransformJacobian;
+    typedef JacobianType & JacobianReferenceType;
+    JacobianReferenceType jacobian = this->m_GetValueAndDerivativePerThreadVariables[threadId].MovingTransformJacobian;
+    JacobianReferenceType jacobianPositional = this->m_GetValueAndDerivativePerThreadVariables[threadId].MovingTransformJacobianPositional;
 
     /** For dense transforms, this returns identity */
-    this->m_Associate->GetMovingTransform()->ComputeJacobianWithRespectToParameters( scanMem.virtualPoint, jacobian );
+    // OLDthis->m_Associate->GetMovingTransform()->ComputeJacobianWithRespectToParameters( scanMem.virtualPoint, jacobian );
+    this->m_Associate->GetMovingTransform()->ComputeJacobianWithRespectToParametersCachedTemporaries(
+      scanMem.virtualPoint, jacobian, jacobianPositional);
 
     NumberOfParametersType numberOfLocalParameters = this->m_Associate->GetMovingTransform()->GetNumberOfLocalParameters();
 

diff --git a/...ion/Metricsv4/include/itkCorrelationImageToImageMetricv4GetValueAndDerivativeThreader.hxx b/...ion/Metricsv4/include/itkCorrelationImageToImageMetricv4GetValueAndDerivativeThreader.hxx
@@ -292,9 +292,12 @@ CorrelationImageToImageMetricv4GetValueAndDerivativeThreader<TDomainPartitioner,
     /* Use a pre-allocated jacobian object for efficiency */
     typedef typename TImageToImageMetric::JacobianType & JacobianReferenceType;
     JacobianReferenceType jacobian = this->m_GetValueAndDerivativePerThreadVariables[threadID].MovingTransformJacobian;
+    JacobianReferenceType jacobianPositional = this->m_GetValueAndDerivativePerThreadVariables[threadID].MovingTransformJacobianPositional;
 
     /** For dense transforms, this returns identity */
-    this->m_CorrelationAssociate->GetMovingTransform()->ComputeJacobianWithRespectToParameters(virtualPoint, jacobian);
+    //OLD this->m_CorrelationAssociate->GetMovingTransform()->ComputeJacobianWithRespectToParameters(virtualPoint, jacobian);
+    this->m_CorrelationAssociate->GetMovingTransform()->ComputeJacobianWithRespectToParametersCachedTemporaries(
+      virtualPoint, jacobian, jacobianPositional);
 
     for (unsigned int par = 0; par < this->m_CorrelationAssociate->GetNumberOfLocalParameters(); par++)
       {

diff --git a/...Registration/Metricsv4/include/itkImageToImageMetricv4GetValueAndDerivativeThreaderBase.h b/...Registration/Metricsv4/include/itkImageToImageMetricv4GetValueAndDerivativeThreaderBase.h
@@ -177,6 +177,7 @@ class ImageToImageMetricv4GetValueAndDerivativeThreaderBase
     /** Pre-allocated transform jacobian objects, for use as needed by dervied
      * classes for efficiency. */
     JacobianType                 MovingTransformJacobian;
+    JacobianType                 MovingTransformJacobianPositional;
     };
   itkPadStruct( ITK_CACHE_LINE_ALIGNMENT, GetValueAndDerivativePerThreadStruct,
                                             PaddedGetValueAndDerivativePerThreadStruct);

diff --git a/...gistration/Metricsv4/include/itkImageToImageMetricv4GetValueAndDerivativeThreaderBase.hxx b/...gistration/Metricsv4/include/itkImageToImageMetricv4GetValueAndDerivativeThreaderBase.hxx
@@ -47,6 +47,10 @@ ImageToImageMetricv4GetValueAndDerivativeThreaderBase< TDomainPartitioner, TImag
 {
   //---------------------------------------------------------------
   // Resize the per thread memory objects.
+  //-----------------------------------------------------------------
+  // Cache some values
+  this->m_CachedNumberOfParameters      = this->m_Associate->GetNumberOfParameters();
+  this->m_CachedNumberOfLocalParameters = this->m_Associate->GetNumberOfLocalParameters();
 
   /* Per-thread results */
   const ThreadIdType numThreadsUsed = this->GetNumberOfThreadsUsed();
@@ -59,8 +63,11 @@ ImageToImageMetricv4GetValueAndDerivativeThreaderBase< TDomainPartitioner, TImag
       {
       /* Allocate intermediary per-thread storage used to get results from
        * derived classes */
-      this->m_GetValueAndDerivativePerThreadVariables[i].LocalDerivatives.SetSize( this->m_Associate->GetNumberOfLocalParameters() );
-      this->m_GetValueAndDerivativePerThreadVariables[i].MovingTransformJacobian.SetSize( this->m_Associate->VirtualImageDimension, this->m_Associate->GetNumberOfLocalParameters() );
+      this->m_GetValueAndDerivativePerThreadVariables[i].LocalDerivatives.SetSize( this->m_CachedNumberOfLocalParameters );
+      this->m_GetValueAndDerivativePerThreadVariables[i].MovingTransformJacobian.SetSize(
+        this->m_Associate->VirtualImageDimension, this->m_CachedNumberOfLocalParameters );
+      this->m_GetValueAndDerivativePerThreadVariables[i].MovingTransformJacobianPositional.SetSize(
+        this->m_Associate->VirtualImageDimension, this->m_Associate->VirtualImageDimension );
       if ( this->m_Associate->m_MovingTransform->GetTransformCategory() == MovingTransformType::DisplacementField )
         {
         /* For transforms with local support, e.g. displacement field,
@@ -76,7 +83,7 @@ ImageToImageMetricv4GetValueAndDerivativeThreaderBase< TDomainPartitioner, TImag
         {
         itkDebugMacro("ImageToImageMetricv4::Initialize: transform does NOT have local support\n");
         /* This size always comes from the moving image */
-        const NumberOfParametersType globalDerivativeSize = this->m_Associate->GetNumberOfParameters();
+        const NumberOfParametersType globalDerivativeSize = this->m_CachedNumberOfParameters;
         /* Global transforms get a separate derivatives container for each thread
          * that holds the result over a particular image region.
          * Use a CompensatedSummation value to provide for better consistency between
@@ -99,18 +106,13 @@ ImageToImageMetricv4GetValueAndDerivativeThreaderBase< TDomainPartitioner, TImag
         /* Be sure to init to 0 here, because the threader may not use
          * all the threads if the region is better split into fewer
          * subregions. */
-        for( NumberOfParametersType p = 0; p < this->m_Associate->GetNumberOfParameters(); p++ )
+        for( NumberOfParametersType p = 0; p < this->m_CachedNumberOfParameters; p++ )
           {
           this->m_GetValueAndDerivativePerThreadVariables[thread].CompensatedDerivatives[p].ResetToZero();
           }
         }
       }
     }
-
-  //-----------------------------------------------------------------
-  // Cache some values
- this->m_CachedNumberOfParameters      = this->m_Associate->GetNumberOfParameters();
- this->m_CachedNumberOfLocalParameters = this->m_Associate->GetNumberOfLocalParameters();
 }
 
 template< typename TDomainPartitioner, typename TImageToImageMetricv4 >

diff --git a/...ation/Metricsv4/include/itkJointHistogramMutualInformationGetValueAndDerivativeThreader.h b/...ation/Metricsv4/include/itkJointHistogramMutualInformationGetValueAndDerivativeThreader.h
@@ -67,7 +67,11 @@ class JointHistogramMutualInformationGetValueAndDerivativeThreader
   typedef typename JointHistogramMetricType::MarginalPDFInterpolatorType    MarginalPDFInterpolatorType;
   typedef typename JointHistogramMetricType::JointPDFInterpolatorPointer    JointPDFInterpolatorPointer;
   typedef typename JointHistogramMetricType::MarginalPDFInterpolatorPointer MarginalPDFInterpolatorPointer;
+  //TODO:  This should just be JacobianType to match the rest of the threaded GetValueAndDerivative types
+  //Deprecate the following
   typedef typename JointHistogramMetricType::FixedTransformJacobianType     FixedTransformJacobianType;
+  //TODO: In the parent class too
+  typedef typename JointHistogramMetricType::FixedTransformJacobianType     JacobianType;
   typedef typename JointHistogramMetricType::NumberOfParametersType         NumberOfParametersType;
   typedef typename JointHistogramMetricType::JointPDFType                   JointPDFType;
   typedef typename JointHistogramMetricType::MarginalPDFType                MarginalPDFType;

diff --git a/...ion/Metricsv4/include/itkJointHistogramMutualInformationGetValueAndDerivativeThreader.hxx b/...ion/Metricsv4/include/itkJointHistogramMutualInformationGetValueAndDerivativeThreader.hxx
@@ -153,11 +153,17 @@ JointHistogramMutualInformationGetValueAndDerivativeThreader< TDomainPartitioner
     }
 
   /* Use a pre-allocated jacobian object for efficiency */
-  FixedTransformJacobianType & jacobian =
-    const_cast< FixedTransformJacobianType &   >(this->m_GetValueAndDerivativePerThreadVariables[threadId].MovingTransformJacobian);
+  typedef JacobianType & JacobianReferenceType;
+  //OLD JacobianReferenceType jacobian =
+  //OLD  const_cast< FixedTransformJacobianType &   >(this->m_GetValueAndDerivativePerThreadVariables[threadId].MovingTransformJacobian);
+  //TODO:  Why const_cast here?
+  JacobianReferenceType jacobian = this->m_GetValueAndDerivativePerThreadVariables[threadId].MovingTransformJacobian;
+  JacobianReferenceType jacobianPositional = this->m_GetValueAndDerivativePerThreadVariables[threadId].MovingTransformJacobianPositional;
 
   /** For dense transforms, this returns identity */
-  this->m_JointAssociate->m_MovingTransform->ComputeJacobianWithRespectToParameters( virtualPoint, jacobian );
+  //OLD this->m_JointAssociate->m_MovingTransform->ComputeJacobianWithRespectToParameters( virtualPoint, jacobian );
+  this->m_JointAssociate->GetMovingTransform()->ComputeJacobianWithRespectToParametersCachedTemporaries(
+      virtualPoint, jacobian, jacobianPositional);
 
   for ( NumberOfParametersType par = 0; par < this->GetCachedNumberOfLocalParameters(); par++ )
     {

diff --git a/...4/include/itkMattesMutualInformationImageToImageMetricv4GetValueAndDerivativeThreader.hxx b/...4/include/itkMattesMutualInformationImageToImageMetricv4GetValueAndDerivativeThreader.hxx
@@ -284,9 +284,13 @@ MattesMutualInformationImageToImageMetricv4GetValueAndDerivativeThreader< TDomai
   // Compute the transform Jacobian.
   typedef JacobianType & JacobianReferenceType;
   JacobianReferenceType jacobian = this->m_GetValueAndDerivativePerThreadVariables[threadID].MovingTransformJacobian;
+  JacobianReferenceType jacobianPositional = this->m_GetValueAndDerivativePerThreadVariables[threadID].MovingTransformJacobianPositional;
   if( this->m_MattesAssociate->GetComputeDerivative() )
     {
-    this->m_MattesAssociate->GetMovingTransform()->ComputeJacobianWithRespectToParameters( virtualPoint, jacobian);
+    //this->m_MattesAssociate->GetMovingTransform()->ComputeJacobianWithRespectToParameters(
+    //  virtualPoint, jacobian);
+    this->m_MattesAssociate->GetMovingTransform()->ComputeJacobianWithRespectToParametersCachedTemporaries(
+      virtualPoint, jacobian, jacobianPositional);
     }
 
   SizeValueType movingParzenBin = 0;

diff --git a/...ion/Metricsv4/include/itkMeanSquaresImageToImageMetricv4GetValueAndDerivativeThreader.hxx b/...ion/Metricsv4/include/itkMeanSquaresImageToImageMetricv4GetValueAndDerivativeThreader.hxx
@@ -58,9 +58,12 @@ MeanSquaresImageToImageMetricv4GetValueAndDerivativeThreader< TDomainPartitioner
   /* Use a pre-allocated jacobian object for efficiency */
   typedef typename TImageToImageMetric::JacobianType & JacobianReferenceType;
   JacobianReferenceType jacobian = this->m_GetValueAndDerivativePerThreadVariables[threadID].MovingTransformJacobian;
+  JacobianReferenceType jacobianPositional = this->m_GetValueAndDerivativePerThreadVariables[threadID].MovingTransformJacobianPositional;
 
   /** For dense transforms, this returns identity */
-  this->m_Associate->GetMovingTransform()->ComputeJacobianWithRespectToParameters( virtualPoint, jacobian );
+  //OLD this->m_Associate->GetMovingTransform()->ComputeJacobianWithRespectToParameters( virtualPoint, jacobian );
+  this->m_Associate->GetMovingTransform()->ComputeJacobianWithRespectToParametersCachedTemporaries(
+      virtualPoint, jacobian, jacobianPositional);
 
   for ( unsigned int par = 0; par < this->GetCachedNumberOfLocalParameters(); par++ )
     {

diff --git a/Modules/Registration/Metricsv4/include/itkPointSetToPointSetMetricv4.hxx b/Modules/Registration/Metricsv4/include/itkPointSetToPointSetMetricv4.hxx
@@ -235,6 +235,7 @@ PointSetToPointSetMetricv4<TFixedPointSet, TMovingPointSet>
 
   value = NumericTraits<MeasureType>::Zero;
   MovingTransformJacobianType  jacobian( MovingPointDimension, this->GetNumberOfLocalParameters() );
+  MovingTransformJacobianType  jacobianPositional( MovingPointDimension, MovingPointDimension );
 
   DerivativeType localTransformDerivative( this->GetNumberOfLocalParameters() );
   localTransformDerivative.Fill( NumericTraits<DerivativeValueType>::Zero );
@@ -286,7 +287,9 @@ PointSetToPointSetMetricv4<TFixedPointSet, TMovingPointSet>
       // Reset to zero since we're not accumulating in the local-support case.
       localTransformDerivative.Fill( NumericTraits<DerivativeValueType>::Zero );
       }
-    this->GetMovingTransform()->ComputeJacobianWithRespectToParameters( virtualIt.Value(), jacobian );
+    //OLD: this->GetMovingTransform()->ComputeJacobianWithRespectToParameters( virtualIt.Value(), jacobian );
+    this->GetMovingTransform()->ComputeJacobianWithRespectToParametersCachedTemporaries(
+      virtualIt.Value(), jacobian, jacobianPositional);
     for ( NumberOfParametersType par = 0; par < this->GetNumberOfLocalParameters(); par++ )
       {
       for( DimensionType d = 0; d < PointDimension; ++d )