]> granicus.if.org Git - imagemagick/blobdiff - magick/statistic.c
(no commit message)
[imagemagick] / magick / statistic.c
index 1dffd9e43e4846f123ffc2aba82eb33302013327..5f4e743d2630c03c2c15642dea90e25497f3e8eb 100644 (file)
 %        SSSSS    T    A   A    T    IIIII  SSSSS    T    IIIII   CCCC        %
 %                                                                             %
 %                                                                             %
-%                          MagickCore Image Methods                           %
+%                     MagickCore Image Statistical Methods                    %
 %                                                                             %
 %                              Software Design                                %
 %                                John Cristy                                  %
 %                                 July 1992                                   %
 %                                                                             %
 %                                                                             %
-%  Copyright 1999-2009 ImageMagick Studio LLC, a non-profit organization      %
+%  Copyright 1999-2010 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  %
@@ -79,6 +79,7 @@
 #include "magick/profile.h"
 #include "magick/quantize.h"
 #include "magick/random_.h"
+#include "magick/random-private.h"
 #include "magick/segment.h"
 #include "magick/semaphore.h"
 #include "magick/signature-private.h"
 %                                                                             %
 %                                                                             %
 %                                                                             %
-%     A v e r a g e I m a g e s                                               %
+%     E v a l u a t e I m a g e                                               %
 %                                                                             %
 %                                                                             %
 %                                                                             %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-%  AverageImages() takes a set of images and averages them together.  Each
-%  image in the set must have the same width and height.  AverageImages()
-%  returns a single image with each corresponding pixel component of each
-%  image averaged.   On failure, a NULL image is returned and exception
-%  describes the reason for the failure.
+%  EvaluateImage() applies a value to the image with an arithmetic, relational,
+%  or logical operator to an image. Use these operations to lighten or darken
+%  an image, to increase or decrease contrast in an image, or to produce the
+%  "negative" of an image.
 %
-%  The format of the AverageImages method is:
+%  The format of the EvaluateImageChannel method is:
 %
-%      Image *AverageImages(Image *image,ExceptionInfo *exception)
+%      MagickBooleanType EvaluateImage(Image *image,
+%        const MagickEvaluateOperator op,const double value,
+%        ExceptionInfo *exception)
+%      MagickBooleanType EvaluateImages(Image *images,
+%        const MagickEvaluateOperator op,const double value,
+%        ExceptionInfo *exception)
+%      MagickBooleanType EvaluateImageChannel(Image *image,
+%        const ChannelType channel,const MagickEvaluateOperator op,
+%        const double value,ExceptionInfo *exception)
 %
 %  A description of each parameter follows:
 %
