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