PERF: Using stl stable_sort instead of a custom class and sort

The algorithm requires a stable sort of image indexes based
on pixel intensity. The paper recommends counting sort but
that has its own problems when using bigger than 8 bit pixel
types or floats - i.e. need additional scans of the data
to figure out range etc. The original implementation
of this class used a structure with both position and brightness.
This version abandons that structure and uses an offset array
only, saving some RAM and simplifying the code.

There's an opportunity to activate stl parallel sorting in the
future.

Modules/Nonunit/Review/include/itkAttributeMorphologyBaseImageFilter.h

@@ -165,7 +165,6 @@ class ITK_TEMPLATE_EXPORT AttributeMorphologyBaseImageFilter : public ImageToIma
   // some constants used several times in the code
   static constexpr OffsetValueType INACTIVE = -1;
   static constexpr OffsetValueType ACTIVE = -2;
-  static constexpr OffsetValueType ROOT = -3;
 
   // Just used for area/volume openings at the moment
   AttributeType * m_AuxData;
@@ -177,38 +176,33 @@ class ITK_TEMPLATE_EXPORT AttributeMorphologyBaseImageFilter : public ImageToIma
   void
   SetupOffsetVec(OffsetDirectVecType & PosOffsets, OffsetVecType & Offsets);
 
-  class GreyAndPos
-  {
-  public:
-    InputPixelType  Val;
-    OffsetValueType Pos;
-  };
-
-  GreyAndPos *      m_SortPixels;
+  // m_SortPixels contains offsets into the raw image
+  // it is sorted with a stable sort by grey level as the
+  // first step in the algorithm. The sorting step avoids
+  // the need to explicitly locate regional extrema.
+  OffsetValueType * m_SortPixels;
   OffsetValueType * m_Parent;
 
   // This is a bit ugly, but I can't see an easy way around
   InputPixelType * m_Raw;
 
-  class ComparePixStruct
+  class CompareOffsetType
   {
   public:
     TFunction m_TFunction;
+    // buf contains the raw data, which is what
+    // we want to sort by. i.e. the first value in
+    // the sorted buffer will be the location of the
+    // largest or smallest pixel.
+    InputPixelType * buf;
     bool
-    operator()(GreyAndPos const & l, GreyAndPos const & r) const
+    operator()(OffsetValueType const & l, OffsetValueType const & r) const
     {
-      if (m_TFunction(l.Val, r.Val))
-      {
-        return true;
-      }
-      if (l.Val == r.Val)
-      {
-        return (l.Pos < r.Pos);
-      }
-      return false;
+      return (m_TFunction(buf[l], buf[r]));
     }
   };
 
+  CompareOffsetType m_CompareOffset;
   // version from PAMI. Note - using the AuxData array rather than the
   // parent array to store area
   void

Modules/Nonunit/Review/include/itkAttributeMorphologyBaseImageFilter.hxx

@@ -97,7 +97,7 @@ AttributeMorphologyBaseImageFilter<TInputImage, TOutputImage, TAttribute, TFunct
 
   fit = faceList.begin();
 
-  m_SortPixels = new GreyAndPos[buffsize];
+  m_SortPixels = new OffsetValueType[buffsize];
   m_Parent = new OffsetValueType[buffsize];
 
   // This is a bit ugly, but I can't see an easy way around
@@ -112,18 +112,19 @@ AttributeMorphologyBaseImageFilter<TInputImage, TOutputImage, TAttribute, TFunct
 
   for (RegIt.GoToBegin(); !RegIt.IsAtEnd(); ++RegIt, ++pos)
   {
-    GreyAndPos P;
-    P.Val = RegIt.Get();
-    P.Pos = pos;
-    m_SortPixels[pos] = P;
-    m_Raw[pos] = P.Val;
+    // GreyAndPos P;
+    // P.Val = RegIt.Get();
+    // P.Pos = pos;
+    m_SortPixels[pos] = pos;
+    m_Raw[pos] = RegIt.Get();
 
     m_Parent[pos] = INACTIVE;
     m_AuxData[pos] = -1; // invalid value;
     progress.CompletedPixel();
   }
   progress.CompletedPixel();
-  std::sort(&(m_SortPixels[0]), &(m_SortPixels[buffsize - 1]), ComparePixStruct());
+  m_CompareOffset.buf = m_Raw;
+  std::stable_sort(&(m_SortPixels[0]), &(m_SortPixels[buffsize - 1]), m_CompareOffset);
   progress.CompletedPixel();
 
   // set up the offset vector
@@ -133,13 +134,13 @@ AttributeMorphologyBaseImageFilter<TInputImage, TOutputImage, TAttribute, TFunct
 
   // the core algorithm
   // process first pixel
-  MakeSet(m_SortPixels[0].Pos);
+  MakeSet(m_SortPixels[0]);
   // m_Processed[0] = true;
   for (SizeValueType k = 1; k < buffsize; ++k)
   {
-    OffsetValueType ThisPos = m_SortPixels[k].Pos;
+    OffsetValueType ThisPos = m_SortPixels[k];
     IndexType       ThisWhere = input->ComputeIndex(ThisPos);
-    InputPixelType  ThisPix = m_SortPixels[k].Val;
+    InputPixelType  ThisPix = m_Raw[ThisPos];
     MakeSet(ThisPos);
     // Some optimization of bounds check
     if (fit->IsInside(ThisWhere))
@@ -187,7 +188,7 @@ AttributeMorphologyBaseImageFilter<TInputImage, TOutputImage, TAttribute, TFunct
   // result in parent
   for (pos = buffsize - 1; pos >= 0; --pos)
   {
-    OffsetValueType RPos = m_SortPixels[pos].Pos;
+    OffsetValueType RPos = m_SortPixels[pos];
     if (m_Parent[RPos] >= 0)
     {
       m_Raw[RPos] = m_Raw[m_Parent[RPos]];
@@ -200,7 +201,6 @@ AttributeMorphologyBaseImageFilter<TInputImage, TOutputImage, TAttribute, TFunct
     progress.CompletedPixel();
   }
 
-
   delete[] m_Raw;
   delete[] m_SortPixels;
   delete[] m_Parent;