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

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