ENH: Optimize AdaptiveHistogramEqualization with MovingHistogram base

Implement the filter as a subclass of the MovingHistogramImageFilter,
providing more efficient updates to a histogram while moving in a
line, and multi-threading.

There can be a significant (70x+) performance improvement with larger
radius and multi-threading. There is reasonable scaling with
multiple threads. However with just one thread, and a 1 radius there
almost a 2x improvement.

The boundary condition has changed. Prior the zero flux boundary
condition was used. Now only the image pixels are used, and the
smaller histogram elements are uniformly weighted.

The use of a temporary float image has been remove in favor of scaling on demand.

The UseLookupTable option has been removed/deprecated.

Change-Id: Ib737905a26578878696a1da9ee10906ce9db3b45

Modules/Filtering/ImageStatistics/include/itkAdaptiveEqualizationHistogram.h

@@ -0,0 +1,135 @@
+/*=========================================================================
+ *
+ *  Copyright Insight Software Consortium
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+#ifndef itkAdaptiveEqualizationHistogram_h
+#define itkAdaptiveEqualizationHistogram_h
+
+#include "itksys/hash_map.hxx"
+#include "itkStructHashFunction.h"
+#include "vnl/vnl_math.h"
+#include <cmath>
+namespace itk
+{
+namespace Function
+{
+
+/* \class AdaptiveEqualizationHistogram
+ *
+ * Implements the function class for a moving histogram algorithm for
+ * adaptive histogram equalization.
+ *
+ * \sa AdaptiveHistogramEqualizationImageFilter
+ * \sa MovingHistogramImageFilter
+ * \ingroup ITKImageStatistics
+ */
+template< class TInputPixel, class TOutputPixel >
+class AdaptiveEqualizationHistogram
+{
+public:
+
+  typedef float RealType;
+
+  AdaptiveEqualizationHistogram()
+    : m_BoundaryCount(0)
+    {
+    }
+
+  // ~AdaptiveEqualizationHistogram()  {} default is ok
+
+  void AddPixel(const TInputPixel & p)
+  {
+    m_Map[p]++;
+  }
+
+  void RemovePixel(const TInputPixel & p)
+  {
+
+    // insert new item if one doesn't exist
+    typename MapType::iterator it = m_Map.find( p );
+
+    itkAssertInDebugAndIgnoreInReleaseMacro( it != m_Map.end() );
+
+    if ( --(it->second) == 0 )
+      {
+      m_Map.erase( it );
+      }
+
+  }
+
+  TOutputPixel GetValue(const TInputPixel &pixel)
+    {
+      const RealType iscale = (RealType)m_Maximum - m_Minimum;
+      const RealType scale = 1.0 / iscale;
+
+      RealType sum = 0;
+      typename MapType::iterator itMap = m_Map.begin();
+      const RealType u = scale * ( (RealType)pixel - m_Minimum ) - 0.5;
+      while ( itMap != m_Map.end() )
+        {
+        const RealType v =  scale * ( (RealType)itMap->first - m_Minimum ) - 0.5;
+        const RealType kernel = 1.0 / (m_KernelSize - m_BoundaryCount);
+        sum += kernel * itMap->second * CumulativeFunction(u,v);
+
+        ++itMap;
+        }
+
+      return (TOutputPixel)( iscale * ( sum + 0.5 ) + m_Minimum );
+    }
+
+  void AddBoundary() {++m_BoundaryCount;}
+
+  void RemoveBoundary() {--m_BoundaryCount;}
+
+  void SetAlpha( RealType alpha ) {m_Alpha=alpha;}
+  void SetBeta( RealType beta ) {m_Beta=beta;}
+  void SetKernelSize( RealType kernelSize ) {m_KernelSize=kernelSize;}
+
+  void SetMinimum( TInputPixel minimum ) {m_Minimum=minimum;}
+  void SetMaximum( TInputPixel maximum ) {m_Maximum=maximum;}
+
+private:
+  RealType m_Alpha;
+  RealType m_Beta;
+  RealType m_KernelSize;
+
+  TInputPixel m_Minimum;
+  TInputPixel m_Maximum;
+
+  RealType CumulativeFunction(RealType u, RealType v)
+  {
+    // Calculate cumulative function
+    const RealType s = vnl_math_sgn(u - v);
+    const RealType ad = vnl_math_abs( 2.0 * ( u - v ) );
+
+    return 0.5 * s * std::pow(ad, m_Alpha) - m_Beta * 0.5 * s * ad + m_Beta * u;
+  }
+
+private:
+  typedef typename itksys::hash_map< TInputPixel,
+                                     size_t,
+                                     StructHashFunction< TInputPixel > > MapType;
+
+
+  MapType       m_Map;
+  size_t        m_BoundaryCount;
+
+};
+
+} // end namespace Function
+} // end namespace itk
+
+#endif  // itkAdaptiveHistogramHistogram_h

Modules/Filtering/ImageStatistics/include/itkAdaptiveHistogramEqualizationImageFilter.h

@@ -18,7 +18,8 @@
 #ifndef itkAdaptiveHistogramEqualizationImageFilter_h
 #define itkAdaptiveHistogramEqualizationImageFilter_h
 
