Fix CLUT interpolation method

[imagemagick] / MagickCore / gem.c

diff --git a/MagickCore/gem.c b/MagickCore/gem.c
index dd7a0b0e1d111d32e900affe5a44ae0fc521cd49..9a111cb262d95695f734be8d4d1fa4b9c0dfb230 100644 (file)
--- a/MagickCore/gem.c
+++ b/MagickCore/gem.c
@@ -17,7 +17,7 @@
 %                                 August 1996                                 %
 %                                                                             %
 %                                                                             %
-%  Copyright 1999-2011 ImageMagick Studio LLC, a non-profit organization      %
+%  Copyright 1999-2012 ImageMagick Studio LLC, a non-profit organization      %
 %  dedicated to making software imaging solutions freely available.           %
 %                                                                             %
 %  You may not use this file except in compliance with the License.  You may  %
@@ -102,7 +102,7 @@ MagickPrivate void ConvertHSBToRGB(const double hue,const double saturation,
   assert(blue != (double *) NULL);
   if (saturation == 0.0)
     {
-      *red=ClampToQuantum((MagickRealType) QuantumRange*brightness);
+      *red=(double) ClampToQuantum((MagickRealType) QuantumRange*brightness);
       *green=(*red);
       *blue=(*red);
       return;
@@ -117,44 +117,44 @@ MagickPrivate void ConvertHSBToRGB(const double hue,const double saturation,
     case 0:
     default:
     {
-      *red=ClampToQuantum((MagickRealType) QuantumRange*brightness);
-      *green=ClampToQuantum((MagickRealType) QuantumRange*t);
-      *blue=ClampToQuantum((MagickRealType) QuantumRange*p);
+      *red=(double) ClampToQuantum((MagickRealType) QuantumRange*brightness);
+      *green=(double) ClampToQuantum((MagickRealType) QuantumRange*t);
+      *blue=(double) ClampToQuantum((MagickRealType) QuantumRange*p);
       break;
     }
     case 1:
     {
-      *red=ClampToQuantum((MagickRealType) QuantumRange*q);
-      *green=ClampToQuantum((MagickRealType) QuantumRange*brightness);
-      *blue=ClampToQuantum((MagickRealType) QuantumRange*p);
+      *red=(double) ClampToQuantum((MagickRealType) QuantumRange*q);
+      *green=(double) ClampToQuantum((MagickRealType) QuantumRange*brightness);
+      *blue=(double) ClampToQuantum((MagickRealType) QuantumRange*p);
       break;
     }
     case 2:
     {
-      *red=ClampToQuantum((MagickRealType) QuantumRange*p);
-      *green=ClampToQuantum((MagickRealType) QuantumRange*brightness);
-      *blue=ClampToQuantum((MagickRealType) QuantumRange*t);
+      *red=(double) ClampToQuantum((MagickRealType) QuantumRange*p);
+      *green=(double) ClampToQuantum((MagickRealType) QuantumRange*brightness);
+      *blue=(double) ClampToQuantum((MagickRealType) QuantumRange*t);
       break;
     }
     case 3:
     {
-      *red=ClampToQuantum((MagickRealType) QuantumRange*p);
-      *green=ClampToQuantum((MagickRealType) QuantumRange*q);
-      *blue=ClampToQuantum((MagickRealType) QuantumRange*brightness);
+      *red=(double) ClampToQuantum((MagickRealType) QuantumRange*p);
+      *green=(double) ClampToQuantum((MagickRealType) QuantumRange*q);
+      *blue=(double) ClampToQuantum((MagickRealType) QuantumRange*brightness);
       break;
     }
     case 4:
     {
-      *red=ClampToQuantum((MagickRealType) QuantumRange*t);
-      *green=ClampToQuantum((MagickRealType) QuantumRange*p);
-      *blue=ClampToQuantum((MagickRealType) QuantumRange*brightness);
+      *red=(double) ClampToQuantum((MagickRealType) QuantumRange*t);
+      *green=(double) ClampToQuantum((MagickRealType) QuantumRange*p);
+      *blue=(double) ClampToQuantum((MagickRealType) QuantumRange*brightness);
       break;
     }
     case 5:
     {
-      *red=ClampToQuantum((MagickRealType) QuantumRange*brightness);
-      *green=ClampToQuantum((MagickRealType) QuantumRange*p);
-      *blue=ClampToQuantum((MagickRealType) QuantumRange*q);
+      *red=(double) ClampToQuantum((MagickRealType) QuantumRange*brightness);
+      *green=(double) ClampToQuantum((MagickRealType) QuantumRange*p);
+      *blue=(double) ClampToQuantum((MagickRealType) QuantumRange*q);
       break;
     }
   }
@@ -222,7 +222,7 @@ MagickExport void ConvertHSLToRGB(const double hue,const double saturation,
   assert(blue != (double *) NULL);
   if (saturation == 0)
     {
-      *red=ClampToQuantum((MagickRealType) QuantumRange*lightness);
+      *red=(double) ClampToQuantum((MagickRealType) QuantumRange*lightness);
       *green=(*red);
       *blue=(*red);
       return;
@@ -235,9 +235,9 @@ MagickExport void ConvertHSLToRGB(const double hue,const double saturation,
   r=ConvertHueToRGB(m1,m2,hue+1.0/3.0);
   g=ConvertHueToRGB(m1,m2,hue);
   b=ConvertHueToRGB(m1,m2,hue-1.0/3.0);
-  *red=ClampToQuantum((MagickRealType) QuantumRange*r);
-  *green=ClampToQuantum((MagickRealType) QuantumRange*g);
-  *blue=ClampToQuantum((MagickRealType) QuantumRange*b);
+  *red=(double) ClampToQuantum((MagickRealType) QuantumRange*r);
+  *green=(double) ClampToQuantum((MagickRealType) QuantumRange*g);
+  *blue=(double) ClampToQuantum((MagickRealType) QuantumRange*b);
 }
 \f
 /*
@@ -290,9 +290,9 @@ MagickPrivate void ConvertHWBToRGB(const double hue,const double whiteness,
   v=1.0-blackness;
   if (hue == 0.0)
     {
-      *red=ClampToQuantum((MagickRealType) QuantumRange*v);
-      *green=ClampToQuantum((MagickRealType) QuantumRange*v);
-      *blue=ClampToQuantum((MagickRealType) QuantumRange*v);
+      *red=(double) ClampToQuantum((MagickRealType) QuantumRange*v);
+      *green=(double) ClampToQuantum((MagickRealType) QuantumRange*v);
+      *blue=(double) ClampToQuantum((MagickRealType) QuantumRange*v);
       return;
     }
   i=(ssize_t) floor(6.0*hue);
@@ -311,9 +311,9 @@ MagickPrivate void ConvertHWBToRGB(const double hue,const double whiteness,
     case 4: r=n; g=whiteness; b=v; break;
     case 5: r=v; g=whiteness; b=n; break;
   }
-  *red=ClampToQuantum((MagickRealType) QuantumRange*r);
-  *green=ClampToQuantum((MagickRealType) QuantumRange*g);
-  *blue=ClampToQuantum((MagickRealType) QuantumRange*b);
+  *red=(double) ClampToQuantum((MagickRealType) QuantumRange*r);
+  *green=(double) ClampToQuantum((MagickRealType) QuantumRange*g);
+  *blue=(double) ClampToQuantum((MagickRealType) QuantumRange*b);
 }
 \f
 /*
@@ -591,8 +591,7 @@ MagickExport double ExpandAffine(const AffineMatrix *affine)
 %  The format of the GenerateDifferentialNoise method is:
 %
 %      double GenerateDifferentialNoise(RandomInfo *random_info,
-%        const Quantum pixel,const NoiseType noise_type,
-%        const MagickRealType attenuate)
+%        const Quantum pixel,const NoiseType noise_type,const double attenuate)
 %
 %  A description of each parameter follows:
 %
@@ -606,16 +605,15 @@ MagickExport double ExpandAffine(const AffineMatrix *affine)
 %
 */
 MagickPrivate double GenerateDifferentialNoise(RandomInfo *random_info,
-  const Quantum pixel,const NoiseType noise_type,const MagickRealType attenuate)
+  const Quantum pixel,const NoiseType noise_type,const double attenuate)
 {
-#define NoiseEpsilon  (attenuate*1.0e-5)
-#define SigmaUniform  (attenuate*4.0)
-#define SigmaGaussian  (attenuate*4.0)
-#define SigmaImpulse  (attenuate*0.10)
-#define SigmaLaplacian (attenuate*10.0)
-#define SigmaMultiplicativeGaussian  (attenuate*1.0)
-#define SigmaPoisson  (attenuate*0.05)
-#define TauGaussian  (attenuate*20.0)
+#define SigmaUniform  (attenuate*0.015625)
+#define SigmaGaussian  (attenuate*0.015625)
+#define SigmaImpulse  (attenuate*0.1)
+#define SigmaLaplacian (attenuate*0.0390625)
+#define SigmaMultiplicativeGaussian  (attenuate*0.5)
+#define SigmaPoisson  (attenuate*12.5)
+#define TauGaussian  (attenuate*0.078125)
 
   double
     alpha,
@@ -629,7 +627,7 @@ MagickPrivate double GenerateDifferentialNoise(RandomInfo *random_info,
     case UniformNoise:
     default:
     {
-      noise=(double) pixel+QuantumRange*SigmaUniform*(alpha-0.5)/255.0;
+      noise=(double) (pixel+QuantumRange*SigmaUniform*(alpha-0.5));
       break;
     }
     case GaussianNoise:
@@ -644,18 +642,8 @@ MagickPrivate double GenerateDifferentialNoise(RandomInfo *random_info,
       gamma=sqrt(-2.0*log(alpha));
       sigma=gamma*cos((double) (2.0*MagickPI*beta));
       tau=gamma*sin((double) (2.0*MagickPI*beta));
-      noise=(double) pixel+sqrt((double) pixel)*SigmaGaussian*sigma/255.0+
-        QuantumRange*TauGaussian*tau/255.0;
-      break;
-    }
-    case MultiplicativeGaussianNoise:
-    {
-      sigma=1.0;
-      if (alpha > NoiseEpsilon)
-        sigma=sqrt(-2.0*log(alpha));
-      beta=GetPseudoRandomValue(random_info);
-      noise=(double) pixel+pixel*SigmaMultiplicativeGaussian*sigma*
-        cos((double) (2.0*MagickPI*beta))/2.0/255.0;
+      noise=(double) (pixel+sqrt((double) pixel)*SigmaGaussian*sigma+
+        QuantumRange*TauGaussian*tau);
       break;
     }
     case ImpulseNoise:
@@ -673,19 +661,28 @@ MagickPrivate double GenerateDifferentialNoise(RandomInfo *random_info,
     {
       if (alpha <= 0.5)
         {
-          if (alpha <= NoiseEpsilon)
-            noise=(double) pixel-(double) QuantumRange;
+          if (alpha <= MagickEpsilon)
+            noise=(double) (pixel-QuantumRange);
           else
-            noise=(double) pixel+QuantumRange*SigmaLaplacian*
-              log(2.0*alpha)/255.0+0.5;
+            noise=(double) (pixel+QuantumRange*SigmaLaplacian*
+              log(2.0*alpha)+0.5);
           break;
         }
       beta=1.0-alpha;
-      if (beta <= (0.5*NoiseEpsilon))
+      if (beta <= (0.5*MagickEpsilon))
         noise=(double) (pixel+QuantumRange);
       else
-        noise=(double) pixel-QuantumRange*SigmaLaplacian*log(2.0*beta)/255.0+
-          0.5;
+        noise=(double) (pixel-QuantumRange*SigmaLaplacian*log(2.0*beta)+0.5);
+      break;
+    }
+    case MultiplicativeGaussianNoise:
+    {
+      sigma=1.0;
+      if (alpha > MagickEpsilon)
+        sigma=sqrt(-2.0*log(alpha));
+      beta=GetPseudoRandomValue(random_info);
+      noise=(double) (pixel+pixel*SigmaMultiplicativeGaussian*sigma*
+        cos((double) (2.0*MagickPI*beta))/2.0);
       break;
     }
     case PoissonNoise:
@@ -696,18 +693,18 @@ MagickPrivate double GenerateDifferentialNoise(RandomInfo *random_info,
       register ssize_t
         i;
 
-      poisson=exp(-SigmaPoisson*255.0*QuantumScale*pixel);
+      poisson=exp(-SigmaPoisson*QuantumScale*pixel);
       for (i=0; alpha > poisson; i++)
       {
         beta=GetPseudoRandomValue(random_info);
         alpha*=beta;
       }
-      noise=(double) QuantumRange*i/SigmaPoisson/255.0;
+      noise=(double) (QuantumRange*i/SigmaPoisson);
       break;
     }
     case RandomNoise:
     {
-      noise=(double) QuantumRange*alpha;
+      noise=(double) (QuantumRange*alpha);
       break;
     }
   }