# HG changeset patch
# User Sebastien Jodogne <s.jodogne@gmail.com>
# Date 1582558019 -3600
# Node ID 5f64c866108ad6d58e4ff588fcaa987fbdd004d6
# Parent  453c0ece560afd86eecab0d5ff7b4020d33086e9
merging implementations of ImageProcessing::ShiftScale() and ApplyWindowing()

diff -r 453c0ece560a -r 5f64c866108a Core/Images/ImageProcessing.cpp
--- a/Core/Images/ImageProcessing.cpp	Thu Feb 20 20:09:24 2020 +0100
+++ b/Core/Images/ImageProcessing.cpp	Mon Feb 24 16:26:59 2020 +0100
@@ -41,11 +41,11 @@
 
 #ifdef __EMSCRIPTEN__
 /* 
-Avoid this error:
------------------
-.../boost/math/special_functions/round.hpp:118:12: warning: implicit conversion from 'std::__2::numeric_limits<long long>::type' (aka 'long long') to 'float' changes value from 9223372036854775807 to 9223372036854775808 [-Wimplicit-int-float-conversion]
-.../mnt/c/osi/dev/orthanc/Core/Images/ImageProcessing.cpp:333:28: note: in instantiation of function template specialization 'boost::math::llround<float>' requested here
-.../mnt/c/osi/dev/orthanc/Core/Images/ImageProcessing.cpp:1006:9: note: in instantiation of function template specialization 'Orthanc::MultiplyConstantInternal<unsigned char, true>' requested here
+   Avoid this error:
+   -----------------
+   .../boost/math/special_functions/round.hpp:118:12: warning: implicit conversion from 'std::__2::numeric_limits<long long>::type' (aka 'long long') to 'float' changes value from 9223372036854775807 to 9223372036854775808 [-Wimplicit-int-float-conversion]
+   .../mnt/c/osi/dev/orthanc/Core/Images/ImageProcessing.cpp:333:28: note: in instantiation of function template specialization 'boost::math::llround<float>' requested here
+   .../mnt/c/osi/dev/orthanc/Core/Images/ImageProcessing.cpp:1006:9: note: in instantiation of function template specialization 'Orthanc::MultiplyConstantInternal<unsigned char, true>' requested here
 */
 #pragma GCC diagnostic ignored "-Wimplicit-int-float-conversion"
 #endif 
@@ -385,28 +385,47 @@
   }
 
 
