BUG: Fix subtle error for dark pixels

Consistently increment pixel intensity prior to computing its
logarithm.  Clarify comments describing the increment and later
decrement.

The code takes logarithms of pixel intensity values, with the intent
of first incrementing each intensity value to avoid log(0).  However,
the code failed to do the increment in one place, which gave small
intensity errors in the output.

Author: Lee Newberg

include/itkStructurePreservingColorNormalizationFilter.hxx

@@ -274,7 +274,6 @@ StructurePreservingColorNormalizationFilter<TImage>::ImageToMatrix(RegionConstIt
 
   SizeValueType numberOfRows = std::min(numberOfPixels, maxNumberOfRows);
 
-  // To avoid zeros, every color intensity is incremented.
   CalcMatrixType matrixV{ numberOfRows, m_NumberOfColors };
   for (iter.GoToBegin(); !iter.IsAtEnd(); ++iter)
   {
@@ -284,6 +283,7 @@ StructurePreservingColorNormalizationFilter<TImage>::ImageToMatrix(RegionConstIt
       PixelType pixelValue = iter.Get();
       for (Eigen::Index color = 0; color < m_NumberOfColors; ++color)
       {
+        // To avoid std::log(0), every color intensity is incremented.
         matrixV(numberOfRows, color) = pixelValue[color] + One;
       }
     }
@@ -575,7 +575,7 @@ StructurePreservingColorNormalizationFilter<TImage>::DistinguishersToColors(Calc
   // m_ColorIndexSuppressedByHematoxylin, and the index of the color
   // suppressed by eosin is indicated by
   // m_ColorIndexSuppressedByEosin.
-  std::vector<int> suppressedIndex;
+  std::vector<int>               suppressedIndex;
   std::vector<CalcRowVectorType> suppressedPixel;
   for (int fromRow = 0; fromRow < distinguishers.rows(); ++fromRow)
   {
@@ -765,7 +765,8 @@ StructurePreservingColorNormalizationFilter<TImage>::NMFsToImage(const CalcMatri
     PixelType pixelValue = inIt.Get();
     for (Eigen::Index color = 0; color < m_NumberOfColors; ++color)
     {
-      matrixV(pixelIndex, color) = static_cast<CalcElementType>(pixelValue[color]);
+      // To avoid std::log(0), every color intensity is incremented.
+      matrixV(pixelIndex, color) = static_cast<CalcElementType>(pixelValue[color]) + One;
     }
   }
 
@@ -800,6 +801,7 @@ StructurePreservingColorNormalizationFilter<TImage>::NMFsToImage(const CalcMatri
     }
     for (Eigen::Index color = 0; color < m_NumberOfColors; ++color)
     {
+      // Every color intensity is decremented to cancel the increment that was applied prior to std::log.
       pixelValue[color] = std::max(std::min(matrixV(pixelIndex, color) - One, upperbound), lowerbound);
     }
     const PixelType inputPixel = inIt.Get();

test/Baseline/itkStructurePreservingColorNormalizationFilterTestOutput.png.sha512

@@ -1 +1 @@
-e9381960df175c97efc56743dd8d9e81eb71bb079186a3788e746ad11f6b95fbc24f5335a8b76573767358961ab16dcd66b614895abf783d967ec82024f368b7
+2a720ed0c8fbe0ee530c88d01de638671fa2380520ec59ccfb2ad126b3c17e5acc5521ce3001604a5c623fcfa14a81e3dad960a5f51be76269d398911e3f8740