granicus.if.org Git - imagemagick/blob - MagickCore/statistic.c

   1 /*
   2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   3 %                                                                             %
   4 %                                                                             %
   5 %                                                                             %
   6 %        SSSSS  TTTTT   AAA   TTTTT  IIIII  SSSSS  TTTTT  IIIII   CCCC        %
   7 %        SS       T    A   A    T      I    SS       T      I    C            %
   8 %         SSS     T    AAAAA    T      I     SSS     T      I    C            %
   9 %           SS    T    A   A    T      I       SS    T      I    C            %
  10 %        SSSSS    T    A   A    T    IIIII  SSSSS    T    IIIII   CCCC        %
  11 %                                                                             %
  12 %                                                                             %
  13 %                     MagickCore Image Statistical Methods                    %
  14 %                                                                             %
  15 %                              Software Design                                %
  16 %                                John Cristy                                  %
  17 %                                 July 1992                                   %
  18 %                                                                             %
  19 %                                                                             %
  20 %  Copyright 1999-2012 ImageMagick Studio LLC, a non-profit organization      %
  21 %  dedicated to making software imaging solutions freely available.           %
  22 %                                                                             %
  23 %  You may not use this file except in compliance with the License.  You may  %
  24 %  obtain a copy of the License at                                            %
  25 %                                                                             %
  26 %    http://www.imagemagick.org/script/license.php                            %
  27 %                                                                             %
  28 %  Unless required by applicable law or agreed to in writing, software        %
  29 %  distributed under the License is distributed on an "AS IS" BASIS,          %
  30 %  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
  31 %  See the License for the specific language governing permissions and        %
  32 %  limitations under the License.                                             %
  33 %                                                                             %
  34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  35 %
  36 %
  37 %
  38 */
  39 \f
  40 /*
  41   Include declarations.
  42 */
  43 #include "MagickCore/studio.h"
  44 #include "MagickCore/property.h"
  45 #include "MagickCore/animate.h"
  46 #include "MagickCore/blob.h"
  47 #include "MagickCore/blob-private.h"
  48 #include "MagickCore/cache.h"
  49 #include "MagickCore/cache-private.h"
  50 #include "MagickCore/cache-view.h"
  51 #include "MagickCore/client.h"
  52 #include "MagickCore/color.h"
  53 #include "MagickCore/color-private.h"
  54 #include "MagickCore/colorspace.h"
  55 #include "MagickCore/colorspace-private.h"
  56 #include "MagickCore/composite.h"
  57 #include "MagickCore/composite-private.h"
  58 #include "MagickCore/compress.h"
  59 #include "MagickCore/constitute.h"
  60 #include "MagickCore/display.h"
  61 #include "MagickCore/draw.h"
  62 #include "MagickCore/enhance.h"
  63 #include "MagickCore/exception.h"
  64 #include "MagickCore/exception-private.h"
  65 #include "MagickCore/gem.h"
  66 #include "MagickCore/gem-private.h"
  67 #include "MagickCore/geometry.h"
  68 #include "MagickCore/list.h"
  69 #include "MagickCore/image-private.h"
  70 #include "MagickCore/magic.h"
  71 #include "MagickCore/magick.h"
  72 #include "MagickCore/memory_.h"
  73 #include "MagickCore/module.h"
  74 #include "MagickCore/monitor.h"
  75 #include "MagickCore/monitor-private.h"
  76 #include "MagickCore/option.h"
  77 #include "MagickCore/paint.h"
  78 #include "MagickCore/pixel-accessor.h"
  79 #include "MagickCore/profile.h"
  80 #include "MagickCore/quantize.h"
  81 #include "MagickCore/quantum-private.h"
  82 #include "MagickCore/random_.h"
  83 #include "MagickCore/random-private.h"
  84 #include "MagickCore/resource_.h"
  85 #include "MagickCore/segment.h"
  86 #include "MagickCore/semaphore.h"
  87 #include "MagickCore/signature-private.h"
  88 #include "MagickCore/statistic.h"
  89 #include "MagickCore/string_.h"
  90 #include "MagickCore/thread-private.h"
  91 #include "MagickCore/timer.h"
  92 #include "MagickCore/utility.h"
  93 #include "MagickCore/version.h"
  94 \f
  95 /*
  96 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  97 %                                                                             %
  98 %                                                                             %
  99 %                                                                             %
 100 %     E v a l u a t e I m a g e                                               %
 101 %                                                                             %
 102 %                                                                             %
 103 %                                                                             %
 104 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 105 %
 106 %  EvaluateImage() applies a value to the image with an arithmetic, relational,
 107 %  or logical operator to an image. Use these operations to lighten or darken
 108 %  an image, to increase or decrease contrast in an image, or to produce the
 109 %  "negative" of an image.
 110 %
 111 %  The format of the EvaluateImage method is:
 112 %
 113 %      MagickBooleanType EvaluateImage(Image *image,
 114 %        const MagickEvaluateOperator op,const double value,
 115 %        ExceptionInfo *exception)
 116 %      MagickBooleanType EvaluateImages(Image *images,
 117 %        const MagickEvaluateOperator op,const double value,
 118 %        ExceptionInfo *exception)
 119 %
 120 %  A description of each parameter follows:
 121 %
 122 %    o image: the image.
 123 %
 124 %    o op: A channel op.
 125 %
 126 %    o value: A value value.
 127 %
 128 %    o exception: return any errors or warnings in this structure.
 129 %
 130 */
 131
 132 typedef struct _PixelChannels
 133 {
 134   MagickRealType
 135     channel[CompositePixelChannel];
 136 } PixelChannels;
 137
 138 static PixelChannels **DestroyPixelThreadSet(PixelChannels **pixels)
 139 {
 140   register ssize_t
 141     i;
 142
 143   assert(pixels != (PixelChannels **) NULL);
 144   for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
 145     if (pixels[i] != (PixelChannels *) NULL)
 146       pixels[i]=(PixelChannels *) RelinquishMagickMemory(pixels[i]);
 147   pixels=(PixelChannels **) RelinquishMagickMemory(pixels);
 148   return(pixels);
 149 }
 150
 151 static PixelChannels **AcquirePixelThreadSet(const Image *image,
 152   const size_t number_images)
 153 {
 154   register ssize_t
 155     i;
 156
 157   PixelChannels
 158     **pixels;
 159
 160   size_t
 161     length,
 162     number_threads;
 163
 164   number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
 165   pixels=(PixelChannels **) AcquireQuantumMemory(number_threads,
 166     sizeof(*pixels));
 167   if (pixels == (PixelChannels **) NULL)
 168     return((PixelChannels **) NULL);
 169   (void) ResetMagickMemory(pixels,0,number_threads*sizeof(*pixels));
 170   for (i=0; i < (ssize_t) number_threads; i++)
 171   {
 172     register ssize_t
 173       j;
 174
 175     length=image->columns;
 176     if (length < number_images)
 177       length=number_images;
 178     pixels[i]=(PixelChannels *) AcquireQuantumMemory(length,sizeof(**pixels));
 179     if (pixels[i] == (PixelChannels *) NULL)
 180       return(DestroyPixelThreadSet(pixels));
 181     for (j=0; j < (ssize_t) length; j++)
 182     {
 183       register ssize_t
 184         k;
 185
 186       for (k=0; k < MaxPixelChannels; k++)
 187         pixels[i][j].channel[k]=0.0;
 188     }
 189   }
 190   return(pixels);
 191 }
 192
 193 static inline double EvaluateMax(const double x,const double y)
 194 {
 195   if (x > y)
 196     return(x);
 197   return(y);
 198 }
 199
 200 #if defined(__cplusplus) || defined(c_plusplus)
 201 extern "C" {
 202 #endif
 203
 204 static int IntensityCompare(const void *x,const void *y)
 205 {
 206   const PixelChannels
 207     *color_1,
 208     *color_2;
 209
 210   MagickRealType
 211     distance;
 212
 213   register ssize_t
 214     i;
 215
 216   color_1=(const PixelChannels *) x;
 217   color_2=(const PixelChannels *) y;
 218   distance=0.0;
 219   for (i=0; i < MaxPixelChannels; i++)
 220     distance+=color_1->channel[i]-(MagickRealType) color_2->channel[i];
 221   return(distance < 0 ? -1 : distance > 0 ? 1 : 0);
 222 }
 223
 224 #if defined(__cplusplus) || defined(c_plusplus)
 225 }
 226 #endif
 227
 228 static inline double MagickMin(const double x,const double y)
 229 {
 230   if (x < y)
 231     return(x);
 232   return(y);
 233 }
 234
 235 static MagickRealType ApplyEvaluateOperator(RandomInfo *random_info,
 236   Quantum pixel,const MagickEvaluateOperator op,const MagickRealType value)
 237 {
 238   MagickRealType
 239     result;
 240
 241   result=0.0;
 242   switch (op)
 243   {
 244     case UndefinedEvaluateOperator:
 245       break;
 246     case AbsEvaluateOperator:
 247     {
 248       result=(MagickRealType) fabs((double) (pixel+value));
 249       break;
 250     }
 251     case AddEvaluateOperator:
 252     {
 253       result=(MagickRealType) (pixel+value);
 254       break;
 255     }
 256     case AddModulusEvaluateOperator:
 257     {
 258       /*
 259         This returns a 'floored modulus' of the addition which is a positive
 260         result.  It differs from % or fmod() that returns a 'truncated modulus'
 261         result, where floor() is replaced by trunc() and could return a
 262         negative result (which is clipped).
 263       */
 264       result=pixel+value;
 265       result-=(QuantumRange+1.0)*floor((double) result/(QuantumRange+1.0));
 266       break;
 267     }
 268     case AndEvaluateOperator:
 269     {
 270       result=(MagickRealType) ((size_t) pixel & (size_t) (value+0.5));
 271       break;
 272     }
 273     case CosineEvaluateOperator:
 274     {
 275       result=(MagickRealType) (QuantumRange*(0.5*cos((double) (2.0*MagickPI*
 276         QuantumScale*pixel*value))+0.5));
 277       break;
 278     }
 279     case DivideEvaluateOperator:
 280     {
 281       result=pixel/(value == 0.0 ? 1.0 : value);
 282       break;
 283     }
 284     case ExponentialEvaluateOperator:
 285     {
 286       result=(MagickRealType) (QuantumRange*exp((double) (value*QuantumScale*
 287         pixel)));
 288       break;
 289     }
 290     case GaussianNoiseEvaluateOperator:
 291     {
 292       result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
 293         GaussianNoise,value);
 294       break;
 295     }
 296     case ImpulseNoiseEvaluateOperator:
 297     {
 298       result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
 299         ImpulseNoise,value);
 300       break;
 301     }
 302     case LaplacianNoiseEvaluateOperator:
 303     {
 304       result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
 305         LaplacianNoise,value);
 306       break;
 307     }
 308     case LeftShiftEvaluateOperator:
 309     {
 310       result=(MagickRealType) ((size_t) pixel << (size_t) (value+0.5));
 311       break;
 312     }
 313     case LogEvaluateOperator:
 314     {
 315       result=(MagickRealType) (QuantumRange*log((double) (QuantumScale*value*
 316         pixel+1.0))/log((double) (value+1.0)));
 317       break;
 318     }
 319     case MaxEvaluateOperator:
 320     {
 321       result=(MagickRealType) EvaluateMax((double) pixel,value);
 322       break;
 323     }
 324     case MeanEvaluateOperator:
 325     {
 326       result=(MagickRealType) (pixel+value);
 327       break;
 328     }
 329     case MedianEvaluateOperator:
 330     {
 331       result=(MagickRealType) (pixel+value);
 332       break;
 333     }
 334     case MinEvaluateOperator:
 335     {
 336       result=(MagickRealType) MagickMin((double) pixel,value);
 337       break;
 338     }
 339     case MultiplicativeNoiseEvaluateOperator:
 340     {
 341       result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
 342         MultiplicativeGaussianNoise,value);
 343       break;
 344     }
 345     case MultiplyEvaluateOperator:
 346     {
 347       result=(MagickRealType) (value*pixel);
 348       break;
 349     }
 350     case OrEvaluateOperator:
 351     {
 352       result=(MagickRealType) ((size_t) pixel | (size_t) (value+0.5));
 353       break;
 354     }
 355     case PoissonNoiseEvaluateOperator:
 356     {
 357       result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
 358         PoissonNoise,value);
 359       break;
 360     }
 361     case PowEvaluateOperator:
 362     {
 363       result=(MagickRealType) (QuantumRange*pow((double) (QuantumScale*pixel),
 364         (double) value));
 365       break;
 366     }
 367     case RightShiftEvaluateOperator:
 368     {
 369       result=(MagickRealType) ((size_t) pixel >> (size_t) (value+0.5));
 370       break;
 371     }
 372     case SetEvaluateOperator:
 373     {
 374       result=value;
 375       break;
 376     }
 377     case SineEvaluateOperator:
 378     {
 379       result=(MagickRealType) (QuantumRange*(0.5*sin((double) (2.0*MagickPI*
 380         QuantumScale*pixel*value))+0.5));
 381       break;
 382     }
 383     case SubtractEvaluateOperator:
 384     {
 385       result=(MagickRealType) (pixel-value);
 386       break;
 387     }
 388     case SumEvaluateOperator:
 389     {
 390       result=(MagickRealType) (pixel+value);
 391       break;
 392     }
 393     case ThresholdEvaluateOperator:
 394     {
 395       result=(MagickRealType) (((MagickRealType) pixel <= value) ? 0 :
 396         QuantumRange);
 397       break;
 398     }
 399     case ThresholdBlackEvaluateOperator:
 400     {
 401       result=(MagickRealType) (((MagickRealType) pixel <= value) ? 0 : pixel);
 402       break;
 403     }
 404     case ThresholdWhiteEvaluateOperator:
 405     {
 406       result=(MagickRealType) (((MagickRealType) pixel > value) ? QuantumRange :
 407         pixel);
 408       break;
 409     }
 410     case UniformNoiseEvaluateOperator:
 411     {
 412       result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
 413         UniformNoise,value);
 414       break;
 415     }
 416     case XorEvaluateOperator:
 417     {
 418       result=(MagickRealType) ((size_t) pixel ^ (size_t) (value+0.5));
 419       break;
 420     }
 421   }
 422   return(result);
 423 }
 424
 425 MagickExport Image *EvaluateImages(const Image *images,
 426   const MagickEvaluateOperator op,ExceptionInfo *exception)
 427 {
 428 #define EvaluateImageTag  "Evaluate/Image"
 429
 430   CacheView
 431     *evaluate_view;
 432
 433   const Image
 434     *next;
 435
 436   Image
 437     *evaluate_image;
 438
 439   MagickBooleanType
 440     status;
 441
 442   MagickOffsetType
 443     progress;
 444
 445   PixelChannels
 446     **restrict evaluate_pixels;
 447
 448   RandomInfo
 449     **restrict random_info;
 450
 451   size_t
 452     number_images;
 453
 454   ssize_t
 455     y;
 456
 457   unsigned long
 458     key;
 459
 460   /*
 461     Ensure the image are the same size.
 462   */
 463   assert(images != (Image *) NULL);
 464   assert(images->signature == MagickSignature);
 465   if (images->debug != MagickFalse)
 466     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
 467   assert(exception != (ExceptionInfo *) NULL);
 468   assert(exception->signature == MagickSignature);
 469   for (next=images; next != (Image *) NULL; next=GetNextImageInList(next))
 470     if ((next->columns != images->columns) || (next->rows != images->rows))
 471       {
 472         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
 473           "ImageWidthsOrHeightsDiffer","'%s'",images->filename);
 474         return((Image *) NULL);
 475       }
 476   /*
 477     Initialize evaluate next attributes.
 478   */
 479   evaluate_image=CloneImage(images,images->columns,images->rows,MagickTrue,
 480     exception);
 481   if (evaluate_image == (Image *) NULL)
 482     return((Image *) NULL);
 483   if (SetImageStorageClass(evaluate_image,DirectClass,exception) == MagickFalse)
 484     {
 485       evaluate_image=DestroyImage(evaluate_image);
 486       return((Image *) NULL);
 487     }
 488   number_images=GetImageListLength(images);
 489   evaluate_pixels=AcquirePixelThreadSet(images,number_images);
 490   if (evaluate_pixels == (PixelChannels **) NULL)
 491     {
 492       evaluate_image=DestroyImage(evaluate_image);
 493       (void) ThrowMagickException(exception,GetMagickModule(),
 494         ResourceLimitError,"MemoryAllocationFailed","'%s'",images->filename);
 495       return((Image *) NULL);
 496     }
 497   /*
 498     Evaluate image pixels.
 499   */
 500   status=MagickTrue;
 501   progress=0;
 502   random_info=AcquireRandomInfoThreadSet();
 503   key=GetRandomSecretKey(random_info[0]);
 504   evaluate_view=AcquireAuthenticCacheView(evaluate_image,exception);
 505   if (op == MedianEvaluateOperator)
 506     {
 507 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 508       #pragma omp parallel for schedule(static) shared(progress,status) \
 509         if (((evaluate_image->rows*evaluate_image->columns) > 8192) && (key == ~0UL)) \
 510           num_threads(GetMagickResourceLimit(ThreadResource))
 511 #endif
 512       for (y=0; y < (ssize_t) evaluate_image->rows; y++)
 513       {
 514         CacheView
 515           *image_view;
 516
 517         const Image
 518           *next;
 519
 520         const int
 521           id = GetOpenMPThreadId();
 522
 523         register PixelChannels
 524           *evaluate_pixel;
 525
 526         register Quantum
 527           *restrict q;
 528
 529         register ssize_t
 530           x;
 531
 532         if (status == MagickFalse)
 533           continue;
 534         q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,
 535           evaluate_image->columns,1,exception);
 536         if (q == (Quantum *) NULL)
 537           {
 538             status=MagickFalse;
 539             continue;
 540           }
 541         evaluate_pixel=evaluate_pixels[id];
 542         for (x=0; x < (ssize_t) evaluate_image->columns; x++)
 543         {
 544           register ssize_t
 545             j,
 546             k;
 547
 548           for (j=0; j < (ssize_t) number_images; j++)
 549             for (k=0; k < MaxPixelChannels; k++)
 550               evaluate_pixel[j].channel[k]=0.0;
 551           next=images;
 552           for (j=0; j < (ssize_t) number_images; j++)
 553           {
 554             register const Quantum
 555               *p;
 556
 557             register ssize_t
 558               i;
 559
 560             image_view=AcquireVirtualCacheView(next,exception);
 561             p=GetCacheViewVirtualPixels(image_view,x,y,1,1,exception);
 562             if (p == (const Quantum *) NULL)
 563               {
 564                 image_view=DestroyCacheView(image_view);
 565                 break;
 566               }
 567             for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++)
 568             {
 569               PixelChannel
 570                 channel;
 571
 572               PixelTrait
 573                 evaluate_traits,
 574                 traits;
 575
 576               channel=GetPixelChannelMapChannel(evaluate_image,i);
 577               evaluate_traits=GetPixelChannelMapTraits(evaluate_image,channel);
 578               traits=GetPixelChannelMapTraits(next,channel);
 579               if ((traits == UndefinedPixelTrait) ||
 580                   (evaluate_traits == UndefinedPixelTrait))
 581                 continue;
 582               if ((evaluate_traits & UpdatePixelTrait) == 0)
 583                 continue;
 584               evaluate_pixel[j].channel[i]=ApplyEvaluateOperator(
 585                 random_info[id],GetPixelChannel(evaluate_image,channel,p),op,
 586                 evaluate_pixel[j].channel[i]);
 587             }
 588             image_view=DestroyCacheView(image_view);
 589             next=GetNextImageInList(next);
 590           }
 591           qsort((void *) evaluate_pixel,number_images,sizeof(*evaluate_pixel),
 592             IntensityCompare);
 593           for (k=0; k < (ssize_t) GetPixelChannels(evaluate_image); k++)
 594             q[k]=ClampToQuantum(evaluate_pixel[j/2].channel[k]);
 595           q+=GetPixelChannels(evaluate_image);
 596         }
 597         if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
 598           status=MagickFalse;
 599         if (images->progress_monitor != (MagickProgressMonitor) NULL)
 600           {
 601             MagickBooleanType
 602               proceed;
 603
 604 #if   defined(MAGICKCORE_OPENMP_SUPPORT)
 605             #pragma omp critical (MagickCore_EvaluateImages)
 606 #endif
 607             proceed=SetImageProgress(images,EvaluateImageTag,progress++,
 608               evaluate_image->rows);
 609             if (proceed == MagickFalse)
 610               status=MagickFalse;
 611           }
 612       }
 613     }
 614   else
 615     {
 616 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 617       #pragma omp parallel for schedule(static) shared(progress,status) \
 618         if (((evaluate_image->rows*evaluate_image->columns) > 8192) && (key == ~0UL)) \
 619           num_threads(GetMagickResourceLimit(ThreadResource))
 620 #endif
 621       for (y=0; y < (ssize_t) evaluate_image->rows; y++)
 622       {
 623         CacheView
 624           *image_view;
 625
 626         const Image
 627           *next;
 628
 629         const int
 630           id = GetOpenMPThreadId();
 631
 632         register ssize_t
 633           i,
 634           x;
 635
 636         register PixelChannels
 637           *evaluate_pixel;
 638
 639         register Quantum
 640           *restrict q;
 641
 642         ssize_t
 643           j;
 644
 645         if (status == MagickFalse)
 646           continue;
 647         q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,
 648           evaluate_image->columns,1,exception);
 649         if (q == (Quantum *) NULL)
 650           {
 651             status=MagickFalse;
 652             continue;
 653           }
 654         evaluate_pixel=evaluate_pixels[id];
 655         for (j=0; j < (ssize_t) evaluate_image->columns; j++)
 656           for (i=0; i < MaxPixelChannels; i++)
 657             evaluate_pixel[j].channel[i]=0.0;
 658         next=images;
 659         for (j=0; j < (ssize_t) number_images; j++)
 660         {
 661           register const Quantum
 662             *p;
 663
 664           image_view=AcquireVirtualCacheView(next,exception);
 665           p=GetCacheViewVirtualPixels(image_view,0,y,next->columns,1,exception);
 666           if (p == (const Quantum *) NULL)
 667             {
 668               image_view=DestroyCacheView(image_view);
 669               break;
 670             }
 671           for (x=0; x < (ssize_t) next->columns; x++)
 672           {
 673             register ssize_t
 674               i;
 675
 676             if (GetPixelMask(next,p) != 0)
 677               {
 678                 p+=GetPixelChannels(next);
 679                 continue;
 680               }
 681             for (i=0; i < (ssize_t) GetPixelChannels(next); i++)
 682             {
 683               PixelChannel
 684                 channel;
 685
 686               PixelTrait
 687                 evaluate_traits,
 688                 traits;
 689
 690               channel=GetPixelChannelMapChannel(evaluate_image,i);
 691               traits=GetPixelChannelMapTraits(next,channel);
 692               evaluate_traits=GetPixelChannelMapTraits(evaluate_image,channel);
 693               if ((traits == UndefinedPixelTrait) ||
 694                   (evaluate_traits == UndefinedPixelTrait))
 695                 continue;
 696               if ((traits & UpdatePixelTrait) == 0)
 697                 continue;
 698               evaluate_pixel[x].channel[i]=ApplyEvaluateOperator(
 699                 random_info[id],GetPixelChannel(evaluate_image,channel,p),j ==
 700                 0 ? AddEvaluateOperator : op,evaluate_pixel[x].channel[i]);
 701             }
 702             p+=GetPixelChannels(next);
 703           }
 704           image_view=DestroyCacheView(image_view);
 705           next=GetNextImageInList(next);
 706         }
 707         for (x=0; x < (ssize_t) evaluate_image->columns; x++)
 708         {
 709           register ssize_t
 710              i;
 711
 712           switch (op)
 713           {
 714             case MeanEvaluateOperator:
 715             {
 716               for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++)
 717                 evaluate_pixel[x].channel[i]/=(MagickRealType) number_images;
 718               break;
 719             }
 720             case MultiplyEvaluateOperator:
 721             {
 722               for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++)
 723               {
 724                 register ssize_t
 725                   j;
 726
 727                 for (j=0; j < (ssize_t) (number_images-1); j++)
 728                   evaluate_pixel[x].channel[i]*=QuantumScale;
 729               }
 730               break;
 731             }
 732             default:
 733               break;
 734           }
 735         }
 736         for (x=0; x < (ssize_t) evaluate_image->columns; x++)
 737         {
 738           register ssize_t
 739             i;
 740
 741           if (GetPixelMask(evaluate_image,q) != 0)
 742             {
 743               q+=GetPixelChannels(evaluate_image);
 744               continue;
 745             }
 746           for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++)
 747           {
 748             PixelChannel
 749               channel;
 750
 751             PixelTrait
 752               traits;
 753
 754             channel=GetPixelChannelMapChannel(evaluate_image,i);
 755             traits=GetPixelChannelMapTraits(evaluate_image,channel);
 756             if (traits == UndefinedPixelTrait)
 757               continue;
 758             if ((traits & UpdatePixelTrait) == 0)
 759               continue;
 760             q[i]=ClampToQuantum(evaluate_pixel[x].channel[i]);
 761           }
 762           q+=GetPixelChannels(evaluate_image);
 763         }
 764         if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
 765           status=MagickFalse;
 766         if (images->progress_monitor != (MagickProgressMonitor) NULL)
 767           {
 768             MagickBooleanType
 769               proceed;
 770
 771 #if   defined(MAGICKCORE_OPENMP_SUPPORT)
 772             #pragma omp critical (MagickCore_EvaluateImages)
 773 #endif
 774             proceed=SetImageProgress(images,EvaluateImageTag,progress++,
 775               evaluate_image->rows);
 776             if (proceed == MagickFalse)
 777               status=MagickFalse;
 778           }
 779       }
 780     }
 781   evaluate_view=DestroyCacheView(evaluate_view);
 782   evaluate_pixels=DestroyPixelThreadSet(evaluate_pixels);
 783   random_info=DestroyRandomInfoThreadSet(random_info);
 784   if (status == MagickFalse)
 785     evaluate_image=DestroyImage(evaluate_image);
 786   return(evaluate_image);
 787 }
 788
 789 MagickExport MagickBooleanType EvaluateImage(Image *image,
 790   const MagickEvaluateOperator op,const double value,ExceptionInfo *exception)
 791 {
 792   CacheView
 793     *image_view;
 794
 795   MagickBooleanType
 796     status;
 797
 798   MagickOffsetType
 799     progress;
 800
 801   RandomInfo
 802     **restrict random_info;
 803
 804   ssize_t
 805     y;
 806
 807   unsigned long
 808     key;
 809
 810   assert(image != (Image *) NULL);
 811   assert(image->signature == MagickSignature);
 812   if (image->debug != MagickFalse)
 813     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 814   assert(exception != (ExceptionInfo *) NULL);
 815   assert(exception->signature == MagickSignature);
 816   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
 817     return(MagickFalse);
 818   status=MagickTrue;
 819   progress=0;
 820   random_info=AcquireRandomInfoThreadSet();
 821   key=GetRandomSecretKey(random_info[0]);
 822   image_view=AcquireAuthenticCacheView(image,exception);
 823 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 824   #pragma omp parallel for schedule(static,4) shared(progress,status) \
 825     if (((image->rows*image->columns) > 8192) && (key == ~0UL)) \
 826       num_threads(GetMagickResourceLimit(ThreadResource))
 827 #endif
 828   for (y=0; y < (ssize_t) image->rows; y++)
 829   {
 830     const int
 831       id = GetOpenMPThreadId();
 832
 833     register Quantum
 834       *restrict q;
 835
 836     register ssize_t
 837       x;
 838
 839     if (status == MagickFalse)
 840       continue;
 841     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
 842     if (q == (Quantum *) NULL)
 843       {
 844         status=MagickFalse;
 845         continue;
 846       }
 847     for (x=0; x < (ssize_t) image->columns; x++)
 848     {
 849       register ssize_t
 850         i;
 851
 852       if (GetPixelMask(image,q) != 0)
 853         {
 854           q+=GetPixelChannels(image);
 855           continue;
 856         }
 857       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 858       {
 859         PixelChannel
 860           channel;
 861
 862         PixelTrait
 863           traits;
 864
 865         channel=GetPixelChannelMapChannel(image,i);
 866         traits=GetPixelChannelMapTraits(image,channel);
 867         if (traits == UndefinedPixelTrait)
 868           continue;
 869         if ((traits & CopyPixelTrait) != 0)
 870           continue;
 871         q[i]=ClampToQuantum(ApplyEvaluateOperator(random_info[id],q[i],op,
 872           value));
 873       }
 874       q+=GetPixelChannels(image);
 875     }
 876     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
 877       status=MagickFalse;
 878     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 879       {
 880         MagickBooleanType
 881           proceed;
 882
 883 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 884         #pragma omp critical (MagickCore_EvaluateImage)
 885 #endif
 886         proceed=SetImageProgress(image,EvaluateImageTag,progress++,image->rows);
 887         if (proceed == MagickFalse)
 888           status=MagickFalse;
 889       }
 890   }
 891   image_view=DestroyCacheView(image_view);
 892   random_info=DestroyRandomInfoThreadSet(random_info);
 893   return(status);
 894 }
 895 \f
 896 /*
 897 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 898 %                                                                             %
 899 %                                                                             %
 900 %                                                                             %
 901 %     F u n c t i o n I m a g e                                               %
 902 %                                                                             %
 903 %                                                                             %
 904 %                                                                             %
 905 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 906 %
 907 %  FunctionImage() applies a value to the image with an arithmetic, relational,
 908 %  or logical operator to an image. Use these operations to lighten or darken
 909 %  an image, to increase or decrease contrast in an image, or to produce the
 910 %  "negative" of an image.
 911 %
 912 %  The format of the FunctionImage method is:
 913 %
 914 %      MagickBooleanType FunctionImage(Image *image,
 915 %        const MagickFunction function,const ssize_t number_parameters,
 916 %        const double *parameters,ExceptionInfo *exception)
 917 %
 918 %  A description of each parameter follows:
 919 %
 920 %    o image: the image.
 921 %
 922 %    o function: A channel function.
 923 %
 924 %    o parameters: one or more parameters.
 925 %
 926 %    o exception: return any errors or warnings in this structure.
 927 %
 928 */
 929
 930 static Quantum ApplyFunction(Quantum pixel,const MagickFunction function,
 931   const size_t number_parameters,const double *parameters,
 932   ExceptionInfo *exception)
 933 {
 934   MagickRealType
 935     result;
 936
 937   register ssize_t
 938     i;
 939
 940   (void) exception;
 941   result=0.0;
 942   switch (function)
 943   {
 944     case PolynomialFunction:
 945     {
 946       /*
 947         Polynomial: polynomial constants, highest to lowest order (e.g. c0*x^3+
 948         c1*x^2+c2*x+c3).
 949       */
 950       result=0.0;
 951       for (i=0; i < (ssize_t) number_parameters; i++)
 952         result=result*QuantumScale*pixel+parameters[i];
 953       result*=QuantumRange;
 954       break;
 955     }
 956     case SinusoidFunction:
 957     {
 958       MagickRealType
 959         amplitude,
 960         bias,
 961         frequency,
 962         phase;
 963
 964       /*
 965         Sinusoid: frequency, phase, amplitude, bias.
 966       */
 967       frequency=(number_parameters >= 1) ? parameters[0] : 1.0;
 968       phase=(number_parameters >= 2) ? parameters[1] : 0.0;
 969       amplitude=(number_parameters >= 3) ? parameters[2] : 0.5;
 970       bias=(number_parameters >= 4) ? parameters[3] : 0.5;
 971       result=(MagickRealType) (QuantumRange*(amplitude*sin((double) (2.0*
 972         MagickPI*(frequency*QuantumScale*pixel+phase/360.0)))+bias));
 973       break;
 974     }
 975     case ArcsinFunction:
 976     {
 977       MagickRealType
 978         bias,
 979         center,
 980         range,
 981         width;
 982
 983       /*
 984         Arcsin (peged at range limits for invalid results): width, center,
 985         range, and bias.
 986       */
 987       width=(number_parameters >= 1) ? parameters[0] : 1.0;
 988       center=(number_parameters >= 2) ? parameters[1] : 0.5;
 989       range=(number_parameters >= 3) ? parameters[2] : 1.0;
 990       bias=(number_parameters >= 4) ? parameters[3] : 0.5;
 991       result=2.0/width*(QuantumScale*pixel-center);
 992       if ( result <= -1.0 )
 993         result=bias-range/2.0;
 994       else
 995         if (result >= 1.0)
 996           result=bias+range/2.0;
 997         else
 998           result=(MagickRealType) (range/MagickPI*asin((double) result)+bias);
 999       result*=QuantumRange;
1000       break;
1001     }
1002     case ArctanFunction:
1003     {
1004       MagickRealType
1005         center,
1006         bias,
1007         range,
1008         slope;
1009
1010       /*
1011         Arctan: slope, center, range, and bias.
1012       */
1013       slope=(number_parameters >= 1) ? parameters[0] : 1.0;
1014       center=(number_parameters >= 2) ? parameters[1] : 0.5;
1015       range=(number_parameters >= 3) ? parameters[2] : 1.0;
1016       bias=(number_parameters >= 4) ? parameters[3] : 0.5;
1017       result=(MagickRealType) (MagickPI*slope*(QuantumScale*pixel-center));
1018       result=(MagickRealType) (QuantumRange*(range/MagickPI*atan((double)
1019         result)+bias));
1020       break;
1021     }
1022     case UndefinedFunction:
1023       break;
1024   }
1025   return(ClampToQuantum(result));
1026 }
1027
1028 MagickExport MagickBooleanType FunctionImage(Image *image,
1029   const MagickFunction function,const size_t number_parameters,
1030   const double *parameters,ExceptionInfo *exception)
1031 {
1032 #define FunctionImageTag  "Function/Image "
1033
1034   CacheView
1035     *image_view;
1036
1037   MagickBooleanType
1038     status;
1039
1040   MagickOffsetType
1041     progress;
1042
1043   ssize_t
1044     y;
1045
1046   assert(image != (Image *) NULL);
1047   assert(image->signature == MagickSignature);
1048   if (image->debug != MagickFalse)
1049     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1050   assert(exception != (ExceptionInfo *) NULL);
1051   assert(exception->signature == MagickSignature);
1052   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
1053     return(MagickFalse);
1054   status=MagickTrue;
1055   progress=0;
1056   image_view=AcquireAuthenticCacheView(image,exception);
1057 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1058   #pragma omp parallel for schedule(static,4) shared(progress,status) \
1059     if ((image->rows*image->columns) > 8192) \
1060       num_threads(GetMagickResourceLimit(ThreadResource))
1061 #endif
1062   for (y=0; y < (ssize_t) image->rows; y++)
1063   {
1064     register Quantum
1065       *restrict q;
1066
1067     register ssize_t
1068       x;
1069
1070     if (status == MagickFalse)
1071       continue;
1072     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
1073     if (q == (Quantum *) NULL)
1074       {
1075         status=MagickFalse;
1076         continue;
1077       }
1078     for (x=0; x < (ssize_t) image->columns; x++)
1079     {
1080       register ssize_t
1081         i;
1082
1083       if (GetPixelMask(image,q) != 0)
1084         {
1085           q+=GetPixelChannels(image);
1086           continue;
1087         }
1088       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1089       {
1090         PixelChannel
1091           channel;
1092
1093         PixelTrait
1094           traits;
1095
1096         channel=GetPixelChannelMapChannel(image,i);
1097         traits=GetPixelChannelMapTraits(image,channel);
1098         if (traits == UndefinedPixelTrait)
1099           continue;
1100         if ((traits & UpdatePixelTrait) == 0)
1101           continue;
1102         q[i]=ApplyFunction(q[i],function,number_parameters,parameters,
1103           exception);
1104       }
1105       q+=GetPixelChannels(image);
1106     }
1107     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1108       status=MagickFalse;
1109     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1110       {
1111         MagickBooleanType
1112           proceed;
1113
1114 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1115         #pragma omp critical (MagickCore_FunctionImage)
1116 #endif
1117         proceed=SetImageProgress(image,FunctionImageTag,progress++,image->rows);
1118         if (proceed == MagickFalse)
1119           status=MagickFalse;
1120       }
1121   }
1122   image_view=DestroyCacheView(image_view);
1123   return(status);
1124 }
1125 \f
1126 /*
1127 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1128 %                                                                             %
1129 %                                                                             %
1130 %                                                                             %
1131 %   G e t I m a g e E x t r e m a                                             %
1132 %                                                                             %
1133 %                                                                             %
1134 %                                                                             %
1135 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1136 %
1137 %  GetImageExtrema() returns the extrema of one or more image channels.
1138 %
1139 %  The format of the GetImageExtrema method is:
1140 %
1141 %      MagickBooleanType GetImageExtrema(const Image *image,size_t *minima,
1142 %        size_t *maxima,ExceptionInfo *exception)
1143 %
1144 %  A description of each parameter follows:
1145 %
1146 %    o image: the image.
1147 %
1148 %    o minima: the minimum value in the channel.
1149 %
1150 %    o maxima: the maximum value in the channel.
1151 %
1152 %    o exception: return any errors or warnings in this structure.
1153 %
1154 */
1155 MagickExport MagickBooleanType GetImageExtrema(const Image *image,
1156   size_t *minima,size_t *maxima,ExceptionInfo *exception)
1157 {
1158   double
1159     max,
1160     min;
1161
1162   MagickBooleanType
1163     status;
1164
1165   assert(image != (Image *) NULL);
1166   assert(image->signature == MagickSignature);
1167   if (image->debug != MagickFalse)
1168     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1169   status=GetImageRange(image,&min,&max,exception);
1170   *minima=(size_t) ceil(min-0.5);
1171   *maxima=(size_t) floor(max+0.5);
1172   return(status);
1173 }
1174 \f
1175 /*
1176 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1177 %                                                                             %
1178 %                                                                             %
1179 %                                                                             %
1180 %   G e t I m a g e M e a n                                                   %
1181 %                                                                             %
1182 %                                                                             %
1183 %                                                                             %
1184 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1185 %
1186 %  GetImageMean() returns the mean and standard deviation of one or more
1187 %  image channels.
1188 %
1189 %  The format of the GetImageMean method is:
1190 %
1191 %      MagickBooleanType GetImageMean(const Image *image,double *mean,
1192 %        double *standard_deviation,ExceptionInfo *exception)
1193 %
1194 %  A description of each parameter follows:
1195 %
1196 %    o image: the image.
1197 %
1198 %    o mean: the average value in the channel.
1199 %
1200 %    o standard_deviation: the standard deviation of the channel.
1201 %
1202 %    o exception: return any errors or warnings in this structure.
1203 %
1204 */
1205 MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean,
1206   double *standard_deviation,ExceptionInfo *exception)
1207 {
1208   ChannelStatistics
1209     *channel_statistics;
1210
1211   register ssize_t
1212     i;
1213
1214   size_t
1215     area;
1216
1217   assert(image != (Image *) NULL);
1218   assert(image->signature == MagickSignature);
1219   if (image->debug != MagickFalse)
1220     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1221   channel_statistics=GetImageStatistics(image,exception);
1222   if (channel_statistics == (ChannelStatistics *) NULL)
1223     return(MagickFalse);
1224   area=0;
1225   channel_statistics[CompositePixelChannel].mean=0.0;
1226   channel_statistics[CompositePixelChannel].standard_deviation=0.0;
1227   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1228   {
1229     PixelChannel
1230       channel;
1231
1232     PixelTrait
1233       traits;
1234
1235     channel=GetPixelChannelMapChannel(image,i);
1236     traits=GetPixelChannelMapTraits(image,channel);
1237     if (traits == UndefinedPixelTrait)
1238       continue;
1239     if ((traits & UpdatePixelTrait) == 0)
1240       continue;
1241     channel_statistics[CompositePixelChannel].mean+=channel_statistics[i].mean;
1242     channel_statistics[CompositePixelChannel].standard_deviation+=
1243       channel_statistics[i].variance-channel_statistics[i].mean*
1244       channel_statistics[i].mean;
1245     area++;
1246   }
1247   channel_statistics[CompositePixelChannel].mean/=area;
1248   channel_statistics[CompositePixelChannel].standard_deviation=
1249     sqrt(channel_statistics[CompositePixelChannel].standard_deviation/area);
1250   *mean=channel_statistics[CompositePixelChannel].mean;
1251   *standard_deviation=
1252     channel_statistics[CompositePixelChannel].standard_deviation;
1253   channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1254     channel_statistics);
1255   return(MagickTrue);
1256 }
1257 \f
1258 /*
1259 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1260 %                                                                             %
1261 %                                                                             %
1262 %                                                                             %
1263 %   G e t I m a g e K u r t o s i s                                           %
1264 %                                                                             %
1265 %                                                                             %
1266 %                                                                             %
1267 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1268 %
1269 %  GetImageKurtosis() returns the kurtosis and skewness of one or more
1270 %  image channels.
1271 %
1272 %  The format of the GetImageKurtosis method is:
1273 %
1274 %      MagickBooleanType GetImageKurtosis(const Image *image,double *kurtosis,
1275 %        double *skewness,ExceptionInfo *exception)
1276 %
1277 %  A description of each parameter follows:
1278 %
1279 %    o image: the image.
1280 %
1281 %    o kurtosis: the kurtosis of the channel.
1282 %
1283 %    o skewness: the skewness of the channel.
1284 %
1285 %    o exception: return any errors or warnings in this structure.
1286 %
1287 */
1288 MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
1289   double *kurtosis,double *skewness,ExceptionInfo *exception)
1290 {
1291   CacheView
1292     *image_view;
1293
1294   double
1295     area,
1296     mean,
1297     standard_deviation,
1298     sum_squares,
1299     sum_cubes,
1300     sum_fourth_power;
1301
1302   MagickBooleanType
1303     status;
1304
1305   ssize_t
1306     y;
1307
1308   assert(image != (Image *) NULL);
1309   assert(image->signature == MagickSignature);
1310   if (image->debug != MagickFalse)
1311     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1312   status=MagickTrue;
1313   *kurtosis=0.0;
1314   *skewness=0.0;
1315   area=0.0;
1316   mean=0.0;
1317   standard_deviation=0.0;
1318   sum_squares=0.0;
1319   sum_cubes=0.0;
1320   sum_fourth_power=0.0;
1321   image_view=AcquireVirtualCacheView(image,exception);
1322 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1323   #pragma omp parallel for schedule(static) shared(status) \
1324     if ((image->rows*image->columns) > 8192) \
1325       num_threads(GetMagickResourceLimit(ThreadResource))
1326 #endif
1327   for (y=0; y < (ssize_t) image->rows; y++)
1328   {
1329     register const Quantum
1330       *restrict p;
1331
1332     register ssize_t
1333       x;
1334
1335     if (status == MagickFalse)
1336       continue;
1337     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1338     if (p == (const Quantum *) NULL)
1339       {
1340         status=MagickFalse;
1341         continue;
1342       }
1343     for (x=0; x < (ssize_t) image->columns; x++)
1344     {
1345       register ssize_t
1346         i;
1347
1348       if (GetPixelMask(image,p) != 0)
1349         {
1350           p+=GetPixelChannels(image);
1351           continue;
1352         }
1353       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1354       {
1355         PixelChannel
1356           channel;
1357
1358         PixelTrait
1359           traits;
1360
1361         channel=GetPixelChannelMapChannel(image,i);
1362         traits=GetPixelChannelMapTraits(image,channel);
1363         if (traits == UndefinedPixelTrait)
1364           continue;
1365         if ((traits & UpdatePixelTrait) == 0)
1366           continue;
1367 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1368         #pragma omp critical (MagickCore_GetImageKurtosis)
1369 #endif
1370         {
1371           mean+=p[i];
1372           sum_squares+=(double) p[i]*p[i];
1373           sum_cubes+=(double) p[i]*p[i]*p[i];
1374           sum_fourth_power+=(double) p[i]*p[i]*p[i]*p[i];
1375           area++;
1376         }
1377       }
1378       p+=GetPixelChannels(image);
1379     }
1380   }
1381   image_view=DestroyCacheView(image_view);
1382   if (area != 0.0)
1383     {
1384       mean/=area;
1385       sum_squares/=area;
1386       sum_cubes/=area;
1387       sum_fourth_power/=area;
1388     }
1389   standard_deviation=sqrt(sum_squares-(mean*mean));
1390   if (standard_deviation != 0.0)
1391     {
1392       *kurtosis=sum_fourth_power-4.0*mean*sum_cubes+6.0*mean*mean*sum_squares-
1393         3.0*mean*mean*mean*mean;
1394       *kurtosis/=standard_deviation*standard_deviation*standard_deviation*
1395         standard_deviation;
1396       *kurtosis-=3.0;
1397       *skewness=sum_cubes-3.0*mean*sum_squares+2.0*mean*mean*mean;
1398       *skewness/=standard_deviation*standard_deviation*standard_deviation;
1399     }
1400   return(status);
1401 }
1402 \f
1403 /*
1404 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1405 %                                                                             %
1406 %                                                                             %
1407 %                                                                             %
1408 %   G e t I m a g e R a n g e                                                 %
1409 %                                                                             %
1410 %                                                                             %
1411 %                                                                             %
1412 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1413 %
1414 %  GetImageRange() returns the range of one or more image channels.
1415 %
1416 %  The format of the GetImageRange method is:
1417 %
1418 %      MagickBooleanType GetImageRange(const Image *image,double *minima,
1419 %        double *maxima,ExceptionInfo *exception)
1420 %
1421 %  A description of each parameter follows:
1422 %
1423 %    o image: the image.
1424 %
1425 %    o minima: the minimum value in the channel.
1426 %
1427 %    o maxima: the maximum value in the channel.
1428 %
1429 %    o exception: return any errors or warnings in this structure.
1430 %
1431 */
1432 MagickExport MagickBooleanType GetImageRange(const Image *image,double *minima,
1433   double *maxima,ExceptionInfo *exception)
1434 {
1435   CacheView
1436     *image_view;
1437
1438   MagickBooleanType
1439     status;
1440
1441   ssize_t
1442     y;
1443
1444   assert(image != (Image *) NULL);
1445   assert(image->signature == MagickSignature);
1446   if (image->debug != MagickFalse)
1447     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1448   status=MagickTrue;
1449   *maxima=(-MagickHuge);
1450   *minima=MagickHuge;
1451   image_view=AcquireVirtualCacheView(image,exception);
1452 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1453   #pragma omp parallel for schedule(static) shared(status) \
1454     if ((image->rows*image->columns) > 8192) \
1455       num_threads(GetMagickResourceLimit(ThreadResource))
1456 #endif
1457   for (y=0; y < (ssize_t) image->rows; y++)
1458   {
1459     register const Quantum
1460       *restrict p;
1461
1462     register ssize_t
1463       x;
1464
1465     if (status == MagickFalse)
1466       continue;
1467     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1468     if (p == (const Quantum *) NULL)
1469       {
1470         status=MagickFalse;
1471         continue;
1472       }
1473     for (x=0; x < (ssize_t) image->columns; x++)
1474     {
1475       register ssize_t
1476         i;
1477
1478       if (GetPixelMask(image,p) != 0)
1479         {
1480           p+=GetPixelChannels(image);
1481           continue;
1482         }
1483       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1484       {
1485         PixelChannel
1486           channel;
1487
1488         PixelTrait
1489           traits;
1490
1491         channel=GetPixelChannelMapChannel(image,i);
1492         traits=GetPixelChannelMapTraits(image,channel);
1493         if (traits == UndefinedPixelTrait)
1494           continue;
1495         if ((traits & UpdatePixelTrait) == 0)
1496           continue;
1497 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1498         #pragma omp critical (MagickCore_GetImageRange)
1499 #endif
1500         {
1501           if ((double) p[i] < *minima)
1502             *minima=(double) p[i];
1503           if ((double) p[i] > *maxima)
1504             *maxima=(double) p[i];
1505         }
1506       }
1507       p+=GetPixelChannels(image);
1508     }
1509   }
1510   image_view=DestroyCacheView(image_view);
1511   return(status);
1512 }
1513 \f
1514 /*
1515 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1516 %                                                                             %
1517 %                                                                             %
1518 %                                                                             %
1519 %   G e t I m a g e S t a t i s t i c s                                       %
1520 %                                                                             %
1521 %                                                                             %
1522 %                                                                             %
1523 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1524 %
1525 %  GetImageStatistics() returns statistics for each channel in the
1526 %  image.  The statistics include the channel depth, its minima, maxima, mean,
1527 %  standard deviation, kurtosis and skewness.  You can access the red channel
1528 %  mean, for example, like this:
1529 %
1530 %      channel_statistics=GetImageStatistics(image,exception);
1531 %      red_mean=channel_statistics[RedPixelChannel].mean;
1532 %
1533 %  Use MagickRelinquishMemory() to free the statistics buffer.
1534 %
1535 %  The format of the GetImageStatistics method is:
1536 %
1537 %      ChannelStatistics *GetImageStatistics(const Image *image,
1538 %        ExceptionInfo *exception)
1539 %
1540 %  A description of each parameter follows:
1541 %
1542 %    o image: the image.
1543 %
1544 %    o exception: return any errors or warnings in this structure.
1545 %
1546 */
1547
1548 static size_t GetImageChannels(const Image *image)
1549 {
1550   register ssize_t
1551     i;
1552
1553   size_t
1554     channels;
1555
1556   channels=0;
1557   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1558   {
1559     PixelChannel
1560       channel;
1561
1562     PixelTrait
1563       traits;
1564
1565     channel=GetPixelChannelMapChannel(image,i);
1566     traits=GetPixelChannelMapTraits(image,channel);
1567     if ((traits & UpdatePixelTrait) != 0)
1568       channels++;
1569   }
1570   return(channels);
1571 }
1572
1573 MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
1574   ExceptionInfo *exception)
1575 {
1576   ChannelStatistics
1577     *channel_statistics;
1578
1579   double
1580     area;
1581
1582   MagickStatusType
1583     status;
1584
1585   QuantumAny
1586     range;
1587
1588   register ssize_t
1589     i;
1590
1591   size_t
1592     channels,
1593     depth;
1594
1595   ssize_t
1596     y;
1597
1598   assert(image != (Image *) NULL);
1599   assert(image->signature == MagickSignature);
1600   if (image->debug != MagickFalse)
1601     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1602   channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(
1603     MaxPixelChannels+1,sizeof(*channel_statistics));
1604   if (channel_statistics == (ChannelStatistics *) NULL)
1605     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1606   (void) ResetMagickMemory(channel_statistics,0,(MaxPixelChannels+1)*
1607     sizeof(*channel_statistics));
1608   for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
1609   {
1610     channel_statistics[i].depth=1;
1611     channel_statistics[i].maxima=(-MagickHuge);
1612     channel_statistics[i].minima=MagickHuge;
1613   }
1614   for (y=0; y < (ssize_t) image->rows; y++)
1615   {
1616     register const Quantum
1617       *restrict p;
1618
1619     register ssize_t
1620       x;
1621
1622     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1623     if (p == (const Quantum *) NULL)
1624       break;
1625     for (x=0; x < (ssize_t) image->columns; x++)
1626     {
1627       register ssize_t
1628         i;
1629
1630       if (GetPixelMask(image,p) != 0)
1631         {
1632           p+=GetPixelChannels(image);
1633           continue;
1634         }
1635       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1636       {
1637         PixelChannel
1638           channel;
1639
1640         PixelTrait
1641           traits;
1642
1643         channel=GetPixelChannelMapChannel(image,i);
1644         traits=GetPixelChannelMapTraits(image,channel);
1645         if (traits == UndefinedPixelTrait)
1646           continue;
1647         if (channel_statistics[channel].depth != MAGICKCORE_QUANTUM_DEPTH)
1648           {
1649             depth=channel_statistics[channel].depth;
1650             range=GetQuantumRange(depth);
1651             status=p[i] != ScaleAnyToQuantum(ScaleQuantumToAny(p[i],range),
1652               range) ? MagickTrue : MagickFalse;
1653             if (status != MagickFalse)
1654               {
1655                 channel_statistics[channel].depth++;
1656                 continue;
1657               }
1658           }
1659         if ((double) p[i] < channel_statistics[channel].minima)
1660           channel_statistics[channel].minima=(double) p[i];
1661         if ((double) p[i] > channel_statistics[channel].maxima)
1662           channel_statistics[channel].maxima=(double) p[i];
1663         channel_statistics[channel].sum+=p[i];
1664         channel_statistics[channel].sum_squared+=(double) p[i]*p[i];
1665         channel_statistics[channel].sum_cubed+=(double) p[i]*p[i]*p[i];
1666         channel_statistics[channel].sum_fourth_power+=(double) p[i]*p[i]*p[i]*
1667           p[i];
1668       }
1669       p+=GetPixelChannels(image);
1670     }
1671   }
1672   area=(double) image->columns*image->rows;
1673   for (i=0; i < (ssize_t) MaxPixelChannels; i++)
1674   {
1675     channel_statistics[i].sum/=area;
1676     channel_statistics[i].sum_squared/=area;
1677     channel_statistics[i].sum_cubed/=area;
1678     channel_statistics[i].sum_fourth_power/=area;
1679     channel_statistics[i].mean=channel_statistics[i].sum;
1680     channel_statistics[i].variance=channel_statistics[i].sum_squared;
1681     channel_statistics[i].standard_deviation=sqrt(
1682       channel_statistics[i].variance-(channel_statistics[i].mean*
1683       channel_statistics[i].mean));
1684   }
1685   for (i=0; i < (ssize_t) MaxPixelChannels; i++)
1686   {
1687     channel_statistics[CompositePixelChannel].depth=(size_t) EvaluateMax(
1688       (double) channel_statistics[CompositePixelChannel].depth,(double)
1689       channel_statistics[i].depth);
1690     channel_statistics[CompositePixelChannel].minima=MagickMin(
1691       channel_statistics[CompositePixelChannel].minima,
1692       channel_statistics[i].minima);
1693     channel_statistics[CompositePixelChannel].maxima=EvaluateMax(
1694       channel_statistics[CompositePixelChannel].maxima,
1695       channel_statistics[i].maxima);
1696     channel_statistics[CompositePixelChannel].sum+=channel_statistics[i].sum;
1697     channel_statistics[CompositePixelChannel].sum_squared+=
1698       channel_statistics[i].sum_squared;
1699     channel_statistics[CompositePixelChannel].sum_cubed+=
1700       channel_statistics[i].sum_cubed;
1701     channel_statistics[CompositePixelChannel].sum_fourth_power+=
1702       channel_statistics[i].sum_fourth_power;
1703     channel_statistics[CompositePixelChannel].mean+=channel_statistics[i].mean;
1704     channel_statistics[CompositePixelChannel].variance+=
1705       channel_statistics[i].variance-channel_statistics[i].mean*
1706       channel_statistics[i].mean;
1707     channel_statistics[CompositePixelChannel].standard_deviation+=
1708       channel_statistics[i].variance-channel_statistics[i].mean*
1709       channel_statistics[i].mean;
1710   }
1711   channels=GetImageChannels(image);
1712   channel_statistics[CompositePixelChannel].sum/=channels;
1713   channel_statistics[CompositePixelChannel].sum_squared/=channels;
1714   channel_statistics[CompositePixelChannel].sum_cubed/=channels;
1715   channel_statistics[CompositePixelChannel].sum_fourth_power/=channels;
1716   channel_statistics[CompositePixelChannel].mean/=channels;
1717   channel_statistics[CompositePixelChannel].variance/=channels;
1718   channel_statistics[CompositePixelChannel].standard_deviation=
1719     sqrt(channel_statistics[CompositePixelChannel].standard_deviation/channels);
1720   channel_statistics[CompositePixelChannel].kurtosis/=channels;
1721   channel_statistics[CompositePixelChannel].skewness/=channels;
1722   for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
1723   {
1724     if (channel_statistics[i].standard_deviation == 0.0)
1725       continue;
1726     channel_statistics[i].skewness=(channel_statistics[i].sum_cubed-3.0*
1727       channel_statistics[i].mean*channel_statistics[i].sum_squared+2.0*
1728       channel_statistics[i].mean*channel_statistics[i].mean*
1729       channel_statistics[i].mean)/(channel_statistics[i].standard_deviation*
1730       channel_statistics[i].standard_deviation*
1731       channel_statistics[i].standard_deviation);
1732     channel_statistics[i].kurtosis=(channel_statistics[i].sum_fourth_power-4.0*
1733       channel_statistics[i].mean*channel_statistics[i].sum_cubed+6.0*
1734       channel_statistics[i].mean*channel_statistics[i].mean*
1735       channel_statistics[i].sum_squared-3.0*channel_statistics[i].mean*
1736       channel_statistics[i].mean*1.0*channel_statistics[i].mean*
1737       channel_statistics[i].mean)/(channel_statistics[i].standard_deviation*
1738       channel_statistics[i].standard_deviation*
1739       channel_statistics[i].standard_deviation*
1740       channel_statistics[i].standard_deviation)-3.0;
1741   }
1742   return(channel_statistics);
1743 }
1744 \f
1745 /*
1746 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1747 %                                                                             %
1748 %                                                                             %
1749 %                                                                             %
1750 %     S t a t i s t i c I m a g e                                             %
1751 %                                                                             %
1752 %                                                                             %
1753 %                                                                             %
1754 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1755 %
1756 %  StatisticImage() makes each pixel the min / max / median / mode / etc. of
1757 %  the neighborhood of the specified width and height.
1758 %
1759 %  The format of the StatisticImage method is:
1760 %
1761 %      Image *StatisticImage(const Image *image,const StatisticType type,
1762 %        const size_t width,const size_t height,ExceptionInfo *exception)
1763 %
1764 %  A description of each parameter follows:
1765 %
1766 %    o image: the image.
1767 %
1768 %    o type: the statistic type (median, mode, etc.).
1769 %
1770 %    o width: the width of the pixel neighborhood.
1771 %
1772 %    o height: the height of the pixel neighborhood.
1773 %
1774 %    o exception: return any errors or warnings in this structure.
1775 %
1776 */
1777
1778 typedef struct _SkipNode
1779 {
1780   size_t
1781     next[9],
1782     count,
1783     signature;
1784 } SkipNode;
1785
1786 typedef struct _SkipList
1787 {
1788   ssize_t
1789     level;
1790
1791   SkipNode
1792     *nodes;
1793 } SkipList;
1794
1795 typedef struct _PixelList
1796 {
1797   size_t
1798     length,
1799     seed;
1800
1801   SkipList
1802     skip_list;
1803
1804   size_t
1805     signature;
1806 } PixelList;
1807
1808 static PixelList *DestroyPixelList(PixelList *pixel_list)
1809 {
1810   if (pixel_list == (PixelList *) NULL)
1811     return((PixelList *) NULL);
1812   if (pixel_list->skip_list.nodes != (SkipNode *) NULL)
1813     pixel_list->skip_list.nodes=(SkipNode *) RelinquishMagickMemory(
1814       pixel_list->skip_list.nodes);
1815   pixel_list=(PixelList *) RelinquishMagickMemory(pixel_list);
1816   return(pixel_list);
1817 }
1818
1819 static PixelList **DestroyPixelListThreadSet(PixelList **pixel_list)
1820 {
1821   register ssize_t
1822     i;
1823
1824   assert(pixel_list != (PixelList **) NULL);
1825   for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
1826     if (pixel_list[i] != (PixelList *) NULL)
1827       pixel_list[i]=DestroyPixelList(pixel_list[i]);
1828   pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list);
1829   return(pixel_list);
1830 }
1831
1832 static PixelList *AcquirePixelList(const size_t width,const size_t height)
1833 {
1834   PixelList
1835     *pixel_list;
1836
1837   pixel_list=(PixelList *) AcquireMagickMemory(sizeof(*pixel_list));
1838   if (pixel_list == (PixelList *) NULL)
1839     return(pixel_list);
1840   (void) ResetMagickMemory((void *) pixel_list,0,sizeof(*pixel_list));
1841   pixel_list->length=width*height;
1842   pixel_list->skip_list.nodes=(SkipNode *) AcquireQuantumMemory(65537UL,
1843     sizeof(*pixel_list->skip_list.nodes));
1844   if (pixel_list->skip_list.nodes == (SkipNode *) NULL)
1845     return(DestroyPixelList(pixel_list));
1846   (void) ResetMagickMemory(pixel_list->skip_list.nodes,0,65537UL*
1847     sizeof(*pixel_list->skip_list.nodes));
1848   pixel_list->signature=MagickSignature;
1849   return(pixel_list);
1850 }
1851
1852 static PixelList **AcquirePixelListThreadSet(const size_t width,
1853   const size_t height)
1854 {
1855   PixelList
1856     **pixel_list;
1857
1858   register ssize_t
1859     i;
1860
1861   size_t
1862     number_threads;
1863
1864   number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
1865   pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
1866     sizeof(*pixel_list));
1867   if (pixel_list == (PixelList **) NULL)
1868     return((PixelList **) NULL);
1869   (void) ResetMagickMemory(pixel_list,0,number_threads*sizeof(*pixel_list));
1870   for (i=0; i < (ssize_t) number_threads; i++)
1871   {
1872     pixel_list[i]=AcquirePixelList(width,height);
1873     if (pixel_list[i] == (PixelList *) NULL)
1874       return(DestroyPixelListThreadSet(pixel_list));
1875   }
1876   return(pixel_list);
1877 }
1878
1879 static void AddNodePixelList(PixelList *pixel_list,const size_t color)
1880 {
1881   register SkipList
1882     *p;
1883
1884   register ssize_t
1885     level;
1886
1887   size_t
1888     search,
1889     update[9];
1890
1891   /*
1892     Initialize the node.
1893   */
1894   p=(&pixel_list->skip_list);
1895   p->nodes[color].signature=pixel_list->signature;
1896   p->nodes[color].count=1;
1897   /*
1898     Determine where it belongs in the list.
1899   */
1900   search=65536UL;
1901   for (level=p->level; level >= 0; level--)
1902   {
1903     while (p->nodes[search].next[level] < color)
1904       search=p->nodes[search].next[level];
1905     update[level]=search;
1906   }
1907   /*
1908     Generate a pseudo-random level for this node.
1909   */
1910   for (level=0; ; level++)
1911   {
1912     pixel_list->seed=(pixel_list->seed*42893621L)+1L;
1913     if ((pixel_list->seed & 0x300) != 0x300)
1914       break;
1915   }
1916   if (level > 8)
1917     level=8;
1918   if (level > (p->level+2))
1919     level=p->level+2;
1920   /*
1921     If we're raising the list's level, link back to the root node.
1922   */
1923   while (level > p->level)
1924   {
1925     p->level++;
1926     update[p->level]=65536UL;
1927   }
1928   /*
1929     Link the node into the skip-list.
1930   */
1931   do
1932   {
1933     p->nodes[color].next[level]=p->nodes[update[level]].next[level];
1934     p->nodes[update[level]].next[level]=color;
1935   } while (level-- > 0);
1936 }
1937
1938 static inline void GetMaximumPixelList(PixelList *pixel_list,Quantum *pixel)
1939 {
1940   register SkipList
1941     *p;
1942
1943   size_t
1944     color,
1945     maximum;
1946
1947   ssize_t
1948     count;
1949
1950   /*
1951     Find the maximum value for each of the color.
1952   */
1953   p=(&pixel_list->skip_list);
1954   color=65536L;
1955   count=0;
1956   maximum=p->nodes[color].next[0];
1957   do
1958   {
1959     color=p->nodes[color].next[0];
1960     if (color > maximum)
1961       maximum=color;
1962     count+=p->nodes[color].count;
1963   } while (count < (ssize_t) pixel_list->length);
1964   *pixel=ScaleShortToQuantum((unsigned short) maximum);
1965 }
1966
1967 static inline void GetMeanPixelList(PixelList *pixel_list,Quantum *pixel)
1968 {
1969   MagickRealType
1970     sum;
1971
1972   register SkipList
1973     *p;
1974
1975   size_t
1976     color;
1977
1978   ssize_t
1979     count;
1980
1981   /*
1982     Find the mean value for each of the color.
1983   */
1984   p=(&pixel_list->skip_list);
1985   color=65536L;
1986   count=0;
1987   sum=0.0;
1988   do
1989   {
1990     color=p->nodes[color].next[0];
1991     sum+=(MagickRealType) p->nodes[color].count*color;
1992     count+=p->nodes[color].count;
1993   } while (count < (ssize_t) pixel_list->length);
1994   sum/=pixel_list->length;
1995   *pixel=ScaleShortToQuantum((unsigned short) sum);
1996 }
1997
1998 static inline void GetMedianPixelList(PixelList *pixel_list,Quantum *pixel)
1999 {
2000   register SkipList
2001     *p;
2002
2003   size_t
2004     color;
2005
2006   ssize_t
2007     count;
2008
2009   /*
2010     Find the median value for each of the color.
2011   */
2012   p=(&pixel_list->skip_list);
2013   color=65536L;
2014   count=0;
2015   do
2016   {
2017     color=p->nodes[color].next[0];
2018     count+=p->nodes[color].count;
2019   } while (count <= (ssize_t) (pixel_list->length >> 1));
2020   *pixel=ScaleShortToQuantum((unsigned short) color);
2021 }
2022
2023 static inline void GetMinimumPixelList(PixelList *pixel_list,Quantum *pixel)
2024 {
2025   register SkipList
2026     *p;
2027
2028   size_t
2029     color,
2030     minimum;
2031
2032   ssize_t
2033     count;
2034
2035   /*
2036     Find the minimum value for each of the color.
2037   */
2038   p=(&pixel_list->skip_list);
2039   count=0;
2040   color=65536UL;
2041   minimum=p->nodes[color].next[0];
2042   do
2043   {
2044     color=p->nodes[color].next[0];
2045     if (color < minimum)
2046       minimum=color;
2047     count+=p->nodes[color].count;
2048   } while (count < (ssize_t) pixel_list->length);
2049   *pixel=ScaleShortToQuantum((unsigned short) minimum);
2050 }
2051
2052 static inline void GetModePixelList(PixelList *pixel_list,Quantum *pixel)
2053 {
2054   register SkipList
2055     *p;
2056
2057   size_t
2058     color,
2059     max_count,
2060     mode;
2061
2062   ssize_t
2063     count;
2064
2065   /*
2066     Make each pixel the 'predominant color' of the specified neighborhood.
2067   */
2068   p=(&pixel_list->skip_list);
2069   color=65536L;
2070   mode=color;
2071   max_count=p->nodes[mode].count;
2072   count=0;
2073   do
2074   {
2075     color=p->nodes[color].next[0];
2076     if (p->nodes[color].count > max_count)
2077       {
2078         mode=color;
2079         max_count=p->nodes[mode].count;
2080       }
2081     count+=p->nodes[color].count;
2082   } while (count < (ssize_t) pixel_list->length);
2083   *pixel=ScaleShortToQuantum((unsigned short) mode);
2084 }
2085
2086 static inline void GetNonpeakPixelList(PixelList *pixel_list,Quantum *pixel)
2087 {
2088   register SkipList
2089     *p;
2090
2091   size_t
2092     color,
2093     next,
2094     previous;
2095
2096   ssize_t
2097     count;
2098
2099   /*
2100     Finds the non peak value for each of the colors.
2101   */
2102   p=(&pixel_list->skip_list);
2103   color=65536L;
2104   next=p->nodes[color].next[0];
2105   count=0;
2106   do
2107   {
2108     previous=color;
2109     color=next;
2110     next=p->nodes[color].next[0];
2111     count+=p->nodes[color].count;
2112   } while (count <= (ssize_t) (pixel_list->length >> 1));
2113   if ((previous == 65536UL) && (next != 65536UL))
2114     color=next;
2115   else
2116     if ((previous != 65536UL) && (next == 65536UL))
2117       color=previous;
2118   *pixel=ScaleShortToQuantum((unsigned short) color);
2119 }
2120
2121 static inline void GetStandardDeviationPixelList(PixelList *pixel_list,
2122   Quantum *pixel)
2123 {
2124   MagickRealType
2125     sum,
2126     sum_squared;
2127
2128   register SkipList
2129     *p;
2130
2131   size_t
2132     color;
2133
2134   ssize_t
2135     count;
2136
2137   /*
2138     Find the standard-deviation value for each of the color.
2139   */
2140   p=(&pixel_list->skip_list);
2141   color=65536L;
2142   count=0;
2143   sum=0.0;
2144   sum_squared=0.0;
2145   do
2146   {
2147     register ssize_t
2148       i;
2149
2150     color=p->nodes[color].next[0];
2151     sum+=(MagickRealType) p->nodes[color].count*color;
2152     for (i=0; i < (ssize_t) p->nodes[color].count; i++)
2153       sum_squared+=((MagickRealType) color)*((MagickRealType) color);
2154     count+=p->nodes[color].count;
2155   } while (count < (ssize_t) pixel_list->length);
2156   sum/=pixel_list->length;
2157   sum_squared/=pixel_list->length;
2158   *pixel=ScaleShortToQuantum((unsigned short) sqrt(sum_squared-(sum*sum)));
2159 }
2160
2161 static inline void InsertPixelList(const Image *image,const Quantum pixel,
2162   PixelList *pixel_list)
2163 {
2164   size_t
2165     signature;
2166
2167   unsigned short
2168     index;
2169
2170   index=ScaleQuantumToShort(pixel);
2171   signature=pixel_list->skip_list.nodes[index].signature;
2172   if (signature == pixel_list->signature)
2173     {
2174       pixel_list->skip_list.nodes[index].count++;
2175       return;
2176     }
2177   AddNodePixelList(pixel_list,index);
2178 }
2179
2180 static inline MagickRealType MagickAbsoluteValue(const MagickRealType x)
2181 {
2182   if (x < 0)
2183     return(-x);
2184   return(x);
2185 }
2186
2187 static inline size_t MagickMax(const size_t x,const size_t y)
2188 {
2189   if (x > y)
2190     return(x);
2191   return(y);
2192 }
2193
2194 static void ResetPixelList(PixelList *pixel_list)
2195 {
2196   int
2197     level;
2198
2199   register SkipNode
2200     *root;
2201
2202   register SkipList
2203     *p;
2204
2205   /*
2206     Reset the skip-list.
2207   */
2208   p=(&pixel_list->skip_list);
2209   root=p->nodes+65536UL;
2210   p->level=0;
2211   for (level=0; level < 9; level++)
2212     root->next[level]=65536UL;
2213   pixel_list->seed=pixel_list->signature++;
2214 }
2215
2216 MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
2217   const size_t width,const size_t height,ExceptionInfo *exception)
2218 {
2219 #define StatisticImageTag  "Statistic/Image"
2220
2221   CacheView
2222     *image_view,
2223     *statistic_view;
2224
2225   Image
2226     *statistic_image;
2227
2228   MagickBooleanType
2229     status;
2230
2231   MagickOffsetType
2232     progress;
2233
2234   PixelList
2235     **restrict pixel_list;
2236
2237   ssize_t
2238     center,
2239     y;
2240
2241   /*
2242     Initialize statistics image attributes.
2243   */
2244   assert(image != (Image *) NULL);
2245   assert(image->signature == MagickSignature);
2246   if (image->debug != MagickFalse)
2247     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2248   assert(exception != (ExceptionInfo *) NULL);
2249   assert(exception->signature == MagickSignature);
2250   statistic_image=CloneImage(image,image->columns,image->rows,MagickTrue,
2251     exception);
2252   if (statistic_image == (Image *) NULL)
2253     return((Image *) NULL);
2254   status=SetImageStorageClass(statistic_image,DirectClass,exception);
2255   if (status == MagickFalse)
2256     {
2257       statistic_image=DestroyImage(statistic_image);
2258       return((Image *) NULL);
2259     }
2260   pixel_list=AcquirePixelListThreadSet(MagickMax(width,1),MagickMax(height,1));
2261   if (pixel_list == (PixelList **) NULL)
2262     {
2263       statistic_image=DestroyImage(statistic_image);
2264       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2265     }
2266   /*
2267     Make each pixel the min / max / median / mode / etc. of the neighborhood.
2268   */
2269   center=(ssize_t) GetPixelChannels(image)*(image->columns+MagickMax(width,1))*
2270     (MagickMax(height,1)/2L)+GetPixelChannels(image)*(MagickMax(width,1)/2L);
2271   status=MagickTrue;
2272   progress=0;
2273   image_view=AcquireVirtualCacheView(image,exception);
2274   statistic_view=AcquireAuthenticCacheView(statistic_image,exception);
2275 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2276   #pragma omp parallel for schedule(static,4) shared(progress,status) \
2277     if ((statistic_image->rows*statistic_image->columns) > 8192) \
2278       num_threads(GetMagickResourceLimit(ThreadResource))
2279 #endif
2280   for (y=0; y < (ssize_t) statistic_image->rows; y++)
2281   {
2282     const int
2283       id = GetOpenMPThreadId();
2284
2285     register const Quantum
2286       *restrict p;
2287
2288     register Quantum
2289       *restrict q;
2290
2291     register ssize_t
2292       x;
2293
2294     if (status == MagickFalse)
2295       continue;
2296     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) MagickMax(width,1)/2L),y-
2297       (ssize_t) (MagickMax(height,1)/2L),image->columns+MagickMax(width,1),
2298       MagickMax(height,1),exception);
2299     q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns,      1,exception);
2300     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2301       {
2302         status=MagickFalse;
2303         continue;
2304       }
2305     for (x=0; x < (ssize_t) statistic_image->columns; x++)
2306     {
2307       register ssize_t
2308         i;
2309
2310       if (GetPixelMask(image,p) != 0)
2311         {
2312           p+=GetPixelChannels(image);
2313           q+=GetPixelChannels(statistic_image);
2314           continue;
2315         }
2316       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2317       {
2318         PixelChannel
2319           channel;
2320
2321         PixelTrait
2322           statistic_traits,
2323           traits;
2324
2325         Quantum
2326           pixel;
2327
2328         register const Quantum
2329           *restrict pixels;
2330
2331         register ssize_t
2332           u;
2333
2334         ssize_t
2335           v;
2336
2337         channel=GetPixelChannelMapChannel(image,i);
2338         traits=GetPixelChannelMapTraits(image,channel);
2339         statistic_traits=GetPixelChannelMapTraits(statistic_image,channel);
2340         if ((traits == UndefinedPixelTrait) ||
2341             (statistic_traits == UndefinedPixelTrait))
2342           continue;
2343         if ((statistic_traits & CopyPixelTrait) != 0)
2344           {
2345             SetPixelChannel(statistic_image,channel,p[center+i],q);
2346             continue;
2347           }
2348         pixels=p;
2349         ResetPixelList(pixel_list[id]);
2350         for (v=0; v < (ssize_t) MagickMax(height,1); v++)
2351         {
2352           for (u=0; u < (ssize_t) MagickMax(width,1); u++)
2353           {
2354             InsertPixelList(image,pixels[i],pixel_list[id]);
2355             pixels+=GetPixelChannels(image);
2356           }
2357           pixels+=image->columns*GetPixelChannels(image);
2358         }
2359         switch (type)
2360         {
2361           case GradientStatistic:
2362           {
2363             MagickRealType
2364               maximum,
2365               minimum;
2366
2367             GetMinimumPixelList(pixel_list[id],&pixel);
2368             minimum=(MagickRealType) pixel;
2369             GetMaximumPixelList(pixel_list[id],&pixel);
2370             maximum=(MagickRealType) pixel;
2371             pixel=ClampToQuantum(MagickAbsoluteValue(maximum-minimum));
2372             break;
2373           }
2374           case MaximumStatistic:
2375           {
2376             GetMaximumPixelList(pixel_list[id],&pixel);
2377             break;
2378           }
2379           case MeanStatistic:
2380           {
2381             GetMeanPixelList(pixel_list[id],&pixel);
2382             break;
2383           }
2384           case MedianStatistic:
2385           default:
2386           {
2387             GetMedianPixelList(pixel_list[id],&pixel);
2388             break;
2389           }
2390           case MinimumStatistic:
2391           {
2392             GetMinimumPixelList(pixel_list[id],&pixel);
2393             break;
2394           }
2395           case ModeStatistic:
2396           {
2397             GetModePixelList(pixel_list[id],&pixel);
2398             break;
2399           }
2400           case NonpeakStatistic:
2401           {
2402             GetNonpeakPixelList(pixel_list[id],&pixel);
2403             break;
2404           }
2405           case StandardDeviationStatistic:
2406           {
2407             GetStandardDeviationPixelList(pixel_list[id],&pixel);
2408             break;
2409           }
2410         }
2411         SetPixelChannel(statistic_image,channel,pixel,q);
2412       }
2413       p+=GetPixelChannels(image);
2414       q+=GetPixelChannels(statistic_image);
2415     }
2416     if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
2417       status=MagickFalse;
2418     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2419       {
2420         MagickBooleanType
2421           proceed;
2422
2423 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2424         #pragma omp critical (MagickCore_StatisticImage)
2425 #endif
2426         proceed=SetImageProgress(image,StatisticImageTag,progress++,
2427           image->rows);
2428         if (proceed == MagickFalse)
2429           status=MagickFalse;
2430       }
2431   }
2432   statistic_view=DestroyCacheView(statistic_view);
2433   image_view=DestroyCacheView(image_view);
2434   pixel_list=DestroyPixelListThreadSet(pixel_list);
2435   return(statistic_image);
2436 }