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