PERF: itk::Image::GetPixel() performance improvment.

itk::Image<>::GetPixel() calls itk::ImageBase<>::ComputeOffset() which
in turn calls itk::ImageBase<>::GetBufferredRegion() which is ... a
virtual function.

AFAIK, this function is overridden only in itk::ImageAdaptator<>, it
doesn't seem to be overridden anywhere else (like in otb::Image<>).  So,
unless the variation point needs to stay in case a child of itk::Image
or otb::Image will require to specialize the behaviour of
GetBufferredRegion(), it'll be best that itk::Image<>::GetPixel() relies
on a non virtual version of GetBufferredRegion().

Why not calling a virtual function from a function which is already
known as inefficient will matters ? Well, first the function will be
less inefficient. Then, some algorithms (like
BSplineInterpolateImageFunction) need to use this GetPixel function. Any
improvement on GetPixel will benefit to these algorithms.

We have benchmarked, with callgrind, the number of cycles spent in the
function BSplineInterpolateImageFunction::EvaluateAtContinousIndex. The
function is called 117290070 times. Before the patch, the function took
212546 million cycles. After the patch, the function takes 133417
million cycles.

EDIT:
Before going further in this direction, we (you?) need to be sure that
there is no plausible scenario where GetBufferredRegion() would be
overridden in a class that inherits from itk::Image.

Then, please note this is a quick and dirty patch. A better
refactorization may be preferred.

Change-Id: Iea7b8da0560e59f9a021fc0a50bfaa57f12939a7

Modules/Core/Common/include/itkImage.h

@@ -160,10 +160,10 @@ class Image:public ImageBase< VImageDimension >
    * typedef typename ImageType::template Rebind< float >::Type OutputImageType;
    *
    */
-  template <typename UPixelType, unsigned int UImageDimension = VImageDimension>
+  template <typename UPixelType, unsigned int NUImageDimension = VImageDimension>
   struct Rebind
     {
-      typedef itk::Image<UPixelType, UImageDimension>  Type;
+      typedef itk::Image<UPixelType, NUImageDimension>  Type;
     };
 
 
@@ -186,7 +186,7 @@ class Image:public ImageBase< VImageDimension >
    * allocated yet. */
   void SetPixel(const IndexType & index, const TPixel & value)
   {
-    OffsetValueType offset = this->ComputeOffset(index);
+    OffsetValueType offset = this->FastComputeOffset(index);
     ( *m_Buffer )[offset] = value;
   }
 
@@ -196,7 +196,7 @@ class Image:public ImageBase< VImageDimension >
    * image has actually been allocated yet. */
   const TPixel & GetPixel(const IndexType & index) const
   {
-    OffsetValueType offset = this->ComputeOffset(index);
+    OffsetValueType offset = this->FastComputeOffset(index);
     return ( ( *m_Buffer )[offset] );
   }
 
@@ -206,7 +206,7 @@ class Image:public ImageBase< VImageDimension >
    * image has actually been allocated yet. */
   TPixel & GetPixel(const IndexType & index)
   {
-    OffsetValueType offset = this->ComputeOffset(index);
+    OffsetValueType offset = this->FastComputeOffset(index);
     return ( ( *m_Buffer )[offset] );
   }
 

Modules/Core/Common/include/itkImageBase.h

@@ -707,6 +707,48 @@ class ImageBase:public DataObject
    * \sa  ReleaseData, Initialize, SetBufferedRegion */
   virtual void InitializeBufferedRegion();
 
+  /** Directly computes an offset from the beginning of the buffer for a pixel
+   * at the specified index.
+   * The index is not checked as to whether it is inside the current buffer, so
+   * the computed offset could conceivably be outside the buffer. If bounds
+   * checking is needed, one can call \c ImageRegion::IsInside(ind) on the
+   * BufferedRegion prior to calling ComputeOffset.
+   * \warning unlike \c ComputeOffset(), this version does not incur a
+   * virtual call. It's meant to be used only for \c itk::Image<>, \c
+   * itk::VectorImage<> and \c itk::SpecialCoordinatesImage<>.
+   */
+  OffsetValueType FastComputeOffset(const IndexType &ind) const
+    {
+    OffsetValueType offset = 0;
+    ImageHelper<VImageDimension,VImageDimension>::ComputeOffset(Self::GetBufferedRegion().GetIndex(),
+                                                                ind,
+                                                                m_OffsetTable,
+                                                                offset);
+    return offset;
+    }
+
+  /** Directly computes the index of the pixel at a specified offset from the
+   * beginning of the buffered region.
+   * Bounds checking is not performed. Thus, the computed index could be
+   * outside the BufferedRegion. To ensure a valid index, the parameter
+   * \c offset should be between 0 and the number of pixels in the
+   * BufferedRegion (the latter can be found using
+   * \c ImageRegion::GetNumberOfPixels()).
+   * \warning unlike \c ComputeOffset(), this version does not incur a
+   * virtual call. It's meant to be used only for \c itk::Image<>, \c
+   * itk::VectorImage<> and \c itk::SpecialCoordinatesImage<>.
+   */
+  IndexType FastComputeIndex(OffsetValueType offset) const
+    {
+    IndexType index;
+    const IndexType &bufferedRegionIndex = Self::GetBufferedRegion().GetIndex();
+    ImageHelper<VImageDimension,VImageDimension>::ComputeIndex(bufferedRegionIndex,
+                                                               offset,
+                                                               m_OffsetTable,
+                                                               index);
+    return index;
+    }
+
 private:
   ImageBase(const Self &);      //purposely not implemented
   void operator=(const Self &); //purposely not implemented