-#include "itkBoxImageFilter.h"
+#include "itkMovingHistogramImageFilter.h"
+#include "itkAdaptiveEqualizationHistogram.h"
 #include "itkImage.h"
 
 namespace itk
@@ -51,6 +52,10 @@ namespace itk
  * By altering alpha, beta and window, a host of equalization and unsharp
  * masking filters is available.
  *
+ * The boundary condition ignores the part of the neighborhood
+ * outside the image, and over-weights the valid part of the
+ * neighborhood.
+ *
  * For detail description, reference "Adaptive Image Contrast
  * Enhancement using Generalizations of Histogram Equalization."
  * J.Alex Stark. IEEE Transactions on Image Processing, May 2000.
@@ -62,16 +67,26 @@ namespace itk
  * \wikiexample{NeedDemo/ImageProcessing/AdaptiveHistogramEqualizationImageFilter,Adaptive histogram equalization}
  * \endwiki
  */
-template< typename TImageType >
+template< typename TImageType , typename TKernel = Neighborhood<bool, TImageType::ImageDimension> >
 class AdaptiveHistogramEqualizationImageFilter:
-  public BoxImageFilter< TImageType, TImageType >
+  public MovingHistogramImageFilter< TImageType,
+                                     TImageType,
+                                     TKernel,
+                                     typename Function::AdaptiveEqualizationHistogram< typename TImageType::PixelType,
+                                                                                       typename TImageType::PixelType > >
+
 {
 public:
   /**
    * Standard class typedefs
    */
   typedef AdaptiveHistogramEqualizationImageFilter     Self;
-  typedef ImageToImageFilter< TImageType, TImageType > Superclass;
+  typedef MovingHistogramImageFilter< TImageType,
+                                     TImageType,
+                                     TKernel,
+                                     typename Function::AdaptiveEqualizationHistogram< typename TImageType::PixelType,
+                                                                                       typename TImageType::PixelType > >
+                                                       Superclass;
   typedef SmartPointer< Self >                         Pointer;
   typedef SmartPointer< const Self >                   ConstPointer;
 
@@ -85,8 +100,9 @@ class AdaptiveHistogramEqualizationImageFilter:
   itkTypeMacro(AdaptiveHistogramEqualizationImageFilter, ImageToImageFilter);
 
   /** Image type typedef support. */
-  typedef TImageType                   ImageType;
-  typedef typename ImageType::SizeType ImageSizeType;
+  typedef TImageType                    ImageType;
+  typedef typename ImageType::PixelType InputPixelType;
+  typedef typename ImageType::SizeType  ImageSizeType;
 
   /** Set/Get the value of alpha.  Alpha=0 produces the adaptive
    * histogram equalization (provided beta=0).  Alpha=1 produces an
@@ -101,11 +117,39 @@ class AdaptiveHistogramEqualizationImageFilter:
   itkSetMacro(Beta, float);
   itkGetConstMacro(Beta, float);
 
+#ifndef ITK_FUTURE_LEGACY_REMOVE
   /** Set/Get whether an optimized lookup table for the intensity
-   * mapping function is used.  Default is off. */
-  itkSetMacro(UseLookupTable, bool);
+   * mapping function is used.  Default is off.
+   *  @deprecated
+   */
+  virtual void SetUseLookupTable( const bool _arg )
+    {
+  itkDebugMacro("setting UseLookupTable to " << _arg );
+  itkGenericLegacyReplaceBodyMacro( "UseLookupTable", "", "nothing" );
+  if (this->m_UseLookupTable != _arg)
+    {
+    this->m_UseLookupTable = _arg;
+  this->Modified();
+    }
+  }
   itkGetConstMacro(UseLookupTable, bool);
   itkBooleanMacro(UseLookupTable);
+#endif
+
+  virtual void ConfigureHistogram( typename Superclass::HistogramType &h)
+    {
+      h.SetAlpha( this->m_Alpha );
+      h.SetBeta( this->m_Beta );
+      h.SetMinimum( this->m_InputMinimum );
+      h.SetMaximum( this->m_InputMaximum );
+
+      typename Superclass::HistogramType::RealType kernelSize = 1;
+      for ( unsigned int i = 0; i < ImageDimension; i++ )
+        {
+        kernelSize *= ( 2 * this->GetRadius()[i] + 1 );
+        }
+      h.SetKernelSize(kernelSize);
+    }
 
 protected:
   AdaptiveHistogramEqualizationImageFilter()
@@ -122,7 +166,7 @@ class AdaptiveHistogramEqualizationImageFilter:
   /**
    * Standard pipeline method
    */
-  void GenerateData() ITK_OVERRIDE;
+  void BeforeThreadedGenerateData() ITK_OVERRIDE;
 
 private:
   AdaptiveHistogramEqualizationImageFilter(const Self &); //purposely not
@@ -137,12 +181,13 @@ class AdaptiveHistogramEqualizationImageFilter:
   /** The alpha parameter of the AdaptiveHistogramEqualization. */
   float m_Beta;
 
+  InputPixelType m_InputMinimum;
+  InputPixelType m_InputMaximum;
+
   /** Should we use a lookup table to optimize the use of the
    * intensity mapping function? */
   bool m_UseLookupTable;
 
-  /** A function which is used in GenerateData(). */
-  float CumulativeFunction(float u, float v);
 };
 } // end namespace itk