]> granicus.if.org Git - imagemagick/blobdiff - MagickCore/statistic.c
Add RobidouxSharp filter depreciate Bessel Filter and Static Gravity
[imagemagick] / MagickCore / statistic.c
index 9e987e8e77053f3b5cb8fae00b15870ec313f477..99f194d9a550cd6cbce778bb14534f7096dc6bde 100644 (file)
@@ -17,7 +17,7 @@
 %                                 July 1992                                   %
 %                                                                             %
 %                                                                             %
-%  Copyright 1999-2011 ImageMagick Studio LLC, a non-profit organization      %
+%  Copyright 1999-2012 ImageMagick Studio LLC, a non-profit organization      %
 %  dedicated to making software imaging solutions freely available.           %
 %                                                                             %
 %  You may not use this file except in compliance with the License.  You may  %
@@ -63,6 +63,7 @@
 #include "MagickCore/exception.h"
 #include "MagickCore/exception-private.h"
 #include "MagickCore/gem.h"
+#include "MagickCore/gem-private.h"
 #include "MagickCore/geometry.h"
 #include "MagickCore/list.h"
 #include "MagickCore/image-private.h"
 %  an image, to increase or decrease contrast in an image, or to produce the
 %  "negative" of an image.
 %
-%  The format of the EvaluateImageChannel method is:
+%  The format of the EvaluateImage method is:
 %
 %      MagickBooleanType EvaluateImage(Image *image,
 %        const MagickEvaluateOperator op,const double value,
 %      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.
 %
-%    o channel: the channel.
-%
 %    o op: A channel op.
 %
 %    o value: A value value.
 %
 */
 
-static PixelInfo **DestroyPixelThreadSet(PixelInfo **pixels)
+typedef struct _PixelChannels
+{
+  MagickRealType
+    channel[CompositePixelChannel];
+} PixelChannels;
+
+static PixelChannels **DestroyPixelThreadSet(PixelChannels **pixels)
 {
   register ssize_t
     i;
 
-  assert(pixels != (PixelInfo **) NULL);
+  assert(pixels != (PixelChannels **) NULL);
   for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
-    if (pixels[i] != (PixelInfo *) NULL)
-      pixels[i]=(PixelInfo *) RelinquishMagickMemory(pixels[i]);
-  pixels=(PixelInfo **) RelinquishMagickMemory(pixels);
+    if (pixels[i] != (PixelChannels *) NULL)
+      pixels[i]=(PixelChannels *) RelinquishMagickMemory(pixels[i]);
+  pixels=(PixelChannels **) RelinquishMagickMemory(pixels);
   return(pixels);
 }
 
-static PixelInfo **AcquirePixelThreadSet(const Image *image,
+static PixelChannels **AcquirePixelThreadSet(const Image *image,
   const size_t number_images)
 {
   register ssize_t
-    i,
-    j;
+    i;
 
-  PixelInfo
+  PixelChannels
     **pixels;
 
   size_t
@@ -160,27 +161,35 @@ static PixelInfo **AcquirePixelThreadSet(const Image *image,
     number_threads;
 
   number_threads=GetOpenMPMaximumThreads();
-  pixels=(PixelInfo **) AcquireQuantumMemory(number_threads,
+  pixels=(PixelChannels **) AcquireQuantumMemory(number_threads,
     sizeof(*pixels));
-  if (pixels == (PixelInfo **) NULL)
-    return((PixelInfo **) NULL);
+  if (pixels == (PixelChannels **) NULL)
+    return((PixelChannels **) NULL);
   (void) ResetMagickMemory(pixels,0,number_threads*sizeof(*pixels));
   for (i=0; i < (ssize_t) number_threads; i++)
   {
+    register ssize_t
+      j;
+
     length=image->columns;
     if (length < number_images)
       length=number_images;
-    pixels[i]=(PixelInfo *) AcquireQuantumMemory(length,
-      sizeof(**pixels));
-    if (pixels[i] == (PixelInfo *) NULL)
+    pixels[i]=(PixelChannels *) AcquireQuantumMemory(length,sizeof(**pixels));
+    if (pixels[i] == (PixelChannels *) NULL)
       return(DestroyPixelThreadSet(pixels));
     for (j=0; j < (ssize_t) length; j++)
-      GetPixelInfo(image,&pixels[i][j]);
+    {
+      register ssize_t
+        k;
+
+      for (k=0; k < MaxPixelChannels; k++)
+        pixels[i][j].channel[k]=0.0;
+    }
   }
   return(pixels);
 }
 
-static inline double MagickMax(const double x,const double y)
+static inline double EvaluateMax(const double x,const double y)
 {
   if (x > y)
     return(x);
@@ -193,18 +202,22 @@ extern "C" {
 
 static int IntensityCompare(const void *x,const void *y)
 {
-  const PixelInfo
+  const PixelChannels
     *color_1,
     *color_2;
 
-  int
-    intensity;
+  MagickRealType
+    distance;
+
+  register ssize_t
+    i;
 
-  color_1=(const PixelInfo *) x;
-  color_2=(const PixelInfo *) y;
-  intensity=(int) GetPixelInfoIntensity(color_2)-(int)
-    GetPixelInfoIntensity(color_1);
-  return(intensity);
+  color_1=(const PixelChannels *) x;
+  color_2=(const PixelChannels *) y;
+  distance=0.0;
+  for (i=0; i < MaxPixelChannels; i++)
+    distance+=color_1->channel[i]-(MagickRealType) color_2->channel[i];
+  return(distance < 0 ? -1 : distance > 0 ? 1 : 0);
 }
 
 #if defined(__cplusplus) || defined(c_plusplus)
@@ -242,10 +255,10 @@ static MagickRealType ApplyEvaluateOperator(RandomInfo *random_info,
     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).
+        This returns a 'floored modulus' of the addition which is a positive
+        result.  It differs from % or fmod() that 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));
@@ -304,7 +317,7 @@ static MagickRealType ApplyEvaluateOperator(RandomInfo *random_info,
     }
     case MaxEvaluateOperator:
     {
-      result=(MagickRealType) MagickMax((double) pixel,value);
+      result=(MagickRealType) EvaluateMax((double) pixel,value);
       break;
     }
     case MeanEvaluateOperator:
@@ -371,6 +384,11 @@ static MagickRealType ApplyEvaluateOperator(RandomInfo *random_info,
       result=(MagickRealType) (pixel-value);
       break;
     }
+    case SumEvaluateOperator:
+    {
+      result=(MagickRealType) (pixel+value);
+      break;
+    }
     case ThresholdEvaluateOperator:
     {
       result=(MagickRealType) (((MagickRealType) pixel <= value) ? 0 :
@@ -403,16 +421,6 @@ static MagickRealType ApplyEvaluateOperator(RandomInfo *random_info,
   return(result);
 }
 
-MagickExport MagickBooleanType EvaluateImage(Image *image,
-  const MagickEvaluateOperator op,const double value,ExceptionInfo *exception)
-{
-  MagickBooleanType
-    status;
-
-  status=EvaluateImageChannel(image,CompositeChannels,op,value,exception);
-  return(status);
-}
-
 MagickExport Image *EvaluateImages(const Image *images,
   const MagickEvaluateOperator op,ExceptionInfo *exception)
 {
@@ -428,14 +436,14 @@ MagickExport Image *EvaluateImages(const Image *images,
     *evaluate_image;
 
   MagickBooleanType
+    concurrent,
     status;
 
   MagickOffsetType
     progress;
 
-  PixelInfo
-    **restrict evaluate_pixels,
-    zero;
+  PixelChannels
+    **restrict evaluate_pixels;
 
   RandomInfo
     **restrict random_info;
@@ -459,7 +467,7 @@ MagickExport Image *EvaluateImages(const Image *images,
     if ((next->columns != images->columns) || (next->rows != images->rows))
       {
         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
-          "ImageWidthsOrHeightsDiffer","`%s'",images->filename);
+          "ImageWidthsOrHeightsDiffer","'%s'",images->filename);
         return((Image *) NULL);
       }
   /*
@@ -469,19 +477,18 @@ MagickExport Image *EvaluateImages(const Image *images,
     exception);
   if (evaluate_image == (Image *) NULL)
     return((Image *) NULL);
-  if (SetImageStorageClass(evaluate_image,DirectClass) == MagickFalse)
+  if (SetImageStorageClass(evaluate_image,DirectClass,exception) == MagickFalse)
     {
-      InheritException(exception,&evaluate_image->exception);
       evaluate_image=DestroyImage(evaluate_image);
       return((Image *) NULL);
     }
   number_images=GetImageListLength(images);
   evaluate_pixels=AcquirePixelThreadSet(images,number_images);
-  if (evaluate_pixels == (PixelInfo **) NULL)
+  if (evaluate_pixels == (PixelChannels **) NULL)
     {
       evaluate_image=DestroyImage(evaluate_image);
       (void) ThrowMagickException(exception,GetMagickModule(),
-        ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
+        ResourceLimitError,"MemoryAllocationFailed","'%s'",images->filename);
       return((Image *) NULL);
     }
   /*
@@ -489,225 +496,281 @@ MagickExport Image *EvaluateImages(const Image *images,
   */
   status=MagickTrue;
   progress=0;
-  GetPixelInfo(images,&zero);
   random_info=AcquireRandomInfoThreadSet();
-  evaluate_view=AcquireCacheView(evaluate_image);
+  concurrent=GetRandomSecretKey(random_info[0]) == ~0UL ? MagickTrue :
+    MagickFalse;
+  evaluate_view=AcquireAuthenticCacheView(evaluate_image,exception);
   if (op == MedianEvaluateOperator)
+    {
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(dynamic) shared(progress,status)
+      #pragma omp parallel for schedule(static) shared(progress,status) omp_concurrent(concurrent)
 #endif
-    for (y=0; y < (ssize_t) evaluate_image->rows; y++)
-    {
-      CacheView
-        *image_view;
-
-      const Image
-        *next;
+      for (y=0; y < (ssize_t) evaluate_image->rows; y++)
+      {
+        CacheView
+          *image_view;
 
-      const int
-        id = GetOpenMPThreadId();
+        const Image
+          *next;
 
-      register PixelInfo
-        *evaluate_pixel;
+        const int
+          id = GetOpenMPThreadId();
 
-      register Quantum
-        *restrict q;
+        register PixelChannels
+          *evaluate_pixel;
 
-      register ssize_t
-        x;
+        register Quantum
+          *restrict q;
 
-      if (status == MagickFalse)
-        continue;
-      q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,evaluate_image->columns,
-        1,exception);
-      if (q == (const Quantum *) NULL)
-        {
-          status=MagickFalse;
-          continue;
-        }
-      evaluate_pixel=evaluate_pixels[id];
-      for (x=0; x < (ssize_t) evaluate_image->columns; x++)
-      {
         register ssize_t
-          i;
+          x;
 
-        for (i=0; i < (ssize_t) number_images; i++)
-          evaluate_pixel[i]=zero;
-        next=images;
-        for (i=0; i < (ssize_t) number_images; i++)
+        if (status == MagickFalse)
+          continue;
+        q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,
+          evaluate_image->columns,1,exception);
+        if (q == (Quantum *) NULL)
+          {
+            status=MagickFalse;
+            continue;
+          }
+        evaluate_pixel=evaluate_pixels[id];
+        for (x=0; x < (ssize_t) evaluate_image->columns; x++)
         {
-          register const Quantum
-            *p;
-
-          image_view=AcquireCacheView(next);
-          p=GetCacheViewVirtualPixels(image_view,x,y,1,1,exception);
-          if (p == (const Quantum *) NULL)
+          register ssize_t
+            j,
+            k;
+
+          for (j=0; j < (ssize_t) number_images; j++)
+            for (k=0; k < MaxPixelChannels; k++)
+              evaluate_pixel[j].channel[k]=0.0;
+          next=images;
+          for (j=0; j < (ssize_t) number_images; j++)
+          {
+            register const Quantum
+              *p;
+
+            register ssize_t
+              i;
+
+            image_view=AcquireVirtualCacheView(next,exception);
+            p=GetCacheViewVirtualPixels(image_view,x,y,1,1,exception);
+            if (p == (const Quantum *) NULL)
+              {
+                image_view=DestroyCacheView(image_view);
+                break;
+              }
+            for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++)
             {
-              image_view=DestroyCacheView(image_view);
-              break;
+              PixelChannel
+                channel;
+
+              PixelTrait
+                evaluate_traits,
+                traits;
+
+              channel=GetPixelChannelMapChannel(evaluate_image,i);
+              evaluate_traits=GetPixelChannelMapTraits(evaluate_image,channel);
+              traits=GetPixelChannelMapTraits(next,channel);
+              if ((traits == UndefinedPixelTrait) ||
+                  (evaluate_traits == UndefinedPixelTrait))
+                continue;
+              if ((evaluate_traits & UpdatePixelTrait) == 0)
+                continue;
+              evaluate_pixel[j].channel[i]=ApplyEvaluateOperator(
+                random_info[id],GetPixelChannel(evaluate_image,channel,p),op,
+                evaluate_pixel[j].channel[i]);
             }
-          evaluate_pixel[i].red=ApplyEvaluateOperator(random_info[id],
-            GetPixelRed(next,p),op,evaluate_pixel[i].red);
-          evaluate_pixel[i].green=ApplyEvaluateOperator(random_info[id],
-            GetPixelGreen(next,p),op,evaluate_pixel[i].green);
-          evaluate_pixel[i].blue=ApplyEvaluateOperator(random_info[id],
-            GetPixelBlue(next,p),op,evaluate_pixel[i].blue);
-          if (evaluate_image->colorspace == CMYKColorspace)
-            evaluate_pixel[i].black=ApplyEvaluateOperator(random_info[id],
-              GetPixelBlack(next,p),op,evaluate_pixel[i].black);
-          evaluate_pixel[i].alpha=ApplyEvaluateOperator(random_info[id],
-            GetPixelAlpha(next,p),op,evaluate_pixel[i].alpha);
-          image_view=DestroyCacheView(image_view);
-          next=GetNextImageInList(next);
+            image_view=DestroyCacheView(image_view);
+            next=GetNextImageInList(next);
+          }
+          qsort((void *) evaluate_pixel,number_images,sizeof(*evaluate_pixel),
+            IntensityCompare);
+          for (k=0; k < (ssize_t) GetPixelChannels(evaluate_image); k++)
+            q[k]=ClampToQuantum(evaluate_pixel[j/2].channel[k]);
+          q+=GetPixelChannels(evaluate_image);
         }
-        qsort((void *) evaluate_pixel,number_images,sizeof(*evaluate_pixel),
-          IntensityCompare);
-        SetPixelRed(evaluate_image,
-          ClampToQuantum(evaluate_pixel[i/2].red),q);
-        SetPixelGreen(evaluate_image,
-          ClampToQuantum(evaluate_pixel[i/2].green),q);
-        SetPixelBlue(evaluate_image,
-          ClampToQuantum(evaluate_pixel[i/2].blue),q);
-        if (evaluate_image->colorspace == CMYKColorspace)
-          SetPixelBlack(evaluate_image,
-          ClampToQuantum(evaluate_pixel[i/2].black),q);
-        if (evaluate_image->matte == MagickFalse)
-          SetPixelAlpha(evaluate_image,
-            ClampToQuantum(evaluate_pixel[i/2].alpha),q);
-        else
-          SetPixelAlpha(evaluate_image,
-            ClampToQuantum(evaluate_pixel[i/2].alpha),q);
-        q+=GetPixelChannels(evaluate_image);
-      }
-      if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
-        status=MagickFalse;
-      if (images->progress_monitor != (MagickProgressMonitor) NULL)
-        {
-          MagickBooleanType
-            proceed;
+        if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
+          status=MagickFalse;
+        if (images->progress_monitor != (MagickProgressMonitor) NULL)
+          {
+            MagickBooleanType
+              proceed;
 
-#if defined(MAGICKCORE_OPENMP_SUPPORT)
-          #pragma omp critical (MagickCore_EvaluateImages)
+#if   defined(MAGICKCORE_OPENMP_SUPPORT)
+            #pragma omp critical (MagickCore_EvaluateImages)
 #endif
-          proceed=SetImageProgress(images,EvaluateImageTag,progress++,
-            evaluate_image->rows);
-          if (proceed == MagickFalse)
-            status=MagickFalse;
-        }
+            proceed=SetImageProgress(images,EvaluateImageTag,progress++,
+              evaluate_image->rows);
+            if (proceed == MagickFalse)
+              status=MagickFalse;
+          }
+      }
     }
   else
+    {
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(dynamic) shared(progress,status)
+      #pragma omp parallel for schedule(static) shared(progress,status) omp_concurrent(concurrent)
 #endif
-    for (y=0; y < (ssize_t) evaluate_image->rows; y++)
-    {
-      CacheView
-        *image_view;
+      for (y=0; y < (ssize_t) evaluate_image->rows; y++)
+      {
+        CacheView
+          *image_view;
 
-      const Image
-        *next;
+        const Image
+          *next;
 
-      const int
-        id = GetOpenMPThreadId();
+        const int
+          id = GetOpenMPThreadId();
 
-      register ssize_t
-        i,
-        x;
+        register ssize_t
+          i,
+          x;
 
-      register PixelInfo
-        *evaluate_pixel;
+        register PixelChannels
+          *evaluate_pixel;
 
-      register Quantum
-        *restrict q;
+        register Quantum
+          *restrict q;
 
-      if (status == MagickFalse)
-        continue;
-      q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,evaluate_image->columns,
-        1,exception);
-      if (q == (const Quantum *) NULL)
-        {
-          status=MagickFalse;
-          continue;
-        }
-      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 Quantum
-          *p;
+        ssize_t
+          j;
 
-        image_view=AcquireCacheView(next);
-        p=GetCacheViewVirtualPixels(image_view,0,y,next->columns,1,exception);
-        if (p == (const Quantum *) NULL)
+        if (status == MagickFalse)
+          continue;
+        q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,
+          evaluate_image->columns,1,exception);
+        if (q == (Quantum *) NULL)
           {
-            image_view=DestroyCacheView(image_view);
-            break;
+            status=MagickFalse;
+            continue;
           }
-        for (x=0; x < (ssize_t) next->columns; x++)
+        evaluate_pixel=evaluate_pixels[id];
+        for (j=0; j < (ssize_t) evaluate_image->columns; j++)
+          for (i=0; i < MaxPixelChannels; i++)
+            evaluate_pixel[j].channel[i]=0.0;
+        next=images;
+        for (j=0; j < (ssize_t) number_images; j++)
         {
-          evaluate_pixel[x].red=ApplyEvaluateOperator(random_info[id],
-            GetPixelRed(next,p),i == 0 ? AddEvaluateOperator : op,
-            evaluate_pixel[x].red);
-          evaluate_pixel[x].green=ApplyEvaluateOperator(random_info[id],
-            GetPixelGreen(next,p),i == 0 ? AddEvaluateOperator : op,
-              evaluate_pixel[x].green);
-          evaluate_pixel[x].blue=ApplyEvaluateOperator(random_info[id],
-            GetPixelBlue(next,p),i == 0 ? AddEvaluateOperator : op,
-              evaluate_pixel[x].blue);
-          if (evaluate_image->colorspace == CMYKColorspace)
-            evaluate_pixel[x].black=ApplyEvaluateOperator(random_info[id],
-              GetPixelBlack(next,p),i == 0 ? AddEvaluateOperator : op,
-              evaluate_pixel[x].black);
-          evaluate_pixel[x].alpha=ApplyEvaluateOperator(random_info[id],
-            GetPixelAlpha(next,p),i == 0 ? AddEvaluateOperator : op,
-            evaluate_pixel[x].alpha);
-          p+=GetPixelChannels(next);
+          register const Quantum
+            *p;
+
+          image_view=AcquireVirtualCacheView(next,exception);
+          p=GetCacheViewVirtualPixels(image_view,0,y,next->columns,1,exception);
+          if (p == (const Quantum *) NULL)
+            {
+              image_view=DestroyCacheView(image_view);
+              break;
+            }
+          for (x=0; x < (ssize_t) next->columns; x++)
+          {
+            register ssize_t
+              i;
+
+            if (GetPixelMask(next,p) != 0)
+              {
+                p+=GetPixelChannels(next);
+                continue;
+              }
+            for (i=0; i < (ssize_t) GetPixelChannels(next); i++)
+            {
+              PixelChannel
+                channel;
+
+              PixelTrait
+                evaluate_traits,
+                traits;
+
+              channel=GetPixelChannelMapChannel(evaluate_image,i);
+              traits=GetPixelChannelMapTraits(next,channel);
+              evaluate_traits=GetPixelChannelMapTraits(evaluate_image,channel);
+              if ((traits == UndefinedPixelTrait) ||
+                  (evaluate_traits == UndefinedPixelTrait))
+                continue;
+              if ((traits & UpdatePixelTrait) == 0)
+                continue;
+              evaluate_pixel[x].channel[i]=ApplyEvaluateOperator(
+                random_info[id],GetPixelChannel(evaluate_image,channel,p),j ==
+                0 ? AddEvaluateOperator : op,evaluate_pixel[x].channel[i]);
+            }
+            p+=GetPixelChannels(next);
+          }
+          image_view=DestroyCacheView(image_view);
+          next=GetNextImageInList(next);
         }
-        image_view=DestroyCacheView(image_view);
-        next=GetNextImageInList(next);
-      }
-      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].black/=number_images;
-          evaluate_pixel[x].alpha/=number_images;
+          register ssize_t
+             i;
+
+          switch (op)
+          {
+            case MeanEvaluateOperator:
+            {
+              for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++)
+                evaluate_pixel[x].channel[i]/=(MagickRealType) number_images;
+              break;
+            }
+            case MultiplyEvaluateOperator:
+            {
+              for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++)
+              {
+                register ssize_t
+                  j;
+
+                for (j=0; j < (ssize_t) (number_images-1); j++)
+                  evaluate_pixel[x].channel[i]*=QuantumScale;
+              }
+              break;
+            }
+            default:
+              break;
+          }
         }
-      for (x=0; x < (ssize_t) evaluate_image->columns; x++)
-      {
-        SetPixelRed(evaluate_image,ClampToQuantum(evaluate_pixel[x].red),q);
-        SetPixelGreen(evaluate_image,ClampToQuantum(evaluate_pixel[x].green),q);
-        SetPixelBlue(evaluate_image,ClampToQuantum(evaluate_pixel[x].blue),q);
-        if (evaluate_image->colorspace == CMYKColorspace)
-          SetPixelBlack(evaluate_image,ClampToQuantum(evaluate_pixel[x].black),
-            q);
-        if (evaluate_image->matte == MagickFalse)
-          SetPixelAlpha(evaluate_image,ClampToQuantum(evaluate_pixel[x].alpha),
-            q);
-        else
-          SetPixelAlpha(evaluate_image,ClampToQuantum(evaluate_pixel[x].alpha),
-            q);
-        q+=GetPixelChannels(evaluate_image);
-      }
-      if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
-        status=MagickFalse;
-      if (images->progress_monitor != (MagickProgressMonitor) NULL)
+        for (x=0; x < (ssize_t) evaluate_image->columns; x++)
         {
-          MagickBooleanType
-            proceed;
+          register ssize_t
+            i;
 
-#if defined(MAGICKCORE_OPENMP_SUPPORT)
-          #pragma omp critical (MagickCore_EvaluateImages)
-#endif
-          proceed=SetImageProgress(images,EvaluateImageTag,progress++,
-            evaluate_image->rows);
-          if (proceed == MagickFalse)
-            status=MagickFalse;
+          if (GetPixelMask(evaluate_image,q) != 0)
+            {
+              q+=GetPixelChannels(evaluate_image);
+              continue;
+            }
+          for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++)
+          {
+            PixelChannel
+              channel;
+
+            PixelTrait
+              traits;
+
+            channel=GetPixelChannelMapChannel(evaluate_image,i);
+            traits=GetPixelChannelMapTraits(evaluate_image,channel);
+            if (traits == UndefinedPixelTrait)
+              continue;
+            if ((traits & UpdatePixelTrait) == 0)
+              continue;
+            q[i]=ClampToQuantum(evaluate_pixel[x].channel[i]);
+          }
+          q+=GetPixelChannels(evaluate_image);
         }
