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