-  template <typename PixelType,
-            bool UseRound>
-  static void ShiftScaleInternal(ImageAccessor& image,
-                                 float offset,
-                                 float scaling,
-                                 const PixelType LowestValue = std::numeric_limits<PixelType>::min())
+  // Computes "a * x + b" at each pixel => Note that this is not the
+  // same convention as in "ShiftScale()"
+  template <typename TargetType,
+            typename SourceType,
+            bool UseRound,
+            bool Invert>
+  static void ShiftScaleInternal(ImageAccessor& target,
+                                 const ImageAccessor& source,
+                                 float a,
+                                 float b,
+                                 const TargetType LowestValue)
+  // This function can be applied inplace (source == target)
   {
-    const PixelType minPixelValue = LowestValue;
-    const PixelType maxPixelValue = std::numeric_limits<PixelType>::max();
+    if (source.GetWidth() != target.GetWidth() ||
+        source.GetHeight() != target.GetHeight())
+    {
+      throw OrthancException(ErrorCode_IncompatibleImageSize);
+    }
+
+    if (&source == &target &&
+        source.GetFormat() != target.GetFormat())
+    {
+      throw OrthancException(ErrorCode_IncompatibleImageFormat);
+    }
+    
+    const TargetType minPixelValue = LowestValue;
+    const TargetType maxPixelValue = std::numeric_limits<TargetType>::max();
     const float minFloatValue = static_cast<float>(LowestValue);
     const float maxFloatValue = static_cast<float>(maxPixelValue);
 
-    const unsigned int height = image.GetHeight();
-    const unsigned int width = image.GetWidth();
+    const unsigned int height = target.GetHeight();
+    const unsigned int width = target.GetWidth();
     
     for (unsigned int y = 0; y < height; y++)
     {
-      PixelType* p = reinterpret_cast<PixelType*>(image.GetRow(y));
+      TargetType* p = reinterpret_cast<TargetType*>(target.GetRow(y));
+      const SourceType* q = reinterpret_cast<const SourceType*>(source.GetRow(y));
 
-      for (unsigned int x = 0; x < width; x++, p++)
+      for (unsigned int x = 0; x < width; x++, p++, q++)
       {
-        float v = (static_cast<float>(*p) + offset) * scaling;
+        float v = a * static_cast<float>(*q) + b;
 
         if (v >= maxFloatValue)
         {
@@ -419,11 +438,16 @@
         else if (UseRound)
         {
           // The "round" operation is very costly
-          *p = static_cast<PixelType>(boost::math::iround(v));
+          *p = static_cast<TargetType>(boost::math::iround(v));
         }
         else
         {
-          *p = static_cast<PixelType>(v);
+          *p = static_cast<TargetType>(std::floor(v));
+        }
+
+        if (Invert)
+        {
+          *p = maxPixelValue - *p;
         }
       }
     }
@@ -498,47 +522,29 @@
                                      float rescaleIntercept,
                                      bool invert)
   {
+    assert(sizeof(SourceType) == source.GetBytesPerPixel() &&
+           sizeof(TargetType) == target.GetBytesPerPixel());
+    
     // WARNING - "::min()" should be replaced by "::lowest()" if
     // dealing with float or double (which is not the case so far)
     assert(sizeof(TargetType) <= 2);  // Safeguard to remember about "float/double"
     const TargetType minTargetValue = std::numeric_limits<TargetType>::min();
     const TargetType maxTargetValue = std::numeric_limits<TargetType>::max();
-    const float minFloatValue = static_cast<float>(minTargetValue);
     const float maxFloatValue = static_cast<float>(maxTargetValue);
-
+    
     const float windowIntercept = windowCenter - windowWidth / 2.0f;
     const float windowSlope = (maxFloatValue + 1.0f) / windowWidth;
 
-    const unsigned int width = source.GetWidth();
-    const unsigned int height = source.GetHeight();
-
-    for (unsigned int y = 0; y < height; y++)
-    {
-      TargetType* t = reinterpret_cast<TargetType*>(target.GetRow(y));
-      const SourceType* s = reinterpret_cast<const SourceType*>(source.GetConstRow(y));
+    const float a = rescaleSlope * windowSlope;
+    const float b = (rescaleIntercept - windowIntercept) * windowSlope;
 
-      for (unsigned int x = 0; x < width; x++, t++, s++)
-      {
-        float rescaledValue = *s * rescaleSlope + rescaleIntercept;
-        float v = (rescaledValue - windowIntercept) * windowSlope;
-        if (v <= minFloatValue)
-        {
-          v = minFloatValue;
-        }
-        else if (v >= maxFloatValue)
-        {
-          v = maxFloatValue;
-        }
-
-        TargetType vv = static_cast<TargetType>(v);
-
-        if (invert)
-        {
-          vv = maxTargetValue - vv;
-        }
-
-        *t = static_cast<TargetType>(vv);
-      }
+    if (invert)
+    {
+      ShiftScaleInternal<TargetType, SourceType, false, true>(target, source, a, b, minTargetValue);
+    }
+    else
+    {
+      ShiftScaleInternal<TargetType, SourceType, false, false>(target, source, a, b, minTargetValue);
     }
   }
 
@@ -562,42 +568,42 @@
       {
         switch (target.GetFormat())
         {
-        case Orthanc::PixelFormat_Grayscale8:
-          ApplyWindowingInternal<uint8_t, float>(target, source, windowCenter, windowWidth, rescaleSlope, rescaleIntercept, invert);
-          break;
-        case Orthanc::PixelFormat_Grayscale16:
-          ApplyWindowingInternal<uint16_t, float>(target, source, windowCenter, windowWidth, rescaleSlope, rescaleIntercept, invert);
-          break;
-        default:
-          throw OrthancException(ErrorCode_NotImplemented);
+          case Orthanc::PixelFormat_Grayscale8:
+            ApplyWindowingInternal<uint8_t, float>(target, source, windowCenter, windowWidth, rescaleSlope, rescaleIntercept, invert);
+            break;
+          case Orthanc::PixelFormat_Grayscale16:
+            ApplyWindowingInternal<uint16_t, float>(target, source, windowCenter, windowWidth, rescaleSlope, rescaleIntercept, invert);
+            break;
+          default:
+            throw OrthancException(ErrorCode_NotImplemented);
         }
       };break;
       case Orthanc::PixelFormat_Grayscale8:
       {
         switch (target.GetFormat())
         {
-        case Orthanc::PixelFormat_Grayscale8:
-          ApplyWindowingInternal<uint8_t, uint8_t>(target, source, windowCenter, windowWidth, rescaleSlope, rescaleIntercept, invert);
-          break;
-        case Orthanc::PixelFormat_Grayscale16:
-          ApplyWindowingInternal<uint16_t, uint8_t>(target, source, windowCenter, windowWidth, rescaleSlope, rescaleIntercept, invert);
-          break;
-        default:
-          throw OrthancException(ErrorCode_NotImplemented);
+          case Orthanc::PixelFormat_Grayscale8:
+            ApplyWindowingInternal<uint8_t, uint8_t>(target, source, windowCenter, windowWidth, rescaleSlope, rescaleIntercept, invert);
+            break;
+          case Orthanc::PixelFormat_Grayscale16:
+            ApplyWindowingInternal<uint16_t, uint8_t>(target, source, windowCenter, windowWidth, rescaleSlope, rescaleIntercept, invert);
+            break;
+          default:
+            throw OrthancException(ErrorCode_NotImplemented);
         }
       };break;
       case Orthanc::PixelFormat_Grayscale16:
       {
         switch (target.GetFormat())
         {
-        case Orthanc::PixelFormat_Grayscale8:
-          ApplyWindowingInternal<uint8_t, uint16_t>(target, source, windowCenter, windowWidth, rescaleSlope, rescaleIntercept, invert);
-          break;
-        case Orthanc::PixelFormat_Grayscale16:
-          ApplyWindowingInternal<uint16_t, uint16_t>(target, source, windowCenter, windowWidth, rescaleSlope, rescaleIntercept, invert);
-          break;
-        default:
-          throw OrthancException(ErrorCode_NotImplemented);
+          case Orthanc::PixelFormat_Grayscale8:
+            ApplyWindowingInternal<uint8_t, uint16_t>(target, source, windowCenter, windowWidth, rescaleSlope, rescaleIntercept, invert);
+            break;
+          case Orthanc::PixelFormat_Grayscale16:
+            ApplyWindowingInternal<uint16_t, uint16_t>(target, source, windowCenter, windowWidth, rescaleSlope, rescaleIntercept, invert);
+            break;
+          default:
+            throw OrthancException(ErrorCode_NotImplemented);
         }
       };break;
       default:
@@ -1181,7 +1187,7 @@
         ShiftRightInternal<uint16_t>(image, shift);
         break;
       }
-    default:
+      default:
         throw OrthancException(ErrorCode_NotImplemented);
     }
   }