+        if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
+          status=MagickFalse;
+        if (images->progress_monitor != (MagickProgressMonitor) NULL)
+          {
+            MagickBooleanType
+              proceed;
+
+#if   defined(MAGICKCORE_OPENMP_SUPPORT)
+            #pragma omp critical (MagickCore_EvaluateImages)
+#endif
+            proceed=SetImageProgress(images,EvaluateImageTag,progress++,
+              evaluate_image->rows);
+            if (proceed == MagickFalse)
+              status=MagickFalse;
+          }
+      }
     }
   evaluate_view=DestroyCacheView(evaluate_view);
   evaluate_pixels=DestroyPixelThreadSet(evaluate_pixels);
@@ -717,14 +780,14 @@ MagickExport Image *EvaluateImages(const Image *images,
   return(evaluate_image);
 }
 
-MagickExport MagickBooleanType EvaluateImageChannel(Image *image,
-  const ChannelType channel,const MagickEvaluateOperator op,const double value,
-  ExceptionInfo *exception)
+MagickExport MagickBooleanType EvaluateImage(Image *image,
+  const MagickEvaluateOperator op,const double value,ExceptionInfo *exception)
 {
   CacheView
     *image_view;
 
   MagickBooleanType
+    concurrent,
     status;
 
   MagickOffsetType
@@ -742,17 +805,16 @@ MagickExport MagickBooleanType EvaluateImageChannel(Image *image,
     (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);
-    }
+  if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
+    return(MagickFalse);
   status=MagickTrue;
   progress=0;
   random_info=AcquireRandomInfoThreadSet();
-  image_view=AcquireCacheView(image);
+  concurrent=GetRandomSecretKey(random_info[0]) == ~0UL ? MagickTrue :
+    MagickFalse;
+  image_view=AcquireAuthenticCacheView(image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
+  #pragma omp parallel for schedule(static,4) shared(progress,status) omp_concurrent(concurrent)
 #endif
   for (y=0; y < (ssize_t) image->rows; y++)
   {
@@ -768,35 +830,38 @@ MagickExport MagickBooleanType EvaluateImageChannel(Image *image,
     if (status == MagickFalse)
       continue;
     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
-    if (q == (const Quantum *) NULL)
+    if (q == (Quantum *) NULL)
       {
         status=MagickFalse;
         continue;
       }
     for (x=0; x < (ssize_t) image->columns; x++)
     {
-      if ((GetPixelRedTraits(image) & ActivePixelTrait) != 0)
-        SetPixelRed(image,ClampToQuantum(ApplyEvaluateOperator(
-          random_info[id],GetPixelRed(image,q),op,value)),q);
-      if ((GetPixelGreenTraits(image) & ActivePixelTrait) != 0)
-        SetPixelGreen(image,ClampToQuantum(ApplyEvaluateOperator(
-          random_info[id],GetPixelGreen(image,q),op,value)),q);
-      if ((GetPixelBlueTraits(image) & ActivePixelTrait) != 0)
-        SetPixelBlue(image,ClampToQuantum(ApplyEvaluateOperator(
-          random_info[id],GetPixelBlue(image,q),op,value)),q);
-      if (((GetPixelBlackTraits(image) & ActivePixelTrait) != 0) &&
-          (image->colorspace == CMYKColorspace))
-        SetPixelBlack(image,ClampToQuantum(ApplyEvaluateOperator(
-          random_info[id],GetPixelBlack(image,q),op,value)),q);
-      if ((GetPixelAlphaTraits(image) & ActivePixelTrait) != 0)
+      register ssize_t
+        i;
+
+      if (GetPixelMask(image,q) != 0)
         {
-          if (image->matte == MagickFalse)
-            SetPixelAlpha(image,ClampToQuantum(ApplyEvaluateOperator(
-              random_info[id],GetPixelAlpha(image,q),op,value)),q);
-          else
-            SetPixelAlpha(image,ClampToQuantum(ApplyEvaluateOperator(
-              random_info[id],GetPixelAlpha(image,q),op,value)),q);
+          q+=GetPixelChannels(image);
+          continue;
         }
+      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
+      {
+        PixelChannel
+          channel;
+
+        PixelTrait
+          traits;
+
+        channel=GetPixelChannelMapChannel(image,i);
+        traits=GetPixelChannelMapTraits(image,channel);
+        if (traits == UndefinedPixelTrait)
+          continue;
+        if ((traits & CopyPixelTrait) != 0)
+          continue;
+        q[i]=ClampToQuantum(ApplyEvaluateOperator(random_info[id],q[i],op,
+          value));
+      }
       q+=GetPixelChannels(image);
     }
     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
@@ -807,7 +872,7 @@ MagickExport MagickBooleanType EvaluateImageChannel(Image *image,
           proceed;
 
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_EvaluateImageChannel)
+  #pragma omp critical (MagickCore_EvaluateImage)
 #endif
         proceed=SetImageProgress(image,EvaluateImageTag,progress++,image->rows);
         if (proceed == MagickFalse)
@@ -835,22 +900,16 @@ MagickExport MagickBooleanType EvaluateImageChannel(Image *image,
 %  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:
+%  The format of the FunctionImage 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.
@@ -876,63 +935,79 @@ static Quantum ApplyFunction(Quantum pixel,const MagickFunction function,
     case PolynomialFunction:
     {
       /*
-       * Polynomial
-       * Parameters:   polynomial constants,  highest to lowest order
-       *   For example:      c0*x^3 + c1*x^2 + c2*x  + c3
-       */
+        Polynomial: polynomial constants, highest to lowest order (e.g. 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;
+        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 ) );
+      MagickRealType
+        amplitude,
+        bias,
+        frequency,
+        phase;
+
+      /*
+        Sinusoid: frequency, phase, amplitude, bias.
+      */
+      frequency=(number_parameters >= 1) ? parameters[0] : 1.0;
+      phase=(number_parameters >= 2) ? parameters[1] : 0.0;
+      amplitude=(number_parameters >= 3) ? parameters[2] : 0.5;
+      bias=(number_parameters >= 4) ? parameters[3] : 0.5;
+      result=(MagickRealType) (QuantumRange*(amplitude*sin((double) (2.0*
+        MagickPI*(frequency*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);
+      MagickRealType
+        bias,
+        center,
+        range,
+        width;
+
+      /*
+        Arcsin (peged at range limits for invalid results): width, center,
+        range, and 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;
+        result=bias-range/2.0;
       else
-        result=(MagickRealType) (range/MagickPI*asin((double) result)+bias);
-      result *= QuantumRange;
+        if (result >= 1.0)
+          result=bias+range/2.0;
+        else
+          result=(MagickRealType) (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;
+      MagickRealType
+        center,
+        bias,
+        range,
+        slope;
+
+      /*
+        Arctan: slope, center, range, and 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=(MagickRealType) (MagickPI*slope*(QuantumScale*pixel-center));
       result=(MagickRealType) (QuantumRange*(range/MagickPI*atan((double)
-                  result) + bias ) );
+        result)+bias));
       break;
     }
     case UndefinedFunction:
@@ -944,19 +1019,6 @@ static Quantum ApplyFunction(Quantum pixel,const MagickFunction function,
 MagickExport MagickBooleanType FunctionImage(Image *image,
   const MagickFunction function,const size_t number_parameters,
   const double *parameters,ExceptionInfo *exception)
-{
-  MagickBooleanType
-    status;
-
-  status=FunctionImageChannel(image,CompositeChannels,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 "
 
@@ -978,57 +1040,57 @@ MagickExport MagickBooleanType FunctionImageChannel(Image *image,
     (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);
-    }
+  if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
+    return(MagickFalse);
   status=MagickTrue;
   progress=0;
-  image_view=AcquireCacheView(image);
+  image_view=AcquireAuthenticCacheView(image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
+  #pragma omp parallel for schedule(static,4) shared(progress,status)
 #endif
   for (y=0; y < (ssize_t) image->rows; y++)
   {
-    register ssize_t
-      x;
-
     register Quantum
       *restrict q;
 
+    register ssize_t
+      x;
+
     if (status == MagickFalse)
       continue;
     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
-    if (q == (const Quantum *) NULL)
+    if (q == (Quantum *) NULL)
       {
         status=MagickFalse;
         continue;
       }
     for (x=0; x < (ssize_t) image->columns; x++)
     {
-      if ((GetPixelRedTraits(image) & ActivePixelTrait) != 0)
-        SetPixelRed(image,ApplyFunction(GetPixelRed(image,q),function,
-          number_parameters,parameters,exception),q);
-      if ((GetPixelGreenTraits(image) & ActivePixelTrait) != 0)
-        SetPixelGreen(image,ApplyFunction(GetPixelGreen(image,q),function,
-          number_parameters,parameters,exception),q);
-      if ((GetPixelBlueTraits(image) & ActivePixelTrait) != 0)
-        SetPixelBlue(image,ApplyFunction(GetPixelBlue(image,q),function,
-          number_parameters,parameters,exception),q);
-      if (((GetPixelBlackTraits(image) & ActivePixelTrait) != 0) &&
-          (image->colorspace == CMYKColorspace))
-        SetPixelBlack(image,ApplyFunction(GetPixelBlack(image,q),function,
-          number_parameters,parameters,exception),q);
-      if ((GetPixelAlphaTraits(image) & ActivePixelTrait) != 0)
+      register ssize_t
+        i;
+
+      if (GetPixelMask(image,q) != 0)
         {
-          if (image->matte == MagickFalse)
-            SetPixelAlpha(image,ApplyFunction(GetPixelAlpha(image,q),function,
-              number_parameters,parameters,exception),q);
-          else
-            SetPixelAlpha(image,ApplyFunction(GetPixelAlpha(image,q),function,
-              number_parameters,parameters,exception),q);
+          q+=GetPixelChannels(image);
+          continue;
         }
+      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
+      {
+        PixelChannel
+          channel;
+
+        PixelTrait
+          traits;
+
+        channel=GetPixelChannelMapChannel(image,i);
+        traits=GetPixelChannelMapTraits(image,channel);
+        if (traits == UndefinedPixelTrait)
+          continue;
+        if ((traits & UpdatePixelTrait) == 0)
+          continue;
+        q[i]=ApplyFunction(q[i],function,number_parameters,parameters,
+          exception);
+      }
       q+=GetPixelChannels(image);
     }
     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
@@ -1039,7 +1101,7 @@ MagickExport MagickBooleanType FunctionImageChannel(Image *image,
           proceed;
 
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_FunctionImageChannel)
+  #pragma omp critical (MagickCore_FunctionImage)
 #endif
         proceed=SetImageProgress(image,FunctionImageTag,progress++,image->rows);
         if (proceed == MagickFalse)
@@ -1055,26 +1117,23 @@ MagickExport MagickBooleanType FunctionImageChannel(Image *image,
 %                                                                             %
 %                                                                             %
 %                                                                             %
-+   G e t I m a g e C h a n n e l E x t r e m a                               %
+%   G e t I m a g e E x t r e m a                                             %
 %                                                                             %
 %                                                                             %
 %                                                                             %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-%  GetImageChannelExtrema() returns the extrema of one or more image channels.
+%  GetImageExtrema() returns the extrema of one or more image channels.
 %
-%  The format of the GetImageChannelExtrema method is:
+%  The format of the GetImageExtrema method is:
 %
-%      MagickBooleanType GetImageChannelExtrema(const Image *image,
-%        const ChannelType channel,size_t *minima,size_t *maxima,
-%        ExceptionInfo *exception)
+%      MagickBooleanType GetImageExtrema(const Image *image,size_t *minima,
+%        size_t *maxima,ExceptionInfo *exception)
 %
 %  A description of each parameter follows:
 %
 %    o image: the image.
 %
-%    o channel: the channel.
-%
 %    o minima: the minimum value in the channel.
 %
 %    o maxima: the maximum value in the channel.
@@ -1082,16 +1141,8 @@ MagickExport MagickBooleanType FunctionImageChannel(Image *image,
 %    o exception: return any errors or warnings in this structure.
 %
 */
-
 MagickExport MagickBooleanType GetImageExtrema(const Image *image,
   size_t *minima,size_t *maxima,ExceptionInfo *exception)
-{
-  return(GetImageChannelExtrema(image,CompositeChannels,minima,maxima,exception));
-}
-
-MagickExport MagickBooleanType GetImageChannelExtrema(const Image *image,
-  const ChannelType channel,size_t *minima,size_t *maxima,
-  ExceptionInfo *exception)
 {
   double
     max,
@@ -1104,7 +1155,7 @@ MagickExport MagickBooleanType GetImageChannelExtrema(const Image *image,
   assert(image->signature == MagickSignature);
   if (image->debug != MagickFalse)
     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
-  status=GetImageChannelRange(image,channel,&min,&max,exception);
+  status=GetImageRange(image,&min,&max,exception);
   *minima=(size_t) ceil(min-0.5);
   *maxima=(size_t) floor(max+0.5);
   return(status);
@@ -1115,27 +1166,24 @@ MagickExport MagickBooleanType GetImageChannelExtrema(const Image *image,
 %                                                                             %
 %                                                                             %
 %                                                                             %
-%   G e t I m a g e C h a n n e l M e a n                                     %
+%   G e t I m a g e M e a n                                                   %
 %                                                                             %
 %                                                                             %
 %                                                                             %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-%  GetImageChannelMean() returns the mean and standard deviation of one or more
+%  GetImageMean() returns the mean and standard deviation of one or more
 %  image channels.
 %
-%  The format of the GetImageChannelMean method is:
+%  The format of the GetImageMean method is:
 %
-%      MagickBooleanType GetImageChannelMean(const Image *image,
-%        const ChannelType channel,double *mean,double *standard_deviation,
-%        ExceptionInfo *exception)
+%      MagickBooleanType GetImageMean(const Image *image,double *mean,
+%        double *standard_deviation,ExceptionInfo *exception)
 %
 %  A description of each parameter follows:
 %
 %    o image: the image.
 %
-%    o channel: the channel.
-%
 %    o mean: the average value in the channel.
 %
 %    o standard_deviation: the standard deviation of the channel.
@@ -1143,95 +1191,54 @@ MagickExport MagickBooleanType GetImageChannelExtrema(const Image *image,
 %    o exception: return any errors or warnings in this structure.
 %
 */
-
 MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean,
   double *standard_deviation,ExceptionInfo *exception)
-{
-  MagickBooleanType
-    status;
-
-  status=GetImageChannelMean(image,CompositeChannels,mean,standard_deviation,
-    exception);
-  return(status);
-}
-
-MagickExport MagickBooleanType GetImageChannelMean(const Image *image,
-  const ChannelType channel,double *mean,double *standard_deviation,
-  ExceptionInfo *exception)
 {
   ChannelStatistics
     *channel_statistics;
 
+  register ssize_t
+    i;
+
   size_t
-    channels;
+    area;
 
   assert(image != (Image *) NULL);
   assert(image->signature == MagickSignature);
   if (image->debug != MagickFalse)
     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
-  channel_statistics=GetImageChannelStatistics(image,exception);
+  channel_statistics=GetImageStatistics(image,exception);
   if (channel_statistics == (ChannelStatistics *) NULL)
     return(MagickFalse);
-  channels=0;
-  channel_statistics[CompositeChannels].mean=0.0;
-  channel_statistics[CompositeChannels].standard_deviation=0.0;
-  if ((GetPixelRedTraits(image) & ActivePixelTrait) != 0)
-    {
-      channel_statistics[CompositeChannels].mean+=
-        channel_statistics[RedChannel].mean;
-      channel_statistics[CompositeChannels].standard_deviation+=
-        channel_statistics[RedChannel].variance-
-        channel_statistics[RedChannel].mean*
-        channel_statistics[RedChannel].mean;
-      channels++;
-    }
-  if ((GetPixelGreenTraits(image) & ActivePixelTrait) != 0)
-    {
-      channel_statistics[CompositeChannels].mean+=
-        channel_statistics[GreenChannel].mean;
-      channel_statistics[CompositeChannels].standard_deviation+=
-        channel_statistics[GreenChannel].variance-
-        channel_statistics[GreenChannel].mean*
-        channel_statistics[GreenChannel].mean;
-      channels++;
-    }
-  if ((GetPixelBlueTraits(image) & ActivePixelTrait) != 0)
-    {
-      channel_statistics[CompositeChannels].mean+=
-        channel_statistics[BlueChannel].mean;
-      channel_statistics[CompositeChannels].standard_deviation+=
-        channel_statistics[BlueChannel].variance-
-        channel_statistics[BlueChannel].mean*
-        channel_statistics[BlueChannel].mean;
-      channels++;
-    }
-  if (((GetPixelBlackTraits(image) & ActivePixelTrait) != 0) &&
-      (image->colorspace == CMYKColorspace))
-    {
-      channel_statistics[CompositeChannels].mean+=
-        channel_statistics[BlackChannel].mean;
-      channel_statistics[CompositeChannels].standard_deviation+=
-        channel_statistics[BlackChannel].variance-
-        channel_statistics[BlackChannel].mean*
-        channel_statistics[BlackChannel].mean;
-      channels++;
-    }
-  if (((GetPixelAlphaTraits(image) & ActivePixelTrait) != 0) &&
-      (image->matte != MagickFalse))
-    {
-      channel_statistics[CompositeChannels].mean+=
-        channel_statistics[AlphaChannel].mean;
-      channel_statistics[CompositeChannels].standard_deviation+=
-        channel_statistics[AlphaChannel].variance-
-        channel_statistics[AlphaChannel].mean*
-        channel_statistics[AlphaChannel].mean;
-      channels++;
-    }
-  channel_statistics[CompositeChannels].mean/=channels;
-  channel_statistics[CompositeChannels].standard_deviation=
-    sqrt(channel_statistics[CompositeChannels].standard_deviation/channels);
-  *mean=channel_statistics[CompositeChannels].mean;
-  *standard_deviation=channel_statistics[CompositeChannels].standard_deviation;
+  area=0;
+  channel_statistics[CompositePixelChannel].mean=0.0;
+  channel_statistics[CompositePixelChannel].standard_deviation=0.0;
+  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
+  {
+    PixelChannel
+      channel;
+
+    PixelTrait
+      traits;
+
+    channel=GetPixelChannelMapChannel(image,i);
+    traits=GetPixelChannelMapTraits(image,channel);
+    if (traits == UndefinedPixelTrait)
+      continue;
+    if ((traits & UpdatePixelTrait) == 0)
+      continue;
+    channel_statistics[CompositePixelChannel].mean+=channel_statistics[i].mean;
+    channel_statistics[CompositePixelChannel].standard_deviation+=
+      channel_statistics[i].variance-channel_statistics[i].mean*
+      channel_statistics[i].mean;
+    area++;
+  }
+  channel_statistics[CompositePixelChannel].mean/=area;
+  channel_statistics[CompositePixelChannel].standard_deviation=
+    sqrt(channel_statistics[CompositePixelChannel].standard_deviation/area);
+  *mean=channel_statistics[CompositePixelChannel].mean;
+  *standard_deviation=
+    channel_statistics[CompositePixelChannel].standard_deviation;
   channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
     channel_statistics);
   return(MagickTrue);
@@ -1242,27 +1249,24 @@ MagickExport MagickBooleanType GetImageChannelMean(const Image *image,
 %                                                                             %
 %                                                                             %
 %                                                                             %
-%   G e t I m a g e C h a n n e l K u r t o s i s                             %
+%   G e t I m a g e K u r t o s i s                                           %
 %                                                                             %
 %                                                                             %
 %                                                                             %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-%  GetImageChannelKurtosis() returns the kurtosis and skewness of one or more
+%  GetImageKurtosis() returns the kurtosis and skewness of one or more
 %  image channels.
 %
-%  The format of the GetImageChannelKurtosis method is:
+%  The format of the GetImageKurtosis method is:
 %
-%      MagickBooleanType GetImageChannelKurtosis(const Image *image,
-%        const ChannelType channel,double *kurtosis,double *skewness,
-%        ExceptionInfo *exception)
+%      MagickBooleanType GetImageKurtosis(const Image *image,double *kurtosis,
+%        double *skewness,ExceptionInfo *exception)
 %
 %  A description of each parameter follows:
 %
 %    o image: the image.
 %
-%    o channel: the channel.
-%
 %    o kurtosis: the kurtosis of the channel.
 %
 %    o skewness: the skewness of the channel.
@@ -1270,22 +1274,12 @@ MagickExport MagickBooleanType GetImageChannelMean(const Image *image,
 %    o exception: return any errors or warnings in this structure.
 %
 */
-
 MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
   double *kurtosis,double *skewness,ExceptionInfo *exception)
 {
-  MagickBooleanType
-    status;
-
-  status=GetImageChannelKurtosis(image,CompositeChannels,kurtosis,skewness,
-    exception);
-  return(status);
-}
+  CacheView
+    *image_view;
 
-MagickExport MagickBooleanType GetImageChannelKurtosis(const Image *image,
-  const ChannelType channel,double *kurtosis,double *skewness,
-  ExceptionInfo *exception)
-{
   double
     area,
     mean,
@@ -1294,6 +1288,9 @@ MagickExport MagickBooleanType GetImageChannelKurtosis(const Image *image,
     sum_cubes,
     sum_fourth_power;
 
+  MagickBooleanType
+    status;
+
   ssize_t
     y;
 
@@ -1301,6 +1298,7 @@ MagickExport MagickBooleanType GetImageChannelKurtosis(const Image *image,
   assert(image->signature == MagickSignature);
   if (image->debug != MagickFalse)
     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
+  status=MagickTrue;
   *kurtosis=0.0;
   *skewness=0.0;
   area=0.0;
@@ -1309,6 +1307,10 @@ MagickExport MagickBooleanType GetImageChannelKurtosis(const Image *image,
   sum_squares=0.0;
   sum_cubes=0.0;
   sum_fourth_power=0.0;
+  image_view=AcquireVirtualCacheView(image,exception);
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+  #pragma omp parallel for schedule(static) shared(status)
+#endif
   for (y=0; y < (ssize_t) image->rows; y++)
   {
     register const Quantum
@@ -1317,71 +1319,53 @@ MagickExport MagickBooleanType GetImageChannelKurtosis(const Image *image,
     register ssize_t
       x;
 
-    p=GetVirtualPixels(image,0,y,image->columns,1,exception);
+    if (status == MagickFalse)
+      continue;
+    p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
     if (p == (const Quantum *) NULL)
-      break;
+      {
+        status=MagickFalse;
+        continue;
+      }
     for (x=0; x < (ssize_t) image->columns; x++)
     {
-      if ((GetPixelRedTraits(image) & ActivePixelTrait) != 0)
-        {
-          mean+=GetPixelRed(image,p);
-          sum_squares+=(double) GetPixelRed(image,p)*GetPixelRed(image,p);
-          sum_cubes+=(double) GetPixelRed(image,p)*GetPixelRed(image,p)*
-            GetPixelRed(image,p);
-          sum_fourth_power+=(double) GetPixelRed(image,p)*
-            GetPixelRed(image,p)*GetPixelRed(image,p)*GetPixelRed(image,p);
-          area++;
-        }
-      if ((GetPixelGreenTraits(image) & ActivePixelTrait) != 0)
-        {
-          mean+=GetPixelGreen(image,p);
-          sum_squares+=(double) GetPixelGreen(image,p)*GetPixelGreen(image,p);
-          sum_cubes+=(double) GetPixelGreen(image,p)*GetPixelGreen(image,p)*
-            GetPixelGreen(image,p);
-          sum_fourth_power+=(double) GetPixelGreen(image,p)*
-            GetPixelGreen(image,p)*GetPixelGreen(image,p)*
-            GetPixelGreen(image,p);
-          area++;
-        }
-      if ((GetPixelBlueTraits(image) & ActivePixelTrait) != 0)
-        {
-          mean+=GetPixelBlue(image,p);
-          sum_squares+=(double) GetPixelBlue(image,p)*GetPixelBlue(image,p);
-          sum_cubes+=(double) GetPixelBlue(image,p)*GetPixelBlue(image,p)*
-            GetPixelBlue(image,p);
-          sum_fourth_power+=(double) GetPixelBlue(image,p)*
-            GetPixelBlue(image,p)*GetPixelBlue(image,p)*
-            GetPixelBlue(image,p);
-          area++;
-        }
-      if (((GetPixelBlackTraits(image) & ActivePixelTrait) != 0) &&
-          (image->colorspace == CMYKColorspace))
+      register ssize_t
+        i;
+
+      if (GetPixelMask(image,p) != 0)
         {
-          mean+=GetPixelBlack(image,p);
-          sum_squares+=(double) GetPixelBlack(image,p)*GetPixelBlack(image,p);
-          sum_cubes+=(double) GetPixelBlack(image,p)*GetPixelBlack(image,p)*
-            GetPixelBlack(image,p);
-          sum_fourth_power+=(double) GetPixelBlack(image,p)*
-            GetPixelBlack(image,p)*GetPixelBlack(image,p)*
-            GetPixelBlack(image,p);
-          area++;
+          p+=GetPixelChannels(image);
+          continue;
         }
-      if ((GetPixelAlphaTraits(image) & ActivePixelTrait) != 0)
+      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
+      {
+        PixelChannel
+          channel;
+
+        PixelTrait
+          traits;
+
+        channel=GetPixelChannelMapChannel(image,i);
+        traits=GetPixelChannelMapTraits(image,channel);
+        if (traits == UndefinedPixelTrait)
+          continue;
+        if ((traits & UpdatePixelTrait) == 0)
+          continue;
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+        #pragma omp critical (MagickCore_GetImageKurtosis)
+#endif
         {
-          mean+=GetPixelAlpha(image,p);
-          sum_squares+=(double) GetPixelAlpha(image,p)*GetPixelAlpha(image,p);
-          sum_cubes+=(double) GetPixelAlpha(image,p)*GetPixelAlpha(image,p)*
-            GetPixelAlpha(image,p);
-          sum_fourth_power+=(double) GetPixelAlpha(image,p)*
-            GetPixelAlpha(image,p)*GetPixelAlpha(image,p)*
-            GetPixelAlpha(image,p);
+          mean+=p[i];
+          sum_squares+=(double) p[i]*p[i];
+          sum_cubes+=(double) p[i]*p[i]*p[i];
+          sum_fourth_power+=(double) p[i]*p[i]*p[i]*p[i];
           area++;
         }
+      }
       p+=GetPixelChannels(image);
     }
   }
-  if (y < (ssize_t) image->rows)
-    return(MagickFalse);
+  image_view=DestroyCacheView(image_view);
   if (area != 0.0)
     {
       mean/=area;
@@ -1400,7 +1384,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 == (ssize_t) image->rows ? MagickTrue : MagickFalse);
+  return(status);
 }
 \f
 /*
@@ -1408,26 +1392,23 @@ MagickExport MagickBooleanType GetImageChannelKurtosis(const Image *image,
 %                                                                             %
 %                                                                             %
 %                                                                             %
-%   G e t I m a g e C h a n n e l R a n g e                                   %
+%   G e t I m a g e R a n g e                                                 %
 %                                                                             %
 %                                                                             %
 %                                                                             %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-%  GetImageChannelRange() returns the range of one or more image channels.
+%  GetImageRange() returns the range of one or more image channels.
 %
-%  The format of the GetImageChannelRange method is:
+%  The format of the GetImageRange method is:
 %
-%      MagickBooleanType GetImageChannelRange(const Image *image,
-%        const ChannelType channel,double *minima,double *maxima,
-%        ExceptionInfo *exception)
+%      MagickBooleanType GetImageRange(const Image *image,double *minima,
+%        double *maxima,ExceptionInfo *exception)
 %
 %  A description of each parameter follows:
 %
 %    o image: the image.
 %
-%    o channel: the channel.
-%
 %    o minima: the minimum value in the channel.
 %
 %    o maxima: the maximum value in the channel.
@@ -1435,19 +1416,14 @@ MagickExport MagickBooleanType GetImageChannelKurtosis(const Image *image,
 %    o exception: return any errors or warnings in this structure.
 %
 */
-
-MagickExport MagickBooleanType GetImageRange(const Image *image,
-  double *minima,double *maxima,ExceptionInfo *exception)
+MagickExport MagickBooleanType GetImageRange(const Image *image,double *minima,
+  double *maxima,ExceptionInfo *exception)
 {
-  return(GetImageChannelRange(image,CompositeChannels,minima,maxima,exception));
-}
+  CacheView
+    *image_view;
 
-MagickExport MagickBooleanType GetImageChannelRange(const Image *image,
-  const ChannelType channel,double *minima,double *maxima,
-  ExceptionInfo *exception)
-{
-  PixelInfo
-    pixel;
+  MagickBooleanType
+    status;
 
   ssize_t
     y;
@@ -1456,9 +1432,13 @@ MagickExport MagickBooleanType GetImageChannelRange(const Image *image,
   assert(image->signature == MagickSignature);
   if (image->debug != MagickFalse)
     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
-  *maxima=(-1.0E-37);
-  *minima=1.0E+37;
-  GetPixelInfo(image,&pixel);
+  status=MagickTrue;
+  *maxima=(-MagickHuge);
+  *minima=MagickHuge;
+  image_view=AcquireVirtualCacheView(image,exception);
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+  #pragma omp parallel for schedule(static) shared(status)
+#endif
   for (y=0; y < (ssize_t) image->rows; y++)
   {
     register const Quantum
@@ -1467,52 +1447,53 @@ MagickExport MagickBooleanType GetImageChannelRange(const Image *image,
     register ssize_t
       x;
 
-    p=GetVirtualPixels(image,0,y,image->columns,1,exception);
+    if (status == MagickFalse)
+      continue;
+    p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
     if (p == (const Quantum *) NULL)
-      break;
+      {
+        status=MagickFalse;
+        continue;
+      }
     for (x=0; x < (ssize_t) image->columns; x++)
     {
-      SetPixelInfo(image,p,&pixel);
-      if ((GetPixelRedTraits(image) & ActivePixelTrait) != 0)
-        {
-          if (pixel.red < *minima)
-            *minima=(double) pixel.red;
-          if (pixel.red > *maxima)
-            *maxima=(double) pixel.red;
-        }
-      if ((GetPixelGreenTraits(image) & ActivePixelTrait) != 0)
-        {
-          if (pixel.green < *minima)
-            *minima=(double) pixel.green;
-          if (pixel.green > *maxima)
-            *maxima=(double) pixel.green;
-        }
-      if ((GetPixelBlueTraits(image) & ActivePixelTrait) != 0)
-        {
-          if (pixel.blue < *minima)
-            *minima=(double) pixel.blue;
-          if (pixel.blue > *maxima)
-            *maxima=(double) pixel.blue;
-        }
-      if (((GetPixelBlackTraits(image) & ActivePixelTrait) != 0) &&
-          (image->colorspace == CMYKColorspace))
+      register ssize_t
+        i;
+
+      if (GetPixelMask(image,p) != 0)
         {
-          if (pixel.black < *minima)
-            *minima=(double) pixel.black;
-          if (pixel.black > *maxima)
-            *maxima=(double) pixel.black;
+          p+=GetPixelChannels(image);
+          continue;
         }
-      if ((GetPixelAlphaTraits(image) & ActivePixelTrait) != 0)
+      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
+      {
+        PixelChannel
+          channel;
+
+        PixelTrait
+          traits;
+
+        channel=GetPixelChannelMapChannel(image,i);
+        traits=GetPixelChannelMapTraits(image,channel);
+        if (traits == UndefinedPixelTrait)
+          continue;
+        if ((traits & UpdatePixelTrait) == 0)
+          continue;
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+        #pragma omp critical (MagickCore_GetImageRange)
+#endif
         {
-          if (pixel.alpha < *minima)
-            *minima=(double) pixel.alpha;
-          if (pixel.alpha > *maxima)
-            *maxima=(double) pixel.alpha;
+          if ((double) p[i] < *minima)
+            *minima=(double) p[i];
+          if ((double) p[i] > *maxima)
+            *maxima=(double) p[i];
         }
+      }
       p+=GetPixelChannels(image);
     }
   }
-  return(y == (ssize_t) image->rows ? MagickTrue : MagickFalse);
+  image_view=DestroyCacheView(image_view);
+  return(status);
 }
 \f
 /*
@@ -1520,25 +1501,25 @@ MagickExport MagickBooleanType GetImageChannelRange(const Image *image,
 %                                                                             %
 %                                                                             %
 %                                                                             %
-%   G e t I m a g e C h a n n e l S t a t i s t i c s                         %
+%   G e t I m a g e S t a t i s t i c s                                       %
 %                                                                             %
 %                                                                             %
 %                                                                             %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-%  GetImageChannelStatistics() returns statistics for each channel in the
+%  GetImageStatistics() returns statistics for each channel in the
 %  image.  The statistics include the channel depth, its minima, maxima, mean,
 %  standard deviation, kurtosis and skewness.  You can access the red channel
 %  mean, for example, like this:
 %
-%      channel_statistics=GetImageChannelStatistics(image,exception);
-%      red_mean=channel_statistics[RedChannel].mean;
+%      channel_statistics=GetImageStatistics(image,exception);
+%      red_mean=channel_statistics[RedPixelChannel].mean;
 %
 %  Use MagickRelinquishMemory() to free the statistics buffer.
 %
-%  The format of the GetImageChannelStatistics method is:
+%  The format of the GetImageStatistics method is:
 %
-%      ChannelStatistics *GetImageChannelStatistics(const Image *image,
+%      ChannelStatistics *GetImageStatistics(const Image *image,
 %        ExceptionInfo *exception)
 %
 %  A description of each parameter follows:
@@ -1548,7 +1529,33 @@ MagickExport MagickBooleanType GetImageChannelRange(const Image *image,
 %    o exception: return any errors or warnings in this structure.
 %
 */
-MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
+
+static size_t GetImageChannels(const Image *image)
+{
+  register ssize_t
+    i;
+
+  size_t
+    channels;
+
+  channels=0;
+  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
+  {
+    PixelChannel
+      channel;
+
+    PixelTrait
+      traits;
+
+    channel=GetPixelChannelMapChannel(image,i);
+    traits=GetPixelChannelMapTraits(image,channel);
+    if ((traits & UpdatePixelTrait) != 0)
+      channels++;
+  }
+  return(channels);
+}
+
+MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
   ExceptionInfo *exception)
 {
   ChannelStatistics
@@ -1568,8 +1575,7 @@ MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
 
   size_t
     channels,
-    depth,
-    length;
+    depth;
 
   ssize_t
     y;
@@ -1578,18 +1584,17 @@ MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
   assert(image->signature == MagickSignature);
   if (image->debug != MagickFalse)
     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
-  length=CompositeChannels+1UL;
-  channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(length,
-    sizeof(*channel_statistics));
+  channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(
+    MaxPixelChannels+1,sizeof(*channel_statistics));
   if (channel_statistics == (ChannelStatistics *) NULL)
     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
-  (void) ResetMagickMemory(channel_statistics,0,length*
+  (void) ResetMagickMemory(channel_statistics,0,(MaxPixelChannels+1)*
     sizeof(*channel_statistics));
-  for (i=0; i <= (ssize_t) CompositeChannels; i++)
+  for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
   {
     channel_statistics[i].depth=1;
-    channel_statistics[i].maxima=(-1.0E-37);
-    channel_statistics[i].minima=1.0E+37;
+    channel_statistics[i].maxima=(-MagickHuge);
+    channel_statistics[i].minima=MagickHuge;
   }
   for (y=0; y < (ssize_t) image->rows; y++)
   {
@@ -1602,155 +1607,55 @@ MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
     if (p == (const Quantum *) NULL)
       break;
-    for (x=0; x < (ssize_t) image->columns; )
+    for (x=0; x < (ssize_t) image->columns; x++)
     {
-      if (channel_statistics[RedChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
-        {
-          depth=channel_statistics[RedChannel].depth;
-          range=GetQuantumRange(depth);
-          status=GetPixelRed(image,p) != ScaleAnyToQuantum(ScaleQuantumToAny(
-            GetPixelRed(image,p),range),range) ? MagickTrue : MagickFalse;
-          if (status != MagickFalse)
-            {
-              channel_statistics[RedChannel].depth++;
-              continue;
-            }
-        }
-      if (channel_statistics[GreenChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
-        {
-          depth=channel_statistics[GreenChannel].depth;
-          range=GetQuantumRange(depth);
-          status=GetPixelGreen(image,p) != ScaleAnyToQuantum(ScaleQuantumToAny(
-            GetPixelGreen(image,p),range),range) ? MagickTrue : MagickFalse;
-          if (status != MagickFalse)
-            {
-              channel_statistics[GreenChannel].depth++;
-              continue;
-            }
-        }
-      if (channel_statistics[BlueChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
+      register ssize_t
+        i;
+
+      if (GetPixelMask(image,p) != 0)
         {
-          depth=channel_statistics[BlueChannel].depth;
-          range=GetQuantumRange(depth);
-          status=GetPixelBlue(image,p) != ScaleAnyToQuantum(ScaleQuantumToAny(
-            GetPixelBlue(image,p),range),range) ? MagickTrue : MagickFalse;
-          if (status != MagickFalse)
-            {
-              channel_statistics[BlueChannel].depth++;
-              continue;
-            }
+          p+=GetPixelChannels(image);
+          continue;
         }
-      if (image->matte != MagickFalse)
-        {
-          if (channel_statistics[AlphaChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
-            {
-              depth=channel_statistics[AlphaChannel].depth;
-              range=GetQuantumRange(depth);
-              status=GetPixelAlpha(image,p) != ScaleAnyToQuantum(
-                ScaleQuantumToAny(GetPixelAlpha(image,p),range),range) ?
-                MagickTrue : MagickFalse;
-              if (status != MagickFalse)
-                {
-                  channel_statistics[AlphaChannel].depth++;
-                  continue;
-                }
-            }
+      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
+      {
+        PixelChannel
+          channel;
+
+        PixelTrait
+          traits;
+
+        channel=GetPixelChannelMapChannel(image,i);
+        traits=GetPixelChannelMapTraits(image,channel);
+        if (traits == UndefinedPixelTrait)
+          continue;
+        if (channel_statistics[channel].depth != MAGICKCORE_QUANTUM_DEPTH)
+          {
+            depth=channel_statistics[channel].depth;
+            range=GetQuantumRange(depth);
+            status=p[i] != ScaleAnyToQuantum(ScaleQuantumToAny(p[i],range),
+              range) ? MagickTrue : MagickFalse;
+            if (status != MagickFalse)
+              {
+                channel_statistics[channel].depth++;
+                continue;
+              }
           }
-      if (image->colorspace == CMYKColorspace)
-        {
-          if (channel_statistics[BlackChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
-            {
-              depth=channel_statistics[BlackChannel].depth;
-              range=GetQuantumRange(depth);
-              status=GetPixelBlack(image,p) != ScaleAnyToQuantum(
-                ScaleQuantumToAny(GetPixelBlack(image,p),range),range) ?
-                MagickTrue : MagickFalse;
-              if (status != MagickFalse)
-                {
-                  channel_statistics[BlackChannel].depth++;
-                  continue;
-                }
-            }
-        }
-      if ((double) GetPixelRed(image,p) < channel_statistics[RedChannel].minima)
-        channel_statistics[RedChannel].minima=(double) GetPixelRed(image,p);
-      if ((double) GetPixelRed(image,p) > channel_statistics[RedChannel].maxima)
-        channel_statistics[RedChannel].maxima=(double) GetPixelRed(image,p);
-      channel_statistics[RedChannel].sum+=GetPixelRed(image,p);
-      channel_statistics[RedChannel].sum_squared+=(double)
-        GetPixelRed(image,p)*GetPixelRed(image,p);
-      channel_statistics[RedChannel].sum_cubed+=(double)
-        GetPixelRed(image,p)*GetPixelRed(image,p)*GetPixelRed(image,p);
-      channel_statistics[RedChannel].sum_fourth_power+=(double)
-        GetPixelRed(image,p)*GetPixelRed(image,p)*GetPixelRed(image,p)*
-        GetPixelRed(image,p);
-      if ((double) GetPixelGreen(image,p) < channel_statistics[GreenChannel].minima)
-        channel_statistics[GreenChannel].minima=(double)
-          GetPixelGreen(image,p);
-      if ((double) GetPixelGreen(image,p) > channel_statistics[GreenChannel].maxima)
-        channel_statistics[GreenChannel].maxima=(double) GetPixelGreen(image,p);
-      channel_statistics[GreenChannel].sum+=GetPixelGreen(image,p);
-      channel_statistics[GreenChannel].sum_squared+=(double)
-        GetPixelGreen(image,p)*GetPixelGreen(image,p);
-      channel_statistics[GreenChannel].sum_cubed+=(double)
-        GetPixelGreen(image,p)*GetPixelGreen(image,p)*GetPixelGreen(image,p);
-      channel_statistics[GreenChannel].sum_fourth_power+=(double)
-        GetPixelGreen(image,p)*GetPixelGreen(image,p)*GetPixelGreen(image,p)*
-        GetPixelGreen(image,p);
-      if ((double) GetPixelBlue(image,p) < channel_statistics[BlueChannel].minima)
-        channel_statistics[BlueChannel].minima=(double) GetPixelBlue(image,p);
-      if ((double) GetPixelBlue(image,p) > channel_statistics[BlueChannel].maxima)
-        channel_statistics[BlueChannel].maxima=(double) GetPixelBlue(image,p);
-      channel_statistics[BlueChannel].sum+=GetPixelBlue(image,p);
-      channel_statistics[BlueChannel].sum_squared+=(double)
-        GetPixelBlue(image,p)*GetPixelBlue(image,p);
-      channel_statistics[BlueChannel].sum_cubed+=(double)
-        GetPixelBlue(image,p)*GetPixelBlue(image,p)*GetPixelBlue(image,p);
-      channel_statistics[BlueChannel].sum_fourth_power+=(double)
-        GetPixelBlue(image,p)*GetPixelBlue(image,p)*GetPixelBlue(image,p)*
-        GetPixelBlue(image,p);
-      if (image->colorspace == CMYKColorspace)
-        {
-          if ((double) GetPixelBlack(image,p) < channel_statistics[BlackChannel].minima)
-            channel_statistics[BlackChannel].minima=(double)
-              GetPixelBlack(image,p);
-          if ((double) GetPixelBlack(image,p) > channel_statistics[BlackChannel].maxima)
-            channel_statistics[BlackChannel].maxima=(double)
-              GetPixelBlack(image,p);
-          channel_statistics[BlackChannel].sum+=GetPixelBlack(image,p);
-          channel_statistics[BlackChannel].sum_squared+=(double)
-            GetPixelBlack(image,p)*GetPixelBlack(image,p);
-          channel_statistics[BlackChannel].sum_cubed+=(double)
-            GetPixelBlack(image,p)*GetPixelBlack(image,p)*
-            GetPixelBlack(image,p);
-          channel_statistics[BlackChannel].sum_fourth_power+=(double)
-            GetPixelBlack(image,p)*GetPixelBlack(image,p)*
-            GetPixelBlack(image,p)*GetPixelBlack(image,p);
-        }
-      if (image->matte != MagickFalse)
-        {
-          if ((double) GetPixelAlpha(image,p) < channel_statistics[AlphaChannel].minima)
-            channel_statistics[AlphaChannel].minima=(double)
-              GetPixelAlpha(image,p);
-          if ((double) GetPixelAlpha(image,p) > channel_statistics[AlphaChannel].maxima)
-            channel_statistics[AlphaChannel].maxima=(double)
-              GetPixelAlpha(image,p);
-          channel_statistics[AlphaChannel].sum+=GetPixelAlpha(image,p);
-          channel_statistics[AlphaChannel].sum_squared+=(double)
-            GetPixelAlpha(image,p)*GetPixelAlpha(image,p);
-          channel_statistics[AlphaChannel].sum_cubed+=(double)
-            GetPixelAlpha(image,p)*GetPixelAlpha(image,p)*
-            GetPixelAlpha(image,p);
-          channel_statistics[AlphaChannel].sum_fourth_power+=(double)
-            GetPixelAlpha(image,p)*GetPixelAlpha(image,p)*
-            GetPixelAlpha(image,p)*GetPixelAlpha(image,p);
-        }
-      x++;
+        if ((double) p[i] < channel_statistics[channel].minima)
+          channel_statistics[channel].minima=(double) p[i];
+        if ((double) p[i] > channel_statistics[channel].maxima)
+          channel_statistics[channel].maxima=(double) p[i];
+        channel_statistics[channel].sum+=p[i];
+        channel_statistics[channel].sum_squared+=(double) p[i]*p[i];
+        channel_statistics[channel].sum_cubed+=(double) p[i]*p[i]*p[i];
+        channel_statistics[channel].sum_fourth_power+=(double) p[i]*p[i]*p[i]*
+          p[i];
+      }
       p+=GetPixelChannels(image);
     }
   }
   area=(double) image->columns*image->rows;
-  for (i=0; i < (ssize_t) CompositeChannels; i++)
+  for (i=0; i < (ssize_t) MaxPixelChannels; i++)
   {
     channel_statistics[i].sum/=area;
     channel_statistics[i].sum_squared/=area;
@@ -1762,60 +1667,56 @@ MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
       channel_statistics[i].variance-(channel_statistics[i].mean*
       channel_statistics[i].mean));
   }
-  for (i=0; i < (ssize_t) CompositeChannels; i++)
+  for (i=0; i < (ssize_t) MaxPixelChannels; i++)
   {
-    channel_statistics[CompositeChannels].depth=(size_t) MagickMax((double)
-      channel_statistics[CompositeChannels].depth,(double)
+    channel_statistics[CompositePixelChannel].depth=(size_t) EvaluateMax(
+      (double) channel_statistics[CompositePixelChannel].depth,(double)
       channel_statistics[i].depth);
-    channel_statistics[CompositeChannels].minima=MagickMin(
-      channel_statistics[CompositeChannels].minima,
+    channel_statistics[CompositePixelChannel].minima=MagickMin(
+      channel_statistics[CompositePixelChannel].minima,
       channel_statistics[i].minima);
-    channel_statistics[CompositeChannels].maxima=MagickMax(
-      channel_statistics[CompositeChannels].maxima,
+    channel_statistics[CompositePixelChannel].maxima=EvaluateMax(
+      channel_statistics[CompositePixelChannel].maxima,
       channel_statistics[i].maxima);
-    channel_statistics[CompositeChannels].sum+=channel_statistics[i].sum;
-    channel_statistics[CompositeChannels].sum_squared+=
+    channel_statistics[CompositePixelChannel].sum+=channel_statistics[i].sum;
+    channel_statistics[CompositePixelChannel].sum_squared+=
       channel_statistics[i].sum_squared;
-    channel_statistics[CompositeChannels].sum_cubed+=
+    channel_statistics[CompositePixelChannel].sum_cubed+=
       channel_statistics[i].sum_cubed;
-    channel_statistics[CompositeChannels].sum_fourth_power+=
+    channel_statistics[CompositePixelChannel].sum_fourth_power+=
       channel_statistics[i].sum_fourth_power;
-    channel_statistics[CompositeChannels].mean+=channel_statistics[i].mean;
-    channel_statistics[CompositeChannels].variance+=
+    channel_statistics[CompositePixelChannel].mean+=channel_statistics[i].mean;
+    channel_statistics[CompositePixelChannel].variance+=
       channel_statistics[i].variance-channel_statistics[i].mean*
       channel_statistics[i].mean;
-    channel_statistics[CompositeChannels].standard_deviation+=
+    channel_statistics[CompositePixelChannel].standard_deviation+=
       channel_statistics[i].variance-channel_statistics[i].mean*
       channel_statistics[i].mean;
   }
-  channels=3;
-  if (image->matte != MagickFalse)
-    channels++;
-  if (image->colorspace == CMYKColorspace)
-    channels++;
-  channel_statistics[CompositeChannels].sum/=channels;
-  channel_statistics[CompositeChannels].sum_squared/=channels;
-  channel_statistics[CompositeChannels].sum_cubed/=channels;
-  channel_statistics[CompositeChannels].sum_fourth_power/=channels;
-  channel_statistics[CompositeChannels].mean/=channels;
-  channel_statistics[CompositeChannels].variance/=channels;
-  channel_statistics[CompositeChannels].standard_deviation=
-    sqrt(channel_statistics[CompositeChannels].standard_deviation/channels);
-  channel_statistics[CompositeChannels].kurtosis/=channels;
-  channel_statistics[CompositeChannels].skewness/=channels;
-  for (i=0; i <= (ssize_t) CompositeChannels; i++)
+  channels=GetImageChannels(image);
+  channel_statistics[CompositePixelChannel].sum/=channels;
+  channel_statistics[CompositePixelChannel].sum_squared/=channels;
+  channel_statistics[CompositePixelChannel].sum_cubed/=channels;
+  channel_statistics[CompositePixelChannel].sum_fourth_power/=channels;
+  channel_statistics[CompositePixelChannel].mean/=channels;
+  channel_statistics[CompositePixelChannel].variance/=channels;
+  channel_statistics[CompositePixelChannel].standard_deviation=
+    sqrt(channel_statistics[CompositePixelChannel].standard_deviation/channels);
+  channel_statistics[CompositePixelChannel].kurtosis/=channels;
+  channel_statistics[CompositePixelChannel].skewness/=channels;
+  for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
   {
     if (channel_statistics[i].standard_deviation == 0.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].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].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*
@@ -1825,3 +1726,694 @@ MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
   }
   return(channel_statistics);
 }
+\f
+/*
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                             %
+%                                                                             %
+%                                                                             %
+%     S t a t i s t i c I m a g e                                             %
+%                                                                             %
+%                                                                             %
+%                                                                             %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  StatisticImage() makes each pixel the min / max / median / mode / etc. of
+%  the neighborhood of the specified width and height.
+%
+%  The format of the StatisticImage method is:
+%
+%      Image *StatisticImage(const Image *image,const StatisticType type,
+%        const size_t width,const size_t height,ExceptionInfo *exception)
+%
+%  A description of each parameter follows:
+%
+%    o image: the image.
+%
+%    o type: the statistic type (median, mode, etc.).
+%
+%    o width: the width of the pixel neighborhood.
+%
+%    o height: the height of the pixel neighborhood.
+%
+%    o exception: return any errors or warnings in this structure.
+%
+*/
+
+typedef struct _SkipNode
+{
+  size_t
+    next[9],
+    count,
+    signature;
+} SkipNode;
+
+typedef struct _SkipList
+{
+  ssize_t
+    level;
+
+  SkipNode
+    *nodes;
+} SkipList;
+
+typedef struct _PixelList
+{
+  size_t
+    length,
+    seed;
+
+  SkipList
+    skip_list;
+
+  size_t
+    signature;
+} PixelList;
+
+static PixelList *DestroyPixelList(PixelList *pixel_list)
+{
+  if (pixel_list == (PixelList *) NULL)
+    return((PixelList *) NULL);
+  if (pixel_list->skip_list.nodes != (SkipNode *) NULL)
+    pixel_list->skip_list.nodes=(SkipNode *) RelinquishMagickMemory(
+      pixel_list->skip_list.nodes);
+  pixel_list=(PixelList *) RelinquishMagickMemory(pixel_list);
+  return(pixel_list);
+}
+
+static PixelList **DestroyPixelListThreadSet(PixelList **pixel_list)
+{
+  register ssize_t
+    i;
+
+  assert(pixel_list != (PixelList **) NULL);
+  for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
+    if (pixel_list[i] != (PixelList *) NULL)
+      pixel_list[i]=DestroyPixelList(pixel_list[i]);
+  pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list);
+  return(pixel_list);
+}
+
+static PixelList *AcquirePixelList(const size_t width,const size_t height)
+{
+  PixelList
+    *pixel_list;
+
+  pixel_list=(PixelList *) AcquireMagickMemory(sizeof(*pixel_list));
+  if (pixel_list == (PixelList *) NULL)
+    return(pixel_list);
+  (void) ResetMagickMemory((void *) pixel_list,0,sizeof(*pixel_list));
+  pixel_list->length=width*height;
+  pixel_list->skip_list.nodes=(SkipNode *) AcquireQuantumMemory(65537UL,
+    sizeof(*pixel_list->skip_list.nodes));
+  if (pixel_list->skip_list.nodes == (SkipNode *) NULL)
+    return(DestroyPixelList(pixel_list));
+  (void) ResetMagickMemory(pixel_list->skip_list.nodes,0,65537UL*
+    sizeof(*pixel_list->skip_list.nodes));
+  pixel_list->signature=MagickSignature;
+  return(pixel_list);
+}
+
+static PixelList **AcquirePixelListThreadSet(const size_t width,
+  const size_t height)
+{
+  PixelList
+    **pixel_list;
+
+  register ssize_t
+    i;
+
+  size_t
+    number_threads;
+
+  number_threads=GetOpenMPMaximumThreads();
+  pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
+    sizeof(*pixel_list));
+  if (pixel_list == (PixelList **) NULL)
+    return((PixelList **) NULL);
+  (void) ResetMagickMemory(pixel_list,0,number_threads*sizeof(*pixel_list));
+  for (i=0; i < (ssize_t) number_threads; i++)
+  {
+    pixel_list[i]=AcquirePixelList(width,height);
+    if (pixel_list[i] == (PixelList *) NULL)
+      return(DestroyPixelListThreadSet(pixel_list));
+  }
+  return(pixel_list);
+}
+
+static void AddNodePixelList(PixelList *pixel_list,const size_t color)
+{
+  register SkipList
+    *p;
+
+  register ssize_t
+    level;
+
+  size_t
+    search,
+    update[9];
+
+  /*
+    Initialize the node.
+  */
+  p=(&pixel_list->skip_list);
+  p->nodes[color].signature=pixel_list->signature;
+  p->nodes[color].count=1;
+  /*
+    Determine where it belongs in the list.
+  */
+  search=65536UL;
+  for (level=p->level; level >= 0; level--)
+  {
+    while (p->nodes[search].next[level] < color)
+      search=p->nodes[search].next[level];
+    update[level]=search;
+  }
+  /*
+    Generate a pseudo-random level for this node.
+  */
+  for (level=0; ; level++)
+  {
+    pixel_list->seed=(pixel_list->seed*42893621L)+1L;
+    if ((pixel_list->seed & 0x300) != 0x300)
+      break;
+  }
+  if (level > 8)
+    level=8;
+  if (level > (p->level+2))
+    level=p->level+2;
+  /*
+    If we're raising the list's level, link back to the root node.
+  */
+  while (level > p->level)
+  {
+    p->level++;
+    update[p->level]=65536UL;
+  }
+  /*
+    Link the node into the skip-list.
+  */
+  do
+  {
+    p->nodes[color].next[level]=p->nodes[update[level]].next[level];
+    p->nodes[update[level]].next[level]=color;
+  } while (level-- > 0);
+}
+
+static inline void GetMaximumPixelList(PixelList *pixel_list,Quantum *pixel)
+{
+  register SkipList
+    *p;
+
+  size_t
+    color,
+    maximum;
+
+  ssize_t
+    count;
+
+  /*
+    Find the maximum value for each of the color.
+  */
+  p=(&pixel_list->skip_list);
+  color=65536L;
+  count=0;
+  maximum=p->nodes[color].next[0];
+  do
+  {
+    color=p->nodes[color].next[0];
+    if (color > maximum)
+      maximum=color;
+    count+=p->nodes[color].count;
+  } while (count < (ssize_t) pixel_list->length);
+  *pixel=ScaleShortToQuantum((unsigned short) maximum);
+}
+
+static inline void GetMeanPixelList(PixelList *pixel_list,Quantum *pixel)
+{
+  MagickRealType
+    sum;
+
+  register SkipList
+    *p;
+
+  size_t
+    color;
+
+  ssize_t
+    count;
+
+  /*
+    Find the mean value for each of the color.
+  */
+  p=(&pixel_list->skip_list);
+  color=65536L;
+  count=0;
+  sum=0.0;
+  do
+  {
+    color=p->nodes[color].next[0];
+    sum+=(MagickRealType) p->nodes[color].count*color;
+    count+=p->nodes[color].count;
+  } while (count < (ssize_t) pixel_list->length);
+  sum/=pixel_list->length;
+  *pixel=ScaleShortToQuantum((unsigned short) sum);
+}
+
+static inline void GetMedianPixelList(PixelList *pixel_list,Quantum *pixel)
+{
+  register SkipList
+    *p;
+
+  size_t
+    color;
+
+  ssize_t
+    count;
+
+  /*
+    Find the median value for each of the color.
+  */
+  p=(&pixel_list->skip_list);
+  color=65536L;
+  count=0;
+  do
+  {
+    color=p->nodes[color].next[0];
+    count+=p->nodes[color].count;
+  } while (count <= (ssize_t) (pixel_list->length >> 1));
+  *pixel=ScaleShortToQuantum((unsigned short) color);
+}
+
+static inline void GetMinimumPixelList(PixelList *pixel_list,Quantum *pixel)
+{
+  register SkipList
+    *p;
+
+  size_t
+    color,
+    minimum;
+
+  ssize_t
+    count;
+
+  /*
+    Find the minimum value for each of the color.
+  */
+  p=(&pixel_list->skip_list);
+  count=0;
+  color=65536UL;
+  minimum=p->nodes[color].next[0];
+  do
+  {
+    color=p->nodes[color].next[0];
+    if (color < minimum)
+      minimum=color;
+    count+=p->nodes[color].count;
+  } while (count < (ssize_t) pixel_list->length);
+  *pixel=ScaleShortToQuantum((unsigned short) minimum);
+}
+
+static inline void GetModePixelList(PixelList *pixel_list,Quantum *pixel)
+{
+  register SkipList
+    *p;
+
+  size_t
+    color,
+    max_count,
+    mode;
+
+  ssize_t
+    count;
+
+  /*
+    Make each pixel the 'predominant color' of the specified neighborhood.
+  */
+  p=(&pixel_list->skip_list);
+  color=65536L;
+  mode=color;
+  max_count=p->nodes[mode].count;
+  count=0;
+  do
+  {
+    color=p->nodes[color].next[0];
+    if (p->nodes[color].count > max_count)
+      {
+        mode=color;
+        max_count=p->nodes[mode].count;
+      }
+    count+=p->nodes[color].count;
+  } while (count < (ssize_t) pixel_list->length);
+  *pixel=ScaleShortToQuantum((unsigned short) mode);
+}
+
+static inline void GetNonpeakPixelList(PixelList *pixel_list,Quantum *pixel)
+{
+  register SkipList
+    *p;
+
+  size_t
+    color,
+    next,
+    previous;
+
+  ssize_t
+    count;
+
+  /*
+    Finds the non peak value for each of the colors.
+  */
+  p=(&pixel_list->skip_list);
+  color=65536L;
+  next=p->nodes[color].next[0];
+  count=0;
+  do
+  {
+    previous=color;
+    color=next;
+    next=p->nodes[color].next[0];
+    count+=p->nodes[color].count;
+  } while (count <= (ssize_t) (pixel_list->length >> 1));
+  if ((previous == 65536UL) && (next != 65536UL))
+    color=next;
+  else
+    if ((previous != 65536UL) && (next == 65536UL))
+      color=previous;
+  *pixel=ScaleShortToQuantum((unsigned short) color);
+}
+
+static inline void GetStandardDeviationPixelList(PixelList *pixel_list,
+  Quantum *pixel)
+{
+  MagickRealType
+    sum,
+    sum_squared;
+
+  register SkipList
+    *p;
+
+  size_t
+    color;
+
+  ssize_t
+    count;
+
+  /*
+    Find the standard-deviation value for each of the color.
+  */
+  p=(&pixel_list->skip_list);
+  color=65536L;
+  count=0;
+  sum=0.0;
+  sum_squared=0.0;
+  do
+  {
+    register ssize_t
+      i;
+
+    color=p->nodes[color].next[0];
+    sum+=(MagickRealType) p->nodes[color].count*color;
+    for (i=0; i < (ssize_t) p->nodes[color].count; i++)
+      sum_squared+=((MagickRealType) color)*((MagickRealType) color);
+    count+=p->nodes[color].count;
+  } while (count < (ssize_t) pixel_list->length);
+  sum/=pixel_list->length;
+  sum_squared/=pixel_list->length;
+  *pixel=ScaleShortToQuantum((unsigned short) sqrt(sum_squared-(sum*sum)));
+}
+
+static inline void InsertPixelList(const Image *image,const Quantum pixel,
+  PixelList *pixel_list)
+{
+  size_t
+    signature;
+
+  unsigned short
+    index;
+
+  index=ScaleQuantumToShort(pixel);
+  signature=pixel_list->skip_list.nodes[index].signature;
+  if (signature == pixel_list->signature)
+    {
+      pixel_list->skip_list.nodes[index].count++;
+      return;
+    }
+  AddNodePixelList(pixel_list,index);
+}
+
+static inline MagickRealType MagickAbsoluteValue(const MagickRealType x)
+{
+  if (x < 0)
+    return(-x);
+  return(x);
+}
+
+static inline size_t MagickMax(const size_t x,const size_t y)
+{
+  if (x > y)
+    return(x);
+  return(y);
+}
+
+static void ResetPixelList(PixelList *pixel_list)
+{
+  int
+    level;
+
+  register SkipNode
+    *root;
+
+  register SkipList
+    *p;
+
+  /*
+    Reset the skip-list.
+  */
+  p=(&pixel_list->skip_list);
+  root=p->nodes+65536UL;
+  p->level=0;
+  for (level=0; level < 9; level++)
+    root->next[level]=65536UL;
+  pixel_list->seed=pixel_list->signature++;
+}
+
+MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
+  const size_t width,const size_t height,ExceptionInfo *exception)
+{
+#define StatisticImageTag  "Statistic/Image"
+
+  CacheView
+    *image_view,
+    *statistic_view;
+
+  Image
+    *statistic_image;
+
+  MagickBooleanType
+    status;
+
+  MagickOffsetType
+    progress;
+
+  PixelList
+    **restrict pixel_list;
+
+  ssize_t
+    center,
+    y;
+
+  /*
+    Initialize statistics image attributes.
+  */
+  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);
+  statistic_image=CloneImage(image,image->columns,image->rows,MagickTrue,
+    exception);
+  if (statistic_image == (Image *) NULL)
+    return((Image *) NULL);
+  status=SetImageStorageClass(statistic_image,DirectClass,exception);
+  if (status == MagickFalse)
+    {
+      statistic_image=DestroyImage(statistic_image);
+      return((Image *) NULL);
+    }
+  pixel_list=AcquirePixelListThreadSet(MagickMax(width,1),MagickMax(height,1));
+  if (pixel_list == (PixelList **) NULL)
+    {
+      statistic_image=DestroyImage(statistic_image);
+      ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
+    }
+  /*
+    Make each pixel the min / max / median / mode / etc. of the neighborhood.
+  */
+  center=(ssize_t) GetPixelChannels(image)*(image->columns+MagickMax(width,1))*
+    (MagickMax(height,1)/2L)+GetPixelChannels(image)*(MagickMax(width,1)/2L);
+  status=MagickTrue;
+  progress=0;
+  image_view=AcquireVirtualCacheView(image,exception);
+  statistic_view=AcquireAuthenticCacheView(statistic_image,exception);
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+  #pragma omp parallel for schedule(static,4) shared(progress,status)
+#endif
+  for (y=0; y < (ssize_t) statistic_image->rows; y++)
+  {
+    const int
+      id = GetOpenMPThreadId();
+
+    register const Quantum
+      *restrict p;
+
+    register Quantum
+      *restrict q;
+
+    register ssize_t
+      x;
+
+    if (status == MagickFalse)
+      continue;
+    p=GetCacheViewVirtualPixels(image_view,-((ssize_t) MagickMax(width,1)/2L),y-
+      (ssize_t) (MagickMax(height,1)/2L),image->columns+MagickMax(width,1),
+      MagickMax(height,1),exception);
+    q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns,      1,exception);
+    if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
+      {
+        status=MagickFalse;
+        continue;
+      }
+    for (x=0; x < (ssize_t) statistic_image->columns; x++)
+    {
+      register ssize_t
+        i;
+
+      if (GetPixelMask(image,p) != 0)
+        {
+          p+=GetPixelChannels(image);
+          q+=GetPixelChannels(statistic_image);
+          continue;
+        }
+      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
+      {
+        PixelChannel
+          channel;
+
+        PixelTrait
+          statistic_traits,
+          traits;
+
+        Quantum
+          pixel;
+
+        register const Quantum
+          *restrict pixels;
+
+        register ssize_t
+          u;
+
+        ssize_t
+          v;
+
+        channel=GetPixelChannelMapChannel(image,i);
+        traits=GetPixelChannelMapTraits(image,channel);
+        statistic_traits=GetPixelChannelMapTraits(statistic_image,channel);
+        if ((traits == UndefinedPixelTrait) ||
+            (statistic_traits == UndefinedPixelTrait))
+          continue;
+        if ((statistic_traits & CopyPixelTrait) != 0)
+          {
+            SetPixelChannel(statistic_image,channel,p[center+i],q);
+            continue;
+          }
+        pixels=p;
+        ResetPixelList(pixel_list[id]);
+        for (v=0; v < (ssize_t) MagickMax(height,1); v++)
+        {
+          for (u=0; u < (ssize_t) MagickMax(width,1); u++)
+          {
+            InsertPixelList(image,pixels[i],pixel_list[id]);
+            pixels+=GetPixelChannels(image);
+          }
+          pixels+=image->columns*GetPixelChannels(image);
+        }
+        switch (type)
+        {
+          case GradientStatistic:
+          {
+            MagickRealType
+              maximum,
+              minimum;
+
+            GetMinimumPixelList(pixel_list[id],&pixel);
+            minimum=(MagickRealType) pixel;
+            GetMaximumPixelList(pixel_list[id],&pixel);
+            maximum=(MagickRealType) pixel;
+            pixel=ClampToQuantum(MagickAbsoluteValue(maximum-minimum));
+            break;
+          }
+          case MaximumStatistic:
+          {
+            GetMaximumPixelList(pixel_list[id],&pixel);
+            break;
+          }
+          case MeanStatistic:
+          {
+            GetMeanPixelList(pixel_list[id],&pixel);
+            break;
+          }
+          case MedianStatistic:
+          default:
+          {
+            GetMedianPixelList(pixel_list[id],&pixel);
+            break;
+          }
+          case MinimumStatistic:
+          {
+            GetMinimumPixelList(pixel_list[id],&pixel);
+            break;
+          }
+          case ModeStatistic:
+          {
+            GetModePixelList(pixel_list[id],&pixel);
+            break;
+          }
+          case NonpeakStatistic:
+          {
+            GetNonpeakPixelList(pixel_list[id],&pixel);
+            break;
+          }
+          case StandardDeviationStatistic:
+          {
+            GetStandardDeviationPixelList(pixel_list[id],&pixel);
+            break;
+          }
+        }
+        SetPixelChannel(statistic_image,channel,pixel,q);
+      }
+      p+=GetPixelChannels(image);
+      q+=GetPixelChannels(statistic_image);
+    }
+    if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
+      status=MagickFalse;
+    if (image->progress_monitor != (MagickProgressMonitor) NULL)
+      {
+        MagickBooleanType
+          proceed;
+
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+  #pragma omp critical (MagickCore_StatisticImage)
+#endif
+        proceed=SetImageProgress(image,StatisticImageTag,progress++,
+          image->rows);
+        if (proceed == MagickFalse)
+          status=MagickFalse;
+      }
+  }
+  statistic_view=DestroyCacheView(statistic_view);
+  image_view=DestroyCacheView(image_view);
+  pixel_list=DestroyPixelListThreadSet(pixel_list);
+  return(statistic_image);
+}