]> granicus.if.org Git - imagemagick/blobdiff - MagickCore/statistic.c
sigmoidal-contrast: prevent argument out of range and remove unnecessary ClampToQuantum
[imagemagick] / MagickCore / statistic.c
index 087cb494a325c7ead5e5f6e0f8db8ee4fc989004..85fbcb3f906d3cbdcf6536ff712fb75ad946f970 100644 (file)
@@ -81,6 +81,7 @@
 #include "MagickCore/quantum-private.h"
 #include "MagickCore/random_.h"
 #include "MagickCore/random-private.h"
+#include "MagickCore/resource_.h"
 #include "MagickCore/segment.h"
 #include "MagickCore/semaphore.h"
 #include "MagickCore/signature-private.h"
 
 typedef struct _PixelChannels
 {
-  MagickRealType
+  double
     channel[CompositePixelChannel];
 } PixelChannels;
 
@@ -140,7 +141,7 @@ static PixelChannels **DestroyPixelThreadSet(PixelChannels **pixels)
     i;
 
   assert(pixels != (PixelChannels **) NULL);
-  for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
+  for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
     if (pixels[i] != (PixelChannels *) NULL)
       pixels[i]=(PixelChannels *) RelinquishMagickMemory(pixels[i]);
   pixels=(PixelChannels **) RelinquishMagickMemory(pixels);
@@ -160,7 +161,7 @@ static PixelChannels **AcquirePixelThreadSet(const Image *image,
     length,
     number_threads;
 
-  number_threads=GetOpenMPMaximumThreads();
+  number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
   pixels=(PixelChannels **) AcquireQuantumMemory(number_threads,
     sizeof(*pixels));
   if (pixels == (PixelChannels **) NULL)
@@ -206,7 +207,7 @@ static int IntensityCompare(const void *x,const void *y)
     *color_1,
     *color_2;
 
-  MagickRealType
+  double
     distance;
 
   register ssize_t
@@ -216,7 +217,7 @@ static int IntensityCompare(const void *x,const void *y)
   color_2=(const PixelChannels *) y;
   distance=0.0;
   for (i=0; i < MaxPixelChannels; i++)
-    distance+=color_1->channel[i]-(MagickRealType) color_2->channel[i];
+    distance+=color_1->channel[i]-(double) color_2->channel[i];
   return(distance < 0 ? -1 : distance > 0 ? 1 : 0);
 }
 
@@ -231,10 +232,10 @@ static inline double MagickMin(const double x,const double y)
   return(y);
 }
 
