Simplify and clean up the horrible sigmoidal code

diff --git a/MagickCore/enhance.c b/MagickCore/enhance.c
index 461ce27c70088ffc0c4ad48ac17b8e822ab49760..0104c8a7134d36932bcde1d4fcf7eab478116b6b 100644 (file)
--- a/MagickCore/enhance.c
+++ b/MagickCore/enhance.c
@@ -3227,8 +3227,9 @@ MagickExport MagickBooleanType NegateImage(Image *image,
 %                                                                             %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-%  NormalizeImage() enhances the contrast of a color image by mapping the
-%  darkest 2 percent of all pixel to black and the brightest 1 percent to white.
+%  The NormalizeImage() method enhances the contrast of a color image by
+%  mapping the darkest 2 percent of all pixel to black and the brightest
+%  1 percent to white.
 %
 %  The format of the NormalizeImage method is:
 %
@@ -3268,9 +3269,10 @@ MagickExport MagickBooleanType NormalizeImage(Image *image,
 %  sigmoidal contrast algorithm.  Increase the contrast of the image using a
 %  sigmoidal transfer function without saturating highlights or shadows.
 %  Contrast indicates how much to increase the contrast (0 is none; 3 is
-%  typical; 20 is pushing it); mid-point indicates where threshold 'knee' of
-%  the curve falls (typical 50% for mid-gray). Set sharpen to MagickTrue to
-%  increase the image contrast otherwise the contrast is reduced.
+%  typical; 20 is pushing it); mid-point indicates where midtones fall in the
+%  resultant image (0 is white; 50% is middle-gray; 100% is black).  Set
+%  sharpen to MagickTrue to increase the image contrast otherwise the contrast
+%  is reduced.
 %
 %  The format of the SigmoidalContrastImage method is:
 %
@@ -3336,22 +3338,56 @@ MagickExport MagickBooleanType SigmoidalContrastImage(Image *image,
   {
     if (sharpen != MagickFalse)
       {
+#define sigmoidal(a,b,x)  (1/(1+exp((a)*((b)-(x)))))
+#if 1
+        /* Simpilified function scaling,
+         * with better 'contrast=0' or 'flatline' handling (greyscale)
+         */
+        double
+          u0 = sigmoidal(contrast,QuantumScale*midpoint,0.0),
+          u1 = sigmoidal(contrast,QuantumScale*midpoint,1.0);
+        sigmoidal_map[i]=(MagickRealType) ScaleMapToQuantum(
+           (MagickRealType)(MaxMap*(
+               (sigmoidal(contrast,QuantumScale*midpoint,(double)i/MaxMap)
+                  -(u0+u1)/2.0)/(u1-u0+MagickEpsilon)+0.5)   ));
+#else
+        /* Scaled sigmoidal formula...
+             (1/(1+exp(a*(b-u))) - 1/(1+exp(a))) /
+                     (1/(1+exp(a*(b-1)))/(1+exp(a)))) */
         sigmoidal_map[i]=(MagickRealType) ScaleMapToQuantum((MagickRealType)
           (MaxMap*((1.0/(1.0+exp(contrast*(midpoint/(double) QuantumRange-
-          (double) i/MaxMap))))-(1.0/(1.0+exp(contrast*(midpoint/(double)
-          QuantumRange)))))/((1.0/(1.0+exp(contrast*(midpoint/(double)
-          QuantumRange-1.0))))-(1.0/(1.0+exp(contrast*(midpoint/(double)
-          QuantumRange)))))+0.5));
+          (double) i/MaxMap))))-(1.0/(1.0+exp(contrast*(midpoint/
+          (double) QuantumRange)))))/((1.0/(1.0+exp(contrast*(midpoint/
+          (double) QuantumRange-1.0))))-(1.0/(1.0+exp(contrast*(midpoint/
+          (double) QuantumRange)))))+0.5));
+#endif
         continue;
       }
+#if 1
+    {
+      /* Inverse -- See
+         http://osdir.com/ml/video.image-magick.devel/2005-04/msg00006.html
+      */
+      double
+        min = sigmoidal(contrast,1.0,0.0),
+        max = sigmoidal(contrast,QuantumScale*midpoint,1.0),
+        xi  = min+(double)i/MaxMap*(max-min);
+      sigmoidal_map[i]=(MagickRealType) ScaleMapToQuantum(
+         (MagickRealType)(MaxMap*(
+             QuantumScale*midpoint-log((1-xi)/xi)/contrast) ));
+    }
+#else
+    /* expanded form of the above */
     sigmoidal_map[i]=(MagickRealType) ScaleMapToQuantum((MagickRealType)
-      (MaxMap*(QuantumScale*midpoint-log((1.0-(1.0/(1.0+exp(midpoint/(double)
-      QuantumRange*contrast))+((double) i/MaxMap)*((1.0/(1.0+exp(contrast*(
-      midpoint/(double) QuantumRange-1.0))))-(1.0/(1.0+exp(midpoint/(double)
-      QuantumRange*contrast))))))/(1.0/(1.0+exp(midpoint/(double) QuantumRange*
-      contrast))+((double) i/MaxMap)*((1.0/(1.0+exp(contrast*(midpoint/(double)
-      QuantumRange-1.0))))-(1.0/(1.0+exp(midpoint/(double) QuantumRange*
-      contrast))))))/contrast)));
+      (MaxMap*(QuantumScale*midpoint-log((1.0-(1.0/(1.0+exp(midpoint/
+      (double) QuantumRange*contrast))+((double) i/MaxMap)*((1.0/
+      (1.0+exp(contrast*(midpoint/(double) QuantumRange-1.0))))-(1.0/
+      (1.0+exp(midpoint/(double) QuantumRange*contrast))))))/
+      (1.0/(1.0+exp(midpoint/(double) QuantumRange*contrast))+
+      ((double) i/MaxMap)*((1.0/(1.0+exp(contrast*(midpoint/
+      (double) QuantumRange-1.0))))-(1.0/(1.0+exp(midpoint/
+      (double) QuantumRange*contrast))))))/contrast)));
+#endif
   }
   if (image->storage_class == PseudoClass)
     {