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

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