[Insight-developers] 2-fold speed-up of itkBSplineDerivativeKernelFunction.h

[wenjia]

Dear Developers,

Thank you for your effort in providing ITK, such an excellent library.

I have a suggestion about 2-fold speeding-up of itkBSplineDerivativeKernelFunction.h, which is used to calculate the derivative of a B-spline function. It uses the following equation for computation, calling itkBSplineKernelFunction twice.

  /** Evaluate the function. */
  inline double Evaluate( const double & u ) const
    {
    return ( m_KernelFunction->Evaluate( u + 0.5 ) - 
      m_KernelFunction->Evaluate( u - 0.5 ) );
    }

The code is very clear for reading. However, m_KernelFunction->Evaluate( u + 0.5) and m_KernelFunction->Evaluate( u - 0.5 ) ) waste time in computing similar equations twice. Therefore, expanding the above equation, as follows, can give a two-fold speed up.

  /** Zeroth order spline. */
  inline double Evaluate (const Dispatch<0>&, const double & u) const
    {
      return 0;
    }

  /** First order spline */
  inline double Evaluate ( const Dispatch<1>&, const double& u) const
    {
      if( u == -1 )
      {
        return 0.5;
	}
      else if( (u > -1) && (u < 0) )
      {
        return 1;
	}
      else if( u == 0 )
      {
        return 0;
	}
      else if( (u > 0) && (u < 1) )
      {
        return -1;
	}
      else if( u == 1 )
      {
        return -0.5;
	}
      else
      {
        return 0;
	}
    }

  /** Second order spline. */
  inline double Evaluate ( const Dispatch<2>&, const double& u) const
    {
      if( (u > -0.5) && (u < 0.5) )
      {
        return ( -2 * u );
	}
      else if( (u >= 0.5) && (u < 1.5) )
      {
        return ( -1.5 + u );
	}
      else if( (u > -1.5) && (u <= -0.5))
      {
        return ( 1.5 + u );
	}
      else
      {
        return 0;
	}
    }

  /**  Third order spline. */
  inline double Evaluate ( const Dispatch<3>&, const double& u) const
    {
      if( (u >= 0) && (u < 1) )
      {
        return ( -2 * u + 1.5 * u * u );
	}
      else if( (u > -1 ) && (u < 0) )
      {
        return ( -2 * u - 1.5 * u * u );
	}
      else if( (u >= 1) && (u < 2) )
      {
        return ( -2 + 2 * u - 0.5 * u * u);
	}
      else if( (u > -2) && (u <= -1) )
      {
        return ( 2 + 2 * u + 0.5 * u * u);
	}
      else
      {
        return 0;
	}
    }

  /** Unimplemented spline order */
  inline double Evaluate ( const DispatchBase&, const double& ) const
    {
    itkExceptionMacro("Evaluate not implemented for spline order " << SplineOrder);
    return 0.0; // This is to avoid compiler warning about missing 
                // return statement.  It should never be evaluated.
    }

Considering that itkBSplineDerivativeKernelFunction is usually called inside several levels of loops for time-consuming numerical compuation, speed is more important than brievity of the code. I think expanding the formula is a good choice.

Best regards,
Wenjia Bai