-%    o image: the image sequence.
+%    o image: the image.
+%
+%    o channel: the channel.
+%
+%    o op: A channel op.
+%
+%    o value: A value value.
 %
 %    o exception: return any errors or warnings in this structure.
 %
 
 static MagickPixelPacket **DestroyPixelThreadSet(MagickPixelPacket **pixels)
 {
-  register long
+  register ssize_t
     i;
 
   assert(pixels != (MagickPixelPacket **) NULL);
-  for (i=0; i < (long) GetOpenMPMaximumThreads(); i++)
+  for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
     if (pixels[i] != (MagickPixelPacket *) NULL)
       pixels[i]=(MagickPixelPacket *) RelinquishMagickMemory(pixels[i]);
   pixels=(MagickPixelPacket **) RelinquishAlignedMemory(pixels);
@@ -133,14 +147,14 @@ static MagickPixelPacket **DestroyPixelThreadSet(MagickPixelPacket **pixels)
 
 static MagickPixelPacket **AcquirePixelThreadSet(const Image *image)
 {
-  register long
+  register ssize_t
     i,
     j;
 
   MagickPixelPacket
     **pixels;
 
-  unsigned long
+  size_t
     number_threads;
 
   number_threads=GetOpenMPMaximumThreads();
@@ -149,88 +163,304 @@ static MagickPixelPacket **AcquirePixelThreadSet(const Image *image)
   if (pixels == (MagickPixelPacket **) NULL)
     return((MagickPixelPacket **) NULL);
   (void) ResetMagickMemory(pixels,0,number_threads*sizeof(*pixels));
-  for (i=0; i < (long) number_threads; i++)
+  for (i=0; i < (ssize_t) number_threads; i++)
   {
     pixels[i]=(MagickPixelPacket *) AcquireQuantumMemory(image->columns,
       sizeof(**pixels));
     if (pixels[i] == (MagickPixelPacket *) NULL)
       return(DestroyPixelThreadSet(pixels));
-    for (j=0; j < (long) image->columns; j++)
+    for (j=0; j < (ssize_t) image->columns; j++)
       GetMagickPixelPacket(image,&pixels[i][j]);
   }
   return(pixels);
 }
 
-MagickExport Image *AverageImages(const Image *image,ExceptionInfo *exception)
+static inline double MagickMax(const double x,const double y)
+{
+  if (x > y)
+    return(x);
+  return(y);
+}
+
+static inline double MagickMin(const double x,const double y)
+{
+  if (x < y)
+    return(x);
+  return(y);
+}
+
+static MagickRealType ApplyEvaluateOperator(RandomInfo *random_info,
+  Quantum pixel,const MagickEvaluateOperator op,const MagickRealType value)
+{
+  MagickRealType
+    result;
+
+  result=0.0;
+  switch (op)
+  {
+    case UndefinedEvaluateOperator:
+      break;
+    case AbsEvaluateOperator:
+    {
+      result=(MagickRealType) fabs((double) (pixel+value));
+      break;
+    }
+    case AddEvaluateOperator:
+    {
+      result=(MagickRealType) (pixel+value);
+      break;
+    }
+    case AddModulusEvaluateOperator:
+    {
+      /*
+        This returns a 'floored modulus' of the addition which is a
+        positive result.  It differs from  % or fmod() which returns a
+        'truncated modulus' result, where floor() is replaced by trunc()
+        and could return a negative result (which is clipped).
+      */
+      result=pixel+value;
+      result-=(QuantumRange+1.0)*floor((double) result/(QuantumRange+1.0));
+      break;
+    }
+    case AndEvaluateOperator:
+    {
+      result=(MagickRealType) ((size_t) pixel & (size_t)
+        (value+0.5));
+      break;
+    }
+    case CosineEvaluateOperator:
+    {
+      result=(MagickRealType) (QuantumRange*(0.5*cos((double) (2.0*MagickPI*
+        QuantumScale*pixel*value))+0.5));
+      break;
+    }
+    case DivideEvaluateOperator:
+    {
+      result=pixel/(value == 0.0 ? 1.0 : value);
+      break;
+    }
+    case GaussianNoiseEvaluateOperator:
+    {
+      result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
+        GaussianNoise,value);
+      break;
+    }
+    case ImpulseNoiseEvaluateOperator:
+    {
+      result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
+        ImpulseNoise,value);
+      break;
+    }
+    case LaplacianNoiseEvaluateOperator:
+    {
+      result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
+        LaplacianNoise,value);
+      break;
+    }
+    case LeftShiftEvaluateOperator:
+    {
+      result=(MagickRealType) ((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)));
+      break;
+    }
+    case MaxEvaluateOperator:
+    {
+      result=(MagickRealType) MagickMax((double) pixel,value);
+      break;
+    }
+    case MeanEvaluateOperator:
+    {
+      result=(MagickRealType) (pixel+value);
+      break;
+    }
+    case MinEvaluateOperator:
+    {
+      result=(MagickRealType) MagickMin((double) pixel,value);
+      break;
+    }
+    case MultiplicativeNoiseEvaluateOperator:
+    {
+      result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
+        MultiplicativeGaussianNoise,value);
+      break;
+    }
+    case MultiplyEvaluateOperator:
+    {
+      result=(MagickRealType) (value*pixel);
+      break;
+    }
+    case OrEvaluateOperator:
+    {
+      result=(MagickRealType) ((size_t) pixel | (size_t)
+        (value+0.5));
+      break;
+    }
+    case PoissonNoiseEvaluateOperator:
+    {
+      result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
+        PoissonNoise,value);
+      break;
+    }
+    case PowEvaluateOperator:
+    {
+      result=(MagickRealType) (QuantumRange*pow((double) (QuantumScale*pixel),
+        (double) value));
+      break;
+    }
+    case RightShiftEvaluateOperator:
+    {
+      result=(MagickRealType) ((size_t) pixel >> (size_t)
+        (value+0.5));
+      break;
+    }
+    case SetEvaluateOperator:
+    {
+      result=value;
+      break;
+    }
+    case SineEvaluateOperator:
+    {
+      result=(MagickRealType) (QuantumRange*(0.5*sin((double) (2.0*MagickPI*
+        QuantumScale*pixel*value))+0.5));
+      break;
+    }
+    case SubtractEvaluateOperator:
+    {
+      result=(MagickRealType) (pixel-value);
+      break;
+    }
+    case ThresholdEvaluateOperator:
+    {
+      result=(MagickRealType) (((MagickRealType) pixel <= value) ? 0 :
+        QuantumRange);
+      break;
+    }
+    case ThresholdBlackEvaluateOperator:
+    {
+      result=(MagickRealType) (((MagickRealType) pixel <= value) ? 0 : pixel);
+      break;
+    }
+    case ThresholdWhiteEvaluateOperator:
+    {
+      result=(MagickRealType) (((MagickRealType) pixel > value) ? QuantumRange :
+        pixel);
+      break;
+    }
+    case UniformNoiseEvaluateOperator:
+    {
+      result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
+        UniformNoise,value);
+      break;
+    }
+    case XorEvaluateOperator:
+    {
+      result=(MagickRealType) ((size_t) pixel ^ (size_t)
+        (value+0.5));
+      break;
+    }
+  }
+  return(result);
+}
+
+MagickExport MagickBooleanType EvaluateImage(Image *image,
+  const MagickEvaluateOperator op,const double value,ExceptionInfo *exception)
+{
+  MagickBooleanType
+    status;
+
+  status=EvaluateImageChannel(image,AllChannels,op,value,exception);
+  return(status);
+}
+
+MagickExport Image *EvaluateImages(const Image *images,
+  const MagickEvaluateOperator op,ExceptionInfo *exception)
 {
-#define AverageImageTag  "Average/Image"
+#define EvaluateImageTag  "Evaluate/Image"
 
   CacheView
-    *average_view;
+    *evaluate_view;
 
   const Image
     *next;
 
   Image
-    *average_image;
-
-  long
-    progress,
-    y;
+    *evaluate_image;
 
   MagickBooleanType
     status;
 
+  MagickOffsetType
+    progress;
+
   MagickPixelPacket
-    **average_pixels,
+    **restrict evaluate_pixels,
     zero;
 
-  unsigned long
+  RandomInfo
+    **restrict random_info;
+
+  size_t
     number_images;
 
+  ssize_t
+    y;
+
   /*
     Ensure the image are the same size.
   */
-  assert(image != (Image *) NULL);
-  assert(image->signature == MagickSignature);
-  if (image->debug != MagickFalse)
-    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
+  assert(images != (Image *) NULL);
+  assert(images->signature == MagickSignature);
+  if (images->debug != MagickFalse)
+    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
   assert(exception != (ExceptionInfo *) NULL);
   assert(exception->signature == MagickSignature);
-  for (next=image; next != (Image *) NULL; next=GetNextImageInList(next))
-    if ((next->columns != image->columns) || (next->rows != image->rows))
-      ThrowImageException(OptionError,"ImageWidthsOrHeightsDiffer");
+  for (next=images; next != (Image *) NULL; next=GetNextImageInList(next))
+    if ((next->columns != images->columns) || (next->rows != images->rows))
+      {
+        (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
+          "ImageWidthsOrHeightsDiffer","`%s'",images->filename);
+        return((Image *) NULL);
+      }
   /*
-    Initialize average next attributes.
+    Initialize evaluate next attributes.
   */
-  average_image=CloneImage(image,image->columns,image->rows,MagickTrue,
+  evaluate_image=CloneImage(images,images->columns,images->rows,MagickTrue,
     exception);
-  if (average_image == (Image *) NULL)
+  if (evaluate_image == (Image *) NULL)
     return((Image *) NULL);
-  if (SetImageStorageClass(average_image,DirectClass) == MagickFalse)
+  if (SetImageStorageClass(evaluate_image,DirectClass) == MagickFalse)
     {
-      InheritException(exception,&average_image->exception);
-      average_image=DestroyImage(average_image);
+      InheritException(exception,&evaluate_image->exception);
+      evaluate_image=DestroyImage(evaluate_image);
       return((Image *) NULL);
     }
-  average_pixels=AcquirePixelThreadSet(image);
-  if (average_pixels == (MagickPixelPacket **) NULL)
+  evaluate_pixels=AcquirePixelThreadSet(images);
+  if (evaluate_pixels == (MagickPixelPacket **) NULL)
     {
-      average_image=DestroyImage(average_image);
-      ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
+      evaluate_image=DestroyImage(evaluate_image);
+      (void) ThrowMagickException(exception,GetMagickModule(),
+        ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
+      return((Image *) NULL);
     }
   /*
-    Average image pixels.
+    Evaluate image pixels.
   */
   status=MagickTrue;
   progress=0;
-  GetMagickPixelPacket(image,&zero);
-  number_images=GetImageListLength(image);
-  average_view=AcquireCacheView(average_image);
+  GetMagickPixelPacket(images,&zero);
+  random_info=AcquireRandomInfoThreadSet();
+  number_images=GetImageListLength(images);
+  evaluate_view=AcquireCacheView(evaluate_image);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
   #pragma omp parallel for schedule(dynamic) shared(progress,status)
 #endif
-  for (y=0; y < (long) average_image->rows; y++)
+  for (y=0; y < (ssize_t) evaluate_image->rows; y++)
   {
     CacheView
       *image_view;
@@ -238,40 +468,42 @@ MagickExport Image *AverageImages(const Image *image,ExceptionInfo *exception)
     const Image
       *next;
 
+    int
+      id;
+
     MagickPixelPacket
       pixel;
 
     register IndexPacket
-      *restrict average_indexes;
+      *restrict evaluate_indexes;
 
-    register long
+    register ssize_t
       i,
-      id,
       x;
 
     register MagickPixelPacket
-      *average_pixel;
+      *evaluate_pixel;
 
     register PixelPacket
       *restrict q;
 
     if (status == MagickFalse)
       continue;
-    q=QueueCacheViewAuthenticPixels(average_view,0,y,average_image->columns,1,
+    q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,evaluate_image->columns,1,
       exception);
     if (q == (PixelPacket *) NULL)
       {
         status=MagickFalse;
         continue;
       }
-    average_indexes=GetCacheViewAuthenticIndexQueue(average_view);
+    evaluate_indexes=GetCacheViewAuthenticIndexQueue(evaluate_view);
     pixel=zero;
     id=GetOpenMPThreadId();
-    average_pixel=average_pixels[id];
-    for (x=0; x < (long) average_image->columns; x++)
-      average_pixel[x]=zero;
-    next=image;
-    for (i=0; i < (long) number_images; i++)
+    evaluate_pixel=evaluate_pixels[id];
+    for (x=0; x < (ssize_t) evaluate_image->columns; x++)
+      evaluate_pixel[x]=zero;
+    next=images;
+    for (i=0; i < (ssize_t) number_images; i++)
     {
       register const IndexPacket
         *indexes;
@@ -287,57 +519,411 @@ MagickExport Image *AverageImages(const Image *image,ExceptionInfo *exception)
           break;
         }
       indexes=GetCacheViewVirtualIndexQueue(image_view);
-      for (x=0; x < (long) next->columns; x++)
+      for (x=0; x < (ssize_t) next->columns; x++)
       {
-        SetMagickPixelPacket(next,p,indexes+x,&pixel);
-        average_pixel[x].red+=QuantumScale*pixel.red;
-        average_pixel[x].green+=QuantumScale*pixel.green;
-        average_pixel[x].blue+=QuantumScale*pixel.blue;
-        average_pixel[x].opacity+=QuantumScale*pixel.opacity;
-        if (average_image->colorspace == CMYKColorspace)
-          average_pixel[x].index+=QuantumScale*pixel.index;
+        evaluate_pixel[x].red=ApplyEvaluateOperator(random_info[id],p->red,
+          i == 0 ? AddEvaluateOperator : op,evaluate_pixel[x].red);
+        evaluate_pixel[x].green=ApplyEvaluateOperator(random_info[id],p->green,
+          i == 0 ? AddEvaluateOperator : op,evaluate_pixel[x].green);
+        evaluate_pixel[x].blue=ApplyEvaluateOperator(random_info[id],p->blue,
+          i == 0 ? AddEvaluateOperator : op,evaluate_pixel[x].blue);
+        evaluate_pixel[x].opacity=ApplyEvaluateOperator(random_info[id],
+          p->opacity,i == 0 ? AddEvaluateOperator : op,
+          evaluate_pixel[x].opacity);
+        if (evaluate_image->colorspace == CMYKColorspace)
+          evaluate_pixel[x].index=ApplyEvaluateOperator(random_info[id],
+            indexes[x],i == 0 ? AddEvaluateOperator : op,
+            evaluate_pixel[x].index);
         p++;
       }
       image_view=DestroyCacheView(image_view);
       next=GetNextImageInList(next);
     }
-    for (x=0; x < (long) average_image->columns; x++)
+    if (op == MeanEvaluateOperator)
+      for (x=0; x < (ssize_t) evaluate_image->columns; x++)
+      {
+        evaluate_pixel[x].red/=number_images;
+        evaluate_pixel[x].green/=number_images;
+        evaluate_pixel[x].blue/=number_images;
+        evaluate_pixel[x].opacity/=number_images;
+        evaluate_pixel[x].index/=number_images;
+      }
+    for (x=0; x < (ssize_t) evaluate_image->columns; x++)
     {
-      average_pixel[x].red=(MagickRealType) (QuantumRange*
-        average_pixel[x].red/number_images);
-      average_pixel[x].green=(MagickRealType) (QuantumRange*
-        average_pixel[x].green/number_images);
-      average_pixel[x].blue=(MagickRealType) (QuantumRange*
-        average_pixel[x].blue/number_images);
-      average_pixel[x].opacity=(MagickRealType) (QuantumRange*
-        average_pixel[x].opacity/number_images);
-      if (average_image->colorspace == CMYKColorspace)
-        average_pixel[x].index=(MagickRealType) (QuantumRange*
-          average_pixel[x].index/number_images);
-      SetPixelPacket(average_image,&average_pixel[x],q,average_indexes+x);
+      q->red=ClampToQuantum(evaluate_pixel[x].red);
+      q->green=ClampToQuantum(evaluate_pixel[x].green);
+      q->blue=ClampToQuantum(evaluate_pixel[x].blue);
+      if (evaluate_image->matte == MagickFalse)
+        q->opacity=ClampToQuantum(evaluate_pixel[x].opacity);
+      else
+        q->opacity=ClampToQuantum(QuantumRange-evaluate_pixel[x].opacity);
+      if (evaluate_image->colorspace == CMYKColorspace)
+        evaluate_indexes[x]=ClampToQuantum(evaluate_pixel[x].index);
       q++;
     }
-    if (SyncCacheViewAuthenticPixels(average_view,exception) == MagickFalse)
+    if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
       status=MagickFalse;
-    if (image->progress_monitor != (MagickProgressMonitor) NULL)
+    if (images->progress_monitor != (MagickProgressMonitor) NULL)
       {
         MagickBooleanType
           proceed;
 
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-        #pragma omp critical (MagickCore_AverageImages)
+        #pragma omp critical (MagickCore_EvaluateImages)
 #endif
-        proceed=SetImageProgress(image,AverageImageTag,progress++,
-          average_image->rows);
+        proceed=SetImageProgress(images,EvaluateImageTag,progress++,
+          evaluate_image->rows);
         if (proceed == MagickFalse)
           status=MagickFalse;
       }
   }
-  average_view=DestroyCacheView(average_view);
-  average_pixels=DestroyPixelThreadSet(average_pixels);
+  evaluate_view=DestroyCacheView(evaluate_view);
+  evaluate_pixels=DestroyPixelThreadSet(evaluate_pixels);
+  random_info=DestroyRandomInfoThreadSet(random_info);
   if (status == MagickFalse)
-    average_image=DestroyImage(average_image);
-  return(average_image);
+    evaluate_image=DestroyImage(evaluate_image);
+  return(evaluate_image);
+}
+
+MagickExport MagickBooleanType EvaluateImageChannel(Image *image,
+  const ChannelType channel,const MagickEvaluateOperator op,const double value,
+  ExceptionInfo *exception)
+{
+  CacheView
+    *image_view;
+
+  MagickBooleanType
+    status;
+
+  MagickOffsetType
+    progress;
+
+  RandomInfo
+    **restrict random_info;
+
+  ssize_t
+    y;
+
+  assert(image != (Image *) NULL);
+  assert(image->signature == MagickSignature);
+  if (image->debug != MagickFalse)
+    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
+  assert(exception != (ExceptionInfo *) NULL);
+  assert(exception->signature == MagickSignature);
+  if (SetImageStorageClass(image,DirectClass) == MagickFalse)
+    {
+      InheritException(exception,&image->exception);
+      return(MagickFalse);
+    }
+  status=MagickTrue;
+  progress=0;
+  random_info=AcquireRandomInfoThreadSet();
+  image_view=AcquireCacheView(image);
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+  #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
+#endif
+  for (y=0; y < (ssize_t) image->rows; y++)
+  {
+    int
+      id;
+
+    register IndexPacket
+      *restrict indexes;
+
+    register ssize_t
+      x;
+
+    register PixelPacket
+      *restrict q;
+
+    if (status == MagickFalse)
+      continue;
+    q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
+    if (q == (PixelPacket *) NULL)
+      {
+        status=MagickFalse;
+        continue;
+      }
+    indexes=GetCacheViewAuthenticIndexQueue(image_view);
+    id=GetOpenMPThreadId();
+    for (x=0; x < (ssize_t) image->columns; x++)
+    {
+      if ((channel & RedChannel) != 0)
+        q->red=ClampToQuantum(ApplyEvaluateOperator(random_info[id],q->red,op,
+          value));
+      if ((channel & GreenChannel) != 0)
+        q->green=ClampToQuantum(ApplyEvaluateOperator(random_info[id],q->green,
+          op,value));
+      if ((channel & BlueChannel) != 0)
+        q->blue=ClampToQuantum(ApplyEvaluateOperator(random_info[id],q->blue,op,
+          value));
+      if ((channel & OpacityChannel) != 0)
+        {
+          if (image->matte == MagickFalse)
+            q->opacity=ClampToQuantum(ApplyEvaluateOperator(random_info[id],
+              q->opacity,op,value));
+          else
+            q->opacity=ClampToQuantum(QuantumRange-ApplyEvaluateOperator(
+              random_info[id],(Quantum) GetAlphaPixelComponent(q),op,value));
+        }
+      if (((channel & IndexChannel) != 0) && (indexes != (IndexPacket *) NULL))
+        indexes[x]=(IndexPacket) ClampToQuantum(ApplyEvaluateOperator(
+          random_info[id],indexes[x],op,value));
+      q++;
+    }
+    if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
+      status=MagickFalse;
+    if (image->progress_monitor != (MagickProgressMonitor) NULL)
+      {
+        MagickBooleanType
+          proceed;
+
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+  #pragma omp critical (MagickCore_EvaluateImageChannel)
+#endif
+        proceed=SetImageProgress(image,EvaluateImageTag,progress++,image->rows);
+        if (proceed == MagickFalse)
+          status=MagickFalse;
+      }
+  }
+  image_view=DestroyCacheView(image_view);
+  random_info=DestroyRandomInfoThreadSet(random_info);
+  return(status);
+}
+\f
+/*
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                             %
+%                                                                             %
+%                                                                             %
+%     F u n c t i o n I m a g e                                               %
+%                                                                             %
+%                                                                             %
+%                                                                             %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  FunctionImage() applies a value to the image with an arithmetic, relational,
+%  or logical operator to an image. Use these operations to lighten or darken
+%  an image, to increase or decrease contrast in an image, or to produce the
+%  "negative" of an image.
+%
+%  The format of the FunctionImageChannel method is:
+%
+%      MagickBooleanType FunctionImage(Image *image,
+%        const MagickFunction function,const ssize_t number_parameters,
+%        const double *parameters,ExceptionInfo *exception)
+%      MagickBooleanType FunctionImageChannel(Image *image,
+%        const ChannelType channel,const MagickFunction function,
+%        const ssize_t number_parameters,const double *argument,
+%        ExceptionInfo *exception)
+%
+%  A description of each parameter follows:
+%
+%    o image: the image.
+%
+%    o channel: the channel.
+%
+%    o function: A channel function.
+%
+%    o parameters: one or more parameters.
+%
+%    o exception: return any errors or warnings in this structure.
+%
+*/
+
+static Quantum ApplyFunction(Quantum pixel,const MagickFunction function,
+  const size_t number_parameters,const double *parameters,
+  ExceptionInfo *exception)
+{
+  MagickRealType
+    result;
+
+  register ssize_t
+    i;
+
+  (void) exception;
+  result=0.0;
+  switch (function)
+  {
+    case PolynomialFunction:
+    {
+      /*
+       * Polynomial
+       * Parameters:   polynomial constants,  highest to lowest order
+       *   For example:      c0*x^3 + c1*x^2 + c2*x  + c3
+       */
+      result=0.0;
+      for (i=0; i < (ssize_t) number_parameters; i++)
+        result = result*QuantumScale*pixel + parameters[i];
+      result *= QuantumRange;
+      break;
+    }
+    case SinusoidFunction:
+    {
+      /* Sinusoid Function
+       * Parameters:   Freq, Phase, Ampl, bias
+       */
+      double  freq,phase,ampl,bias;
+      freq  = ( number_parameters >= 1 ) ? parameters[0] : 1.0;
+      phase = ( number_parameters >= 2 ) ? parameters[1] : 0.0;
+      ampl  = ( number_parameters >= 3 ) ? parameters[2] : 0.5;
+      bias  = ( number_parameters >= 4 ) ? parameters[3] : 0.5;
+      result=(MagickRealType) (QuantumRange*(ampl*sin((double) (2.0*MagickPI*
+        (freq*QuantumScale*pixel + phase/360.0) )) + bias ) );
+      break;
+    }
+    case ArcsinFunction:
+    {
+      /* Arcsin Function  (peged at range limits for invalid results)
+       * Parameters:   Width, Center, Range, Bias
+       */
+      double  width,range,center,bias;
+      width  = ( number_parameters >= 1 ) ? parameters[0] : 1.0;
+      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 = 2.0/width*(QuantumScale*pixel - center);
+      if ( result <= -1.0 )
+        result = bias - range/2.0;
+      else if ( result >= 1.0 )
+        result = bias + range/2.0;
+      else
+        result=range/MagickPI*asin((double)result) + bias;
+      result *= QuantumRange;
+      break;
+    }
+    case ArctanFunction:
+    {
+      /* Arctan Function
+       * Parameters:   Slope, Center, Range, Bias
+       */
+      double  slope,range,center,bias;
+      slope  = ( number_parameters >= 1 ) ? parameters[0] : 1.0;
+      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 = MagickPI*slope*(QuantumScale*pixel - center);
+      result=(MagickRealType) (QuantumRange*(range/MagickPI*atan((double)
+                  result) + bias ) );
+      break;
+    }
+    case UndefinedFunction:
+      break;
+  }
+  return(ClampToQuantum(result));
+}
+
+MagickExport MagickBooleanType FunctionImage(Image *image,
+  const MagickFunction function,const size_t number_parameters,
+  const double *parameters,ExceptionInfo *exception)
+{
+  MagickBooleanType
+    status;
+
+  status=FunctionImageChannel(image,AllChannels,function,number_parameters,
+    parameters,exception);
+  return(status);
+}
+
+MagickExport MagickBooleanType FunctionImageChannel(Image *image,
+  const ChannelType channel,const MagickFunction function,
+  const size_t number_parameters,const double *parameters,
+  ExceptionInfo *exception)
+{
+#define FunctionImageTag  "Function/Image "
+
+  CacheView
+    *image_view;
+
+  MagickBooleanType
+    status;
+
+  MagickOffsetType
+    progress;
+
+  ssize_t
+    y;
+
+  assert(image != (Image *) NULL);
+  assert(image->signature == MagickSignature);
+  if (image->debug != MagickFalse)
+    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
+  assert(exception != (ExceptionInfo *) NULL);
+  assert(exception->signature == MagickSignature);
+  if (SetImageStorageClass(image,DirectClass) == MagickFalse)
+    {
+      InheritException(exception,&image->exception);
+      return(MagickFalse);
+    }
+  status=MagickTrue;
+  progress=0;
+  image_view=AcquireCacheView(image);
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+  #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
+#endif
+  for (y=0; y < (ssize_t) image->rows; y++)
+  {
+    register IndexPacket
+      *restrict indexes;
+
+    register ssize_t
+      x;
+
+    register PixelPacket
+      *restrict q;
+
+    if (status == MagickFalse)
+      continue;
+    q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
+    if (q == (PixelPacket *) NULL)
+      {
+        status=MagickFalse;
+        continue;
+      }
+    indexes=GetCacheViewAuthenticIndexQueue(image_view);
+    for (x=0; x < (ssize_t) image->columns; x++)
+    {
+      if ((channel & RedChannel) != 0)
+        q->red=ApplyFunction(q->red,function,number_parameters,parameters,
+          exception);
+      if ((channel & GreenChannel) != 0)
+        q->green=ApplyFunction(q->green,function,number_parameters,parameters,
+          exception);
+      if ((channel & BlueChannel) != 0)
+        q->blue=ApplyFunction(q->blue,function,number_parameters,parameters,
+          exception);
+      if ((channel & OpacityChannel) != 0)
+        {
+          if (image->matte == MagickFalse)
+            q->opacity=ApplyFunction(q->opacity,function,number_parameters,
+              parameters,exception);
+          else
+            q->opacity=(Quantum) QuantumRange-ApplyFunction((Quantum)
+              GetAlphaPixelComponent(q),function,number_parameters,parameters,
+              exception);
+        }
+      if (((channel & IndexChannel) != 0) && (indexes != (IndexPacket *) NULL))
+        indexes[x]=(IndexPacket) ApplyFunction(indexes[x],function,
+          number_parameters,parameters,exception);
+      q++;
+    }
+    if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
+      status=MagickFalse;
+    if (image->progress_monitor != (MagickProgressMonitor) NULL)
+      {
+        MagickBooleanType
+          proceed;
+
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+  #pragma omp critical (MagickCore_FunctionImageChannel)
+#endif
+        proceed=SetImageProgress(image,FunctionImageTag,progress++,image->rows);
+        if (proceed == MagickFalse)
+          status=MagickFalse;
+      }
+  }
+  image_view=DestroyCacheView(image_view);
+  return(status);
 }
 \f
 /*
@@ -356,7 +942,7 @@ MagickExport Image *AverageImages(const Image *image,ExceptionInfo *exception)
 %  The format of the GetImageChannelExtrema method is:
 %
 %      MagickBooleanType GetImageChannelExtrema(const Image *image,
-%        const ChannelType channel,unsigned long *minima,unsigned long *maxima,
+%        const ChannelType channel,size_t *minima,size_t *maxima,
 %        ExceptionInfo *exception)
 %
 %  A description of each parameter follows:
@@ -374,13 +960,13 @@ MagickExport Image *AverageImages(const Image *image,ExceptionInfo *exception)
 */
 
 MagickExport MagickBooleanType GetImageExtrema(const Image *image,
-  unsigned long *minima,unsigned long *maxima,ExceptionInfo *exception)
+  size_t *minima,size_t *maxima,ExceptionInfo *exception)
 {
   return(GetImageChannelExtrema(image,AllChannels,minima,maxima,exception));
 }
 
 MagickExport MagickBooleanType GetImageChannelExtrema(const Image *image,
-  const ChannelType channel,unsigned long *minima,unsigned long *maxima,
+  const ChannelType channel,size_t *minima,size_t *maxima,
   ExceptionInfo *exception)
 {
   double
@@ -395,8 +981,8 @@ MagickExport MagickBooleanType GetImageChannelExtrema(const Image *image,
   if (image->debug != MagickFalse)
     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
   status=GetImageChannelRange(image,channel,&min,&max,exception);
-  *minima=(unsigned long) (min+0.5);
-  *maxima=(unsigned long) (max+0.5);
+  *minima=(size_t) ceil(min-0.5);
+  *maxima=(size_t) floor(max+0.5);
   return(status);
 }
 \f
@@ -449,79 +1035,82 @@ MagickExport MagickBooleanType GetImageChannelMean(const Image *image,
   const ChannelType channel,double *mean,double *standard_deviation,
   ExceptionInfo *exception)
 {
-  double
-    area;
+  ChannelStatistics
+    *channel_statistics;
 
-  long
-    y;
+  size_t
+    channels;
 
   assert(image != (Image *) NULL);
   assert(image->signature == MagickSignature);
   if (image->debug != MagickFalse)
     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
-  *mean=0.0;
-  *standard_deviation=0.0;
-  area=0.0;
-  for (y=0; y < (long) image->rows; y++)
-  {
-    register const IndexPacket
-      *restrict indexes;
-
-    register const PixelPacket
-      *restrict p;
-
-    register long
-      x;
-
-    p=GetVirtualPixels(image,0,y,image->columns,1,exception);
-    if (p == (const PixelPacket *) NULL)
-      break;
-    indexes=GetVirtualIndexQueue(image);
-    for (x=0; x < (long) image->columns; x++)
+  channel_statistics=GetImageChannelStatistics(image,exception);
+  if (channel_statistics == (ChannelStatistics *) NULL)
+    return(MagickFalse);
+  channels=0;
+  channel_statistics[AllChannels].mean=0.0;
+  channel_statistics[AllChannels].standard_deviation=0.0;
+  if ((channel & RedChannel) != 0)
     {
-      if ((channel & RedChannel) != 0)
-        {
-          *mean+=p->red;
-          *standard_deviation+=(double) p->red*p->red;
-          area++;
-        }
-      if ((channel & GreenChannel) != 0)
-        {
-          *mean+=p->green;
-          *standard_deviation+=(double) p->green*p->green;
-          area++;
-        }
-      if ((channel & BlueChannel) != 0)
-        {
-          *mean+=p->blue;
-          *standard_deviation+=(double) p->blue*p->blue;
-          area++;
-        }
-      if ((channel & OpacityChannel) != 0)
-        {
-          *mean+=p->opacity;
-          *standard_deviation+=(double) p->opacity*p->opacity;
-          area++;
-        }
-      if (((channel & IndexChannel) != 0) &&
-          (image->colorspace == CMYKColorspace))
-        {
-          *mean+=indexes[x];
-          *standard_deviation+=(double) indexes[x]*indexes[x];
-          area++;
-        }
-      p++;
+      channel_statistics[AllChannels].mean+=
+        channel_statistics[RedChannel].mean;
+      channel_statistics[AllChannels].standard_deviation+=
+        channel_statistics[RedChannel].variance-
+        channel_statistics[RedChannel].mean*
+        channel_statistics[RedChannel].mean;
+      channels++;
     }
-  }
-  if (y < (long) image->rows)
-    return(MagickFalse);
-  if (area != 0)
+  if ((channel & GreenChannel) != 0)
+    {
+      channel_statistics[AllChannels].mean+=
+        channel_statistics[GreenChannel].mean;
+      channel_statistics[AllChannels].standard_deviation+=
+        channel_statistics[GreenChannel].variance-
+        channel_statistics[GreenChannel].mean*
+        channel_statistics[GreenChannel].mean;
+      channels++;
+    }
+  if ((channel & BlueChannel) != 0)
+    {
+      channel_statistics[AllChannels].mean+=
+        channel_statistics[BlueChannel].mean;
+      channel_statistics[AllChannels].standard_deviation+=
+        channel_statistics[BlueChannel].variance-
+        channel_statistics[BlueChannel].mean*
+        channel_statistics[BlueChannel].mean;
+      channels++;
+    }
+  if (((channel & OpacityChannel) != 0) &&
+      (image->matte != MagickFalse))
+    {
+      channel_statistics[AllChannels].mean+=
+        channel_statistics[OpacityChannel].mean;
+      channel_statistics[AllChannels].standard_deviation+=
+        channel_statistics[OpacityChannel].variance-
+        channel_statistics[OpacityChannel].mean*
+        channel_statistics[OpacityChannel].mean;
+      channels++;
+    }
+  if (((channel & IndexChannel) != 0) &&
+      (image->colorspace == CMYKColorspace))
     {
-      *mean/=area;
-      *standard_deviation/=area;
+      channel_statistics[AllChannels].mean+=
+        channel_statistics[BlackChannel].mean;
+      channel_statistics[AllChannels].standard_deviation+=
+        channel_statistics[BlackChannel].variance-
+        channel_statistics[BlackChannel].mean*
+        channel_statistics[BlackChannel].mean;
+      channels++;
     }
-  *standard_deviation=sqrt(*standard_deviation-(*mean*(*mean)));
-  return(y == (long) image->rows ? MagickTrue : MagickFalse);
+  channel_statistics[AllChannels].mean/=channels;
+  channel_statistics[AllChannels].standard_deviation=
+    sqrt(channel_statistics[AllChannels].standard_deviation/channels);
+  *mean=channel_statistics[AllChannels].mean;
+  *standard_deviation=channel_statistics[AllChannels].standard_deviation;
+  channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
+    channel_statistics);
+  return(MagickTrue);
 }
 \f
 /*
@@ -581,7 +1170,7 @@ MagickExport MagickBooleanType GetImageChannelKurtosis(const Image *image,
     sum_cubes,
     sum_fourth_power;
 
-  long
+  ssize_t
     y;
 
   assert(image != (Image *) NULL);
@@ -596,7 +1185,7 @@ MagickExport MagickBooleanType GetImageChannelKurtosis(const Image *image,
   sum_squares=0.0;
   sum_cubes=0.0;
   sum_fourth_power=0.0;
-  for (y=0; y < (long) image->rows; y++)
+  for (y=0; y < (ssize_t) image->rows; y++)
   {
     register const IndexPacket
       *restrict indexes;
@@ -604,46 +1193,49 @@ MagickExport MagickBooleanType GetImageChannelKurtosis(const Image *image,
     register const PixelPacket
       *restrict p;
 
-    register long
+    register ssize_t
       x;
 
     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
     if (p == (const PixelPacket *) NULL)
       break;
     indexes=GetVirtualIndexQueue(image);
-    for (x=0; x < (long) image->columns; x++)
+    for (x=0; x < (ssize_t) image->columns; x++)
     {
       if ((channel & RedChannel) != 0)
         {
-          mean+=p->red;
-          sum_squares+=(double) p->red*p->red;
-          sum_cubes+=(double) p->red*p->red*p->red;
-          sum_fourth_power+=(double) p->red*p->red*p->red*p->red;
+          mean+=GetRedPixelComponent(p);
+          sum_squares+=(double) p->red*GetRedPixelComponent(p);
+          sum_cubes+=(double) p->red*p->red*GetRedPixelComponent(p);
+          sum_fourth_power+=(double) p->red*p->red*p->red*
+            GetRedPixelComponent(p);
           area++;
         }
       if ((channel & GreenChannel) != 0)
         {
-          mean+=p->green;
-          sum_squares+=(double) p->green*p->green;
-          sum_cubes+=(double) p->green*p->green*p->green;
-          sum_fourth_power+=(double) p->green*p->green*p->green*p->green;
+          mean+=GetGreenPixelComponent(p);
+          sum_squares+=(double) p->green*GetGreenPixelComponent(p);
+          sum_cubes+=(double) p->green*p->green*GetGreenPixelComponent(p);
+          sum_fourth_power+=(double) p->green*p->green*p->green*
+            GetGreenPixelComponent(p);
           area++;
         }
       if ((channel & BlueChannel) != 0)
         {
-          mean+=p->blue;
-          sum_squares+=(double) p->blue*p->blue;
-          sum_cubes+=(double) p->blue*p->blue*p->blue;
-          sum_fourth_power+=(double) p->blue*p->blue*p->blue*p->blue;
+          mean+=GetBluePixelComponent(p);
+          sum_squares+=(double) p->blue*GetBluePixelComponent(p);
+          sum_cubes+=(double) p->blue*p->blue*GetBluePixelComponent(p);
+          sum_fourth_power+=(double) p->blue*p->blue*p->blue*
+            GetBluePixelComponent(p);
           area++;
         }
       if ((channel & OpacityChannel) != 0)
         {
-          mean+=p->opacity;
-          sum_squares+=(double) p->opacity*p->opacity;
-          sum_cubes+=(double) p->opacity*p->opacity*p->opacity;
+          mean+=GetOpacityPixelComponent(p);
+          sum_squares+=(double) p->opacity*GetOpacityPixelComponent(p);
+          sum_cubes+=(double) p->opacity*p->opacity*GetOpacityPixelComponent(p);
           sum_fourth_power+=(double) p->opacity*p->opacity*p->opacity*
-            p->opacity;
+            GetOpacityPixelComponent(p);
           area++;
         }
       if (((channel & IndexChannel) != 0) &&
@@ -659,7 +1251,7 @@ MagickExport MagickBooleanType GetImageChannelKurtosis(const Image *image,
       p++;
     }
   }
-  if (y < (long) image->rows)
+  if (y < (ssize_t) image->rows)
     return(MagickFalse);
   if (area != 0.0)
     {
@@ -679,9 +1271,9 @@ MagickExport MagickBooleanType GetImageChannelKurtosis(const Image *image,
       *skewness=sum_cubes-3.0*mean*sum_squares+2.0*mean*mean*mean;
       *skewness/=standard_deviation*standard_deviation*standard_deviation;
     }
-  return(y == (long) image->rows ? MagickTrue : MagickFalse);
+  return(y == (ssize_t) image->rows ? MagickTrue : MagickFalse);
 }
-
+\f
 /*
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %                                                                             %
@@ -725,7 +1317,7 @@ MagickExport MagickBooleanType GetImageChannelRange(const Image *image,
   const ChannelType channel,double *minima,double *maxima,
   ExceptionInfo *exception)
 {
-  long
+  ssize_t
     y;
 
   MagickPixelPacket
@@ -738,7 +1330,7 @@ MagickExport MagickBooleanType GetImageChannelRange(const Image *image,
   *maxima=(-1.0E-37);
   *minima=1.0E+37;
   GetMagickPixelPacket(image,&pixel);
-  for (y=0; y < (long) image->rows; y++)
+  for (y=0; y < (ssize_t) image->rows; y++)
   {
     register const IndexPacket
       *restrict indexes;
@@ -746,14 +1338,14 @@ MagickExport MagickBooleanType GetImageChannelRange(const Image *image,
     register const PixelPacket
       *restrict p;
 
-    register long
+    register ssize_t
       x;
 
     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
     if (p == (const PixelPacket *) NULL)
       break;
     indexes=GetVirtualIndexQueue(image);
-    for (x=0; x < (long) image->columns; x++)
+    for (x=0; x < (ssize_t) image->columns; x++)
     {
       SetMagickPixelPacket(image,p,indexes+x,&pixel);
       if ((channel & RedChannel) != 0)
@@ -795,7 +1387,7 @@ MagickExport MagickBooleanType GetImageChannelRange(const Image *image,
       p++;
     }
   }
-  return(y == (long) image->rows ? MagickTrue : MagickFalse);
+  return(y == (ssize_t) image->rows ? MagickTrue : MagickFalse);
 }
 \f
 /*
@@ -831,21 +1423,6 @@ MagickExport MagickBooleanType GetImageChannelRange(const Image *image,
 %    o exception: return any errors or warnings in this structure.
 %
 */