@@ -1210,7 +1216,7 @@
         ShiftLeftInternal<uint16_t>(image, shift);
         break;
       }
-    default:
+      default:
         throw OrthancException(ErrorCode_NotImplemented);
     }
   }
@@ -1366,49 +1372,55 @@
                                    float scaling,
                                    bool useRound)
   {
+    // Rewrite "(x + offset) * scaling" as "a * x + b"
+
+    const float a = scaling;
+    const float b = offset * scaling;
+    
     switch (image.GetFormat())
     {
       case PixelFormat_Grayscale8:
         if (useRound)
         {
-          ShiftScaleInternal<uint8_t, true>(image, offset, scaling);
+          ShiftScaleInternal<uint8_t, uint8_t, true, false>(image, image, a, b, std::numeric_limits<uint8_t>::min());
         }
         else
         {
-          ShiftScaleInternal<uint8_t, false>(image, offset, scaling);
+          ShiftScaleInternal<uint8_t, uint8_t, false, false>(image, image, a, b, std::numeric_limits<uint8_t>::min());
         }
         return;
 
       case PixelFormat_Grayscale16:
         if (useRound)
         {
-          ShiftScaleInternal<uint16_t, true>(image, offset, scaling);
+          ShiftScaleInternal<uint16_t, uint16_t, true, false>(image, image, a, b, std::numeric_limits<uint16_t>::min());
         }
         else
         {
-          ShiftScaleInternal<uint16_t, false>(image, offset, scaling);
+          ShiftScaleInternal<uint16_t, uint16_t, false, false>(image, image, a, b, std::numeric_limits<uint16_t>::min());
         }
         return;
 
       case PixelFormat_SignedGrayscale16:
         if (useRound)
         {
-          ShiftScaleInternal<int16_t, true>(image, offset, scaling);
+          ShiftScaleInternal<int16_t, int16_t, true, false>(image, image, a, b, std::numeric_limits<int16_t>::min());
         }
         else
         {
-          ShiftScaleInternal<int16_t, false>(image, offset, scaling);
+          ShiftScaleInternal<int16_t, int16_t, false, false>(image, image, a, b, std::numeric_limits<int16_t>::min());
         }
         return;
 
       case PixelFormat_Float32:
+        // "::min()" must be replaced by "::lowest()" if dealing with float or double
         if (useRound)
         {
-          ShiftScaleInternal<float, true>(image, offset, scaling, -std::numeric_limits<float>::max());
+          ShiftScaleInternal<float, float, true, false>(image, image, a, b, std::numeric_limits<float>::lowest());
         }
         else
         {
-          ShiftScaleInternal<float, false>(image, offset, scaling, -std::numeric_limits<float>::max());
+          ShiftScaleInternal<float, float, false, false>(image, image, a, b, std::numeric_limits<float>::lowest());
         }
         return;
 
@@ -2216,7 +2228,7 @@
 
         // Deal with the right border
         for (unsigned int x = static_cast<unsigned int>(
-          horizontalAnchor + width - horizontal.size() + 1); x < width; x++)
+               horizontalAnchor + width - horizontal.size() + 1); x < width; x++)
         {
           for (unsigned int c = 0; c < ChannelsCount; c++, p++)
           {
diff -r 453c0ece560a -r 5f64c866108a Core/Images/ImageProcessing.h
--- a/Core/Images/ImageProcessing.h	Thu Feb 20 20:09:24 2020 +0100
+++ b/Core/Images/ImageProcessing.h	Mon Feb 24 16:26:59 2020 +0100
@@ -127,7 +127,7 @@
                           float factor,
                           bool useRound);
 
-    // "useRound" is expensive
+    // Computes "(x + offset) * scaling" inplace. "useRound" is expensive.
     void ShiftScale(ImageAccessor& image,
                     float offset,
                     float scaling,
diff -r 453c0ece560a -r 5f64c866108a UnitTestsSources/ImageProcessingTests.cpp
--- a/UnitTestsSources/ImageProcessingTests.cpp	Thu Feb 20 20:09:24 2020 +0100
+++ b/UnitTestsSources/ImageProcessingTests.cpp	Mon Feb 24 16:26:59 2020 +0100
@@ -285,7 +285,7 @@
 static void SetGrayscale16Pixel(ImageAccessor& image,
                                 unsigned int x,
                                 unsigned int y,
-                                uint8_t value)
+                                uint16_t value)
 {
   ImageTraits<PixelFormat_Grayscale16>::SetPixel(image, value, x, y);
 }
@@ -301,6 +301,25 @@
   return p == value;
 }
 