Modules/Core/Common/include/itkSpecialCoordinatesImage.h

@@ -187,7 +187,7 @@ class SpecialCoordinatesImage:public ImageBase< VImageDimension >
    * allocated yet. */
   void SetPixel(const IndexType & index, const TPixel & value)
   {
-    OffsetValueType offset = this->ComputeOffset(index);
+    OffsetValueType offset = this->FastComputeOffset(index);
     ( *m_Buffer )[offset] = value;
   }
 
@@ -197,7 +197,7 @@ class SpecialCoordinatesImage:public ImageBase< VImageDimension >
    * image has actually been allocated yet. */
   const TPixel & GetPixel(const IndexType & index) const
   {
-    OffsetValueType offset = this->ComputeOffset(index);
+    OffsetValueType offset = this->FastComputeOffset(index);
     return ( ( *m_Buffer )[offset] );
   }
 
@@ -207,7 +207,7 @@ class SpecialCoordinatesImage:public ImageBase< VImageDimension >
    * image has actually been allocated yet. */
   TPixel & GetPixel(const IndexType & index)
   {
-    OffsetValueType offset = this->ComputeOffset(index);
+    OffsetValueType offset = this->FastComputeOffset(index);
     return ( ( *m_Buffer )[offset] );
   }
 

Modules/Core/Common/include/itkVectorImage.h

@@ -182,17 +182,17 @@ class VectorImage:
    *
    * \sa Image::Rebind
    */
-  template <typename UPixelType, unsigned int UImageDimension = VImageDimension>
+  template <typename UPixelType, unsigned int NUImageDimension = VImageDimension>
   struct Rebind
   {
-    typedef itk::VectorImage<UPixelType, UImageDimension>  Type;
+    typedef itk::VectorImage<UPixelType, NUImageDimension>  Type;
   };
 
   /** \cond HIDE_SPECIALIZATION_DOCUMENTATION */
-  template <typename UElementType, unsigned int UImageDimension>
-  struct Rebind< VariableLengthVector< UElementType >, UImageDimension>
+  template <typename UElementType, unsigned int NUImageDimension>
+  struct Rebind< VariableLengthVector< UElementType >, NUImageDimension>
   {
-    typedef itk::VectorImage<UElementType, UImageDimension>  Type;
+    typedef itk::VectorImage<UElementType, NUImageDimension>  Type;
   };
   /** \endcond */
 
@@ -215,7 +215,7 @@ class VectorImage:
    * allocated yet. */
   void SetPixel(const IndexType & index, const PixelType & value)
   {
-    OffsetValueType offset = m_VectorLength * this->ComputeOffset(index);
+    OffsetValueType offset = m_VectorLength * this->FastComputeOffset(index);
 
     for ( VectorLengthType i = 0; i < m_VectorLength; i++ )
       {
@@ -230,7 +230,7 @@ class VectorImage:
    * pixel on the stack. */
   const PixelType GetPixel(const IndexType & index) const
   {
-    OffsetValueType offset = m_VectorLength * this->ComputeOffset(index);
+    OffsetValueType offset = m_VectorLength * this->FastComputeOffset(index);
 
     // Do not create a local for this method, to use return value
     // optimization.
@@ -248,7 +248,7 @@ class VectorImage:
    * image has actually been allocated yet. */
   PixelType  GetPixel(const IndexType & index)
   {
-    OffsetValueType offset = m_VectorLength * this->ComputeOffset(index);
+    OffsetValueType offset = m_VectorLength * this->FastComputeOffset(index);
 
     // Correctness of this method relies of return value optimization, do
     // not create a local for the value.