-static MagickRealType ApplyEvaluateOperator(RandomInfo *random_info,
-  Quantum pixel,const MagickEvaluateOperator op,const MagickRealType value)
+static double ApplyEvaluateOperator(RandomInfo *random_info,
+  Quantum pixel,const MagickEvaluateOperator op,const double value)
 {
-  MagickRealType
+  double
     result;
 
   result=0.0;
@@ -244,12 +245,12 @@ static MagickRealType ApplyEvaluateOperator(RandomInfo *random_info,
       break;
     case AbsEvaluateOperator:
     {
-      result=(MagickRealType) fabs((double) (pixel+value));
+      result=(double) fabs((double) (pixel+value));
       break;
     }
     case AddEvaluateOperator:
     {
-      result=(MagickRealType) (pixel+value);
+      result=(double) (pixel+value);
       break;
     }
     case AddModulusEvaluateOperator:
@@ -266,12 +267,12 @@ static MagickRealType ApplyEvaluateOperator(RandomInfo *random_info,
     }
     case AndEvaluateOperator:
     {
-      result=(MagickRealType) ((size_t) pixel & (size_t) (value+0.5));
+      result=(double) ((size_t) pixel & (size_t) (value+0.5));
       break;
     }
     case CosineEvaluateOperator:
     {
-      result=(MagickRealType) (QuantumRange*(0.5*cos((double) (2.0*MagickPI*
+      result=(double) (QuantumRange*(0.5*cos((double) (2.0*MagickPI*
         QuantumScale*pixel*value))+0.5));
       break;
     }
@@ -282,90 +283,91 @@ static MagickRealType ApplyEvaluateOperator(RandomInfo *random_info,
     }
     case ExponentialEvaluateOperator:
     {
-      result=(MagickRealType) (QuantumRange*exp((double) (value*QuantumScale*
+      result=(double) (QuantumRange*exp((double) (value*QuantumScale*
         pixel)));
       break;
     }
     case GaussianNoiseEvaluateOperator:
     {
-      result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
+      result=(double) GenerateDifferentialNoise(random_info,pixel,
         GaussianNoise,value);
       break;
     }
     case ImpulseNoiseEvaluateOperator:
     {
-      result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
+      result=(double) GenerateDifferentialNoise(random_info,pixel,
         ImpulseNoise,value);
       break;
     }
     case LaplacianNoiseEvaluateOperator:
     {
-      result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
+      result=(double) GenerateDifferentialNoise(random_info,pixel,
         LaplacianNoise,value);
       break;
     }
     case LeftShiftEvaluateOperator:
     {
-      result=(MagickRealType) ((size_t) pixel << (size_t) (value+0.5));
+      result=(double) ((size_t) pixel << (size_t) (value+0.5));
       break;
     }
     case LogEvaluateOperator:
     {
-      result=(MagickRealType) (QuantumRange*log((double) (QuantumScale*value*
-        pixel+1.0))/log((double) (value+1.0)));
+      if ((QuantumScale*pixel) >= MagickEpsilon)
+        result=(double) (QuantumRange*log((double) (QuantumScale*value*
+          pixel+1.0))/log((double) (value+1.0)));
       break;
     }
     case MaxEvaluateOperator:
     {
-      result=(MagickRealType) EvaluateMax((double) pixel,value);
+      result=(double) EvaluateMax((double) pixel,value);
       break;
     }
     case MeanEvaluateOperator:
     {
-      result=(MagickRealType) (pixel+value);
+      result=(double) (pixel+value);
       break;
     }
     case MedianEvaluateOperator:
     {
-      result=(MagickRealType) (pixel+value);
+      result=(double) (pixel+value);
       break;
     }
     case MinEvaluateOperator:
     {
-      result=(MagickRealType) MagickMin((double) pixel,value);
+      result=(double) MagickMin((double) pixel,value);
       break;
     }
     case MultiplicativeNoiseEvaluateOperator:
     {
-      result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
+      result=(double) GenerateDifferentialNoise(random_info,pixel,
         MultiplicativeGaussianNoise,value);
       break;
     }
     case MultiplyEvaluateOperator:
     {
-      result=(MagickRealType) (value*pixel);
+      result=(double) (value*pixel);
       break;
     }
     case OrEvaluateOperator:
     {
-      result=(MagickRealType) ((size_t) pixel | (size_t) (value+0.5));
+      result=(double) ((size_t) pixel | (size_t) (value+0.5));
       break;
     }
     case PoissonNoiseEvaluateOperator:
     {
-      result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
+      result=(double) GenerateDifferentialNoise(random_info,pixel,
         PoissonNoise,value);
       break;
     }
     case PowEvaluateOperator:
     {
-      result=(MagickRealType) (QuantumRange*pow((double) (QuantumScale*pixel),
+      result=(double) (QuantumRange*pow((double) (QuantumScale*pixel),
         (double) value));
       break;
     }
     case RightShiftEvaluateOperator:
     {
-      result=(MagickRealType) ((size_t) pixel >> (size_t) (value+0.5));
+      result=(double) ((size_t) pixel >> (size_t) (value+0.5));
       break;
     }
     case SetEvaluateOperator:
@@ -375,46 +377,46 @@ static MagickRealType ApplyEvaluateOperator(RandomInfo *random_info,
     }
     case SineEvaluateOperator:
     {
-      result=(MagickRealType) (QuantumRange*(0.5*sin((double) (2.0*MagickPI*
+      result=(double) (QuantumRange*(0.5*sin((double) (2.0*MagickPI*
         QuantumScale*pixel*value))+0.5));
       break;
     }
     case SubtractEvaluateOperator:
     {
-      result=(MagickRealType) (pixel-value);
+      result=(double) (pixel-value);
       break;
     }
     case SumEvaluateOperator:
     {
-      result=(MagickRealType) (pixel+value);
+      result=(double) (pixel+value);
       break;
     }
     case ThresholdEvaluateOperator:
     {
-      result=(MagickRealType) (((MagickRealType) pixel <= value) ? 0 :
+      result=(double) (((double) pixel <= value) ? 0 :
         QuantumRange);
       break;
     }
     case ThresholdBlackEvaluateOperator:
     {
-      result=(MagickRealType) (((MagickRealType) pixel <= value) ? 0 : pixel);
+      result=(double) (((double) pixel <= value) ? 0 : pixel);
       break;
     }
     case ThresholdWhiteEvaluateOperator:
     {
-      result=(MagickRealType) (((MagickRealType) pixel > value) ? QuantumRange :
+      result=(double) (((double) pixel > value) ? QuantumRange :
         pixel);
       break;
     }
     case UniformNoiseEvaluateOperator:
     {
-      result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
+      result=(double) GenerateDifferentialNoise(random_info,pixel,
         UniformNoise,value);
       break;
     }
     case XorEvaluateOperator:
     {
-      result=(MagickRealType) ((size_t) pixel ^ (size_t) (value+0.5));
+      result=(double) ((size_t) pixel ^ (size_t) (value+0.5));
       break;
     }
   }
@@ -433,10 +435,9 @@ MagickExport Image *EvaluateImages(const Image *images,
     *next;
 
   Image
-    *evaluate_image;
+    *image;
 
   MagickBooleanType
-    concurrent,
     status;
 
   MagickOffsetType
@@ -454,6 +455,9 @@ MagickExport Image *EvaluateImages(const Image *images,
   ssize_t
     y;
 
+  unsigned long
+    key;
+
   /*
     Ensure the image are the same size.
   */
@@ -467,28 +471,28 @@ MagickExport Image *EvaluateImages(const Image *images,
     if ((next->columns != images->columns) || (next->rows != images->rows))
       {
         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
-          "ImageWidthsOrHeightsDiffer","`%s'",images->filename);
+          "ImageWidthsOrHeightsDiffer","'%s'",images->filename);
         return((Image *) NULL);
       }
   /*
     Initialize evaluate next attributes.
   */
-  evaluate_image=CloneImage(images,images->columns,images->rows,MagickTrue,
+  image=CloneImage(images,images->columns,images->rows,MagickTrue,
     exception);
-  if (evaluate_image == (Image *) NULL)
+  if (image == (Image *) NULL)
     return((Image *) NULL);
-  if (SetImageStorageClass(evaluate_image,DirectClass,exception) == MagickFalse)
+  if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
     {
-      evaluate_image=DestroyImage(evaluate_image);
+      image=DestroyImage(image);
       return((Image *) NULL);
     }
   number_images=GetImageListLength(images);
   evaluate_pixels=AcquirePixelThreadSet(images,number_images);
   if (evaluate_pixels == (PixelChannels **) NULL)
     {
-      evaluate_image=DestroyImage(evaluate_image);
+      image=DestroyImage(image);
       (void) ThrowMagickException(exception,GetMagickModule(),
-        ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
+        ResourceLimitError,"MemoryAllocationFailed","'%s'",images->filename);
       return((Image *) NULL);
     }
   /*
@@ -497,15 +501,15 @@ MagickExport Image *EvaluateImages(const Image *images,
   status=MagickTrue;
   progress=0;
   random_info=AcquireRandomInfoThreadSet();
-  concurrent=GetRandomSecretKey(random_info[0]) == ~0UL ? MagickTrue :
-    MagickFalse;
-  evaluate_view=AcquireAuthenticCacheView(evaluate_image,exception);
+  key=GetRandomSecretKey(random_info[0]);
+  evaluate_view=AcquireAuthenticCacheView(image,exception);
   if (op == MedianEvaluateOperator)
     {
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-      #pragma omp parallel for schedule(static) shared(progress,status) omp_concurrent(concurrent)
+      #pragma omp parallel for schedule(static) shared(progress,status) \
+        dynamic_number_threads(image,image->columns,image->rows,key == ~0UL)
 #endif
-      for (y=0; y < (ssize_t) evaluate_image->rows; y++)
+      for (y=0; y < (ssize_t) image->rows; y++)
       {
         CacheView
           *image_view;
@@ -528,14 +532,14 @@ MagickExport Image *EvaluateImages(const Image *images,
         if (status == MagickFalse)
           continue;
         q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,
-          evaluate_image->columns,1,exception);
+          image->columns,1,exception);
         if (q == (Quantum *) NULL)
           {
             status=MagickFalse;
             continue;
           }
         evaluate_pixel=evaluate_pixels[id];
-        for (x=0; x < (ssize_t) evaluate_image->columns; x++)
+        for (x=0; x < (ssize_t) image->columns; x++)
         {
           register ssize_t
             j,
@@ -560,7 +564,7 @@ MagickExport Image *EvaluateImages(const Image *images,
                 image_view=DestroyCacheView(image_view);
                 break;
               }
-            for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++)
+            for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
             {
               PixelChannel
                 channel;
@@ -569,16 +573,16 @@ MagickExport Image *EvaluateImages(const Image *images,
                 evaluate_traits,
                 traits;
 
-              channel=GetPixelChannelMapChannel(evaluate_image,i);
-              evaluate_traits=GetPixelChannelMapTraits(evaluate_image,channel);
-              traits=GetPixelChannelMapTraits(next,channel);
+              channel=GetPixelChannelChannel(image,i);
+              evaluate_traits=GetPixelChannelTraits(image,channel);
+              traits=GetPixelChannelTraits(next,channel);
               if ((traits == UndefinedPixelTrait) ||
                   (evaluate_traits == UndefinedPixelTrait))
                 continue;
               if ((evaluate_traits & UpdatePixelTrait) == 0)
                 continue;
               evaluate_pixel[j].channel[i]=ApplyEvaluateOperator(
-                random_info[id],GetPixelChannel(evaluate_image,channel,p),op,
+                random_info[id],GetPixelChannel(image,channel,p),op,
                 evaluate_pixel[j].channel[i]);
             }
             image_view=DestroyCacheView(image_view);
@@ -586,9 +590,9 @@ MagickExport Image *EvaluateImages(const Image *images,
           }
           qsort((void *) evaluate_pixel,number_images,sizeof(*evaluate_pixel),
             IntensityCompare);
-          for (k=0; k < (ssize_t) GetPixelChannels(evaluate_image); k++)
+          for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
             q[k]=ClampToQuantum(evaluate_pixel[j/2].channel[k]);
-          q+=GetPixelChannels(evaluate_image);
+          q+=GetPixelChannels(image);
         }
         if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
           status=MagickFalse;
@@ -601,7 +605,7 @@ MagickExport Image *EvaluateImages(const Image *images,
             #pragma omp critical (MagickCore_EvaluateImages)
 #endif
             proceed=SetImageProgress(images,EvaluateImageTag,progress++,
-              evaluate_image->rows);
+              image->rows);
             if (proceed == MagickFalse)
               status=MagickFalse;
           }
@@ -610,9 +614,10 @@ MagickExport Image *EvaluateImages(const Image *images,
   else
     {
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-      #pragma omp parallel for schedule(static) shared(progress,status) omp_concurrent(concurrent)
+      #pragma omp parallel for schedule(static) shared(progress,status) \
+        dynamic_number_threads(image,image->columns,image->rows,key == ~0UL)
 #endif
-      for (y=0; y < (ssize_t) evaluate_image->rows; y++)
+      for (y=0; y < (ssize_t) image->rows; y++)
       {
         CacheView
           *image_view;
@@ -639,14 +644,14 @@ MagickExport Image *EvaluateImages(const Image *images,
         if (status == MagickFalse)
           continue;
         q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,
-          evaluate_image->columns,1,exception);
+          image->columns,1,exception);
         if (q == (Quantum *) NULL)
           {
             status=MagickFalse;
             continue;
           }
         evaluate_pixel=evaluate_pixels[id];
-        for (j=0; j < (ssize_t) evaluate_image->columns; j++)
+        for (j=0; j < (ssize_t) image->columns; j++)
           for (i=0; i < MaxPixelChannels; i++)
             evaluate_pixel[j].channel[i]=0.0;
         next=images;
@@ -681,16 +686,16 @@ MagickExport Image *EvaluateImages(const Image *images,
                 evaluate_traits,
                 traits;
 
-              channel=GetPixelChannelMapChannel(evaluate_image,i);
-              traits=GetPixelChannelMapTraits(next,channel);
-              evaluate_traits=GetPixelChannelMapTraits(evaluate_image,channel);
+              channel=GetPixelChannelChannel(image,i);
+              traits=GetPixelChannelTraits(next,channel);
+              evaluate_traits=GetPixelChannelTraits(image,channel);
               if ((traits == UndefinedPixelTrait) ||
                   (evaluate_traits == UndefinedPixelTrait))
                 continue;
               if ((traits & UpdatePixelTrait) == 0)
                 continue;
               evaluate_pixel[x].channel[i]=ApplyEvaluateOperator(
-                random_info[id],GetPixelChannel(evaluate_image,channel,p),j ==
+                random_info[id],GetPixelChannel(image,channel,p),j ==
                 0 ? AddEvaluateOperator : op,evaluate_pixel[x].channel[i]);
             }
             p+=GetPixelChannels(next);
@@ -698,7 +703,7 @@ MagickExport Image *EvaluateImages(const Image *images,
           image_view=DestroyCacheView(image_view);
           next=GetNextImageInList(next);
         }
-        for (x=0; x < (ssize_t) evaluate_image->columns; x++)
+        for (x=0; x < (ssize_t) image->columns; x++)
         {
           register ssize_t
              i;
@@ -707,13 +712,13 @@ MagickExport Image *EvaluateImages(const Image *images,
           {
             case MeanEvaluateOperator:
             {
-              for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++)
-                evaluate_pixel[x].channel[i]/=(MagickRealType) number_images;
+              for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
+                evaluate_pixel[x].channel[i]/=(double) number_images;
               break;
             }
             case MultiplyEvaluateOperator:
             {
-              for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++)
+              for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
               {
                 register ssize_t
                   j;
@@ -727,17 +732,17 @@ MagickExport Image *EvaluateImages(const Image *images,
               break;
           }
         }
-        for (x=0; x < (ssize_t) evaluate_image->columns; x++)
+        for (x=0; x < (ssize_t) image->columns; x++)
         {
           register ssize_t
             i;
 
-          if (GetPixelMask(evaluate_image,q) != 0)
+          if (GetPixelMask(image,q) != 0)
             {
-              q+=GetPixelChannels(evaluate_image);
+              q+=GetPixelChannels(image);
               continue;
             }
-          for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++)
+          for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
           {
             PixelChannel
               channel;
@@ -745,15 +750,15 @@ MagickExport Image *EvaluateImages(const Image *images,
             PixelTrait
               traits;
 
-            channel=GetPixelChannelMapChannel(evaluate_image,i);
-            traits=GetPixelChannelMapTraits(evaluate_image,channel);
+            channel=GetPixelChannelChannel(image,i);
+            traits=GetPixelChannelTraits(image,channel);
             if (traits == UndefinedPixelTrait)
               continue;
             if ((traits & UpdatePixelTrait) == 0)
               continue;
             q[i]=ClampToQuantum(evaluate_pixel[x].channel[i]);
           }
-          q+=GetPixelChannels(evaluate_image);
+          q+=GetPixelChannels(image);
         }
         if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
           status=MagickFalse;
@@ -766,7 +771,7 @@ MagickExport Image *EvaluateImages(const Image *images,
             #pragma omp critical (MagickCore_EvaluateImages)
 #endif
             proceed=SetImageProgress(images,EvaluateImageTag,progress++,
-              evaluate_image->rows);
+              image->rows);
             if (proceed == MagickFalse)
               status=MagickFalse;
           }
@@ -776,8 +781,8 @@ MagickExport Image *EvaluateImages(const Image *images,
   evaluate_pixels=DestroyPixelThreadSet(evaluate_pixels);
   random_info=DestroyRandomInfoThreadSet(random_info);
   if (status == MagickFalse)
-    evaluate_image=DestroyImage(evaluate_image);
-  return(evaluate_image);
+    image=DestroyImage(image);
+  return(image);
 }
 
 MagickExport MagickBooleanType EvaluateImage(Image *image,
@@ -787,7 +792,6 @@ MagickExport MagickBooleanType EvaluateImage(Image *image,
     *image_view;
 
   MagickBooleanType
-    concurrent,
     status;
 
   MagickOffsetType
@@ -799,6 +803,9 @@ MagickExport MagickBooleanType EvaluateImage(Image *image,
   ssize_t
     y;
 
+  unsigned long
+    key;
+
   assert(image != (Image *) NULL);
   assert(image->signature == MagickSignature);
   if (image->debug != MagickFalse)
@@ -810,11 +817,11 @@ MagickExport MagickBooleanType EvaluateImage(Image *image,
   status=MagickTrue;
   progress=0;
   random_info=AcquireRandomInfoThreadSet();
-  concurrent=GetRandomSecretKey(random_info[0]) == ~0UL ? MagickTrue :
-    MagickFalse;
+  key=GetRandomSecretKey(random_info[0]);
   image_view=AcquireAuthenticCacheView(image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(static,4) shared(progress,status) omp_concurrent(concurrent)
+  #pragma omp parallel for schedule(static,4) shared(progress,status) \
+    dynamic_number_threads(image,image->columns,image->rows,key == ~0UL)
 #endif
   for (y=0; y < (ssize_t) image->rows; y++)
   {
@@ -840,11 +847,6 @@ MagickExport MagickBooleanType EvaluateImage(Image *image,
       register ssize_t
         i;
 
-      if (GetPixelMask(image,q) != 0)
-        {
-          q+=GetPixelChannels(image);
-          continue;
-        }
       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
       {
         PixelChannel
@@ -853,11 +855,12 @@ MagickExport MagickBooleanType EvaluateImage(Image *image,
         PixelTrait
           traits;
 
-        channel=GetPixelChannelMapChannel(image,i);
-        traits=GetPixelChannelMapTraits(image,channel);
+        channel=GetPixelChannelChannel(image,i);
+        traits=GetPixelChannelTraits(image,channel);
         if (traits == UndefinedPixelTrait)
           continue;
-        if ((traits & CopyPixelTrait) != 0)
+        if (((traits & CopyPixelTrait) != 0) ||
+            (GetPixelMask(image,q) != 0))
           continue;
         q[i]=ClampToQuantum(ApplyEvaluateOperator(random_info[id],q[i],op,
           value));
@@ -872,7 +875,7 @@ MagickExport MagickBooleanType EvaluateImage(Image *image,
           proceed;
 
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_EvaluateImage)
+        #pragma omp critical (MagickCore_EvaluateImage)
 #endif
         proceed=SetImageProgress(image,EvaluateImageTag,progress++,image->rows);
         if (proceed == MagickFalse)
@@ -922,7 +925,7 @@ static Quantum ApplyFunction(Quantum pixel,const MagickFunction function,
   const size_t number_parameters,const double *parameters,
   ExceptionInfo *exception)
 {
-  MagickRealType
+  double
     result;
 
   register ssize_t
@@ -946,7 +949,7 @@ static Quantum ApplyFunction(Quantum pixel,const MagickFunction function,
     }
     case SinusoidFunction:
     {
-      MagickRealType
+      double
         amplitude,
         bias,
         frequency,
@@ -959,13 +962,13 @@ static Quantum ApplyFunction(Quantum pixel,const MagickFunction function,
       phase=(number_parameters >= 2) ? parameters[1] : 0.0;
       amplitude=(number_parameters >= 3) ? parameters[2] : 0.5;
       bias=(number_parameters >= 4) ? parameters[3] : 0.5;
-      result=(MagickRealType) (QuantumRange*(amplitude*sin((double) (2.0*
+      result=(double) (QuantumRange*(amplitude*sin((double) (2.0*
         MagickPI*(frequency*QuantumScale*pixel+phase/360.0)))+bias));
       break;
     }
     case ArcsinFunction:
     {
-      MagickRealType
+      double
         bias,
         center,
         range,
@@ -986,13 +989,13 @@ static Quantum ApplyFunction(Quantum pixel,const MagickFunction function,
         if (result >= 1.0)
           result=bias+range/2.0;
         else
-          result=(MagickRealType) (range/MagickPI*asin((double) result)+bias);
+          result=(double) (range/MagickPI*asin((double) result)+bias);
       result*=QuantumRange;
       break;
     }
     case ArctanFunction:
     {
-      MagickRealType
+      double
         center,
         bias,
         range,
@@ -1005,8 +1008,8 @@ static Quantum ApplyFunction(Quantum pixel,const MagickFunction function,
       center=(number_parameters >= 2) ? parameters[1] : 0.5;
       range=(number_parameters >= 3) ? parameters[2] : 1.0;
       bias=(number_parameters >= 4) ? parameters[3] : 0.5;
-      result=(MagickRealType) (MagickPI*slope*(QuantumScale*pixel-center));
-      result=(MagickRealType) (QuantumRange*(range/MagickPI*atan((double)
+      result=(double) (MagickPI*slope*(QuantumScale*pixel-center));
+      result=(double) (QuantumRange*(range/MagickPI*atan((double)
         result)+bias));
       break;
     }
@@ -1046,7 +1049,8 @@ MagickExport MagickBooleanType FunctionImage(Image *image,
   progress=0;
   image_view=AcquireAuthenticCacheView(image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(static,4) shared(progress,status)
+  #pragma omp parallel for schedule(static,4) shared(progress,status) \
+    dynamic_number_threads(image,image->columns,image->rows,1)
 #endif
   for (y=0; y < (ssize_t) image->rows; y++)
   {
@@ -1082,8 +1086,8 @@ MagickExport MagickBooleanType FunctionImage(Image *image,
         PixelTrait
           traits;
 
-        channel=GetPixelChannelMapChannel(image,i);
-        traits=GetPixelChannelMapTraits(image,channel);
+        channel=GetPixelChannelChannel(image,i);
+        traits=GetPixelChannelTraits(image,channel);
         if (traits == UndefinedPixelTrait)
           continue;
         if ((traits & UpdatePixelTrait) == 0)
@@ -1101,7 +1105,7 @@ MagickExport MagickBooleanType FunctionImage(Image *image,
           proceed;
 
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_FunctionImage)
+        #pragma omp critical (MagickCore_FunctionImage)
 #endif
         proceed=SetImageProgress(image,FunctionImageTag,progress++,image->rows);
         if (proceed == MagickFalse)
@@ -1172,8 +1176,8 @@ MagickExport MagickBooleanType GetImageExtrema(const Image *image,
 %                                                                             %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-%  GetImageMean() returns the mean and standard deviation of one or more
-%  image channels.
+%  GetImageMean() returns the mean and standard deviation of one or more image
+%  channels.
 %
 %  The format of the GetImageMean method is:
 %
@@ -1194,15 +1198,15 @@ MagickExport MagickBooleanType GetImageExtrema(const Image *image,
 MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean,
   double *standard_deviation,ExceptionInfo *exception)
 {
+  double
+    area;
+
   ChannelStatistics
     *channel_statistics;
 
   register ssize_t
     i;
 
-  size_t
-    area;
-
   assert(image != (Image *) NULL);
   assert(image->signature == MagickSignature);
   if (image->debug != MagickFalse)
@@ -1210,7 +1214,7 @@ MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean,
   channel_statistics=GetImageStatistics(image,exception);
   if (channel_statistics == (ChannelStatistics *) NULL)
     return(MagickFalse);
-  area=0;
+  area=0.0;
   channel_statistics[CompositePixelChannel].mean=0.0;
   channel_statistics[CompositePixelChannel].standard_deviation=0.0;
   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
@@ -1221,8 +1225,8 @@ MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean,
     PixelTrait
       traits;
 
-    channel=GetPixelChannelMapChannel(image,i);
-    traits=GetPixelChannelMapTraits(image,channel);
+    channel=GetPixelChannelChannel(image,i);
+    traits=GetPixelChannelTraits(image,channel);
     if (traits == UndefinedPixelTrait)
       continue;
     if ((traits & UpdatePixelTrait) == 0)
@@ -1309,7 +1313,8 @@ MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
   sum_fourth_power=0.0;
   image_view=AcquireVirtualCacheView(image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(static) shared(status)
+  #pragma omp parallel for schedule(static) shared(status) \
+    dynamic_number_threads(image,image->columns,image->rows,1)
 #endif
   for (y=0; y < (ssize_t) image->rows; y++)
   {
@@ -1332,11 +1337,6 @@ MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
       register ssize_t
         i;
 
-      if (GetPixelMask(image,p) != 0)
-        {
-          p+=GetPixelChannels(image);
-          continue;
-        }
       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
       {
         PixelChannel
@@ -1345,11 +1345,12 @@ MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
         PixelTrait
           traits;
 
-        channel=GetPixelChannelMapChannel(image,i);
-        traits=GetPixelChannelMapTraits(image,channel);
+        channel=GetPixelChannelChannel(image,i);
+        traits=GetPixelChannelTraits(image,channel);
         if (traits == UndefinedPixelTrait)
           continue;
-        if ((traits & UpdatePixelTrait) == 0)
+        if (((traits & UpdatePixelTrait) == 0) ||
+            (GetPixelMask(image,p) != 0))
           continue;
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
         #pragma omp critical (MagickCore_GetImageKurtosis)
@@ -1423,6 +1424,7 @@ MagickExport MagickBooleanType GetImageRange(const Image *image,double *minima,
     *image_view;
 
   MagickBooleanType
+    initialize,
     status;
 
   ssize_t
@@ -1433,11 +1435,13 @@ MagickExport MagickBooleanType GetImageRange(const Image *image,double *minima,
   if (image->debug != MagickFalse)
     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
   status=MagickTrue;
-  *maxima=(-MagickHuge);
-  *minima=MagickHuge;
+  initialize=MagickTrue;
+  *maxima=0.0;
+  *minima=0.0;
   image_view=AcquireVirtualCacheView(image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(static) shared(status)
+  #pragma omp parallel for schedule(static) shared(status,initialize) \
+    dynamic_number_threads(image,image->columns,image->rows,1)
 #endif
   for (y=0; y < (ssize_t) image->rows; y++)
   {
@@ -1473,8 +1477,8 @@ MagickExport MagickBooleanType GetImageRange(const Image *image,double *minima,
         PixelTrait
           traits;
 
-        channel=GetPixelChannelMapChannel(image,i);
-        traits=GetPixelChannelMapTraits(image,channel);
+        channel=GetPixelChannelChannel(image,i);
+        traits=GetPixelChannelTraits(image,channel);
         if (traits == UndefinedPixelTrait)
           continue;
         if ((traits & UpdatePixelTrait) == 0)
@@ -1483,10 +1487,19 @@ MagickExport MagickBooleanType GetImageRange(const Image *image,double *minima,
         #pragma omp critical (MagickCore_GetImageRange)
 #endif
         {
-          if ((double) p[i] < *minima)
-            *minima=(double) p[i];
-          if ((double) p[i] > *maxima)
-            *maxima=(double) p[i];
+          if (initialize != MagickFalse)
+            {
+              *minima=(double) p[i];
+              *maxima=(double) p[i];
+              initialize=MagickFalse;
+            }
+          else
+            {
+              if ((double) p[i] < *minima)
+                *minima=(double) p[i];
+              if ((double) p[i] > *maxima)
+                *maxima=(double) p[i];
+           }
         }
       }
       p+=GetPixelChannels(image);
@@ -1507,10 +1520,10 @@ MagickExport MagickBooleanType GetImageRange(const Image *image,double *minima,
 %                                                                             %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-%  GetImageStatistics() returns statistics for each channel in the
-%  image.  The statistics include the channel depth, its minima, maxima, mean,
-%  standard deviation, kurtosis and skewness.  You can access the red channel
-%  mean, for example, like this:
+%  GetImageStatistics() returns statistics for each channel in the image.  The
+%  statistics include the channel depth, its minima, maxima, mean, standard
+%  deviation, kurtosis and skewness.  You can access the red channel mean, for
+%  example, like this:
 %
 %      channel_statistics=GetImageStatistics(image,exception);
 %      red_mean=channel_statistics[RedPixelChannel].mean;
@@ -1547,8 +1560,8 @@ static size_t GetImageChannels(const Image *image)
     PixelTrait
       traits;
 
-    channel=GetPixelChannelMapChannel(image,i);
-    traits=GetPixelChannelMapTraits(image,channel);
+    channel=GetPixelChannelChannel(image,i);
+    traits=GetPixelChannelTraits(image,channel);
     if ((traits & UpdatePixelTrait) != 0)
       channels++;
   }
@@ -1561,9 +1574,6 @@ MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
   ChannelStatistics
     *channel_statistics;
 
-  double
-    area;
-
   MagickStatusType
     status;
 
@@ -1625,8 +1635,8 @@ MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
         PixelTrait
           traits;
 
-        channel=GetPixelChannelMapChannel(image,i);
-        traits=GetPixelChannelMapTraits(image,channel);
+        channel=GetPixelChannelChannel(image,i);
+        traits=GetPixelChannelTraits(image,channel);
         if (traits == UndefinedPixelTrait)
           continue;
         if (channel_statistics[channel].depth != MAGICKCORE_QUANTUM_DEPTH)
@@ -1638,6 +1648,7 @@ MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
             if (status != MagickFalse)
               {
                 channel_statistics[channel].depth++;
+                i--;
                 continue;
               }
           }
@@ -1650,17 +1661,21 @@ MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
         channel_statistics[channel].sum_cubed+=(double) p[i]*p[i]*p[i];
         channel_statistics[channel].sum_fourth_power+=(double) p[i]*p[i]*p[i]*
           p[i];
+        channel_statistics[channel].area++;
       }
       p+=GetPixelChannels(image);
     }
   }
-  area=(double) image->columns*image->rows;
   for (i=0; i < (ssize_t) MaxPixelChannels; i++)
   {
-    channel_statistics[i].sum/=area;
-    channel_statistics[i].sum_squared/=area;
-    channel_statistics[i].sum_cubed/=area;
-    channel_statistics[i].sum_fourth_power/=area;
+    double
+      area;
+
+    area=MagickEpsilonReciprocal(channel_statistics[i].area);
+    channel_statistics[i].sum*=area;
+    channel_statistics[i].sum_squared*=area;
+    channel_statistics[i].sum_cubed*=area;
+    channel_statistics[i].sum_fourth_power*=area;
     channel_statistics[i].mean=channel_statistics[i].sum;
     channel_statistics[i].variance=channel_statistics[i].sum_squared;
     channel_statistics[i].standard_deviation=sqrt(
@@ -1706,23 +1721,23 @@ MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
   channel_statistics[CompositePixelChannel].skewness/=channels;
   for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
   {
-    if (channel_statistics[i].standard_deviation == 0.0)
-      continue;
+    double
+      standard_deviation;
+
+    standard_deviation=MagickEpsilonReciprocal(
+      channel_statistics[i].standard_deviation);
     channel_statistics[i].skewness=(channel_statistics[i].sum_cubed-3.0*
       channel_statistics[i].mean*channel_statistics[i].sum_squared+2.0*
       channel_statistics[i].mean*channel_statistics[i].mean*
-      channel_statistics[i].mean)/(channel_statistics[i].standard_deviation*
-      channel_statistics[i].standard_deviation*
-      channel_statistics[i].standard_deviation);
+      channel_statistics[i].mean)*(standard_deviation*standard_deviation*
+      standard_deviation);
     channel_statistics[i].kurtosis=(channel_statistics[i].sum_fourth_power-4.0*
       channel_statistics[i].mean*channel_statistics[i].sum_cubed+6.0*
       channel_statistics[i].mean*channel_statistics[i].mean*
       channel_statistics[i].sum_squared-3.0*channel_statistics[i].mean*
       channel_statistics[i].mean*1.0*channel_statistics[i].mean*
-      channel_statistics[i].mean)/(channel_statistics[i].standard_deviation*
-      channel_statistics[i].standard_deviation*
-      channel_statistics[i].standard_deviation*
-      channel_statistics[i].standard_deviation)-3.0;
+      channel_statistics[i].mean)*(standard_deviation*standard_deviation*
+      standard_deviation*standard_deviation)-3.0;
   }
   return(channel_statistics);
 }
@@ -1807,7 +1822,7 @@ static PixelList **DestroyPixelListThreadSet(PixelList **pixel_list)
     i;
 
   assert(pixel_list != (PixelList **) NULL);
-  for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
+  for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
     if (pixel_list[i] != (PixelList *) NULL)
       pixel_list[i]=DestroyPixelList(pixel_list[i]);
   pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list);
@@ -1846,7 +1861,7 @@ static PixelList **AcquirePixelListThreadSet(const size_t width,
   size_t
     number_threads;
 
-  number_threads=GetOpenMPMaximumThreads();
+  number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
   pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
     sizeof(*pixel_list));
   if (pixel_list == (PixelList **) NULL)
@@ -1951,7 +1966,7 @@ static inline void GetMaximumPixelList(PixelList *pixel_list,Quantum *pixel)
 
 static inline void GetMeanPixelList(PixelList *pixel_list,Quantum *pixel)
 {
-  MagickRealType
+  double
     sum;
 
   register SkipList
@@ -1973,7 +1988,7 @@ static inline void GetMeanPixelList(PixelList *pixel_list,Quantum *pixel)
   do
   {
     color=p->nodes[color].next[0];
-    sum+=(MagickRealType) p->nodes[color].count*color;
+    sum+=(double) p->nodes[color].count*color;
     count+=p->nodes[color].count;
   } while (count < (ssize_t) pixel_list->length);
   sum/=pixel_list->length;
@@ -2106,7 +2121,7 @@ static inline void GetNonpeakPixelList(PixelList *pixel_list,Quantum *pixel)
 static inline void GetStandardDeviationPixelList(PixelList *pixel_list,
   Quantum *pixel)
 {
-  MagickRealType
+  double
     sum,
     sum_squared;
 
@@ -2133,9 +2148,9 @@ static inline void GetStandardDeviationPixelList(PixelList *pixel_list,
       i;
 
     color=p->nodes[color].next[0];
-    sum+=(MagickRealType) p->nodes[color].count*color;
+    sum+=(double) p->nodes[color].count*color;
     for (i=0; i < (ssize_t) p->nodes[color].count; i++)
-      sum_squared+=((MagickRealType) color)*((MagickRealType) color);
+      sum_squared+=((double) color)*((double) color);
     count+=p->nodes[color].count;
   } while (count < (ssize_t) pixel_list->length);
   sum/=pixel_list->length;
@@ -2162,7 +2177,7 @@ static inline void InsertPixelList(const Image *image,const Quantum pixel,
   AddNodePixelList(pixel_list,index);
 }
 
-static inline MagickRealType MagickAbsoluteValue(const MagickRealType x)
+static inline double MagickAbsoluteValue(const double x)
 {
   if (x < 0)
     return(-x);
@@ -2258,7 +2273,8 @@ MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
   image_view=AcquireVirtualCacheView(image,exception);
   statistic_view=AcquireAuthenticCacheView(statistic_image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(static,4) shared(progress,status)
+  #pragma omp parallel for schedule(static,4) shared(progress,status) \
+    dynamic_number_threads(image,image->columns,image->rows,1)
 #endif
   for (y=0; y < (ssize_t) statistic_image->rows; y++)
   {
@@ -2290,12 +2306,6 @@ MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
       register ssize_t
         i;
 
-      if (GetPixelMask(image,p) != 0)
-        {
-          p+=GetPixelChannels(image);
-          q+=GetPixelChannels(statistic_image);
-          continue;
-        }
       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
       {
         PixelChannel
@@ -2317,13 +2327,14 @@ MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
         ssize_t
           v;
 
-        channel=GetPixelChannelMapChannel(image,i);
-        traits=GetPixelChannelMapTraits(image,channel);
-        statistic_traits=GetPixelChannelMapTraits(statistic_image,channel);
+        channel=GetPixelChannelChannel(image,i);
+        traits=GetPixelChannelTraits(image,channel);
+        statistic_traits=GetPixelChannelTraits(statistic_image,channel);
         if ((traits == UndefinedPixelTrait) ||
             (statistic_traits == UndefinedPixelTrait))
           continue;
-        if ((statistic_traits & CopyPixelTrait) != 0)
+        if (((statistic_traits & CopyPixelTrait) != 0) ||
+            (GetPixelMask(image,p) != 0))
           {
             SetPixelChannel(statistic_image,channel,p[center+i],q);
             continue;
@@ -2343,14 +2354,14 @@ MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
         {
           case GradientStatistic:
           {
-            MagickRealType
+            double
               maximum,
               minimum;
 
             GetMinimumPixelList(pixel_list[id],&pixel);
-            minimum=(MagickRealType) pixel;
+            minimum=(double) pixel;
             GetMaximumPixelList(pixel_list[id],&pixel);
-            maximum=(MagickRealType) pixel;
+            maximum=(double) pixel;
             pixel=ClampToQuantum(MagickAbsoluteValue(maximum-minimum));
             break;
           }
@@ -2404,7 +2415,7 @@ MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
           proceed;
 
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_StatisticImage)
+        #pragma omp critical (MagickCore_StatisticImage)
 #endif
         proceed=SetImageProgress(image,StatisticImageTag,progress++,
           image->rows);