+static void SetSignedGrayscale16Pixel(ImageAccessor& image,
+                                unsigned int x,
+                                unsigned int y,
+                                int16_t value)
+{
+  ImageTraits<PixelFormat_SignedGrayscale16>::SetPixel(image, value, x, y);
+}
+
+static bool TestSignedGrayscale16Pixel(const ImageAccessor& image,
+                                       unsigned int x,
+                                       unsigned int y,
+                                       int16_t value)
+{
+  PixelTraits<PixelFormat_SignedGrayscale16>::PixelType p;
+  ImageTraits<PixelFormat_SignedGrayscale16>::GetPixel(p, image, x, y);
+  if (p != value) printf("%d %d\n", p, value);
+  return p == value;
+}
+
 static void SetRGB24Pixel(ImageAccessor& image,
                           unsigned int x,
                           unsigned int y,
@@ -831,7 +850,7 @@
 
     {
       Image target(PixelFormat_Grayscale8, 6, 1, false);
-      ImageProcessing::ApplyWindowing(target, image, 5000.0f, 10000.0f, 1000.0f, 0.0f, false);
+      ImageProcessing::ApplyWindowing(target, image, 5000.0f, 10000.01f, 1000.0f, 0.0f, false);
 
       ASSERT_TRUE(TestGrayscale8Pixel(target, 0, 0, 0));
       ASSERT_TRUE(TestGrayscale8Pixel(target, 1, 0, 0));
@@ -843,7 +862,7 @@
 
     {
       Image target(PixelFormat_Grayscale8, 6, 1, false);
-      ImageProcessing::ApplyWindowing(target, image, 5000.0f, 10000.0f, 1000.0f, 0.0f, true);
+      ImageProcessing::ApplyWindowing(target, image, 5000.0f, 10000.01f, 1000.0f, 0.0f, true);
 
       ASSERT_TRUE(TestGrayscale8Pixel(target, 0, 0, 255));
       ASSERT_TRUE(TestGrayscale8Pixel(target, 1, 0, 255));
@@ -855,7 +874,7 @@
 
     {
       Image target(PixelFormat_Grayscale8, 6, 1, false);
-      ImageProcessing::ApplyWindowing(target, image, 50.0f, 100.0f, 10.0f, 30.0f, false);
+      ImageProcessing::ApplyWindowing(target, image, 50.0f, 100.1f, 10.0f, 30.0f, false);
 
       ASSERT_TRUE(TestGrayscale8Pixel(target, 0, 0, 0));  // (-5 * 10) + 30 => pixel value = -20 => 0
       ASSERT_TRUE(TestGrayscale8Pixel(target, 1, 0, 256*30/100));  // ((0 * 10) + 30 => pixel value = 30 => 30%
@@ -938,3 +957,57 @@
     }
   }
 }
+
+
+TEST(ImageProcessing, ShiftScaleGrayscale8)
+{
+  Image image(PixelFormat_Grayscale8, 5, 1, false);
+  SetGrayscale8Pixel(image, 0, 0, 0);
+  SetGrayscale8Pixel(image, 1, 0, 2);
+  SetGrayscale8Pixel(image, 2, 0, 5);
+  SetGrayscale8Pixel(image, 3, 0, 10);
+  SetGrayscale8Pixel(image, 4, 0, 255);
+
+  ImageProcessing::ShiftScale(image, -1.1, 1.5, true);
+  ASSERT_TRUE(TestGrayscale8Pixel(image, 0, 0, 0));
+  ASSERT_TRUE(TestGrayscale8Pixel(image, 1, 0, 1));
+  ASSERT_TRUE(TestGrayscale8Pixel(image, 2, 0, 6));
+  ASSERT_TRUE(TestGrayscale8Pixel(image, 3, 0, 13));
+  ASSERT_TRUE(TestGrayscale8Pixel(image, 4, 0, 255));
+}
+
+
+TEST(ImageProcessing, ShiftScaleGrayscale16)
+{
+  Image image(PixelFormat_Grayscale16, 5, 1, false);
+  SetGrayscale16Pixel(image, 0, 0, 0);
+  SetGrayscale16Pixel(image, 1, 0, 2);
+  SetGrayscale16Pixel(image, 2, 0, 5);
+  SetGrayscale16Pixel(image, 3, 0, 10);
+  SetGrayscale16Pixel(image, 4, 0, 255);
+
+  ImageProcessing::ShiftScale(image, -1.1, 1.5, true);
+  ASSERT_TRUE(TestGrayscale16Pixel(image, 0, 0, 0));
+  ASSERT_TRUE(TestGrayscale16Pixel(image, 1, 0, 1));
+  ASSERT_TRUE(TestGrayscale16Pixel(image, 2, 0, 6));
+  ASSERT_TRUE(TestGrayscale16Pixel(image, 3, 0, 13));
+  ASSERT_TRUE(TestGrayscale16Pixel(image, 4, 0, 381));
+}
+
+
+TEST(ImageProcessing, ShiftScaleSignedGrayscale16)
+{
+  Image image(PixelFormat_SignedGrayscale16, 5, 1, false);
+  SetSignedGrayscale16Pixel(image, 0, 0, 0);
+  SetSignedGrayscale16Pixel(image, 1, 0, 2);
+  SetSignedGrayscale16Pixel(image, 2, 0, 5);
+  SetSignedGrayscale16Pixel(image, 3, 0, 10);
+  SetSignedGrayscale16Pixel(image, 4, 0, 255);
+
+  ImageProcessing::ShiftScale(image, -17.1, 11.5, true);
+  ASSERT_TRUE(TestSignedGrayscale16Pixel(image, 0, 0, -197));
+  ASSERT_TRUE(TestSignedGrayscale16Pixel(image, 1, 0, -174));
+  ASSERT_TRUE(TestSignedGrayscale16Pixel(image, 2, 0, -139));
+  ASSERT_TRUE(TestSignedGrayscale16Pixel(image, 3, 0, -82));
+  ASSERT_TRUE(TestSignedGrayscale16Pixel(image, 4, 0, 2736));
+}