[Insight-developers] 2-fold speed-up of itkBSplineDerivativeKernelFunction.h
Hans Johnson
hans-johnson at uiowa.edu
Sat Aug 15 19:54:11 EDT 2009
Wenjia,
This is great news! Have you posted this code somewhere that I could get it
and start using it right away (i.e. Before it is introduced into ITK
formally?).
Thanks,
Hans
On 8/15/09 5:38 PM, "Wenjia Bai" <wenjia at robots.ox.ac.uk> wrote:
> Hi Luis,
>
> Sorry for the late reply.
>
> I have just written a program to test the improved
> itkBSplineDerivativeKernelFunction class, which is named
> itkBSplineDerivativeKernelFunction2. The testing program and the
> itkBSplineDerivativeKernelFunction2 class are attached.
>
> Experimental results show that,
>
> (1) The difference between the derivative evaluations of the original
> class and new class is negligible. Therefore, the new class is validated.
>
> B-Spline Order Average Evaluation Difference
> 1 0
> 2 1.79E-017
> 3 2.08E-017
>
> (2) The new class improves the speed by 3-4 times.
>
> Computation Time per Million Evaluations (sec)
> B-Spline Order Original New
> 1 0.1215 0.0393
> 2 0.1210 0.0340
> 3 0.1925 0.0406
>
> Best regards
> Wenjia
>
> Luis Ibanez wrote:
>> Hi Wenjia,
>>
>> This sounds like a great idea.
>>
>> Could you please post this as a patch ?
>> (just send it to the list attached to your email).
>>
>> Also,
>> Have you ran Experimental builds with this change ?
>>
>> Did you find all tests passing ?
>>
>>
>> Please let us know,
>>
>>
>> Luis
>>
>>
>>
>> -------------------
>> wenjia wrote:
>>> 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
>>> _______________________________________________
>>> Powered by www.kitware.com
>>>
>>> Visit other Kitware open-source projects at
>>> http://www.kitware.com/opensource/opensource.html
>>>
>>> Please keep messages on-topic and check the ITK FAQ at:
>>> http://www.itk.org/Wiki/ITK_FAQ
>>>
>>> Follow this link to subscribe/unsubscribe:
>>> http://www.itk.org/mailman/listinfo/insight-developers
>>>
>
> PROJECT(comp_bspline_derivative_kernel)
FIND_PACKAGE(ITK
> REQUIRED)
IF(ITK_FOUND)
> INCLUDE(${ITK_USE_FILE})
ENDIF(ITK_FOUND)
ADD_EXECUTABLE(comp_bspline_derivat
> ive_kernel
> comp_bspline_derivative_kernel.cxx)
TARGET_LINK_LIBRARIES(comp_bspline_derivat
> ive_kernel ITKNumerics)
> #include "itkBSplineDerivativeKernelFunction.h"
> #include "itkBSplineDerivativeKernelFunction2.h"
> #include "itkTimeProbesCollectorBase.h"
>
> int main(int argc, char *argv[])
> {
> // BSpline order
> // Change this manually in testing different B-spline orders
> const unsigned int order = 3;
>
> // Original BSplineDerivativeKernelFunction class
> typedef itk::BSplineDerivativeKernelFunction<order> DerivativeKernelType;
> DerivativeKernelType::Pointer kernel = DerivativeKernelType::New();
>
> // Improved BSplineDerivativeKernelFunction class
> typedef itk::BSplineDerivativeKernelFunction2<order> DerivativeKernelType2;
> DerivativeKernelType2::Pointer kernel2 = DerivativeKernelType2::New();
>
> // Compare the function values in the support region
> double x;
> double xmax = ((double)order + 1) / 2;
> double xmin = - xmax;
> double step = 1E-6;
> double val, val2;
> double diff;
> double err = 0;
> unsigned int n = 0;
>
> for(x=xmin; x<xmax; x+=step)
> {
> val = kernel->Evaluate(x);
> val2 = kernel2->Evaluate(x);
> diff = val2 - val;
> err += fabs(diff);
> n++;
> }
> err /= n;
>
> // Report the difference
> std::cout << "Comparison of BSplineDerivativeKernelFunction and
> BSplineDerivativeKernelFunction2" << std::endl;
> std::cout << " BSpline order = " << order << std::endl;
> std::cout << " Support region = [" << xmin << ", " << xmax << "]" <<
> std::endl;
> std::cout << std::endl;
> std::cout << " The average difference between the evaluations = " << err <<
> std::endl;
> std::cout << std::endl;
>
> // Compare the computation time
> std::cout << " Computation time for evaluating " << n << " times" <<
> std::endl;
>
> // The time probe
> itk::TimeProbesCollectorBase collector;
> collector.Start("Original");
>
> for(x=xmin; x<xmax; x+=step)
> {
> val = kernel->Evaluate(x);
> }
>
> // Report the computation time
> collector.Stop("Original");
> collector.Report();
>
> // The time probe
> itk::TimeProbesCollectorBase collector2;
> collector2.Start("Improved");
>
> for(x=xmin; x<xmax; x+=step)
> {
> val2 = kernel2->Evaluate(x);
> }
>
> // Report the computation time
> collector2.Stop("Improved");
> collector2.Report();
>
> return 0;
> }
> /* itkBSplineDerivativeKernelFunction2.h
> *
> * Use explicit formula to compute the derivative instead of the recursive
> formula
> * so that computation could be faster.
> *
> * Modified from itkBSplineDerivativeKernelFunction.h
> * Modified by: Wenjia Bai
> * Last modified: 2009.08.15
> */
> /*=========================================================================
>
> Program: Insight Segmentation & Registration Toolkit
> Module: $RCSfile: itkBSplineDerivativeKernelFunction2.h,v $
> Language: C++
> Date: $Date: 2008-06-25 11:00:19 $
> Version: $Revision: 1.7 $
>
> Copyright (c) Insight Software Consortium. All rights reserved.
> See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
>
> This software is distributed WITHOUT ANY WARRANTY; without even
> the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
> PURPOSE. See the above copyright notices for more information.
>
> =========================================================================*/
> #ifndef __itkBSplineDerivativeKernelFunction2_h
> #define __itkBSplineDerivativeKernelFunction2_h
>
> #include "itkKernelFunction.h"
> #include "vnl/vnl_math.h"
>
> namespace itk
> {
>
> /** \class BSplineDerivativeKernelFunction
> * \brief Derivative of a BSpline kernel used for density estimation and
> * nonparameteric regression.
> *
> * This class encapsulates the derivative of a BSpline kernel for
> * density estimation or nonparameteric regression.
> * See documentation for KernelFunction for more details.
> *
> * This class is templated over the spline order.
> * \warning Evaluate is only implemented for spline order 0 to 3
> *
> * \sa KernelFunction
> *
> * \ingroup Functions
> */
> template <unsigned int VSplineOrder = 3>
> class ITK_EXPORT BSplineDerivativeKernelFunction2 : public KernelFunction
> {
> public:
> /** Standard class typedefs. */
> typedef BSplineDerivativeKernelFunction2 Self;
> typedef KernelFunction Superclass;
> typedef SmartPointer<Self> Pointer;
>
> /** Method for creation through the object factory. */
> itkNewMacro(Self);
>
> /** Run-time type information (and related methods). */
> itkTypeMacro(BSplineDerivativeKernelFunction2, KernelFunction);
>
> /** Enum of for spline order. */
> itkStaticConstMacro(SplineOrder, unsigned int, VSplineOrder);
>
> /** Evaluate the function. */
> inline double Evaluate( const double & u ) const
> {
> return this->Evaluate( Dispatch<VSplineOrder>(), u );
> }
>
> protected:
>
> BSplineDerivativeKernelFunction2(){};
> ~BSplineDerivativeKernelFunction2(){};
> void PrintSelf(std::ostream& os, Indent indent) const
> {
> Superclass::PrintSelf( os, indent );
> os << indent << "Spline Order: " << SplineOrder << std::endl;
> }
>
> private:
> BSplineDerivativeKernelFunction2(const Self&); //purposely not implemented
> void operator=(const Self&); //purposely not implemented
>
> /** Structures to control overloaded versions of Evaluate */
> struct DispatchBase {};
> template<unsigned int>
> struct Dispatch : DispatchBase{};
>
> /** 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.
> }
>
> };
>
> } // end namespace itk
>
> #endif
> _______________________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-developers
More information about the Insight-developers
mailing list