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

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