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