-
-static inline double MagickMax(const double x,const double y)
-{
-  if (x > y)
-    return(x);
-  return(y);
-}
-
-static inline double MagickMin(const double x,const double y)
-{
-  if (x < y)
-    return(x);
-  return(y);
-}
-
 MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
   ExceptionInfo *exception)
 {
@@ -853,11 +1430,9 @@ MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
     *channel_statistics;
 
   double
-    area,
-    sum_squares,
-    sum_cubes;
+    area;
 
-  long
+  ssize_t
     y;
 
   MagickStatusType
@@ -866,13 +1441,13 @@ MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
   QuantumAny
     range;
 
-  register long
+  register ssize_t
     i;
 
   size_t
     length;
 
-  unsigned long
+  size_t
     channels,
     depth;
 
@@ -892,12 +1467,8 @@ MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
     channel_statistics[i].depth=1;
     channel_statistics[i].maxima=(-1.0E-37);
     channel_statistics[i].minima=1.0E+37;
-    channel_statistics[i].mean=0.0;
-    channel_statistics[i].standard_deviation=0.0;
-    channel_statistics[i].kurtosis=0.0;
-    channel_statistics[i].skewness=0.0;
   }
-  for (y=0; y < (long) image->rows; y++)
+  for (y=0; y < (ssize_t) image->rows; y++)
   {
     register const IndexPacket
       *restrict indexes;
@@ -905,14 +1476,14 @@ MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
     register const PixelPacket
       *restrict p;
 
-    register long
+    register ssize_t
       x;
 
     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
     if (p == (const PixelPacket *) NULL)
       break;
     indexes=GetVirtualIndexQueue(image);
-    for (x=0; x < (long) image->columns; )
+    for (x=0; x < (ssize_t) image->columns; )
     {
       if (channel_statistics[RedChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
         {
@@ -981,49 +1552,57 @@ MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
             }
         }
       if ((double) p->red < channel_statistics[RedChannel].minima)
-        channel_statistics[RedChannel].minima=(double) p->red;
+        channel_statistics[RedChannel].minima=(double) GetRedPixelComponent(p);
       if ((double) p->red > channel_statistics[RedChannel].maxima)
-        channel_statistics[RedChannel].maxima=(double) p->red;
-      channel_statistics[RedChannel].mean+=p->red;
-      channel_statistics[RedChannel].standard_deviation+=(double) p->red*p->red;
-      channel_statistics[RedChannel].kurtosis+=(double) p->red*p->red*
-        p->red*p->red;
-      channel_statistics[RedChannel].skewness+=(double) p->red*p->red*p->red;
+        channel_statistics[RedChannel].maxima=(double) GetRedPixelComponent(p);
+      channel_statistics[RedChannel].sum+=GetRedPixelComponent(p);
+      channel_statistics[RedChannel].sum_squared+=(double) p->red*
+        GetRedPixelComponent(p);
+      channel_statistics[RedChannel].sum_cubed+=(double) p->red*p->red*
+        GetRedPixelComponent(p);
+      channel_statistics[RedChannel].sum_fourth_power+=(double) p->red*p->red*
+        p->red*GetRedPixelComponent(p);
       if ((double) p->green < channel_statistics[GreenChannel].minima)
-        channel_statistics[GreenChannel].minima=(double) p->green;
+        channel_statistics[GreenChannel].minima=(double)
+          GetGreenPixelComponent(p);
       if ((double) p->green > channel_statistics[GreenChannel].maxima)
-        channel_statistics[GreenChannel].maxima=(double) p->green;
-      channel_statistics[GreenChannel].mean+=p->green;
-      channel_statistics[GreenChannel].standard_deviation+=(double) p->green*
-        p->green;
-      channel_statistics[GreenChannel].kurtosis+=(double) p->green*p->green*
-        p->green*p->green;
-      channel_statistics[GreenChannel].skewness+=(double) p->green*p->green*
-        p->green;
+        channel_statistics[GreenChannel].maxima=(double)
+          GetGreenPixelComponent(p);
+      channel_statistics[GreenChannel].sum+=GetGreenPixelComponent(p);
+      channel_statistics[GreenChannel].sum_squared+=(double) p->green*
+        GetGreenPixelComponent(p);
+      channel_statistics[GreenChannel].sum_cubed+=(double) p->green*p->green*
+        GetGreenPixelComponent(p);
+      channel_statistics[GreenChannel].sum_fourth_power+=(double) p->green*
+        p->green*p->green*GetGreenPixelComponent(p);
       if ((double) p->blue < channel_statistics[BlueChannel].minima)
-        channel_statistics[BlueChannel].minima=(double) p->blue;
+        channel_statistics[BlueChannel].minima=(double)
+          GetBluePixelComponent(p);
       if ((double) p->blue > channel_statistics[BlueChannel].maxima)
-        channel_statistics[BlueChannel].maxima=(double) p->blue;
-      channel_statistics[BlueChannel].mean+=p->blue;
-      channel_statistics[BlueChannel].standard_deviation+=(double) p->blue*
-        p->blue;
-      channel_statistics[BlueChannel].kurtosis+=(double) p->blue*p->blue*
-        p->blue*p->blue;
-      channel_statistics[BlueChannel].skewness+=(double) p->blue*p->blue*
-        p->blue;
+        channel_statistics[BlueChannel].maxima=(double)
+          GetBluePixelComponent(p);
+      channel_statistics[BlueChannel].sum+=GetBluePixelComponent(p);
+      channel_statistics[BlueChannel].sum_squared+=(double) p->blue*
+        GetBluePixelComponent(p);
+      channel_statistics[BlueChannel].sum_cubed+=(double) p->blue*p->blue*
+        GetBluePixelComponent(p);
+      channel_statistics[BlueChannel].sum_fourth_power+=(double) p->blue*
+        p->blue*p->blue*GetBluePixelComponent(p);
       if (image->matte != MagickFalse)
         {
           if ((double) p->opacity < channel_statistics[OpacityChannel].minima)
-            channel_statistics[OpacityChannel].minima=(double) p->opacity;
+            channel_statistics[OpacityChannel].minima=(double)
+              GetOpacityPixelComponent(p);
           if ((double) p->opacity > channel_statistics[OpacityChannel].maxima)
-            channel_statistics[OpacityChannel].maxima=(double) p->opacity;
-          channel_statistics[OpacityChannel].mean+=p->opacity;
-          channel_statistics[OpacityChannel].standard_deviation+=(double)
-            p->opacity*p->opacity;
-          channel_statistics[OpacityChannel].kurtosis+=(double) p->opacity*
-            p->opacity*p->opacity*p->opacity;
-          channel_statistics[OpacityChannel].skewness+=(double) p->opacity*
-            p->opacity*p->opacity;
+            channel_statistics[OpacityChannel].maxima=(double)
+              GetOpacityPixelComponent(p);
+          channel_statistics[OpacityChannel].sum+=GetOpacityPixelComponent(p);
+          channel_statistics[OpacityChannel].sum_squared+=(double)
+            p->opacity*GetOpacityPixelComponent(p);
+          channel_statistics[OpacityChannel].sum_cubed+=(double) p->opacity*
+            p->opacity*GetOpacityPixelComponent(p);
+          channel_statistics[OpacityChannel].sum_fourth_power+=(double)
+            p->opacity*p->opacity*p->opacity*GetOpacityPixelComponent(p);
         }
       if (image->colorspace == CMYKColorspace)
         {
@@ -1031,13 +1610,13 @@ MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
             channel_statistics[BlackChannel].minima=(double) indexes[x];
           if ((double) indexes[x] > channel_statistics[BlackChannel].maxima)
             channel_statistics[BlackChannel].maxima=(double) indexes[x];
-          channel_statistics[BlackChannel].mean+=indexes[x];
-          channel_statistics[BlackChannel].standard_deviation+=(double)
+          channel_statistics[BlackChannel].sum+=indexes[x];
+          channel_statistics[BlackChannel].sum_squared+=(double)
             indexes[x]*indexes[x];
-          channel_statistics[BlackChannel].kurtosis+=(double) indexes[x]*
-            indexes[x]*indexes[x]*indexes[x];
-          channel_statistics[BlackChannel].skewness+=(double) indexes[x]*
+          channel_statistics[BlackChannel].sum_cubed+=(double) indexes[x]*
             indexes[x]*indexes[x];
+          channel_statistics[BlackChannel].sum_fourth_power+=(double)
+            indexes[x]*indexes[x]*indexes[x]*indexes[x];
         }
       x++;
       p++;
@@ -1046,66 +1625,72 @@ MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
   area=(double) image->columns*image->rows;
   for (i=0; i < AllChannels; i++)
   {
-    channel_statistics[i].mean/=area;
-    channel_statistics[i].standard_deviation/=area;
-    channel_statistics[i].kurtosis/=area;
-    channel_statistics[i].skewness/=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(
+      channel_statistics[i].variance-(channel_statistics[i].mean*
+      channel_statistics[i].mean));
   }
   for (i=0; i < AllChannels; i++)
   {
-    channel_statistics[AllChannels].depth=(unsigned long) MagickMax((double)
+    channel_statistics[AllChannels].depth=(size_t) MagickMax((double)
       channel_statistics[AllChannels].depth,(double)
       channel_statistics[i].depth);
     channel_statistics[AllChannels].minima=MagickMin(
       channel_statistics[AllChannels].minima,channel_statistics[i].minima);
     channel_statistics[AllChannels].maxima=MagickMax(
       channel_statistics[AllChannels].maxima,channel_statistics[i].maxima);
+    channel_statistics[AllChannels].sum+=channel_statistics[i].sum;
+    channel_statistics[AllChannels].sum_squared+=
+      channel_statistics[i].sum_squared;
+    channel_statistics[AllChannels].sum_cubed+=channel_statistics[i].sum_cubed;
+    channel_statistics[AllChannels].sum_fourth_power+=
+      channel_statistics[i].sum_fourth_power;
     channel_statistics[AllChannels].mean+=channel_statistics[i].mean;
+    channel_statistics[AllChannels].variance+=channel_statistics[i].variance-
+      channel_statistics[i].mean*channel_statistics[i].mean;
     channel_statistics[AllChannels].standard_deviation+=
-      channel_statistics[i].standard_deviation;
-    channel_statistics[AllChannels].kurtosis+=channel_statistics[i].kurtosis;
-    channel_statistics[AllChannels].skewness+=channel_statistics[i].skewness;
+      channel_statistics[i].variance-channel_statistics[i].mean*
+      channel_statistics[i].mean;
   }
-  channels=4;
+  channels=3;
+  if (image->matte != MagickFalse)
+    channels++;
   if (image->colorspace == CMYKColorspace)
     channels++;
+  channel_statistics[AllChannels].sum/=channels;
+  channel_statistics[AllChannels].sum_squared/=channels;
+  channel_statistics[AllChannels].sum_cubed/=channels;
+  channel_statistics[AllChannels].sum_fourth_power/=channels;
   channel_statistics[AllChannels].mean/=channels;
-  channel_statistics[AllChannels].standard_deviation/=channels;
+  channel_statistics[AllChannels].variance/=channels;
+  channel_statistics[AllChannels].standard_deviation=
+    sqrt(channel_statistics[AllChannels].standard_deviation/channels);
   channel_statistics[AllChannels].kurtosis/=channels;
   channel_statistics[AllChannels].skewness/=channels;
   for (i=0; i <= AllChannels; i++)
   {
-    sum_squares=0.0;
-    sum_squares=channel_statistics[i].standard_deviation;
-    sum_cubes=0.0;
-    sum_cubes=channel_statistics[i].skewness;
-    channel_statistics[i].standard_deviation=sqrt(
-      channel_statistics[i].standard_deviation-
-       (channel_statistics[i].mean*channel_statistics[i].mean));
     if (channel_statistics[i].standard_deviation == 0.0)
-      {
-        channel_statistics[i].kurtosis=0.0;
-        channel_statistics[i].skewness=0.0;
-      }
-    else
-      {
-        channel_statistics[i].skewness=(channel_statistics[i].skewness-
-          3.0*channel_statistics[i].mean*sum_squares+
-          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].kurtosis=(channel_statistics[i].kurtosis-
-          4.0*channel_statistics[i].mean*sum_cubes+
-          6.0*channel_statistics[i].mean*channel_statistics[i].mean*sum_squares-
-          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;
-      }
+      continue;
+    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].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;
   }
   return(channel_statistics);
 }