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