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=AcquireCacheView(evaluate_image);
 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=AcquireCacheView(next);
 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=AcquireCacheView(next);
 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=AcquireCacheView(image);
 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         q[i]=ClampToQuantum(ApplyEvaluateOperator(random_info[id],q[i],op,
 861           value));
 862       }
 863       q+=GetPixelChannels(image);
 864     }
 865     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
 866       status=MagickFalse;
 867     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 868       {
 869         MagickBooleanType
 870           proceed;
 871
 872 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 873   #pragma omp critical (MagickCore_EvaluateImage)
 874 #endif
 875         proceed=SetImageProgress(image,EvaluateImageTag,progress++,image->rows);
 876         if (proceed == MagickFalse)
 877           status=MagickFalse;
 878       }
 879   }
 880   image_view=DestroyCacheView(image_view);
 881   random_info=DestroyRandomInfoThreadSet(random_info);
 882   return(status);
 883 }
 884 \f
 885 /*
 886 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 887 %                                                                             %
 888 %                                                                             %
 889 %                                                                             %
 890 %     F u n c t i o n I m a g e                                               %
 891 %                                                                             %
 892 %                                                                             %
 893 %                                                                             %
 894 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 895 %
 896 %  FunctionImage() applies a value to the image with an arithmetic, relational,
 897 %  or logical operator to an image. Use these operations to lighten or darken
 898 %  an image, to increase or decrease contrast in an image, or to produce the
 899 %  "negative" of an image.
 900 %
 901 %  The format of the FunctionImage method is:
 902 %
 903 %      MagickBooleanType FunctionImage(Image *image,
 904 %        const MagickFunction function,const ssize_t number_parameters,
 905 %        const double *parameters,ExceptionInfo *exception)
 906 %
 907 %  A description of each parameter follows:
 908 %
 909 %    o image: the image.
 910 %
 911 %    o function: A channel function.
 912 %
 913 %    o parameters: one or more parameters.
 914 %
 915 %    o exception: return any errors or warnings in this structure.
 916 %
 917 */
 918
 919 static Quantum ApplyFunction(Quantum pixel,const MagickFunction function,
 920   const size_t number_parameters,const double *parameters,
 921   ExceptionInfo *exception)
 922 {
 923   MagickRealType
 924     result;
 925
 926   register ssize_t
 927     i;
 928
 929   (void) exception;
 930   result=0.0;
 931   switch (function)
 932   {
 933     case PolynomialFunction:
 934     {
 935       /*
 936         Polynomial: polynomial constants, highest to lowest order (e.g. c0*x^3+
 937         c1*x^2+c2*x+c3).
 938       */
 939       result=0.0;
 940       for (i=0; i < (ssize_t) number_parameters; i++)
 941         result=result*QuantumScale*pixel+parameters[i];
 942       result*=QuantumRange;
 943       break;
 944     }
 945     case SinusoidFunction:
 946     {
 947       MagickRealType
 948         amplitude,
 949         bias,
 950         frequency,
 951         phase;
 952
 953       /*
 954         Sinusoid: frequency, phase, amplitude, bias.
 955       */
 956       frequency=(number_parameters >= 1) ? parameters[0] : 1.0;
 957       phase=(number_parameters >= 2) ? parameters[1] : 0.0;
 958       amplitude=(number_parameters >= 3) ? parameters[2] : 0.5;
 959       bias=(number_parameters >= 4) ? parameters[3] : 0.5;
 960       result=(MagickRealType) (QuantumRange*(amplitude*sin((double) (2.0*
 961         MagickPI*(frequency*QuantumScale*pixel+phase/360.0)))+bias));
 962       break;
 963     }
 964     case ArcsinFunction:
 965     {
 966       MagickRealType
 967         bias,
 968         center,
 969         range,
 970         width;
 971
 972       /*
 973         Arcsin (peged at range limits for invalid results): width, center,
 974         range, and bias.
 975       */
 976       width=(number_parameters >= 1) ? parameters[0] : 1.0;
 977       center=(number_parameters >= 2) ? parameters[1] : 0.5;
 978       range=(number_parameters >= 3) ? parameters[2] : 1.0;
 979       bias=(number_parameters >= 4) ? parameters[3] : 0.5;
 980       result=2.0/width*(QuantumScale*pixel-center);
 981       if ( result <= -1.0 )
 982         result=bias-range/2.0;
 983       else
 984         if (result >= 1.0)
 985           result=bias+range/2.0;
 986         else
 987           result=(MagickRealType) (range/MagickPI*asin((double) result)+bias);
 988       result*=QuantumRange;
 989       break;
 990     }
 991     case ArctanFunction:
 992     {
 993       MagickRealType
 994         center,
 995         bias,
 996         range,
 997         slope;
 998
 999       /*
1000         Arctan: slope, center, range, and bias.
1001       */
1002       slope=(number_parameters >= 1) ? parameters[0] : 1.0;
1003       center=(number_parameters >= 2) ? parameters[1] : 0.5;
1004       range=(number_parameters >= 3) ? parameters[2] : 1.0;
1005       bias=(number_parameters >= 4) ? parameters[3] : 0.5;
1006       result=(MagickRealType) (MagickPI*slope*(QuantumScale*pixel-center));
1007       result=(MagickRealType) (QuantumRange*(range/MagickPI*atan((double)
1008         result)+bias));
1009       break;
1010     }
1011     case UndefinedFunction:
1012       break;
1013   }
1014   return(ClampToQuantum(result));
1015 }
1016
1017 MagickExport MagickBooleanType FunctionImage(Image *image,
1018   const MagickFunction function,const size_t number_parameters,
1019   const double *parameters,ExceptionInfo *exception)
1020 {
1021 #define FunctionImageTag  "Function/Image "
1022
1023   CacheView
1024     *image_view;
1025
1026   MagickBooleanType
1027     status;
1028
1029   MagickOffsetType
1030     progress;
1031
1032   ssize_t
1033     y;
1034
1035   assert(image != (Image *) NULL);
1036   assert(image->signature == MagickSignature);
1037   if (image->debug != MagickFalse)
1038     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1039   assert(exception != (ExceptionInfo *) NULL);
1040   assert(exception->signature == MagickSignature);
1041   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
1042     return(MagickFalse);
1043   status=MagickTrue;
1044   progress=0;
1045   image_view=AcquireCacheView(image);
1046 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1047   #pragma omp parallel for schedule(static,4) shared(progress,status)
1048 #endif
1049   for (y=0; y < (ssize_t) image->rows; y++)
1050   {
1051     register Quantum
1052       *restrict q;
1053
1054     register ssize_t
1055       x;
1056
1057     if (status == MagickFalse)
1058       continue;
1059     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
1060     if (q == (Quantum *) NULL)
1061       {
1062         status=MagickFalse;
1063         continue;
1064       }
1065     for (x=0; x < (ssize_t) image->columns; x++)
1066     {
1067       register ssize_t
1068         i;
1069
1070       if (GetPixelMask(image,q) != 0)
1071         {
1072           q+=GetPixelChannels(image);
1073           continue;
1074         }
1075       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1076       {
1077         PixelChannel
1078           channel;
1079
1080         PixelTrait
1081           traits;
1082
1083         channel=GetPixelChannelMapChannel(image,i);
1084         traits=GetPixelChannelMapTraits(image,channel);
1085         if (traits == UndefinedPixelTrait)
1086           continue;
1087         if ((traits & UpdatePixelTrait) == 0)
1088           continue;
1089         q[i]=ApplyFunction(q[i],function,number_parameters,parameters,
1090           exception);
1091       }
1092       q+=GetPixelChannels(image);
1093     }
1094     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1095       status=MagickFalse;
1096     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1097       {
1098         MagickBooleanType
1099           proceed;
1100
1101 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1102   #pragma omp critical (MagickCore_FunctionImage)
1103 #endif
1104         proceed=SetImageProgress(image,FunctionImageTag,progress++,image->rows);
1105         if (proceed == MagickFalse)
1106           status=MagickFalse;
1107       }
1108   }
1109   image_view=DestroyCacheView(image_view);
1110   return(status);
1111 }
1112 \f
1113 /*
1114 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1115 %                                                                             %
1116 %                                                                             %
1117 %                                                                             %
1118 %   G e t I m a g e E x t r e m a                                             %
1119 %                                                                             %
1120 %                                                                             %
1121 %                                                                             %
1122 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1123 %
1124 %  GetImageExtrema() returns the extrema of one or more image channels.
1125 %
1126 %  The format of the GetImageExtrema method is:
1127 %
1128 %      MagickBooleanType GetImageExtrema(const Image *image,size_t *minima,
1129 %        size_t *maxima,ExceptionInfo *exception)
1130 %
1131 %  A description of each parameter follows:
1132 %
1133 %    o image: the image.
1134 %
1135 %    o minima: the minimum value in the channel.
1136 %
1137 %    o maxima: the maximum value in the channel.
1138 %
1139 %    o exception: return any errors or warnings in this structure.
1140 %
1141 */
1142 MagickExport MagickBooleanType GetImageExtrema(const Image *image,
1143   size_t *minima,size_t *maxima,ExceptionInfo *exception)
1144 {
1145   double
1146     max,
1147     min;
1148
1149   MagickBooleanType
1150     status;
1151
1152   assert(image != (Image *) NULL);
1153   assert(image->signature == MagickSignature);
1154   if (image->debug != MagickFalse)
1155     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1156   status=GetImageRange(image,&min,&max,exception);
1157   *minima=(size_t) ceil(min-0.5);
1158   *maxima=(size_t) floor(max+0.5);
1159   return(status);
1160 }
1161 \f
1162 /*
1163 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1164 %                                                                             %
1165 %                                                                             %
1166 %                                                                             %
1167 %   G e t I m a g e M e a n                                                   %
1168 %                                                                             %
1169 %                                                                             %
1170 %                                                                             %
1171 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1172 %
1173 %  GetImageMean() returns the mean and standard deviation of one or more
1174 %  image channels.
1175 %
1176 %  The format of the GetImageMean method is:
1177 %
1178 %      MagickBooleanType GetImageMean(const Image *image,double *mean,
1179 %        double *standard_deviation,ExceptionInfo *exception)
1180 %
1181 %  A description of each parameter follows:
1182 %
1183 %    o image: the image.
1184 %
1185 %    o mean: the average value in the channel.
1186 %
1187 %    o standard_deviation: the standard deviation of the channel.
1188 %
1189 %    o exception: return any errors or warnings in this structure.
1190 %
1191 */
1192 MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean,
1193   double *standard_deviation,ExceptionInfo *exception)
1194 {
1195   ChannelStatistics
1196     *channel_statistics;
1197
1198   register ssize_t
1199     i;
1200
1201   size_t
1202     area;
1203
1204   assert(image != (Image *) NULL);
1205   assert(image->signature == MagickSignature);
1206   if (image->debug != MagickFalse)
1207     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1208   channel_statistics=GetImageStatistics(image,exception);
1209   if (channel_statistics == (ChannelStatistics *) NULL)
1210     return(MagickFalse);
1211   area=0;
1212   channel_statistics[CompositePixelChannel].mean=0.0;
1213   channel_statistics[CompositePixelChannel].standard_deviation=0.0;
1214   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1215   {
1216     PixelChannel
1217       channel;
1218
1219     PixelTrait
1220       traits;
1221
1222     channel=GetPixelChannelMapChannel(image,i);
1223     traits=GetPixelChannelMapTraits(image,channel);
1224     if (traits == UndefinedPixelTrait)
1225       continue;
1226     if ((traits & UpdatePixelTrait) == 0)
1227       continue;
1228     channel_statistics[CompositePixelChannel].mean+=channel_statistics[i].mean;
1229     channel_statistics[CompositePixelChannel].standard_deviation+=
1230       channel_statistics[i].variance-channel_statistics[i].mean*
1231       channel_statistics[i].mean;
1232     area++;
1233   }
1234   channel_statistics[CompositePixelChannel].mean/=area;
1235   channel_statistics[CompositePixelChannel].standard_deviation=
1236     sqrt(channel_statistics[CompositePixelChannel].standard_deviation/area);
1237   *mean=channel_statistics[CompositePixelChannel].mean;
1238   *standard_deviation=
1239     channel_statistics[CompositePixelChannel].standard_deviation;
1240   channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1241     channel_statistics);
1242   return(MagickTrue);
1243 }
1244 \f
1245 /*
1246 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1247 %                                                                             %
1248 %                                                                             %
1249 %                                                                             %
1250 %   G e t I m a g e K u r t o s i s                                           %
1251 %                                                                             %
1252 %                                                                             %
1253 %                                                                             %
1254 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1255 %
1256 %  GetImageKurtosis() returns the kurtosis and skewness of one or more
1257 %  image channels.
1258 %
1259 %  The format of the GetImageKurtosis method is:
1260 %
1261 %      MagickBooleanType GetImageKurtosis(const Image *image,double *kurtosis,
1262 %        double *skewness,ExceptionInfo *exception)
1263 %
1264 %  A description of each parameter follows:
1265 %
1266 %    o image: the image.
1267 %
1268 %    o kurtosis: the kurtosis of the channel.
1269 %
1270 %    o skewness: the skewness of the channel.
1271 %
1272 %    o exception: return any errors or warnings in this structure.
1273 %
1274 */
1275 MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
1276   double *kurtosis,double *skewness,ExceptionInfo *exception)
1277 {
1278   CacheView
1279     *image_view;
1280
1281   double
1282     area,
1283     mean,
1284     standard_deviation,
1285     sum_squares,
1286     sum_cubes,
1287     sum_fourth_power;
1288
1289   MagickBooleanType
1290     status;
1291
1292   ssize_t
1293     y;
1294
1295   assert(image != (Image *) NULL);
1296   assert(image->signature == MagickSignature);
1297   if (image->debug != MagickFalse)
1298     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1299   status=MagickTrue;
1300   *kurtosis=0.0;
1301   *skewness=0.0;
1302   area=0.0;
1303   mean=0.0;
1304   standard_deviation=0.0;
1305   sum_squares=0.0;
1306   sum_cubes=0.0;
1307   sum_fourth_power=0.0;
1308   image_view=AcquireCacheView(image);
1309 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1310   #pragma omp parallel for schedule(static) shared(status)
1311 #endif
1312   for (y=0; y < (ssize_t) image->rows; y++)
1313   {
1314     register const Quantum
1315       *restrict p;
1316
1317     register ssize_t
1318       x;
1319
1320     if (status == MagickFalse)
1321       continue;
1322     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1323     if (p == (const Quantum *) NULL)
1324       {
1325         status=MagickFalse;
1326         continue;
1327       }
1328     for (x=0; x < (ssize_t) image->columns; x++)
1329     {
1330       register ssize_t
1331         i;
1332
1333       if (GetPixelMask(image,p) != 0)
1334         {
1335           p+=GetPixelChannels(image);
1336           continue;
1337         }
1338       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1339       {
1340         PixelChannel
1341           channel;
1342
1343         PixelTrait
1344           traits;
1345
1346         channel=GetPixelChannelMapChannel(image,i);
1347         traits=GetPixelChannelMapTraits(image,channel);
1348         if (traits == UndefinedPixelTrait)
1349           continue;
1350         if ((traits & UpdatePixelTrait) == 0)
1351           continue;
1352 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1353         #pragma omp critical (MagickCore_GetImageKurtosis)
1354 #endif
1355         {
1356           mean+=p[i];
1357           sum_squares+=(double) p[i]*p[i];
1358           sum_cubes+=(double) p[i]*p[i]*p[i];
1359           sum_fourth_power+=(double) p[i]*p[i]*p[i]*p[i];
1360           area++;
1361         }
1362       }
1363       p+=GetPixelChannels(image);
1364     }
1365   }
1366   image_view=DestroyCacheView(image_view);
1367   if (area != 0.0)
1368     {
1369       mean/=area;
1370       sum_squares/=area;
1371       sum_cubes/=area;
1372       sum_fourth_power/=area;
1373     }
1374   standard_deviation=sqrt(sum_squares-(mean*mean));
1375   if (standard_deviation != 0.0)
1376     {
1377       *kurtosis=sum_fourth_power-4.0*mean*sum_cubes+6.0*mean*mean*sum_squares-
1378         3.0*mean*mean*mean*mean;
1379       *kurtosis/=standard_deviation*standard_deviation*standard_deviation*
1380         standard_deviation;
1381       *kurtosis-=3.0;
1382       *skewness=sum_cubes-3.0*mean*sum_squares+2.0*mean*mean*mean;
1383       *skewness/=standard_deviation*standard_deviation*standard_deviation;
1384     }
1385   return(status);
1386 }
1387 \f
1388 /*
1389 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1390 %                                                                             %
1391 %                                                                             %
1392 %                                                                             %
1393 %   G e t I m a g e R a n g e                                                 %
1394 %                                                                             %
1395 %                                                                             %
1396 %                                                                             %
1397 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1398 %
1399 %  GetImageRange() returns the range of one or more image channels.
1400 %
1401 %  The format of the GetImageRange method is:
1402 %
1403 %      MagickBooleanType GetImageRange(const Image *image,double *minima,
1404 %        double *maxima,ExceptionInfo *exception)
1405 %
1406 %  A description of each parameter follows:
1407 %
1408 %    o image: the image.
1409 %
1410 %    o minima: the minimum value in the channel.
1411 %
1412 %    o maxima: the maximum value in the channel.
1413 %
1414 %    o exception: return any errors or warnings in this structure.
1415 %
1416 */
1417 MagickExport MagickBooleanType GetImageRange(const Image *image,double *minima,
1418   double *maxima,ExceptionInfo *exception)
1419 {
1420   CacheView
1421     *image_view;
1422
1423   MagickBooleanType
1424     status;
1425
1426   ssize_t
1427     y;
1428
1429   assert(image != (Image *) NULL);
1430   assert(image->signature == MagickSignature);
1431   if (image->debug != MagickFalse)
1432     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1433   status=MagickTrue;
1434   *maxima=(-MagickHuge);
1435   *minima=MagickHuge;
1436   image_view=AcquireCacheView(image);
1437 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1438   #pragma omp parallel for schedule(static) shared(status)
1439 #endif
1440   for (y=0; y < (ssize_t) image->rows; y++)
1441   {
1442     register const Quantum
1443       *restrict p;
1444
1445     register ssize_t
1446       x;
1447
1448     if (status == MagickFalse)
1449       continue;
1450     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1451     if (p == (const Quantum *) NULL)
1452       {
1453         status=MagickFalse;
1454         continue;
1455       }
1456     for (x=0; x < (ssize_t) image->columns; x++)
1457     {
1458       register ssize_t
1459         i;
1460
1461       if (GetPixelMask(image,p) != 0)
1462         {
1463           p+=GetPixelChannels(image);
1464           continue;
1465         }
1466       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1467       {
1468         PixelChannel
1469           channel;
1470
1471         PixelTrait
1472           traits;
1473
1474         channel=GetPixelChannelMapChannel(image,i);
1475         traits=GetPixelChannelMapTraits(image,channel);
1476         if (traits == UndefinedPixelTrait)
1477           continue;
1478         if ((traits & UpdatePixelTrait) == 0)
1479           continue;
1480 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1481         #pragma omp critical (MagickCore_GetImageRange)
1482 #endif
1483         {
1484           if ((double) p[i] < *minima)
1485             *minima=(double) p[i];
1486           if ((double) p[i] > *maxima)
1487             *maxima=(double) p[i];
1488         }
1489       }
1490       p+=GetPixelChannels(image);
1491     }
1492   }
1493   image_view=DestroyCacheView(image_view);
1494   return(status);
1495 }
1496 \f
1497 /*
1498 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1499 %                                                                             %
1500 %                                                                             %
1501 %                                                                             %
1502 %   G e t I m a g e S t a t i s t i c s                                       %
1503 %                                                                             %
1504 %                                                                             %
1505 %                                                                             %
1506 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1507 %
1508 %  GetImageStatistics() returns statistics for each channel in the
1509 %  image.  The statistics include the channel depth, its minima, maxima, mean,
1510 %  standard deviation, kurtosis and skewness.  You can access the red channel
1511 %  mean, for example, like this:
1512 %
1513 %      channel_statistics=GetImageStatistics(image,exception);
1514 %      red_mean=channel_statistics[RedPixelChannel].mean;
1515 %
1516 %  Use MagickRelinquishMemory() to free the statistics buffer.
1517 %
1518 %  The format of the GetImageStatistics method is:
1519 %
1520 %      ChannelStatistics *GetImageStatistics(const Image *image,
1521 %        ExceptionInfo *exception)
1522 %
1523 %  A description of each parameter follows:
1524 %
1525 %    o image: the image.
1526 %
1527 %    o exception: return any errors or warnings in this structure.
1528 %
1529 */
1530
1531 static size_t GetImageChannels(const Image *image)
1532 {
1533   register ssize_t
1534     i;
1535
1536   size_t
1537     channels;
1538
1539   channels=0;
1540   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1541   {
1542     PixelChannel
1543       channel;
1544
1545     PixelTrait
1546       traits;
1547
1548     channel=GetPixelChannelMapChannel(image,i);
1549     traits=GetPixelChannelMapTraits(image,channel);
1550     if ((traits & UpdatePixelTrait) != 0)
1551       channels++;
1552   }
1553   return(channels);
1554 }
1555
1556 MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
1557   ExceptionInfo *exception)
1558 {
1559   ChannelStatistics
1560     *channel_statistics;
1561
1562   double
1563     area;
1564
1565   MagickStatusType
1566     status;
1567
1568   QuantumAny
1569     range;
1570
1571   register ssize_t
1572     i;
1573
1574   size_t
1575     channels,
1576     depth;
1577
1578   ssize_t
1579     y;
1580
1581   assert(image != (Image *) NULL);
1582   assert(image->signature == MagickSignature);
1583   if (image->debug != MagickFalse)
1584     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1585   channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(
1586     MaxPixelChannels+1,sizeof(*channel_statistics));
1587   if (channel_statistics == (ChannelStatistics *) NULL)
1588     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1589   (void) ResetMagickMemory(channel_statistics,0,(MaxPixelChannels+1)*
1590     sizeof(*channel_statistics));
1591   for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
1592   {
1593     channel_statistics[i].depth=1;
1594     channel_statistics[i].maxima=(-MagickHuge);
1595     channel_statistics[i].minima=MagickHuge;
1596   }
1597   for (y=0; y < (ssize_t) image->rows; y++)
1598   {
1599     register const Quantum
1600       *restrict p;
1601
1602     register ssize_t
1603       x;
1604
1605     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1606     if (p == (const Quantum *) NULL)
1607       break;
1608     for (x=0; x < (ssize_t) image->columns; x++)
1609     {
1610       register ssize_t
1611         i;
1612
1613       if (GetPixelMask(image,p) != 0)
1614         {
1615           p+=GetPixelChannels(image);
1616           continue;
1617         }
1618       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1619       {
1620         PixelChannel
1621           channel;
1622
1623         PixelTrait
1624           traits;
1625
1626         channel=GetPixelChannelMapChannel(image,i);
1627         traits=GetPixelChannelMapTraits(image,channel);
1628         if (traits == UndefinedPixelTrait)
1629           continue;
1630         if (channel_statistics[channel].depth != MAGICKCORE_QUANTUM_DEPTH)
1631           {
1632             depth=channel_statistics[channel].depth;
1633             range=GetQuantumRange(depth);
1634             status=p[i] != ScaleAnyToQuantum(ScaleQuantumToAny(p[i],range),
1635               range) ? MagickTrue : MagickFalse;
1636             if (status != MagickFalse)
1637               {
1638                 channel_statistics[channel].depth++;
1639                 continue;
1640               }
1641           }
1642         if ((double) p[i] < channel_statistics[channel].minima)
1643           channel_statistics[channel].minima=(double) p[i];
1644         if ((double) p[i] > channel_statistics[channel].maxima)
1645           channel_statistics[channel].maxima=(double) p[i];
1646         channel_statistics[channel].sum+=p[i];
1647         channel_statistics[channel].sum_squared+=(double) p[i]*p[i];
1648         channel_statistics[channel].sum_cubed+=(double) p[i]*p[i]*p[i];
1649         channel_statistics[channel].sum_fourth_power+=(double) p[i]*p[i]*p[i]*
1650           p[i];
1651       }
1652       p+=GetPixelChannels(image);
1653     }
1654   }
1655   area=(double) image->columns*image->rows;
1656   for (i=0; i < (ssize_t) MaxPixelChannels; i++)
1657   {
1658     channel_statistics[i].sum/=area;
1659     channel_statistics[i].sum_squared/=area;
1660     channel_statistics[i].sum_cubed/=area;
1661     channel_statistics[i].sum_fourth_power/=area;
1662     channel_statistics[i].mean=channel_statistics[i].sum;
1663     channel_statistics[i].variance=channel_statistics[i].sum_squared;
1664     channel_statistics[i].standard_deviation=sqrt(
1665       channel_statistics[i].variance-(channel_statistics[i].mean*
1666       channel_statistics[i].mean));
1667   }
1668   for (i=0; i < (ssize_t) MaxPixelChannels; i++)
1669   {
1670     channel_statistics[CompositePixelChannel].depth=(size_t) EvaluateMax(
1671       (double) channel_statistics[CompositePixelChannel].depth,(double)
1672       channel_statistics[i].depth);
1673     channel_statistics[CompositePixelChannel].minima=MagickMin(
1674       channel_statistics[CompositePixelChannel].minima,
1675       channel_statistics[i].minima);
1676     channel_statistics[CompositePixelChannel].maxima=EvaluateMax(
1677       channel_statistics[CompositePixelChannel].maxima,
1678       channel_statistics[i].maxima);
1679     channel_statistics[CompositePixelChannel].sum+=channel_statistics[i].sum;
1680     channel_statistics[CompositePixelChannel].sum_squared+=
1681       channel_statistics[i].sum_squared;
1682     channel_statistics[CompositePixelChannel].sum_cubed+=
1683       channel_statistics[i].sum_cubed;
1684     channel_statistics[CompositePixelChannel].sum_fourth_power+=
1685       channel_statistics[i].sum_fourth_power;
1686     channel_statistics[CompositePixelChannel].mean+=channel_statistics[i].mean;
1687     channel_statistics[CompositePixelChannel].variance+=
1688       channel_statistics[i].variance-channel_statistics[i].mean*
1689       channel_statistics[i].mean;
1690     channel_statistics[CompositePixelChannel].standard_deviation+=
1691       channel_statistics[i].variance-channel_statistics[i].mean*
1692       channel_statistics[i].mean;
1693   }
1694   channels=GetImageChannels(image);
1695   channel_statistics[CompositePixelChannel].sum/=channels;
1696   channel_statistics[CompositePixelChannel].sum_squared/=channels;
1697   channel_statistics[CompositePixelChannel].sum_cubed/=channels;
1698   channel_statistics[CompositePixelChannel].sum_fourth_power/=channels;
1699   channel_statistics[CompositePixelChannel].mean/=channels;
1700   channel_statistics[CompositePixelChannel].variance/=channels;
1701   channel_statistics[CompositePixelChannel].standard_deviation=
1702     sqrt(channel_statistics[CompositePixelChannel].standard_deviation/channels);
1703   channel_statistics[CompositePixelChannel].kurtosis/=channels;
1704   channel_statistics[CompositePixelChannel].skewness/=channels;
1705   for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
1706   {
1707     if (channel_statistics[i].standard_deviation == 0.0)
1708       continue;
1709     channel_statistics[i].skewness=(channel_statistics[i].sum_cubed-3.0*
1710       channel_statistics[i].mean*channel_statistics[i].sum_squared+2.0*
1711       channel_statistics[i].mean*channel_statistics[i].mean*
1712       channel_statistics[i].mean)/(channel_statistics[i].standard_deviation*
1713       channel_statistics[i].standard_deviation*
1714       channel_statistics[i].standard_deviation);
1715     channel_statistics[i].kurtosis=(channel_statistics[i].sum_fourth_power-4.0*
1716       channel_statistics[i].mean*channel_statistics[i].sum_cubed+6.0*
1717       channel_statistics[i].mean*channel_statistics[i].mean*
1718       channel_statistics[i].sum_squared-3.0*channel_statistics[i].mean*
1719       channel_statistics[i].mean*1.0*channel_statistics[i].mean*
1720       channel_statistics[i].mean)/(channel_statistics[i].standard_deviation*
1721       channel_statistics[i].standard_deviation*
1722       channel_statistics[i].standard_deviation*
1723       channel_statistics[i].standard_deviation)-3.0;
1724   }
1725   return(channel_statistics);
1726 }
1727 \f
1728 /*
1729 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1730 %                                                                             %
1731 %                                                                             %
1732 %                                                                             %
1733 %     S t a t i s t i c I m a g e                                             %
1734 %                                                                             %
1735 %                                                                             %
1736 %                                                                             %
1737 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1738 %
1739 %  StatisticImage() makes each pixel the min / max / median / mode / etc. of
1740 %  the neighborhood of the specified width and height.
1741 %
1742 %  The format of the StatisticImage method is:
1743 %
1744 %      Image *StatisticImage(const Image *image,const StatisticType type,
1745 %        const size_t width,const size_t height,ExceptionInfo *exception)
1746 %
1747 %  A description of each parameter follows:
1748 %
1749 %    o image: the image.
1750 %
1751 %    o type: the statistic type (median, mode, etc.).
1752 %
1753 %    o width: the width of the pixel neighborhood.
1754 %
1755 %    o height: the height of the pixel neighborhood.
1756 %
1757 %    o exception: return any errors or warnings in this structure.
1758 %
1759 */
1760
1761 typedef struct _SkipNode
1762 {
1763   size_t
1764     next[9],
1765     count,
1766     signature;
1767 } SkipNode;
1768
1769 typedef struct _SkipList
1770 {
1771   ssize_t
1772     level;
1773
1774   SkipNode
1775     *nodes;
1776 } SkipList;
1777
1778 typedef struct _PixelList
1779 {
1780   size_t
1781     length,
1782     seed;
1783
1784   SkipList
1785     skip_list;
1786
1787   size_t
1788     signature;
1789 } PixelList;
1790
1791 static PixelList *DestroyPixelList(PixelList *pixel_list)
1792 {
1793   if (pixel_list == (PixelList *) NULL)
1794     return((PixelList *) NULL);
1795   if (pixel_list->skip_list.nodes != (SkipNode *) NULL)
1796     pixel_list->skip_list.nodes=(SkipNode *) RelinquishMagickMemory(
1797       pixel_list->skip_list.nodes);
1798   pixel_list=(PixelList *) RelinquishMagickMemory(pixel_list);
1799   return(pixel_list);
1800 }
1801
1802 static PixelList **DestroyPixelListThreadSet(PixelList **pixel_list)
1803 {
1804   register ssize_t
1805     i;
1806
1807   assert(pixel_list != (PixelList **) NULL);
1808   for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
1809     if (pixel_list[i] != (PixelList *) NULL)
1810       pixel_list[i]=DestroyPixelList(pixel_list[i]);
1811   pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list);
1812   return(pixel_list);
1813 }
1814
1815 static PixelList *AcquirePixelList(const size_t width,const size_t height)
1816 {
1817   PixelList
1818     *pixel_list;
1819
1820   pixel_list=(PixelList *) AcquireMagickMemory(sizeof(*pixel_list));
1821   if (pixel_list == (PixelList *) NULL)
1822     return(pixel_list);
1823   (void) ResetMagickMemory((void *) pixel_list,0,sizeof(*pixel_list));
1824   pixel_list->length=width*height;
1825   pixel_list->skip_list.nodes=(SkipNode *) AcquireQuantumMemory(65537UL,
1826     sizeof(*pixel_list->skip_list.nodes));
1827   if (pixel_list->skip_list.nodes == (SkipNode *) NULL)
1828     return(DestroyPixelList(pixel_list));
1829   (void) ResetMagickMemory(pixel_list->skip_list.nodes,0,65537UL*
1830     sizeof(*pixel_list->skip_list.nodes));
1831   pixel_list->signature=MagickSignature;
1832   return(pixel_list);
1833 }
1834
1835 static PixelList **AcquirePixelListThreadSet(const size_t width,
1836   const size_t height)
1837 {
1838   PixelList
1839     **pixel_list;
1840
1841   register ssize_t
1842     i;
1843
1844   size_t
1845     number_threads;
1846
1847   number_threads=GetOpenMPMaximumThreads();
1848   pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
1849     sizeof(*pixel_list));
1850   if (pixel_list == (PixelList **) NULL)
1851     return((PixelList **) NULL);
1852   (void) ResetMagickMemory(pixel_list,0,number_threads*sizeof(*pixel_list));
1853   for (i=0; i < (ssize_t) number_threads; i++)
1854   {
1855     pixel_list[i]=AcquirePixelList(width,height);
1856     if (pixel_list[i] == (PixelList *) NULL)
1857       return(DestroyPixelListThreadSet(pixel_list));
1858   }
1859   return(pixel_list);
1860 }
1861
1862 static void AddNodePixelList(PixelList *pixel_list,const size_t color)
1863 {
1864   register SkipList
1865     *p;
1866
1867   register ssize_t
1868     level;
1869
1870   size_t
1871     search,
1872     update[9];
1873
1874   /*
1875     Initialize the node.
1876   */
1877   p=(&pixel_list->skip_list);
1878   p->nodes[color].signature=pixel_list->signature;
1879   p->nodes[color].count=1;
1880   /*
1881     Determine where it belongs in the list.
1882   */
1883   search=65536UL;
1884   for (level=p->level; level >= 0; level--)
1885   {
1886     while (p->nodes[search].next[level] < color)
1887       search=p->nodes[search].next[level];
1888     update[level]=search;
1889   }
1890   /*
1891     Generate a pseudo-random level for this node.
1892   */
1893   for (level=0; ; level++)
1894   {
1895     pixel_list->seed=(pixel_list->seed*42893621L)+1L;
1896     if ((pixel_list->seed & 0x300) != 0x300)
1897       break;
1898   }
1899   if (level > 8)
1900     level=8;
1901   if (level > (p->level+2))
1902     level=p->level+2;
1903   /*
1904     If we're raising the list's level, link back to the root node.
1905   */
1906   while (level > p->level)
1907   {
1908     p->level++;
1909     update[p->level]=65536UL;
1910   }
1911   /*
1912     Link the node into the skip-list.
1913   */
1914   do
1915   {
1916     p->nodes[color].next[level]=p->nodes[update[level]].next[level];
1917     p->nodes[update[level]].next[level]=color;
1918   } while (level-- > 0);
1919 }
1920
1921 static inline void GetMaximumPixelList(PixelList *pixel_list,Quantum *pixel)
1922 {
1923   register SkipList
1924     *p;
1925
1926   size_t
1927     color,
1928     maximum;
1929
1930   ssize_t
1931     count;
1932
1933   /*
1934     Find the maximum value for each of the color.
1935   */
1936   p=(&pixel_list->skip_list);
1937   color=65536L;
1938   count=0;
1939   maximum=p->nodes[color].next[0];
1940   do
1941   {
1942     color=p->nodes[color].next[0];
1943     if (color > maximum)
1944       maximum=color;
1945     count+=p->nodes[color].count;
1946   } while (count < (ssize_t) pixel_list->length);
1947   *pixel=ScaleShortToQuantum((unsigned short) maximum);
1948 }
1949
1950 static inline void GetMeanPixelList(PixelList *pixel_list,Quantum *pixel)
1951 {
1952   MagickRealType
1953     sum;
1954
1955   register SkipList
1956     *p;
1957
1958   size_t
1959     color;
1960
1961   ssize_t
1962     count;
1963
1964   /*
1965     Find the mean value for each of the color.
1966   */
1967   p=(&pixel_list->skip_list);
1968   color=65536L;
1969   count=0;
1970   sum=0.0;
1971   do
1972   {
1973     color=p->nodes[color].next[0];
1974     sum+=(MagickRealType) p->nodes[color].count*color;
1975     count+=p->nodes[color].count;
1976   } while (count < (ssize_t) pixel_list->length);
1977   sum/=pixel_list->length;
1978   *pixel=ScaleShortToQuantum((unsigned short) sum);
1979 }
1980
1981 static inline void GetMedianPixelList(PixelList *pixel_list,Quantum *pixel)
1982 {
1983   register SkipList
1984     *p;
1985
1986   size_t
1987     color;
1988
1989   ssize_t
1990     count;
1991
1992   /*
1993     Find the median value for each of the color.
1994   */
1995   p=(&pixel_list->skip_list);
1996   color=65536L;
1997   count=0;
1998   do
1999   {
2000     color=p->nodes[color].next[0];
2001     count+=p->nodes[color].count;
2002   } while (count <= (ssize_t) (pixel_list->length >> 1));
2003   *pixel=ScaleShortToQuantum((unsigned short) color);
2004 }
2005
2006 static inline void GetMinimumPixelList(PixelList *pixel_list,Quantum *pixel)
2007 {
2008   register SkipList
2009     *p;
2010
2011   size_t
2012     color,
2013     minimum;
2014
2015   ssize_t
2016     count;
2017
2018   /*
2019     Find the minimum value for each of the color.
2020   */
2021   p=(&pixel_list->skip_list);
2022   count=0;
2023   color=65536UL;
2024   minimum=p->nodes[color].next[0];
2025   do
2026   {
2027     color=p->nodes[color].next[0];
2028     if (color < minimum)
2029       minimum=color;
2030     count+=p->nodes[color].count;
2031   } while (count < (ssize_t) pixel_list->length);
2032   *pixel=ScaleShortToQuantum((unsigned short) minimum);
2033 }
2034
2035 static inline void GetModePixelList(PixelList *pixel_list,Quantum *pixel)
2036 {
2037   register SkipList
2038     *p;
2039
2040   size_t
2041     color,
2042     max_count,
2043     mode;
2044
2045   ssize_t
2046     count;
2047
2048   /*
2049     Make each pixel the 'predominant color' of the specified neighborhood.
2050   */
2051   p=(&pixel_list->skip_list);
2052   color=65536L;
2053   mode=color;
2054   max_count=p->nodes[mode].count;
2055   count=0;
2056   do
2057   {
2058     color=p->nodes[color].next[0];
2059     if (p->nodes[color].count > max_count)
2060       {
2061         mode=color;
2062         max_count=p->nodes[mode].count;
2063       }
2064     count+=p->nodes[color].count;
2065   } while (count < (ssize_t) pixel_list->length);
2066   *pixel=ScaleShortToQuantum((unsigned short) mode);
2067 }
2068
2069 static inline void GetNonpeakPixelList(PixelList *pixel_list,Quantum *pixel)
2070 {
2071   register SkipList
2072     *p;
2073
2074   size_t
2075     color,
2076     next,
2077     previous;
2078
2079   ssize_t
2080     count;
2081
2082   /*
2083     Finds the non peak value for each of the colors.
2084   */
2085   p=(&pixel_list->skip_list);
2086   color=65536L;
2087   next=p->nodes[color].next[0];
2088   count=0;
2089   do
2090   {
2091     previous=color;
2092     color=next;
2093     next=p->nodes[color].next[0];
2094     count+=p->nodes[color].count;
2095   } while (count <= (ssize_t) (pixel_list->length >> 1));
2096   if ((previous == 65536UL) && (next != 65536UL))
2097     color=next;
2098   else
2099     if ((previous != 65536UL) && (next == 65536UL))
2100       color=previous;
2101   *pixel=ScaleShortToQuantum((unsigned short) color);
2102 }
2103
2104 static inline void GetStandardDeviationPixelList(PixelList *pixel_list,
2105   Quantum *pixel)
2106 {
2107   MagickRealType
2108     sum,
2109     sum_squared;
2110
2111   register SkipList
2112     *p;
2113
2114   size_t
2115     color;
2116
2117   ssize_t
2118     count;
2119
2120   /*
2121     Find the standard-deviation value for each of the color.
2122   */
2123   p=(&pixel_list->skip_list);
2124   color=65536L;
2125   count=0;
2126   sum=0.0;
2127   sum_squared=0.0;
2128   do
2129   {
2130     register ssize_t
2131       i;
2132
2133     color=p->nodes[color].next[0];
2134     sum+=(MagickRealType) p->nodes[color].count*color;
2135     for (i=0; i < (ssize_t) p->nodes[color].count; i++)
2136       sum_squared+=((MagickRealType) color)*((MagickRealType) color);
2137     count+=p->nodes[color].count;
2138   } while (count < (ssize_t) pixel_list->length);
2139   sum/=pixel_list->length;
2140   sum_squared/=pixel_list->length;
2141   *pixel=ScaleShortToQuantum((unsigned short) sqrt(sum_squared-(sum*sum)));
2142 }
2143
2144 static inline void InsertPixelList(const Image *image,const Quantum pixel,
2145   PixelList *pixel_list)
2146 {
2147   size_t
2148     signature;
2149
2150   unsigned short
2151     index;
2152
2153   index=ScaleQuantumToShort(pixel);
2154   signature=pixel_list->skip_list.nodes[index].signature;
2155   if (signature == pixel_list->signature)
2156     {
2157       pixel_list->skip_list.nodes[index].count++;
2158       return;
2159     }
2160   AddNodePixelList(pixel_list,index);
2161 }
2162
2163 static inline MagickRealType MagickAbsoluteValue(const MagickRealType x)
2164 {
2165   if (x < 0)
2166     return(-x);
2167   return(x);
2168 }
2169
2170 static inline size_t MagickMax(const size_t x,const size_t y)
2171 {
2172   if (x > y)
2173     return(x);
2174   return(y);
2175 }
2176
2177 static void ResetPixelList(PixelList *pixel_list)
2178 {
2179   int
2180     level;
2181
2182   register SkipNode
2183     *root;
2184
2185   register SkipList
2186     *p;
2187
2188   /*
2189     Reset the skip-list.
2190   */
2191   p=(&pixel_list->skip_list);
2192   root=p->nodes+65536UL;
2193   p->level=0;
2194   for (level=0; level < 9; level++)
2195     root->next[level]=65536UL;
2196   pixel_list->seed=pixel_list->signature++;
2197 }
2198
2199 MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
2200   const size_t width,const size_t height,ExceptionInfo *exception)
2201 {
2202 #define StatisticImageTag  "Statistic/Image"
2203
2204   CacheView
2205     *image_view,
2206     *statistic_view;
2207
2208   Image
2209     *statistic_image;
2210
2211   MagickBooleanType
2212     status;
2213
2214   MagickOffsetType
2215     progress;
2216
2217   PixelList
2218     **restrict pixel_list;
2219
2220   ssize_t
2221     center,
2222     y;
2223
2224   /*
2225     Initialize statistics image attributes.
2226   */
2227   assert(image != (Image *) NULL);
2228   assert(image->signature == MagickSignature);
2229   if (image->debug != MagickFalse)
2230     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2231   assert(exception != (ExceptionInfo *) NULL);
2232   assert(exception->signature == MagickSignature);
2233   statistic_image=CloneImage(image,image->columns,image->rows,MagickTrue,
2234     exception);
2235   if (statistic_image == (Image *) NULL)
2236     return((Image *) NULL);
2237   status=SetImageStorageClass(statistic_image,DirectClass,exception);
2238   if (status == MagickFalse)
2239     {
2240       statistic_image=DestroyImage(statistic_image);
2241       return((Image *) NULL);
2242     }
2243   pixel_list=AcquirePixelListThreadSet(MagickMax(width,1),MagickMax(height,1));
2244   if (pixel_list == (PixelList **) NULL)
2245     {
2246       statistic_image=DestroyImage(statistic_image);
2247       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2248     }
2249   /*
2250     Make each pixel the min / max / median / mode / etc. of the neighborhood.
2251   */
2252   center=(ssize_t) GetPixelChannels(image)*(image->columns+MagickMax(width,1))*
2253     (MagickMax(height,1)/2L)+GetPixelChannels(image)*(MagickMax(width,1)/2L);
2254   status=MagickTrue;
2255   progress=0;
2256   image_view=AcquireCacheView(image);
2257   statistic_view=AcquireCacheView(statistic_image);
2258 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2259   #pragma omp parallel for schedule(static,4) shared(progress,status)
2260 #endif
2261   for (y=0; y < (ssize_t) statistic_image->rows; y++)
2262   {
2263     const int
2264       id = GetOpenMPThreadId();
2265
2266     register const Quantum
2267       *restrict p;
2268
2269     register Quantum
2270       *restrict q;
2271
2272     register ssize_t
2273       x;
2274
2275     if (status == MagickFalse)
2276       continue;
2277     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) MagickMax(width,1)/2L),y-
2278       (ssize_t) (MagickMax(height,1)/2L),image->columns+MagickMax(width,1),
2279       MagickMax(height,1),exception);
2280     q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns,      1,exception);
2281     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2282       {
2283         status=MagickFalse;
2284         continue;
2285       }
2286     for (x=0; x < (ssize_t) statistic_image->columns; x++)
2287     {
2288       register ssize_t
2289         i;
2290
2291       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2292       {
2293         PixelChannel
2294           channel;
2295
2296         PixelTrait
2297           statistic_traits,
2298           traits;
2299
2300         Quantum
2301           pixel;
2302
2303         register const Quantum
2304           *restrict pixels;
2305
2306         register ssize_t
2307           u;
2308
2309         ssize_t
2310           v;
2311
2312         channel=GetPixelChannelMapChannel(image,i);
2313         traits=GetPixelChannelMapTraits(image,channel);
2314         statistic_traits=GetPixelChannelMapTraits(statistic_image,channel);
2315         if ((traits == UndefinedPixelTrait) ||
2316             (statistic_traits == UndefinedPixelTrait))
2317           continue;
2318         if (((statistic_traits & CopyPixelTrait) != 0) ||
2319             (GetPixelMask(image,p) != 0))
2320           {
2321             SetPixelChannel(statistic_image,channel,p[center+i],q);
2322             continue;
2323           }
2324         pixels=p;
2325         ResetPixelList(pixel_list[id]);
2326         for (v=0; v < (ssize_t) MagickMax(height,1); v++)
2327         {
2328           for (u=0; u < (ssize_t) MagickMax(width,1); u++)
2329           {
2330             InsertPixelList(image,pixels[i],pixel_list[id]);
2331             pixels+=GetPixelChannels(image);
2332           }
2333           pixels+=image->columns*GetPixelChannels(image);
2334         }
2335         switch (type)
2336         {
2337           case GradientStatistic:
2338           {
2339             MagickRealType
2340               maximum,
2341               minimum;
2342
2343             GetMinimumPixelList(pixel_list[id],&pixel);
2344             minimum=(MagickRealType) pixel;
2345             GetMaximumPixelList(pixel_list[id],&pixel);
2346             maximum=(MagickRealType) pixel;
2347             pixel=ClampToQuantum(MagickAbsoluteValue(maximum-minimum));
2348             break;
2349           }
2350           case MaximumStatistic:
2351           {
2352             GetMaximumPixelList(pixel_list[id],&pixel);
2353             break;
2354           }
2355           case MeanStatistic:
2356           {
2357             GetMeanPixelList(pixel_list[id],&pixel);
2358             break;
2359           }
2360           case MedianStatistic:
2361           default:
2362           {
2363             GetMedianPixelList(pixel_list[id],&pixel);
2364             break;
2365           }
2366           case MinimumStatistic:
2367           {
2368             GetMinimumPixelList(pixel_list[id],&pixel);
2369             break;
2370           }
2371           case ModeStatistic:
2372           {
2373             GetModePixelList(pixel_list[id],&pixel);
2374             break;
2375           }
2376           case NonpeakStatistic:
2377           {
2378             GetNonpeakPixelList(pixel_list[id],&pixel);
2379             break;
2380           }
2381           case StandardDeviationStatistic:
2382           {
2383             GetStandardDeviationPixelList(pixel_list[id],&pixel);
2384             break;
2385           }
2386         }
2387         SetPixelChannel(statistic_image,channel,pixel,q);
2388       }
2389       p+=GetPixelChannels(image);
2390       q+=GetPixelChannels(statistic_image);
2391     }
2392     if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
2393       status=MagickFalse;
2394     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2395       {
2396         MagickBooleanType
2397           proceed;
2398
2399 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2400   #pragma omp critical (MagickCore_StatisticImage)
2401 #endif
2402         proceed=SetImageProgress(image,StatisticImageTag,progress++,
2403           image->rows);
2404         if (proceed == MagickFalse)
2405           status=MagickFalse;
2406       }
2407   }
2408   statistic_view=DestroyCacheView(statistic_view);
2409   image_view=DestroyCacheView(image_view);
2410   pixel_list=DestroyPixelListThreadSet(pixel_list);
2411   return(statistic_image);
2412 }