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 == 0 ?
 700                 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 %     P o l y n o m i a l I m a g e                                           %
1756 %                                                                             %
1757 %                                                                             %
1758 %                                                                             %
1759 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1760 %
1761 %  PolynomialImage() returns a new image where each pixel is the sum of the
1762 %  pixels in the image sequence after applying its corresponding terms
1763 %  (coefficient and degree pairs) and a constant.
1764 %
1765 %  The format of the PolynomialImage method is:
1766 %
1767 %      Image *PolynomialImage(const Image *image,const size_t number_terms,
1768 %        const double *terms,ExceptionInfo *exception)
1769 %
1770 %  A description of each parameter follows:
1771 %
1772 %    o image: the image.
1773 %
1774 %    o number_terms: the number of terms in the list.  The actual list length
1775 %      is 2 x number_terms + 1 (the constant).
1776 %
1777 %    o terms: the list of polynomial coefficients and degree pairs and a
1778 %      constant.
1779 %
1780 %    o exception: return any errors or warnings in this structure.
1781 %
1782 */
1783
1784 MagickExport Image *PolynomialImage(const Image *image,
1785   const size_t number_terms,const double *terms,ExceptionInfo *exception)
1786 {
1787   Image
1788     *polynomial_image;
1789
1790   MagickBooleanType
1791     status;
1792
1793   assert(image != (Image *) NULL);
1794   assert(image->signature == MagickSignature);
1795   if (image->debug != MagickFalse)
1796     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1797   assert(terms != (double *) NULL);
1798   assert(exception != (ExceptionInfo *) NULL);
1799   assert(exception->signature == MagickSignature);
1800   polynomial_image=CloneImage(image,image->columns,image->rows,MagickTrue,
1801     exception);
1802   if (polynomial_image == (Image *) NULL)
1803     return((Image *) NULL);
1804   status=SetImageStorageClass(polynomial_image,DirectClass,exception);
1805   if (status == MagickFalse)
1806     {
1807       polynomial_image=DestroyImage(polynomial_image);
1808       return((Image *) NULL);
1809     }
1810   return(polynomial_image);
1811 }
1812 \f
1813 /*
1814 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1815 %                                                                             %
1816 %                                                                             %
1817 %                                                                             %
1818 %     S t a t i s t i c I m a g e                                             %
1819 %                                                                             %
1820 %                                                                             %
1821 %                                                                             %
1822 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1823 %
1824 %  StatisticImage() makes each pixel the min / max / median / mode / etc. of
1825 %  the neighborhood of the specified width and height.
1826 %
1827 %  The format of the StatisticImage method is:
1828 %
1829 %      Image *StatisticImage(const Image *image,const StatisticType type,
1830 %        const size_t width,const size_t height,ExceptionInfo *exception)
1831 %
1832 %  A description of each parameter follows:
1833 %
1834 %    o image: the image.
1835 %
1836 %    o type: the statistic type (median, mode, etc.).
1837 %
1838 %    o width: the width of the pixel neighborhood.
1839 %
1840 %    o height: the height of the pixel neighborhood.
1841 %
1842 %    o exception: return any errors or warnings in this structure.
1843 %
1844 */
1845
1846 typedef struct _SkipNode
1847 {
1848   size_t
1849     next[9],
1850     count,
1851     signature;
1852 } SkipNode;
1853
1854 typedef struct _SkipList
1855 {
1856   ssize_t
1857     level;
1858
1859   SkipNode
1860     *nodes;
1861 } SkipList;
1862
1863 typedef struct _PixelList
1864 {
1865   size_t
1866     length,
1867     seed;
1868
1869   SkipList
1870     skip_list;
1871
1872   size_t
1873     signature;
1874 } PixelList;
1875
1876 static PixelList *DestroyPixelList(PixelList *pixel_list)
1877 {
1878   if (pixel_list == (PixelList *) NULL)
1879     return((PixelList *) NULL);
1880   if (pixel_list->skip_list.nodes != (SkipNode *) NULL)
1881     pixel_list->skip_list.nodes=(SkipNode *) RelinquishMagickMemory(
1882       pixel_list->skip_list.nodes);
1883   pixel_list=(PixelList *) RelinquishMagickMemory(pixel_list);
1884   return(pixel_list);
1885 }
1886
1887 static PixelList **DestroyPixelListThreadSet(PixelList **pixel_list)
1888 {
1889   register ssize_t
1890     i;
1891
1892   assert(pixel_list != (PixelList **) NULL);
1893   for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
1894     if (pixel_list[i] != (PixelList *) NULL)
1895       pixel_list[i]=DestroyPixelList(pixel_list[i]);
1896   pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list);
1897   return(pixel_list);
1898 }
1899
1900 static PixelList *AcquirePixelList(const size_t width,const size_t height)
1901 {
1902   PixelList
1903     *pixel_list;
1904
1905   pixel_list=(PixelList *) AcquireMagickMemory(sizeof(*pixel_list));
1906   if (pixel_list == (PixelList *) NULL)
1907     return(pixel_list);
1908   (void) ResetMagickMemory((void *) pixel_list,0,sizeof(*pixel_list));
1909   pixel_list->length=width*height;
1910   pixel_list->skip_list.nodes=(SkipNode *) AcquireQuantumMemory(65537UL,
1911     sizeof(*pixel_list->skip_list.nodes));
1912   if (pixel_list->skip_list.nodes == (SkipNode *) NULL)
1913     return(DestroyPixelList(pixel_list));
1914   (void) ResetMagickMemory(pixel_list->skip_list.nodes,0,65537UL*
1915     sizeof(*pixel_list->skip_list.nodes));
1916   pixel_list->signature=MagickSignature;
1917   return(pixel_list);
1918 }
1919
1920 static PixelList **AcquirePixelListThreadSet(const size_t width,
1921   const size_t height)
1922 {
1923   PixelList
1924     **pixel_list;
1925
1926   register ssize_t
1927     i;
1928
1929   size_t
1930     number_threads;
1931
1932   number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
1933   pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
1934     sizeof(*pixel_list));
1935   if (pixel_list == (PixelList **) NULL)
1936     return((PixelList **) NULL);
1937   (void) ResetMagickMemory(pixel_list,0,number_threads*sizeof(*pixel_list));
1938   for (i=0; i < (ssize_t) number_threads; i++)
1939   {
1940     pixel_list[i]=AcquirePixelList(width,height);
1941     if (pixel_list[i] == (PixelList *) NULL)
1942       return(DestroyPixelListThreadSet(pixel_list));
1943   }
1944   return(pixel_list);
1945 }
1946
1947 static void AddNodePixelList(PixelList *pixel_list,const size_t color)
1948 {
1949   register SkipList
1950     *p;
1951
1952   register ssize_t
1953     level;
1954
1955   size_t
1956     search,
1957     update[9];
1958
1959   /*
1960     Initialize the node.
1961   */
1962   p=(&pixel_list->skip_list);
1963   p->nodes[color].signature=pixel_list->signature;
1964   p->nodes[color].count=1;
1965   /*
1966     Determine where it belongs in the list.
1967   */
1968   search=65536UL;
1969   for (level=p->level; level >= 0; level--)
1970   {
1971     while (p->nodes[search].next[level] < color)
1972       search=p->nodes[search].next[level];
1973     update[level]=search;
1974   }
1975   /*
1976     Generate a pseudo-random level for this node.
1977   */
1978   for (level=0; ; level++)
1979   {
1980     pixel_list->seed=(pixel_list->seed*42893621L)+1L;
1981     if ((pixel_list->seed & 0x300) != 0x300)
1982       break;
1983   }
1984   if (level > 8)
1985     level=8;
1986   if (level > (p->level+2))
1987     level=p->level+2;
1988   /*
1989     If we're raising the list's level, link back to the root node.
1990   */
1991   while (level > p->level)
1992   {
1993     p->level++;
1994     update[p->level]=65536UL;
1995   }
1996   /*
1997     Link the node into the skip-list.
1998   */
1999   do
2000   {
2001     p->nodes[color].next[level]=p->nodes[update[level]].next[level];
2002     p->nodes[update[level]].next[level]=color;
2003   } while (level-- > 0);
2004 }
2005
2006 static inline void GetMaximumPixelList(PixelList *pixel_list,Quantum *pixel)
2007 {
2008   register SkipList
2009     *p;
2010
2011   size_t
2012     color,
2013     maximum;
2014
2015   ssize_t
2016     count;
2017
2018   /*
2019     Find the maximum value for each of the color.
2020   */
2021   p=(&pixel_list->skip_list);
2022   color=65536L;
2023   count=0;
2024   maximum=p->nodes[color].next[0];
2025   do
2026   {
2027     color=p->nodes[color].next[0];
2028     if (color > maximum)
2029       maximum=color;
2030     count+=p->nodes[color].count;
2031   } while (count < (ssize_t) pixel_list->length);
2032   *pixel=ScaleShortToQuantum((unsigned short) maximum);
2033 }
2034
2035 static inline void GetMeanPixelList(PixelList *pixel_list,Quantum *pixel)
2036 {
2037   double
2038     sum;
2039
2040   register SkipList
2041     *p;
2042
2043   size_t
2044     color;
2045
2046   ssize_t
2047     count;
2048
2049   /*
2050     Find the mean value for each of the color.
2051   */
2052   p=(&pixel_list->skip_list);
2053   color=65536L;
2054   count=0;
2055   sum=0.0;
2056   do
2057   {
2058     color=p->nodes[color].next[0];
2059     sum+=(double) p->nodes[color].count*color;
2060     count+=p->nodes[color].count;
2061   } while (count < (ssize_t) pixel_list->length);
2062   sum/=pixel_list->length;
2063   *pixel=ScaleShortToQuantum((unsigned short) sum);
2064 }
2065
2066 static inline void GetMedianPixelList(PixelList *pixel_list,Quantum *pixel)
2067 {
2068   register SkipList
2069     *p;
2070
2071   size_t
2072     color;
2073
2074   ssize_t
2075     count;
2076
2077   /*
2078     Find the median value for each of the color.
2079   */
2080   p=(&pixel_list->skip_list);
2081   color=65536L;
2082   count=0;
2083   do
2084   {
2085     color=p->nodes[color].next[0];
2086     count+=p->nodes[color].count;
2087   } while (count <= (ssize_t) (pixel_list->length >> 1));
2088   *pixel=ScaleShortToQuantum((unsigned short) color);
2089 }
2090
2091 static inline void GetMinimumPixelList(PixelList *pixel_list,Quantum *pixel)
2092 {
2093   register SkipList
2094     *p;
2095
2096   size_t
2097     color,
2098     minimum;
2099
2100   ssize_t
2101     count;
2102
2103   /*
2104     Find the minimum value for each of the color.
2105   */
2106   p=(&pixel_list->skip_list);
2107   count=0;
2108   color=65536UL;
2109   minimum=p->nodes[color].next[0];
2110   do
2111   {
2112     color=p->nodes[color].next[0];
2113     if (color < minimum)
2114       minimum=color;
2115     count+=p->nodes[color].count;
2116   } while (count < (ssize_t) pixel_list->length);
2117   *pixel=ScaleShortToQuantum((unsigned short) minimum);
2118 }
2119
2120 static inline void GetModePixelList(PixelList *pixel_list,Quantum *pixel)
2121 {
2122   register SkipList
2123     *p;
2124
2125   size_t
2126     color,
2127     max_count,
2128     mode;
2129
2130   ssize_t
2131     count;
2132
2133   /*
2134     Make each pixel the 'predominant color' of the specified neighborhood.
2135   */
2136   p=(&pixel_list->skip_list);
2137   color=65536L;
2138   mode=color;
2139   max_count=p->nodes[mode].count;
2140   count=0;
2141   do
2142   {
2143     color=p->nodes[color].next[0];
2144     if (p->nodes[color].count > max_count)
2145       {
2146         mode=color;
2147         max_count=p->nodes[mode].count;
2148       }
2149     count+=p->nodes[color].count;
2150   } while (count < (ssize_t) pixel_list->length);
2151   *pixel=ScaleShortToQuantum((unsigned short) mode);
2152 }
2153
2154 static inline void GetNonpeakPixelList(PixelList *pixel_list,Quantum *pixel)
2155 {
2156   register SkipList
2157     *p;
2158
2159   size_t
2160     color,
2161     next,
2162     previous;
2163
2164   ssize_t
2165     count;
2166
2167   /*
2168     Finds the non peak value for each of the colors.
2169   */
2170   p=(&pixel_list->skip_list);
2171   color=65536L;
2172   next=p->nodes[color].next[0];
2173   count=0;
2174   do
2175   {
2176     previous=color;
2177     color=next;
2178     next=p->nodes[color].next[0];
2179     count+=p->nodes[color].count;
2180   } while (count <= (ssize_t) (pixel_list->length >> 1));
2181   if ((previous == 65536UL) && (next != 65536UL))
2182     color=next;
2183   else
2184     if ((previous != 65536UL) && (next == 65536UL))
2185       color=previous;
2186   *pixel=ScaleShortToQuantum((unsigned short) color);
2187 }
2188
2189 static inline void GetStandardDeviationPixelList(PixelList *pixel_list,
2190   Quantum *pixel)
2191 {
2192   double
2193     sum,
2194     sum_squared;
2195
2196   register SkipList
2197     *p;
2198
2199   size_t
2200     color;
2201
2202   ssize_t
2203     count;
2204
2205   /*
2206     Find the standard-deviation value for each of the color.
2207   */
2208   p=(&pixel_list->skip_list);
2209   color=65536L;
2210   count=0;
2211   sum=0.0;
2212   sum_squared=0.0;
2213   do
2214   {
2215     register ssize_t
2216       i;
2217
2218     color=p->nodes[color].next[0];
2219     sum+=(double) p->nodes[color].count*color;
2220     for (i=0; i < (ssize_t) p->nodes[color].count; i++)
2221       sum_squared+=((double) color)*((double) color);
2222     count+=p->nodes[color].count;
2223   } while (count < (ssize_t) pixel_list->length);
2224   sum/=pixel_list->length;
2225   sum_squared/=pixel_list->length;
2226   *pixel=ScaleShortToQuantum((unsigned short) sqrt(sum_squared-(sum*sum)));
2227 }
2228
2229 static inline void InsertPixelList(const Image *image,const Quantum pixel,
2230   PixelList *pixel_list)
2231 {
2232   size_t
2233     signature;
2234
2235   unsigned short
2236     index;
2237
2238   index=ScaleQuantumToShort(pixel);
2239   signature=pixel_list->skip_list.nodes[index].signature;
2240   if (signature == pixel_list->signature)
2241     {
2242       pixel_list->skip_list.nodes[index].count++;
2243       return;
2244     }
2245   AddNodePixelList(pixel_list,index);
2246 }
2247
2248 static inline double MagickAbsoluteValue(const double x)
2249 {
2250   if (x < 0)
2251     return(-x);
2252   return(x);
2253 }
2254
2255 static inline size_t MagickMax(const size_t x,const size_t y)
2256 {
2257   if (x > y)
2258     return(x);
2259   return(y);
2260 }
2261
2262 static void ResetPixelList(PixelList *pixel_list)
2263 {
2264   int
2265     level;
2266
2267   register SkipNode
2268     *root;
2269
2270   register SkipList
2271     *p;
2272
2273   /*
2274     Reset the skip-list.
2275   */
2276   p=(&pixel_list->skip_list);
2277   root=p->nodes+65536UL;
2278   p->level=0;
2279   for (level=0; level < 9; level++)
2280     root->next[level]=65536UL;
2281   pixel_list->seed=pixel_list->signature++;
2282 }
2283
2284 MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
2285   const size_t width,const size_t height,ExceptionInfo *exception)
2286 {
2287 #define StatisticImageTag  "Statistic/Image"
2288
2289   CacheView
2290     *image_view,
2291     *statistic_view;
2292
2293   Image
2294     *statistic_image;
2295
2296   MagickBooleanType
2297     status;
2298
2299   MagickOffsetType
2300     progress;
2301
2302   PixelList
2303     **restrict pixel_list;
2304
2305   ssize_t
2306     center,
2307     y;
2308
2309   /*
2310     Initialize statistics image attributes.
2311   */
2312   assert(image != (Image *) NULL);
2313   assert(image->signature == MagickSignature);
2314   if (image->debug != MagickFalse)
2315     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2316   assert(exception != (ExceptionInfo *) NULL);
2317   assert(exception->signature == MagickSignature);
2318   statistic_image=CloneImage(image,image->columns,image->rows,MagickTrue,
2319     exception);
2320   if (statistic_image == (Image *) NULL)
2321     return((Image *) NULL);
2322   status=SetImageStorageClass(statistic_image,DirectClass,exception);
2323   if (status == MagickFalse)
2324     {
2325       statistic_image=DestroyImage(statistic_image);
2326       return((Image *) NULL);
2327     }
2328   pixel_list=AcquirePixelListThreadSet(MagickMax(width,1),MagickMax(height,1));
2329   if (pixel_list == (PixelList **) NULL)
2330     {
2331       statistic_image=DestroyImage(statistic_image);
2332       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2333     }
2334   /*
2335     Make each pixel the min / max / median / mode / etc. of the neighborhood.
2336   */
2337   center=(ssize_t) GetPixelChannels(image)*(image->columns+MagickMax(width,1))*
2338     (MagickMax(height,1)/2L)+GetPixelChannels(image)*(MagickMax(width,1)/2L);
2339   status=MagickTrue;
2340   progress=0;
2341   image_view=AcquireVirtualCacheView(image,exception);
2342   statistic_view=AcquireAuthenticCacheView(statistic_image,exception);
2343 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2344   #pragma omp parallel for schedule(static,4) shared(progress,status) \
2345     dynamic_number_threads(image,image->columns,image->rows,1)
2346 #endif
2347   for (y=0; y < (ssize_t) statistic_image->rows; y++)
2348   {
2349     const int
2350       id = GetOpenMPThreadId();
2351
2352     register const Quantum
2353       *restrict p;
2354
2355     register Quantum
2356       *restrict q;
2357
2358     register ssize_t
2359       x;
2360
2361     if (status == MagickFalse)
2362       continue;
2363     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) MagickMax(width,1)/2L),y-
2364       (ssize_t) (MagickMax(height,1)/2L),image->columns+MagickMax(width,1),
2365       MagickMax(height,1),exception);
2366     q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns,      1,exception);
2367     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2368       {
2369         status=MagickFalse;
2370         continue;
2371       }
2372     for (x=0; x < (ssize_t) statistic_image->columns; x++)
2373     {
2374       register ssize_t
2375         i;
2376
2377       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2378       {
2379         PixelChannel
2380           channel;
2381
2382         PixelTrait
2383           statistic_traits,
2384           traits;
2385
2386         Quantum
2387           pixel;
2388
2389         register const Quantum
2390           *restrict pixels;
2391
2392         register ssize_t
2393           u;
2394
2395         ssize_t
2396           v;
2397
2398         channel=GetPixelChannelChannel(image,i);
2399         traits=GetPixelChannelTraits(image,channel);
2400         statistic_traits=GetPixelChannelTraits(statistic_image,channel);
2401         if ((traits == UndefinedPixelTrait) ||
2402             (statistic_traits == UndefinedPixelTrait))
2403           continue;
2404         if (((statistic_traits & CopyPixelTrait) != 0) ||
2405             (GetPixelMask(image,p) != 0))
2406           {
2407             SetPixelChannel(statistic_image,channel,p[center+i],q);
2408             continue;
2409           }
2410         pixels=p;
2411         ResetPixelList(pixel_list[id]);
2412         for (v=0; v < (ssize_t) MagickMax(height,1); v++)
2413         {
2414           for (u=0; u < (ssize_t) MagickMax(width,1); u++)
2415           {
2416             InsertPixelList(image,pixels[i],pixel_list[id]);
2417             pixels+=GetPixelChannels(image);
2418           }
2419           pixels+=image->columns*GetPixelChannels(image);
2420         }
2421         switch (type)
2422         {
2423           case GradientStatistic:
2424           {
2425             double
2426               maximum,
2427               minimum;
2428
2429             GetMinimumPixelList(pixel_list[id],&pixel);
2430             minimum=(double) pixel;
2431             GetMaximumPixelList(pixel_list[id],&pixel);
2432             maximum=(double) pixel;
2433             pixel=ClampToQuantum(MagickAbsoluteValue(maximum-minimum));
2434             break;
2435           }
2436           case MaximumStatistic:
2437           {
2438             GetMaximumPixelList(pixel_list[id],&pixel);
2439             break;
2440           }
2441           case MeanStatistic:
2442           {
2443             GetMeanPixelList(pixel_list[id],&pixel);
2444             break;
2445           }
2446           case MedianStatistic:
2447           default:
2448           {
2449             GetMedianPixelList(pixel_list[id],&pixel);
2450             break;
2451           }
2452           case MinimumStatistic:
2453           {
2454             GetMinimumPixelList(pixel_list[id],&pixel);
2455             break;
2456           }
2457           case ModeStatistic:
2458           {
2459             GetModePixelList(pixel_list[id],&pixel);
2460             break;
2461           }
2462           case NonpeakStatistic:
2463           {
2464             GetNonpeakPixelList(pixel_list[id],&pixel);
2465             break;
2466           }
2467           case StandardDeviationStatistic:
2468           {
2469             GetStandardDeviationPixelList(pixel_list[id],&pixel);
2470             break;
2471           }
2472         }
2473         SetPixelChannel(statistic_image,channel,pixel,q);
2474       }
2475       p+=GetPixelChannels(image);
2476       q+=GetPixelChannels(statistic_image);
2477     }
2478     if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
2479       status=MagickFalse;
2480     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2481       {
2482         MagickBooleanType
2483           proceed;
2484
2485 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2486         #pragma omp critical (MagickCore_StatisticImage)
2487 #endif
2488         proceed=SetImageProgress(image,StatisticImageTag,progress++,
2489           image->rows);
2490         if (proceed == MagickFalse)
2491           status=MagickFalse;
2492       }
2493   }
2494   statistic_view=DestroyCacheView(statistic_view);
2495   image_view=DestroyCacheView(image_view);
2496   pixel_list=DestroyPixelListThreadSet(pixel_list);
2497   return(statistic_image);
2498 }