]> granicus.if.org Git - imagemagick/blobdiff - magick/statistic.c
(no commit message)
[imagemagick] / magick / statistic.c
index 200a3f6f0c339d5b6f8f144c4756b8c7a8c029c1..5f4e743d2630c03c2c15642dea90e25497f3e8eb 100644 (file)
 
 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);
@@ -147,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();
@@ -163,13 +163,13 @@ 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);
@@ -200,6 +200,11 @@ static MagickRealType ApplyEvaluateOperator(RandomInfo *random_info,
   {
     case UndefinedEvaluateOperator:
       break;
+    case AbsEvaluateOperator:
+    {
+      result=(MagickRealType) fabs((double) (pixel+value));
+      break;
+    }
     case AddEvaluateOperator:
     {
       result=(MagickRealType) (pixel+value);
@@ -214,12 +219,12 @@ static MagickRealType ApplyEvaluateOperator(RandomInfo *random_info,
         and could return a negative result (which is clipped).
       */
       result=pixel+value;
-      result-=(QuantumRange+1)*floor(result/(QuantumRange+1));
+      result-=(QuantumRange+1.0)*floor((double) result/(QuantumRange+1.0));
       break;
     }
     case AndEvaluateOperator:
     {
-      result=(MagickRealType) ((unsigned long) pixel & (unsigned long)
+      result=(MagickRealType) ((size_t) pixel & (size_t)
         (value+0.5));
       break;
     }
@@ -254,7 +259,7 @@ static MagickRealType ApplyEvaluateOperator(RandomInfo *random_info,
     }
     case LeftShiftEvaluateOperator:
     {
-      result=(MagickRealType) ((unsigned long) pixel << (unsigned long)
+      result=(MagickRealType) ((size_t) pixel << (size_t)
         (value+0.5));
       break;
     }
@@ -292,7 +297,7 @@ static MagickRealType ApplyEvaluateOperator(RandomInfo *random_info,
     }
     case OrEvaluateOperator:
     {
-      result=(MagickRealType) ((unsigned long) pixel | (unsigned long)
+      result=(MagickRealType) ((size_t) pixel | (size_t)
         (value+0.5));
       break;
     }
@@ -310,7 +315,7 @@ static MagickRealType ApplyEvaluateOperator(RandomInfo *random_info,
     }
     case RightShiftEvaluateOperator:
     {
-      result=(MagickRealType) ((unsigned long) pixel >> (unsigned long)
+      result=(MagickRealType) ((size_t) pixel >> (size_t)
         (value+0.5));
       break;
     }
@@ -355,7 +360,7 @@ static MagickRealType ApplyEvaluateOperator(RandomInfo *random_info,
     }
     case XorEvaluateOperator:
     {
-      result=(MagickRealType) ((unsigned long) pixel ^ (unsigned long)
+      result=(MagickRealType) ((size_t) pixel ^ (size_t)
         (value+0.5));
       break;
     }
@@ -387,13 +392,12 @@ MagickExport Image *EvaluateImages(const Image *images,
   Image
     *evaluate_image;
 
-  long
-    progress,
-    y;
-
   MagickBooleanType
     status;
 
+  MagickOffsetType
+    progress;
+
   MagickPixelPacket
     **restrict evaluate_pixels,
     zero;
@@ -401,9 +405,12 @@ MagickExport Image *EvaluateImages(const Image *images,
   RandomInfo
     **restrict random_info;
 
-  unsigned long
+  size_t
     number_images;
 
+  ssize_t
+    y;
+
   /*
     Ensure the image are the same size.
   */
@@ -453,7 +460,7 @@ MagickExport Image *EvaluateImages(const Image *images,
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
   #pragma omp parallel for schedule(dynamic) shared(progress,status)
 #endif
-  for (y=0; y < (long) evaluate_image->rows; y++)
+  for (y=0; y < (ssize_t) evaluate_image->rows; y++)
   {
     CacheView
       *image_view;
@@ -461,15 +468,17 @@ MagickExport Image *EvaluateImages(const Image *images,
     const Image
       *next;
 
+    int
+      id;
+
     MagickPixelPacket
       pixel;
 
     register IndexPacket
       *restrict evaluate_indexes;
 
-    register long
+    register ssize_t
       i,
-      id,
       x;
 
     register MagickPixelPacket
@@ -491,10 +500,10 @@ MagickExport Image *EvaluateImages(const Image *images,
     pixel=zero;
     id=GetOpenMPThreadId();
     evaluate_pixel=evaluate_pixels[id];
-    for (x=0; x < (long) evaluate_image->columns; x++)
+    for (x=0; x < (ssize_t) evaluate_image->columns; x++)
       evaluate_pixel[x]=zero;
     next=images;
-    for (i=0; i < (long) number_images; i++)
+    for (i=0; i < (ssize_t) number_images; i++)
     {
       register const IndexPacket
         *indexes;
@@ -510,7 +519,7 @@ MagickExport Image *EvaluateImages(const Image *images,
           break;
         }
       indexes=GetCacheViewVirtualIndexQueue(image_view);
-      for (x=0; x < (long) next->columns; x++)
+      for (x=0; x < (ssize_t) next->columns; x++)
       {
         evaluate_pixel[x].red=ApplyEvaluateOperator(random_info[id],p->red,
           i == 0 ? AddEvaluateOperator : op,evaluate_pixel[x].red);
@@ -531,7 +540,7 @@ MagickExport Image *EvaluateImages(const Image *images,
       next=GetNextImageInList(next);
     }
     if (op == MeanEvaluateOperator)
-      for (x=0; x < (long) evaluate_image->columns; x++)
+      for (x=0; x < (ssize_t) evaluate_image->columns; x++)
       {
         evaluate_pixel[x].red/=number_images;
         evaluate_pixel[x].green/=number_images;
@@ -539,7 +548,7 @@ MagickExport Image *EvaluateImages(const Image *images,
         evaluate_pixel[x].opacity/=number_images;
         evaluate_pixel[x].index/=number_images;
       }
-    for (x=0; x < (long) evaluate_image->columns; x++)
+    for (x=0; x < (ssize_t) evaluate_image->columns; x++)
     {
       q->red=ClampToQuantum(evaluate_pixel[x].red);
       q->green=ClampToQuantum(evaluate_pixel[x].green);
@@ -583,16 +592,18 @@ MagickExport MagickBooleanType EvaluateImageChannel(Image *image,
   CacheView
     *image_view;
 
-  long
-    progress,
-    y;
-
   MagickBooleanType
     status;
 
+  MagickOffsetType
+    progress;
+
   RandomInfo
     **restrict random_info;
 
+  ssize_t
+    y;
+
   assert(image != (Image *) NULL);
   assert(image->signature == MagickSignature);
   if (image->debug != MagickFalse)
@@ -611,13 +622,15 @@ MagickExport MagickBooleanType EvaluateImageChannel(Image *image,
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
 #endif
-  for (y=0; y < (long) image->rows; y++)
+  for (y=0; y < (ssize_t) image->rows; y++)
   {
+    int
+      id;
+
     register IndexPacket
       *restrict indexes;
 
-    register long
-      id,
+    register ssize_t
       x;
 
     register PixelPacket
@@ -633,7 +646,7 @@ MagickExport MagickBooleanType EvaluateImageChannel(Image *image,
       }
     indexes=GetCacheViewAuthenticIndexQueue(image_view);
     id=GetOpenMPThreadId();
-    for (x=0; x < (long) image->columns; x++)
+    for (x=0; x < (ssize_t) image->columns; x++)
     {
       if ((channel & RedChannel) != 0)
         q->red=ClampToQuantum(ApplyEvaluateOperator(random_info[id],q->red,op,
@@ -697,11 +710,11 @@ MagickExport MagickBooleanType EvaluateImageChannel(Image *image,
 %  The format of the FunctionImageChannel method is:
 %
 %      MagickBooleanType FunctionImage(Image *image,
-%        const MagickFunction function,const long number_parameters,
+%        const MagickFunction function,const ssize_t number_parameters,
 %        const double *parameters,ExceptionInfo *exception)
 %      MagickBooleanType FunctionImageChannel(Image *image,
 %        const ChannelType channel,const MagickFunction function,
-%        const long number_parameters,const double *argument,
+%        const ssize_t number_parameters,const double *argument,
 %        ExceptionInfo *exception)
 %
 %  A description of each parameter follows:
@@ -719,13 +732,13 @@ MagickExport MagickBooleanType EvaluateImageChannel(Image *image,
 */
 
 static Quantum ApplyFunction(Quantum pixel,const MagickFunction function,
-  const unsigned long number_parameters,const double *parameters,
+  const size_t number_parameters,const double *parameters,
   ExceptionInfo *exception)
 {
   MagickRealType
     result;
 
-  register long
+  register ssize_t
     i;
 
   (void) exception;
@@ -740,7 +753,7 @@ static Quantum ApplyFunction(Quantum pixel,const MagickFunction function,
        *   For example:      c0*x^3 + c1*x^2 + c2*x  + c3
        */
       result=0.0;
-      for (i=0; i < (long) number_parameters; i++)
+      for (i=0; i < (ssize_t) number_parameters; i++)
         result = result*QuantumScale*pixel + parameters[i];
       result *= QuantumRange;
       break;
@@ -801,7 +814,7 @@ static Quantum ApplyFunction(Quantum pixel,const MagickFunction function,
 }
 
 MagickExport MagickBooleanType FunctionImage(Image *image,
-  const MagickFunction function,const unsigned long number_parameters,
+  const MagickFunction function,const size_t number_parameters,
   const double *parameters,ExceptionInfo *exception)
 {
   MagickBooleanType
@@ -814,7 +827,7 @@ MagickExport MagickBooleanType FunctionImage(Image *image,
 
 MagickExport MagickBooleanType FunctionImageChannel(Image *image,
   const ChannelType channel,const MagickFunction function,
-  const unsigned long number_parameters,const double *parameters,
+  const size_t number_parameters,const double *parameters,
   ExceptionInfo *exception)
 {
 #define FunctionImageTag  "Function/Image "
@@ -822,13 +835,15 @@ MagickExport MagickBooleanType FunctionImageChannel(Image *image,
   CacheView
     *image_view;
 
-  long
-    progress,
-    y;
-
   MagickBooleanType
     status;
 
+  MagickOffsetType
+    progress;
+
+  ssize_t
+    y;
+
   assert(image != (Image *) NULL);
   assert(image->signature == MagickSignature);
   if (image->debug != MagickFalse)
@@ -846,12 +861,12 @@ MagickExport MagickBooleanType FunctionImageChannel(Image *image,
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
 #endif
-  for (y=0; y < (long) image->rows; y++)
+  for (y=0; y < (ssize_t) image->rows; y++)
   {
     register IndexPacket
       *restrict indexes;
 
-    register long
+    register ssize_t
       x;
 
     register PixelPacket
@@ -866,7 +881,7 @@ MagickExport MagickBooleanType FunctionImageChannel(Image *image,
         continue;
       }
     indexes=GetCacheViewAuthenticIndexQueue(image_view);
-    for (x=0; x < (long) image->columns; x++)
+    for (x=0; x < (ssize_t) image->columns; x++)
     {
       if ((channel & RedChannel) != 0)
         q->red=ApplyFunction(q->red,function,number_parameters,parameters,
@@ -927,7 +942,7 @@ MagickExport MagickBooleanType FunctionImageChannel(Image *image,
 %  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:
@@ -945,13 +960,13 @@ MagickExport MagickBooleanType FunctionImageChannel(Image *image,
 */
 
 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
@@ -966,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) ceil(min-0.5);
-  *maxima=(unsigned long) floor(max+0.5);
+  *minima=(size_t) ceil(min-0.5);
+  *maxima=(size_t) floor(max+0.5);
   return(status);
 }
 \f
@@ -1020,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+=GetRedPixelComponent(p);
-          *standard_deviation+=(double) p->red*GetRedPixelComponent(p);
-          area++;
-        }
-      if ((channel & GreenChannel) != 0)
-        {
-          *mean+=GetGreenPixelComponent(p);
-          *standard_deviation+=(double) p->green*GetGreenPixelComponent(p);
-          area++;
-        }
-      if ((channel & BlueChannel) != 0)
-        {
-          *mean+=GetBluePixelComponent(p);
-          *standard_deviation+=(double) p->blue*GetBluePixelComponent(p);
-          area++;
-        }
-      if ((channel & OpacityChannel) != 0)
-        {
-          *mean+=GetOpacityPixelComponent(p);
-          *standard_deviation+=(double) p->opacity*GetOpacityPixelComponent(p);
-          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)
     {
-      *mean/=area;
-      *standard_deviation/=area;
+      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++;
     }
-  *standard_deviation=sqrt(*standard_deviation-(*mean*(*mean)));
-  return(y == (long) image->rows ? MagickTrue : MagickFalse);
+  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))
+    {
+      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++;
+    }
+  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
 /*
@@ -1152,7 +1170,7 @@ MagickExport MagickBooleanType GetImageChannelKurtosis(const Image *image,
     sum_cubes,
     sum_fourth_power;
 
-  long
+  ssize_t
     y;
 
   assert(image != (Image *) NULL);
@@ -1167,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;
@@ -1175,14 +1193,14 @@ 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)
         {
@@ -1233,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)
     {
@@ -1253,7 +1271,7 @@ 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
 /*
@@ -1299,7 +1317,7 @@ MagickExport MagickBooleanType GetImageChannelRange(const Image *image,
   const ChannelType channel,double *minima,double *maxima,
   ExceptionInfo *exception)
 {
-  long
+  ssize_t
     y;
 
   MagickPixelPacket
@@ -1312,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;
@@ -1320,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)
@@ -1369,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
 /*
@@ -1412,11 +1430,9 @@ MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
     *channel_statistics;
 
   double
-    area,
-    sum_squares,
-    sum_cubes;
+    area;
 
-  long
+  ssize_t
     y;
 
   MagickStatusType
@@ -1425,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;
 
@@ -1451,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;
@@ -1464,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)
         {
@@ -1543,39 +1555,39 @@ MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
         channel_statistics[RedChannel].minima=(double) GetRedPixelComponent(p);
       if ((double) p->red > channel_statistics[RedChannel].maxima)
         channel_statistics[RedChannel].maxima=(double) GetRedPixelComponent(p);
-      channel_statistics[RedChannel].mean+=GetRedPixelComponent(p);
-      channel_statistics[RedChannel].standard_deviation+=(double) p->red*
+      channel_statistics[RedChannel].sum+=GetRedPixelComponent(p);
+      channel_statistics[RedChannel].sum_squared+=(double) p->red*
         GetRedPixelComponent(p);
-      channel_statistics[RedChannel].kurtosis+=(double) p->red*p->red*
-        p->red*GetRedPixelComponent(p);
-      channel_statistics[RedChannel].skewness+=(double) p->red*p->red*
+      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)
           GetGreenPixelComponent(p);
       if ((double) p->green > channel_statistics[GreenChannel].maxima)
         channel_statistics[GreenChannel].maxima=(double)
           GetGreenPixelComponent(p);
-      channel_statistics[GreenChannel].mean+=GetGreenPixelComponent(p);
-      channel_statistics[GreenChannel].standard_deviation+=(double) p->green*
+      channel_statistics[GreenChannel].sum+=GetGreenPixelComponent(p);
+      channel_statistics[GreenChannel].sum_squared+=(double) p->green*
         GetGreenPixelComponent(p);
-      channel_statistics[GreenChannel].kurtosis+=(double) p->green*p->green*
-        p->green*GetGreenPixelComponent(p);
-      channel_statistics[GreenChannel].skewness+=(double) p->green*p->green*
+      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)
           GetBluePixelComponent(p);
       if ((double) p->blue > channel_statistics[BlueChannel].maxima)
         channel_statistics[BlueChannel].maxima=(double)
           GetBluePixelComponent(p);
-      channel_statistics[BlueChannel].mean+=GetBluePixelComponent(p);
-      channel_statistics[BlueChannel].standard_deviation+=(double) p->blue*
+      channel_statistics[BlueChannel].sum+=GetBluePixelComponent(p);
+      channel_statistics[BlueChannel].sum_squared+=(double) p->blue*
         GetBluePixelComponent(p);
-      channel_statistics[BlueChannel].kurtosis+=(double) p->blue*p->blue*
-        p->blue*GetBluePixelComponent(p);
-      channel_statistics[BlueChannel].skewness+=(double) p->blue*p->blue*
+      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)
@@ -1584,13 +1596,13 @@ MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
           if ((double) p->opacity > channel_statistics[OpacityChannel].maxima)
             channel_statistics[OpacityChannel].maxima=(double)
               GetOpacityPixelComponent(p);
-          channel_statistics[OpacityChannel].mean+=GetOpacityPixelComponent(p);
-          channel_statistics[OpacityChannel].standard_deviation+=(double)
+          channel_statistics[OpacityChannel].sum+=GetOpacityPixelComponent(p);
+          channel_statistics[OpacityChannel].sum_squared+=(double)
             p->opacity*GetOpacityPixelComponent(p);
-          channel_statistics[OpacityChannel].kurtosis+=(double) p->opacity*
-            p->opacity*p->opacity*GetOpacityPixelComponent(p);
-          channel_statistics[OpacityChannel].skewness+=(double) p->opacity*
+          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)
         {
@@ -1598,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++;
@@ -1613,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);
 }