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-2011 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/segment.h"
  85 #include "MagickCore/semaphore.h"
  86 #include "MagickCore/signature-private.h"
  87 #include "MagickCore/statistic.h"
  88 #include "MagickCore/string_.h"
  89 #include "MagickCore/thread-private.h"
  90 #include "MagickCore/timer.h"
  91 #include "MagickCore/utility.h"
  92 #include "MagickCore/version.h"
  93 \f
  94 /*
  95 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  96 %                                                                             %
  97 %                                                                             %
  98 %                                                                             %
  99 %     E v a l u a t e I m a g e                                               %
 100 %                                                                             %
 101 %                                                                             %
 102 %                                                                             %
 103 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 104 %
 105 %  EvaluateImage() applies a value to the image with an arithmetic, relational,
 106 %  or logical operator to an image. Use these operations to lighten or darken
 107 %  an image, to increase or decrease contrast in an image, or to produce the
 108 %  "negative" of an image.
 109 %
 110 %  The format of the EvaluateImage method is:
 111 %
 112 %      MagickBooleanType EvaluateImage(Image *image,
 113 %        const MagickEvaluateOperator op,const double value,
 114 %        ExceptionInfo *exception)
 115 %      MagickBooleanType EvaluateImages(Image *images,
 116 %        const MagickEvaluateOperator op,const double value,
 117 %        ExceptionInfo *exception)
 118 %
 119 %  A description of each parameter follows:
 120 %
 121 %    o image: the image.
 122 %
 123 %    o op: A channel op.
 124 %
 125 %    o value: A value value.
 126 %
 127 %    o exception: return any errors or warnings in this structure.
 128 %
 129 */
 130
 131 static PixelInfo **DestroyPixelThreadSet(PixelInfo **pixels)
 132 {
 133   register ssize_t
 134     i;
 135
 136   assert(pixels != (PixelInfo **) NULL);
 137   for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
 138     if (pixels[i] != (PixelInfo *) NULL)
 139       pixels[i]=(PixelInfo *) RelinquishMagickMemory(pixels[i]);
 140   pixels=(PixelInfo **) RelinquishMagickMemory(pixels);
 141   return(pixels);
 142 }
 143
 144 static PixelInfo **AcquirePixelThreadSet(const Image *image,
 145   const size_t number_images)
 146 {
 147   register ssize_t
 148     i,
 149     j;
 150
 151   PixelInfo
 152     **pixels;
 153
 154   size_t
 155     length,
 156     number_threads;
 157
 158   number_threads=GetOpenMPMaximumThreads();
 159   pixels=(PixelInfo **) AcquireQuantumMemory(number_threads,
 160     sizeof(*pixels));
 161   if (pixels == (PixelInfo **) NULL)
 162     return((PixelInfo **) NULL);
 163   (void) ResetMagickMemory(pixels,0,number_threads*sizeof(*pixels));
 164   for (i=0; i < (ssize_t) number_threads; i++)
 165   {
 166     length=image->columns;
 167     if (length < number_images)
 168       length=number_images;
 169     pixels[i]=(PixelInfo *) AcquireQuantumMemory(length,
 170       sizeof(**pixels));
 171     if (pixels[i] == (PixelInfo *) NULL)
 172       return(DestroyPixelThreadSet(pixels));
 173     for (j=0; j < (ssize_t) length; j++)
 174       GetPixelInfo(image,&pixels[i][j]);
 175   }
 176   return(pixels);
 177 }
 178
 179 static inline double MagickMax(const double x,const double y)
 180 {
 181   if (x > y)
 182     return(x);
 183   return(y);
 184 }
 185
 186 #if defined(__cplusplus) || defined(c_plusplus)
 187 extern "C" {
 188 #endif
 189
 190 static int IntensityCompare(const void *x,const void *y)
 191 {
 192   const PixelInfo
 193     *color_1,
 194     *color_2;
 195
 196   int
 197     intensity;
 198
 199   color_1=(const PixelInfo *) x;
 200   color_2=(const PixelInfo *) y;
 201   intensity=(int) GetPixelInfoIntensity(color_2)-(int)
 202     GetPixelInfoIntensity(color_1);
 203   return(intensity);
 204 }
 205
 206 #if defined(__cplusplus) || defined(c_plusplus)
 207 }
 208 #endif
 209
 210 static inline double MagickMin(const double x,const double y)
 211 {
 212   if (x < y)
 213     return(x);
 214   return(y);
 215 }
 216
 217 static MagickRealType ApplyEvaluateOperator(RandomInfo *random_info,
 218   Quantum pixel,const MagickEvaluateOperator op,const MagickRealType value)
 219 {
 220   MagickRealType
 221     result;
 222
 223   result=0.0;
 224   switch (op)
 225   {
 226     case UndefinedEvaluateOperator:
 227       break;
 228     case AbsEvaluateOperator:
 229     {
 230       result=(MagickRealType) fabs((double) (pixel+value));
 231       break;
 232     }
 233     case AddEvaluateOperator:
 234     {
 235       result=(MagickRealType) (pixel+value);
 236       break;
 237     }
 238     case AddModulusEvaluateOperator:
 239     {
 240       /*
 241         This returns a 'floored modulus' of the addition which is a
 242         positive result.  It differs from  % or fmod() which returns a
 243         'truncated modulus' result, where floor() is replaced by trunc()
 244         and could return a negative result (which is clipped).
 245       */
 246       result=pixel+value;
 247       result-=(QuantumRange+1.0)*floor((double) result/(QuantumRange+1.0));
 248       break;
 249     }
 250     case AndEvaluateOperator:
 251     {
 252       result=(MagickRealType) ((size_t) pixel & (size_t) (value+0.5));
 253       break;
 254     }
 255     case CosineEvaluateOperator:
 256     {
 257       result=(MagickRealType) (QuantumRange*(0.5*cos((double) (2.0*MagickPI*
 258         QuantumScale*pixel*value))+0.5));
 259       break;
 260     }
 261     case DivideEvaluateOperator:
 262     {
 263       result=pixel/(value == 0.0 ? 1.0 : value);
 264       break;
 265     }
 266     case ExponentialEvaluateOperator:
 267     {
 268       result=(MagickRealType) (QuantumRange*exp((double) (value*QuantumScale*
 269         pixel)));
 270       break;
 271     }
 272     case GaussianNoiseEvaluateOperator:
 273     {
 274       result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
 275         GaussianNoise,value);
 276       break;
 277     }
 278     case ImpulseNoiseEvaluateOperator:
 279     {
 280       result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
 281         ImpulseNoise,value);
 282       break;
 283     }
 284     case LaplacianNoiseEvaluateOperator:
 285     {
 286       result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
 287         LaplacianNoise,value);
 288       break;
 289     }
 290     case LeftShiftEvaluateOperator:
 291     {
 292       result=(MagickRealType) ((size_t) pixel << (size_t) (value+0.5));
 293       break;
 294     }
 295     case LogEvaluateOperator:
 296     {
 297       result=(MagickRealType) (QuantumRange*log((double) (QuantumScale*value*
 298         pixel+1.0))/log((double) (value+1.0)));
 299       break;
 300     }
 301     case MaxEvaluateOperator:
 302     {
 303       result=(MagickRealType) MagickMax((double) pixel,value);
 304       break;
 305     }
 306     case MeanEvaluateOperator:
 307     {
 308       result=(MagickRealType) (pixel+value);
 309       break;
 310     }
 311     case MedianEvaluateOperator:
 312     {
 313       result=(MagickRealType) (pixel+value);
 314       break;
 315     }
 316     case MinEvaluateOperator:
 317     {
 318       result=(MagickRealType) MagickMin((double) pixel,value);
 319       break;
 320     }
 321     case MultiplicativeNoiseEvaluateOperator:
 322     {
 323       result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
 324         MultiplicativeGaussianNoise,value);
 325       break;
 326     }
 327     case MultiplyEvaluateOperator:
 328     {
 329       result=(MagickRealType) (value*pixel);
 330       break;
 331     }
 332     case OrEvaluateOperator:
 333     {
 334       result=(MagickRealType) ((size_t) pixel | (size_t) (value+0.5));
 335       break;
 336     }
 337     case PoissonNoiseEvaluateOperator:
 338     {
 339       result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
 340         PoissonNoise,value);
 341       break;
 342     }
 343     case PowEvaluateOperator:
 344     {
 345       result=(MagickRealType) (QuantumRange*pow((double) (QuantumScale*pixel),
 346         (double) value));
 347       break;
 348     }
 349     case RightShiftEvaluateOperator:
 350     {
 351       result=(MagickRealType) ((size_t) pixel >> (size_t) (value+0.5));
 352       break;
 353     }
 354     case SetEvaluateOperator:
 355     {
 356       result=value;
 357       break;
 358     }
 359     case SineEvaluateOperator:
 360     {
 361       result=(MagickRealType) (QuantumRange*(0.5*sin((double) (2.0*MagickPI*
 362         QuantumScale*pixel*value))+0.5));
 363       break;
 364     }
 365     case SubtractEvaluateOperator:
 366     {
 367       result=(MagickRealType) (pixel-value);
 368       break;
 369     }
 370     case ThresholdEvaluateOperator:
 371     {
 372       result=(MagickRealType) (((MagickRealType) pixel <= value) ? 0 :
 373         QuantumRange);
 374       break;
 375     }
 376     case ThresholdBlackEvaluateOperator:
 377     {
 378       result=(MagickRealType) (((MagickRealType) pixel <= value) ? 0 : pixel);
 379       break;
 380     }
 381     case ThresholdWhiteEvaluateOperator:
 382     {
 383       result=(MagickRealType) (((MagickRealType) pixel > value) ? QuantumRange :
 384         pixel);
 385       break;
 386     }
 387     case UniformNoiseEvaluateOperator:
 388     {
 389       result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
 390         UniformNoise,value);
 391       break;
 392     }
 393     case XorEvaluateOperator:
 394     {
 395       result=(MagickRealType) ((size_t) pixel ^ (size_t) (value+0.5));
 396       break;
 397     }
 398   }
 399   return(result);
 400 }
 401
 402 MagickExport Image *EvaluateImages(const Image *images,
 403   const MagickEvaluateOperator op,ExceptionInfo *exception)
 404 {
 405 #define EvaluateImageTag  "Evaluate/Image"
 406
 407   CacheView
 408     *evaluate_view;
 409
 410   const Image
 411     *next;
 412
 413   Image
 414     *evaluate_image;
 415
 416   MagickBooleanType
 417     status;
 418
 419   MagickOffsetType
 420     progress;
 421
 422   PixelInfo
 423     **restrict evaluate_pixels,
 424     zero;
 425
 426   RandomInfo
 427     **restrict random_info;
 428
 429   size_t
 430     number_images;
 431
 432   ssize_t
 433     y;
 434
 435   /*
 436     Ensure the image are the same size.
 437   */
 438   assert(images != (Image *) NULL);
 439   assert(images->signature == MagickSignature);
 440   if (images->debug != MagickFalse)
 441     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
 442   assert(exception != (ExceptionInfo *) NULL);
 443   assert(exception->signature == MagickSignature);
 444   for (next=images; next != (Image *) NULL; next=GetNextImageInList(next))
 445     if ((next->columns != images->columns) || (next->rows != images->rows))
 446       {
 447         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
 448           "ImageWidthsOrHeightsDiffer","`%s'",images->filename);
 449         return((Image *) NULL);
 450       }
 451   /*
 452     Initialize evaluate next attributes.
 453   */
 454   evaluate_image=CloneImage(images,images->columns,images->rows,MagickTrue,
 455     exception);
 456   if (evaluate_image == (Image *) NULL)
 457     return((Image *) NULL);
 458   if (SetImageStorageClass(evaluate_image,DirectClass,exception) == MagickFalse)
 459     {
 460       evaluate_image=DestroyImage(evaluate_image);
 461       return((Image *) NULL);
 462     }
 463   number_images=GetImageListLength(images);
 464   evaluate_pixels=AcquirePixelThreadSet(images,number_images);
 465   if (evaluate_pixels == (PixelInfo **) NULL)
 466     {
 467       evaluate_image=DestroyImage(evaluate_image);
 468       (void) ThrowMagickException(exception,GetMagickModule(),
 469         ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
 470       return((Image *) NULL);
 471     }
 472   /*
 473     Evaluate image pixels.
 474   */
 475   status=MagickTrue;
 476   progress=0;
 477   GetPixelInfo(images,&zero);
 478   random_info=AcquireRandomInfoThreadSet();
 479   evaluate_view=AcquireCacheView(evaluate_image);
 480   if (op == MedianEvaluateOperator)
 481 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 482   #pragma omp parallel for schedule(dynamic) shared(progress,status)
 483 #endif
 484     for (y=0; y < (ssize_t) evaluate_image->rows; y++)
 485     {
 486       CacheView
 487         *image_view;
 488
 489       const Image
 490         *next;
 491
 492       const int
 493         id = GetOpenMPThreadId();
 494
 495       register PixelInfo
 496         *evaluate_pixel;
 497
 498       register Quantum
 499         *restrict q;
 500
 501       register ssize_t
 502         x;
 503
 504       if (status == MagickFalse)
 505         continue;
 506       q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,evaluate_image->columns,
 507         1,exception);
 508       if (q == (Quantum *) NULL)
 509         {
 510           status=MagickFalse;
 511           continue;
 512         }
 513       evaluate_pixel=evaluate_pixels[id];
 514       for (x=0; x < (ssize_t) evaluate_image->columns; x++)
 515       {
 516         register ssize_t
 517           i;
 518
 519         for (i=0; i < (ssize_t) number_images; i++)
 520           evaluate_pixel[i]=zero;
 521         next=images;
 522         for (i=0; i < (ssize_t) number_images; i++)
 523         {
 524           register const Quantum
 525             *p;
 526
 527           image_view=AcquireCacheView(next);
 528           p=GetCacheViewVirtualPixels(image_view,x,y,1,1,exception);
 529           if (p == (const Quantum *) NULL)
 530             {
 531               image_view=DestroyCacheView(image_view);
 532               break;
 533             }
 534           evaluate_pixel[i].red=ApplyEvaluateOperator(random_info[id],
 535             GetPixelRed(next,p),op,evaluate_pixel[i].red);
 536           evaluate_pixel[i].green=ApplyEvaluateOperator(random_info[id],
 537             GetPixelGreen(next,p),op,evaluate_pixel[i].green);
 538           evaluate_pixel[i].blue=ApplyEvaluateOperator(random_info[id],
 539             GetPixelBlue(next,p),op,evaluate_pixel[i].blue);
 540           if (evaluate_image->colorspace == CMYKColorspace)
 541             evaluate_pixel[i].black=ApplyEvaluateOperator(random_info[id],
 542               GetPixelBlack(next,p),op,evaluate_pixel[i].black);
 543           evaluate_pixel[i].alpha=ApplyEvaluateOperator(random_info[id],
 544             GetPixelAlpha(next,p),op,evaluate_pixel[i].alpha);
 545           image_view=DestroyCacheView(image_view);
 546           next=GetNextImageInList(next);
 547         }
 548         qsort((void *) evaluate_pixel,number_images,sizeof(*evaluate_pixel),
 549           IntensityCompare);
 550         SetPixelRed(evaluate_image,
 551           ClampToQuantum(evaluate_pixel[i/2].red),q);
 552         SetPixelGreen(evaluate_image,
 553           ClampToQuantum(evaluate_pixel[i/2].green),q);
 554         SetPixelBlue(evaluate_image,
 555           ClampToQuantum(evaluate_pixel[i/2].blue),q);
 556         if (evaluate_image->colorspace == CMYKColorspace)
 557           SetPixelBlack(evaluate_image,
 558           ClampToQuantum(evaluate_pixel[i/2].black),q);
 559         if (evaluate_image->matte == MagickFalse)
 560           SetPixelAlpha(evaluate_image,
 561             ClampToQuantum(evaluate_pixel[i/2].alpha),q);
 562         else
 563           SetPixelAlpha(evaluate_image,
 564             ClampToQuantum(evaluate_pixel[i/2].alpha),q);
 565         q+=GetPixelChannels(evaluate_image);
 566       }
 567       if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
 568         status=MagickFalse;
 569       if (images->progress_monitor != (MagickProgressMonitor) NULL)
 570         {
 571           MagickBooleanType
 572             proceed;
 573
 574 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 575           #pragma omp critical (MagickCore_EvaluateImages)
 576 #endif
 577           proceed=SetImageProgress(images,EvaluateImageTag,progress++,
 578             evaluate_image->rows);
 579           if (proceed == MagickFalse)
 580             status=MagickFalse;
 581         }
 582     }
 583   else
 584 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 585   #pragma omp parallel for schedule(dynamic) shared(progress,status)
 586 #endif
 587     for (y=0; y < (ssize_t) evaluate_image->rows; y++)
 588     {
 589       CacheView
 590         *image_view;
 591
 592       const Image
 593         *next;
 594
 595       const int
 596         id = GetOpenMPThreadId();
 597
 598       register ssize_t
 599         i,
 600         x;
 601
 602       register PixelInfo
 603         *evaluate_pixel;
 604
 605       register Quantum
 606         *restrict q;
 607
 608       if (status == MagickFalse)
 609         continue;
 610       q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,evaluate_image->columns,
 611         1,exception);
 612       if (q == (Quantum *) NULL)
 613         {
 614           status=MagickFalse;
 615           continue;
 616         }
 617       evaluate_pixel=evaluate_pixels[id];
 618       for (x=0; x < (ssize_t) evaluate_image->columns; x++)
 619         evaluate_pixel[x]=zero;
 620       next=images;
 621       for (i=0; i < (ssize_t) number_images; i++)
 622       {
 623         register const Quantum
 624           *p;
 625
 626         image_view=AcquireCacheView(next);
 627         p=GetCacheViewVirtualPixels(image_view,0,y,next->columns,1,exception);
 628         if (p == (const Quantum *) NULL)
 629           {
 630             image_view=DestroyCacheView(image_view);
 631             break;
 632           }
 633         for (x=0; x < (ssize_t) next->columns; x++)
 634         {
 635           evaluate_pixel[x].red=ApplyEvaluateOperator(random_info[id],
 636             GetPixelRed(next,p),i == 0 ? AddEvaluateOperator : op,
 637             evaluate_pixel[x].red);
 638           evaluate_pixel[x].green=ApplyEvaluateOperator(random_info[id],
 639             GetPixelGreen(next,p),i == 0 ? AddEvaluateOperator : op,
 640               evaluate_pixel[x].green);
 641           evaluate_pixel[x].blue=ApplyEvaluateOperator(random_info[id],
 642             GetPixelBlue(next,p),i == 0 ? AddEvaluateOperator : op,
 643               evaluate_pixel[x].blue);
 644           if (evaluate_image->colorspace == CMYKColorspace)
 645             evaluate_pixel[x].black=ApplyEvaluateOperator(random_info[id],
 646               GetPixelBlack(next,p),i == 0 ? AddEvaluateOperator : op,
 647               evaluate_pixel[x].black);
 648           evaluate_pixel[x].alpha=ApplyEvaluateOperator(random_info[id],
 649             GetPixelAlpha(next,p),i == 0 ? AddEvaluateOperator : op,
 650             evaluate_pixel[x].alpha);
 651           p+=GetPixelChannels(next);
 652         }
 653         image_view=DestroyCacheView(image_view);
 654         next=GetNextImageInList(next);
 655       }
 656       if (op == MeanEvaluateOperator)
 657         for (x=0; x < (ssize_t) evaluate_image->columns; x++)
 658         {
 659           evaluate_pixel[x].red/=number_images;
 660           evaluate_pixel[x].green/=number_images;
 661           evaluate_pixel[x].blue/=number_images;
 662           evaluate_pixel[x].black/=number_images;
 663           evaluate_pixel[x].alpha/=number_images;
 664         }
 665       for (x=0; x < (ssize_t) evaluate_image->columns; x++)
 666       {
 667         SetPixelRed(evaluate_image,ClampToQuantum(evaluate_pixel[x].red),q);
 668         SetPixelGreen(evaluate_image,ClampToQuantum(evaluate_pixel[x].green),q);
 669         SetPixelBlue(evaluate_image,ClampToQuantum(evaluate_pixel[x].blue),q);
 670         if (evaluate_image->colorspace == CMYKColorspace)
 671           SetPixelBlack(evaluate_image,ClampToQuantum(evaluate_pixel[x].black),
 672             q);
 673         if (evaluate_image->matte == MagickFalse)
 674           SetPixelAlpha(evaluate_image,ClampToQuantum(evaluate_pixel[x].alpha),
 675             q);
 676         else
 677           SetPixelAlpha(evaluate_image,ClampToQuantum(evaluate_pixel[x].alpha),
 678             q);
 679         q+=GetPixelChannels(evaluate_image);
 680       }
 681       if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
 682         status=MagickFalse;
 683       if (images->progress_monitor != (MagickProgressMonitor) NULL)
 684         {
 685           MagickBooleanType
 686             proceed;
 687
 688 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 689           #pragma omp critical (MagickCore_EvaluateImages)
 690 #endif
 691           proceed=SetImageProgress(images,EvaluateImageTag,progress++,
 692             evaluate_image->rows);
 693           if (proceed == MagickFalse)
 694             status=MagickFalse;
 695         }
 696     }
 697   evaluate_view=DestroyCacheView(evaluate_view);
 698   evaluate_pixels=DestroyPixelThreadSet(evaluate_pixels);
 699   random_info=DestroyRandomInfoThreadSet(random_info);
 700   if (status == MagickFalse)
 701     evaluate_image=DestroyImage(evaluate_image);
 702   return(evaluate_image);
 703 }
 704
 705 MagickExport MagickBooleanType EvaluateImage(Image *image,
 706   const MagickEvaluateOperator op,const double value,ExceptionInfo *exception)
 707 {
 708   CacheView
 709     *image_view;
 710
 711   MagickBooleanType
 712     status;
 713
 714   MagickOffsetType
 715     progress;
 716
 717   RandomInfo
 718     **restrict random_info;
 719
 720   ssize_t
 721     y;
 722
 723   assert(image != (Image *) NULL);
 724   assert(image->signature == MagickSignature);
 725   if (image->debug != MagickFalse)
 726     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 727   assert(exception != (ExceptionInfo *) NULL);
 728   assert(exception->signature == MagickSignature);
 729   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
 730     return(MagickFalse);
 731   status=MagickTrue;
 732   progress=0;
 733   random_info=AcquireRandomInfoThreadSet();
 734   image_view=AcquireCacheView(image);
 735 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 736   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
 737 #endif
 738   for (y=0; y < (ssize_t) image->rows; y++)
 739   {
 740     const int
 741       id = GetOpenMPThreadId();
 742
 743     register Quantum
 744       *restrict q;
 745
 746     register ssize_t
 747       x;
 748
 749     if (status == MagickFalse)
 750       continue;
 751     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
 752     if (q == (Quantum *) NULL)
 753       {
 754         status=MagickFalse;
 755         continue;
 756       }
 757     for (x=0; x < (ssize_t) image->columns; x++)
 758     {
 759       if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
 760         SetPixelRed(image,ClampToQuantum(ApplyEvaluateOperator(
 761           random_info[id],GetPixelRed(image,q),op,value)),q);
 762       if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
 763         SetPixelGreen(image,ClampToQuantum(ApplyEvaluateOperator(
 764           random_info[id],GetPixelGreen(image,q),op,value)),q);
 765       if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
 766         SetPixelBlue(image,ClampToQuantum(ApplyEvaluateOperator(
 767           random_info[id],GetPixelBlue(image,q),op,value)),q);
 768       if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
 769           (image->colorspace == CMYKColorspace))
 770         SetPixelBlack(image,ClampToQuantum(ApplyEvaluateOperator(
 771           random_info[id],GetPixelBlack(image,q),op,value)),q);
 772       if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0)
 773         {
 774           if (image->matte == MagickFalse)
 775             SetPixelAlpha(image,ClampToQuantum(ApplyEvaluateOperator(
 776               random_info[id],GetPixelAlpha(image,q),op,value)),q);
 777           else
 778             SetPixelAlpha(image,ClampToQuantum(ApplyEvaluateOperator(
 779               random_info[id],GetPixelAlpha(image,q),op,value)),q);
 780         }
 781       q+=GetPixelChannels(image);
 782     }
 783     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
 784       status=MagickFalse;
 785     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 786       {
 787         MagickBooleanType
 788           proceed;
 789
 790 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 791   #pragma omp critical (MagickCore_EvaluateImage)
 792 #endif
 793         proceed=SetImageProgress(image,EvaluateImageTag,progress++,image->rows);
 794         if (proceed == MagickFalse)
 795           status=MagickFalse;
 796       }
 797   }
 798   image_view=DestroyCacheView(image_view);
 799   random_info=DestroyRandomInfoThreadSet(random_info);
 800   return(status);
 801 }
 802 \f
 803 /*
 804 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 805 %                                                                             %
 806 %                                                                             %
 807 %                                                                             %
 808 %     F u n c t i o n I m a g e                                               %
 809 %                                                                             %
 810 %                                                                             %
 811 %                                                                             %
 812 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 813 %
 814 %  FunctionImage() applies a value to the image with an arithmetic, relational,
 815 %  or logical operator to an image. Use these operations to lighten or darken
 816 %  an image, to increase or decrease contrast in an image, or to produce the
 817 %  "negative" of an image.
 818 %
 819 %  The format of the FunctionImage method is:
 820 %
 821 %      MagickBooleanType FunctionImage(Image *image,
 822 %        const MagickFunction function,const ssize_t number_parameters,
 823 %        const double *parameters,ExceptionInfo *exception)
 824 %
 825 %  A description of each parameter follows:
 826 %
 827 %    o image: the image.
 828 %
 829 %    o function: A channel function.
 830 %
 831 %    o parameters: one or more parameters.
 832 %
 833 %    o exception: return any errors or warnings in this structure.
 834 %
 835 */
 836
 837 static Quantum ApplyFunction(Quantum pixel,const MagickFunction function,
 838   const size_t number_parameters,const double *parameters,
 839   ExceptionInfo *exception)
 840 {
 841   MagickRealType
 842     result;
 843
 844   register ssize_t
 845     i;
 846
 847   (void) exception;
 848   result=0.0;
 849   switch (function)
 850   {
 851     case PolynomialFunction:
 852     {
 853       /*
 854        * Polynomial
 855        * Parameters:   polynomial constants,  highest to lowest order
 856        *   For example:      c0*x^3 + c1*x^2 + c2*x  + c3
 857        */
 858       result=0.0;
 859       for (i=0; i < (ssize_t) number_parameters; i++)
 860         result = result*QuantumScale*pixel + parameters[i];
 861       result *= QuantumRange;
 862       break;
 863     }
 864     case SinusoidFunction:
 865     {
 866       /* Sinusoid Function
 867        * Parameters:   Freq, Phase, Ampl, bias
 868        */
 869       double  freq,phase,ampl,bias;
 870       freq  = ( number_parameters >= 1 ) ? parameters[0] : 1.0;
 871       phase = ( number_parameters >= 2 ) ? parameters[1] : 0.0;
 872       ampl  = ( number_parameters >= 3 ) ? parameters[2] : 0.5;
 873       bias  = ( number_parameters >= 4 ) ? parameters[3] : 0.5;
 874       result=(MagickRealType) (QuantumRange*(ampl*sin((double) (2.0*MagickPI*
 875         (freq*QuantumScale*pixel + phase/360.0) )) + bias ) );
 876       break;
 877     }
 878     case ArcsinFunction:
 879     {
 880       /* Arcsin Function  (peged at range limits for invalid results)
 881        * Parameters:   Width, Center, Range, Bias
 882        */
 883       double  width,range,center,bias;
 884       width  = ( number_parameters >= 1 ) ? parameters[0] : 1.0;
 885       center = ( number_parameters >= 2 ) ? parameters[1] : 0.5;
 886       range  = ( number_parameters >= 3 ) ? parameters[2] : 1.0;
 887       bias   = ( number_parameters >= 4 ) ? parameters[3] : 0.5;
 888       result = 2.0/width*(QuantumScale*pixel - center);
 889       if ( result <= -1.0 )
 890         result = bias - range/2.0;
 891       else if ( result >= 1.0 )
 892         result = bias + range/2.0;
 893       else
 894         result=(MagickRealType) (range/MagickPI*asin((double) result)+bias);
 895       result *= QuantumRange;
 896       break;
 897     }
 898     case ArctanFunction:
 899     {
 900       /* Arctan Function
 901        * Parameters:   Slope, Center, Range, Bias
 902        */
 903       double  slope,range,center,bias;
 904       slope  = ( number_parameters >= 1 ) ? parameters[0] : 1.0;
 905       center = ( number_parameters >= 2 ) ? parameters[1] : 0.5;
 906       range  = ( number_parameters >= 3 ) ? parameters[2] : 1.0;
 907       bias   = ( number_parameters >= 4 ) ? parameters[3] : 0.5;
 908       result=(MagickRealType) (MagickPI*slope*(QuantumScale*pixel-center));
 909       result=(MagickRealType) (QuantumRange*(range/MagickPI*atan((double)
 910                   result) + bias ) );
 911       break;
 912     }
 913     case UndefinedFunction:
 914       break;
 915   }
 916   return(ClampToQuantum(result));
 917 }
 918
 919 MagickExport MagickBooleanType FunctionImage(Image *image,
 920   const MagickFunction function,const size_t number_parameters,
 921   const double *parameters,ExceptionInfo *exception)
 922 {
 923 #define FunctionImageTag  "Function/Image "
 924
 925   CacheView
 926     *image_view;
 927
 928   MagickBooleanType
 929     status;
 930
 931   MagickOffsetType
 932     progress;
 933
 934   ssize_t
 935     y;
 936
 937   assert(image != (Image *) NULL);
 938   assert(image->signature == MagickSignature);
 939   if (image->debug != MagickFalse)
 940     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 941   assert(exception != (ExceptionInfo *) NULL);
 942   assert(exception->signature == MagickSignature);
 943   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
 944     return(MagickFalse);
 945   status=MagickTrue;
 946   progress=0;
 947   image_view=AcquireCacheView(image);
 948 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 949   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
 950 #endif
 951   for (y=0; y < (ssize_t) image->rows; y++)
 952   {
 953     register ssize_t
 954       x;
 955
 956     register Quantum
 957       *restrict q;
 958
 959     if (status == MagickFalse)
 960       continue;
 961     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
 962     if (q == (Quantum *) NULL)
 963       {
 964         status=MagickFalse;
 965         continue;
 966       }
 967     for (x=0; x < (ssize_t) image->columns; x++)
 968     {
 969       if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
 970         SetPixelRed(image,ApplyFunction(GetPixelRed(image,q),function,
 971           number_parameters,parameters,exception),q);
 972       if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
 973         SetPixelGreen(image,ApplyFunction(GetPixelGreen(image,q),function,
 974           number_parameters,parameters,exception),q);
 975       if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
 976         SetPixelBlue(image,ApplyFunction(GetPixelBlue(image,q),function,
 977           number_parameters,parameters,exception),q);
 978       if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
 979           (image->colorspace == CMYKColorspace))
 980         SetPixelBlack(image,ApplyFunction(GetPixelBlack(image,q),function,
 981           number_parameters,parameters,exception),q);
 982       if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0)
 983         {
 984           if (image->matte == MagickFalse)
 985             SetPixelAlpha(image,ApplyFunction(GetPixelAlpha(image,q),function,
 986               number_parameters,parameters,exception),q);
 987           else
 988             SetPixelAlpha(image,ApplyFunction(GetPixelAlpha(image,q),function,
 989               number_parameters,parameters,exception),q);
 990         }
 991       q+=GetPixelChannels(image);
 992     }
 993     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
 994       status=MagickFalse;
 995     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 996       {
 997         MagickBooleanType
 998           proceed;
 999
1000 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1001   #pragma omp critical (MagickCore_FunctionImage)
1002 #endif
1003         proceed=SetImageProgress(image,FunctionImageTag,progress++,image->rows);
1004         if (proceed == MagickFalse)
1005           status=MagickFalse;
1006       }
1007   }
1008   image_view=DestroyCacheView(image_view);
1009   return(status);
1010 }
1011 \f
1012 /*
1013 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1014 %                                                                             %
1015 %                                                                             %
1016 %                                                                             %
1017 %   G e t I m a g e E x t r e m a                                             %
1018 %                                                                             %
1019 %                                                                             %
1020 %                                                                             %
1021 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1022 %
1023 %  GetImageExtrema() returns the extrema of one or more image channels.
1024 %
1025 %  The format of the GetImageExtrema method is:
1026 %
1027 %      MagickBooleanType GetImageExtrema(const Image *image,size_t *minima,
1028 %        size_t *maxima,ExceptionInfo *exception)
1029 %
1030 %  A description of each parameter follows:
1031 %
1032 %    o image: the image.
1033 %
1034 %    o minima: the minimum value in the channel.
1035 %
1036 %    o maxima: the maximum value in the channel.
1037 %
1038 %    o exception: return any errors or warnings in this structure.
1039 %
1040 */
1041 MagickExport MagickBooleanType GetImageExtrema(const Image *image,
1042   size_t *minima,size_t *maxima,ExceptionInfo *exception)
1043 {
1044   double
1045     max,
1046     min;
1047
1048   MagickBooleanType
1049     status;
1050
1051   assert(image != (Image *) NULL);
1052   assert(image->signature == MagickSignature);
1053   if (image->debug != MagickFalse)
1054     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1055   status=GetImageRange(image,&min,&max,exception);
1056   *minima=(size_t) ceil(min-0.5);
1057   *maxima=(size_t) floor(max+0.5);
1058   return(status);
1059 }
1060 \f
1061 /*
1062 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1063 %                                                                             %
1064 %                                                                             %
1065 %                                                                             %
1066 %   G e t I m a g e M e a n                                                   %
1067 %                                                                             %
1068 %                                                                             %
1069 %                                                                             %
1070 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1071 %
1072 %  GetImageMean() returns the mean and standard deviation of one or more
1073 %  image channels.
1074 %
1075 %  The format of the GetImageMean method is:
1076 %
1077 %      MagickBooleanType GetImageMean(const Image *image,double *mean,
1078 %        double *standard_deviation,ExceptionInfo *exception)
1079 %
1080 %  A description of each parameter follows:
1081 %
1082 %    o image: the image.
1083 %
1084 %    o mean: the average value in the channel.
1085 %
1086 %    o standard_deviation: the standard deviation of the channel.
1087 %
1088 %    o exception: return any errors or warnings in this structure.
1089 %
1090 */
1091 MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean,
1092   double *standard_deviation,ExceptionInfo *exception)
1093 {
1094   ChannelStatistics
1095     *channel_statistics;
1096
1097   size_t
1098     channels;
1099
1100   assert(image != (Image *) NULL);
1101   assert(image->signature == MagickSignature);
1102   if (image->debug != MagickFalse)
1103     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1104   channel_statistics=GetImageStatistics(image,exception);
1105   if (channel_statistics == (ChannelStatistics *) NULL)
1106     return(MagickFalse);
1107   channels=0;
1108   channel_statistics[CompositeChannels].mean=0.0;
1109   channel_statistics[CompositeChannels].standard_deviation=0.0;
1110   if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
1111     {
1112       channel_statistics[CompositeChannels].mean+=
1113         channel_statistics[RedChannel].mean;
1114       channel_statistics[CompositeChannels].standard_deviation+=
1115         channel_statistics[RedChannel].variance-
1116         channel_statistics[RedChannel].mean*
1117         channel_statistics[RedChannel].mean;
1118       channels++;
1119     }
1120   if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
1121     {
1122       channel_statistics[CompositeChannels].mean+=
1123         channel_statistics[GreenChannel].mean;
1124       channel_statistics[CompositeChannels].standard_deviation+=
1125         channel_statistics[GreenChannel].variance-
1126         channel_statistics[GreenChannel].mean*
1127         channel_statistics[GreenChannel].mean;
1128       channels++;
1129     }
1130   if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
1131     {
1132       channel_statistics[CompositeChannels].mean+=
1133         channel_statistics[BlueChannel].mean;
1134       channel_statistics[CompositeChannels].standard_deviation+=
1135         channel_statistics[BlueChannel].variance-
1136         channel_statistics[BlueChannel].mean*
1137         channel_statistics[BlueChannel].mean;
1138       channels++;
1139     }
1140   if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
1141       (image->colorspace == CMYKColorspace))
1142     {
1143       channel_statistics[CompositeChannels].mean+=
1144         channel_statistics[BlackChannel].mean;
1145       channel_statistics[CompositeChannels].standard_deviation+=
1146         channel_statistics[BlackChannel].variance-
1147         channel_statistics[BlackChannel].mean*
1148         channel_statistics[BlackChannel].mean;
1149       channels++;
1150     }
1151   if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
1152       (image->matte != MagickFalse))
1153     {
1154       channel_statistics[CompositeChannels].mean+=
1155         channel_statistics[AlphaChannel].mean;
1156       channel_statistics[CompositeChannels].standard_deviation+=
1157         channel_statistics[AlphaChannel].variance-
1158         channel_statistics[AlphaChannel].mean*
1159         channel_statistics[AlphaChannel].mean;
1160       channels++;
1161     }
1162   channel_statistics[CompositeChannels].mean/=channels;
1163   channel_statistics[CompositeChannels].standard_deviation=
1164     sqrt(channel_statistics[CompositeChannels].standard_deviation/channels);
1165   *mean=channel_statistics[CompositeChannels].mean;
1166   *standard_deviation=channel_statistics[CompositeChannels].standard_deviation;
1167   channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1168     channel_statistics);
1169   return(MagickTrue);
1170 }
1171 \f
1172 /*
1173 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1174 %                                                                             %
1175 %                                                                             %
1176 %                                                                             %
1177 %   G e t I m a g e K u r t o s i s                                           %
1178 %                                                                             %
1179 %                                                                             %
1180 %                                                                             %
1181 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1182 %
1183 %  GetImageKurtosis() returns the kurtosis and skewness of one or more
1184 %  image channels.
1185 %
1186 %  The format of the GetImageKurtosis method is:
1187 %
1188 %      MagickBooleanType GetImageKurtosis(const Image *image,double *kurtosis,
1189 %        double *skewness,ExceptionInfo *exception)
1190 %
1191 %  A description of each parameter follows:
1192 %
1193 %    o image: the image.
1194 %
1195 %    o kurtosis: the kurtosis of the channel.
1196 %
1197 %    o skewness: the skewness of the channel.
1198 %
1199 %    o exception: return any errors or warnings in this structure.
1200 %
1201 */
1202 MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
1203   double *kurtosis,double *skewness,ExceptionInfo *exception)
1204 {
1205   double
1206     area,
1207     mean,
1208     standard_deviation,
1209     sum_squares,
1210     sum_cubes,
1211     sum_fourth_power;
1212
1213   ssize_t
1214     y;
1215
1216   assert(image != (Image *) NULL);
1217   assert(image->signature == MagickSignature);
1218   if (image->debug != MagickFalse)
1219     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1220   *kurtosis=0.0;
1221   *skewness=0.0;
1222   area=0.0;
1223   mean=0.0;
1224   standard_deviation=0.0;
1225   sum_squares=0.0;
1226   sum_cubes=0.0;
1227   sum_fourth_power=0.0;
1228   for (y=0; y < (ssize_t) image->rows; y++)
1229   {
1230     register const Quantum
1231       *restrict p;
1232
1233     register ssize_t
1234       x;
1235
1236     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1237     if (p == (const Quantum *) NULL)
1238       break;
1239     for (x=0; x < (ssize_t) image->columns; x++)
1240     {
1241       if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
1242         {
1243           mean+=GetPixelRed(image,p);
1244           sum_squares+=(double) GetPixelRed(image,p)*GetPixelRed(image,p);
1245           sum_cubes+=(double) GetPixelRed(image,p)*GetPixelRed(image,p)*
1246             GetPixelRed(image,p);
1247           sum_fourth_power+=(double) GetPixelRed(image,p)*
1248             GetPixelRed(image,p)*GetPixelRed(image,p)*GetPixelRed(image,p);
1249           area++;
1250         }
1251       if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
1252         {
1253           mean+=GetPixelGreen(image,p);
1254           sum_squares+=(double) GetPixelGreen(image,p)*GetPixelGreen(image,p);
1255           sum_cubes+=(double) GetPixelGreen(image,p)*GetPixelGreen(image,p)*
1256             GetPixelGreen(image,p);
1257           sum_fourth_power+=(double) GetPixelGreen(image,p)*
1258             GetPixelGreen(image,p)*GetPixelGreen(image,p)*
1259             GetPixelGreen(image,p);
1260           area++;
1261         }
1262       if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
1263         {
1264           mean+=GetPixelBlue(image,p);
1265           sum_squares+=(double) GetPixelBlue(image,p)*GetPixelBlue(image,p);
1266           sum_cubes+=(double) GetPixelBlue(image,p)*GetPixelBlue(image,p)*
1267             GetPixelBlue(image,p);
1268           sum_fourth_power+=(double) GetPixelBlue(image,p)*
1269             GetPixelBlue(image,p)*GetPixelBlue(image,p)*
1270             GetPixelBlue(image,p);
1271           area++;
1272         }
1273       if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
1274           (image->colorspace == CMYKColorspace))
1275         {
1276           mean+=GetPixelBlack(image,p);
1277           sum_squares+=(double) GetPixelBlack(image,p)*GetPixelBlack(image,p);
1278           sum_cubes+=(double) GetPixelBlack(image,p)*GetPixelBlack(image,p)*
1279             GetPixelBlack(image,p);
1280           sum_fourth_power+=(double) GetPixelBlack(image,p)*
1281             GetPixelBlack(image,p)*GetPixelBlack(image,p)*
1282             GetPixelBlack(image,p);
1283           area++;
1284         }
1285       if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0)
1286         {
1287           mean+=GetPixelAlpha(image,p);
1288           sum_squares+=(double) GetPixelAlpha(image,p)*GetPixelAlpha(image,p);
1289           sum_cubes+=(double) GetPixelAlpha(image,p)*GetPixelAlpha(image,p)*
1290             GetPixelAlpha(image,p);
1291           sum_fourth_power+=(double) GetPixelAlpha(image,p)*
1292             GetPixelAlpha(image,p)*GetPixelAlpha(image,p)*
1293             GetPixelAlpha(image,p);
1294           area++;
1295         }
1296       p+=GetPixelChannels(image);
1297     }
1298   }
1299   if (y < (ssize_t) image->rows)
1300     return(MagickFalse);
1301   if (area != 0.0)
1302     {
1303       mean/=area;
1304       sum_squares/=area;
1305       sum_cubes/=area;
1306       sum_fourth_power/=area;
1307     }
1308   standard_deviation=sqrt(sum_squares-(mean*mean));
1309   if (standard_deviation != 0.0)
1310     {
1311       *kurtosis=sum_fourth_power-4.0*mean*sum_cubes+6.0*mean*mean*sum_squares-
1312         3.0*mean*mean*mean*mean;
1313       *kurtosis/=standard_deviation*standard_deviation*standard_deviation*
1314         standard_deviation;
1315       *kurtosis-=3.0;
1316       *skewness=sum_cubes-3.0*mean*sum_squares+2.0*mean*mean*mean;
1317       *skewness/=standard_deviation*standard_deviation*standard_deviation;
1318     }
1319   return(y == (ssize_t) image->rows ? MagickTrue : MagickFalse);
1320 }
1321 \f
1322 /*
1323 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1324 %                                                                             %
1325 %                                                                             %
1326 %                                                                             %
1327 %   G e t I m a g e R a n g e                                                 %
1328 %                                                                             %
1329 %                                                                             %
1330 %                                                                             %
1331 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1332 %
1333 %  GetImageRange() returns the range of one or more image channels.
1334 %
1335 %  The format of the GetImageRange method is:
1336 %
1337 %      MagickBooleanType GetImageRange(const Image *image,double *minima,
1338 %        double *maxima,ExceptionInfo *exception)
1339 %
1340 %  A description of each parameter follows:
1341 %
1342 %    o image: the image.
1343 %
1344 %    o minima: the minimum value in the channel.
1345 %
1346 %    o maxima: the maximum value in the channel.
1347 %
1348 %    o exception: return any errors or warnings in this structure.
1349 %
1350 */
1351 MagickExport MagickBooleanType GetImageRange(const Image *image,double *minima,
1352   double *maxima,ExceptionInfo *exception)
1353 {
1354   PixelInfo
1355     pixel;
1356
1357   ssize_t
1358     y;
1359
1360   assert(image != (Image *) NULL);
1361   assert(image->signature == MagickSignature);
1362   if (image->debug != MagickFalse)
1363     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1364   *maxima=(-1.0E-37);
1365   *minima=1.0E+37;
1366   GetPixelInfo(image,&pixel);
1367   for (y=0; y < (ssize_t) image->rows; y++)
1368   {
1369     register const Quantum
1370       *restrict p;
1371
1372     register ssize_t
1373       x;
1374
1375     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1376     if (p == (const Quantum *) NULL)
1377       break;
1378     for (x=0; x < (ssize_t) image->columns; x++)
1379     {
1380       SetPixelInfo(image,p,&pixel);
1381       if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
1382         {
1383           if (pixel.red < *minima)
1384             *minima=(double) pixel.red;
1385           if (pixel.red > *maxima)
1386             *maxima=(double) pixel.red;
1387         }
1388       if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
1389         {
1390           if (pixel.green < *minima)
1391             *minima=(double) pixel.green;
1392           if (pixel.green > *maxima)
1393             *maxima=(double) pixel.green;
1394         }
1395       if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
1396         {
1397           if (pixel.blue < *minima)
1398             *minima=(double) pixel.blue;
1399           if (pixel.blue > *maxima)
1400             *maxima=(double) pixel.blue;
1401         }
1402       if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
1403           (image->colorspace == CMYKColorspace))
1404         {
1405           if (pixel.black < *minima)
1406             *minima=(double) pixel.black;
1407           if (pixel.black > *maxima)
1408             *maxima=(double) pixel.black;
1409         }
1410       if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0)
1411         {
1412           if (pixel.alpha < *minima)
1413             *minima=(double) pixel.alpha;
1414           if (pixel.alpha > *maxima)
1415             *maxima=(double) pixel.alpha;
1416         }
1417       p+=GetPixelChannels(image);
1418     }
1419   }
1420   return(y == (ssize_t) image->rows ? MagickTrue : MagickFalse);
1421 }
1422 \f
1423 /*
1424 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1425 %                                                                             %
1426 %                                                                             %
1427 %                                                                             %
1428 %   G e t I m a g e S t a t i s t i c s                                       %
1429 %                                                                             %
1430 %                                                                             %
1431 %                                                                             %
1432 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1433 %
1434 %  GetImageStatistics() returns statistics for each channel in the
1435 %  image.  The statistics include the channel depth, its minima, maxima, mean,
1436 %  standard deviation, kurtosis and skewness.  You can access the red channel
1437 %  mean, for example, like this:
1438 %
1439 %      channel_statistics=GetImageStatistics(image,exception);
1440 %      red_mean=channel_statistics[RedChannel].mean;
1441 %
1442 %  Use MagickRelinquishMemory() to free the statistics buffer.
1443 %
1444 %  The format of the GetImageStatistics method is:
1445 %
1446 %      ChannelStatistics *GetImageStatistics(const Image *image,
1447 %        ExceptionInfo *exception)
1448 %
1449 %  A description of each parameter follows:
1450 %
1451 %    o image: the image.
1452 %
1453 %    o exception: return any errors or warnings in this structure.
1454 %
1455 */
1456 MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
1457   ExceptionInfo *exception)
1458 {
1459   ChannelStatistics
1460     *channel_statistics;
1461
1462   double
1463     area;
1464
1465   MagickStatusType
1466     status;
1467
1468   QuantumAny
1469     range;
1470
1471   register ssize_t
1472     i;
1473
1474   size_t
1475     channels,
1476     depth,
1477     length;
1478
1479   ssize_t
1480     y;
1481
1482   assert(image != (Image *) NULL);
1483   assert(image->signature == MagickSignature);
1484   if (image->debug != MagickFalse)
1485     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1486   length=CompositeChannels+1UL;
1487   channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(length,
1488     sizeof(*channel_statistics));
1489   if (channel_statistics == (ChannelStatistics *) NULL)
1490     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1491   (void) ResetMagickMemory(channel_statistics,0,length*
1492     sizeof(*channel_statistics));
1493   for (i=0; i <= (ssize_t) CompositeChannels; i++)
1494   {
1495     channel_statistics[i].depth=1;
1496     channel_statistics[i].maxima=(-1.0E-37);
1497     channel_statistics[i].minima=1.0E+37;
1498   }
1499   for (y=0; y < (ssize_t) image->rows; y++)
1500   {
1501     register const Quantum
1502       *restrict p;
1503
1504     register ssize_t
1505       x;
1506
1507     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1508     if (p == (const Quantum *) NULL)
1509       break;
1510     for (x=0; x < (ssize_t) image->columns; )
1511     {
1512       if (channel_statistics[RedChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
1513         {
1514           depth=channel_statistics[RedChannel].depth;
1515           range=GetQuantumRange(depth);
1516           status=GetPixelRed(image,p) != ScaleAnyToQuantum(ScaleQuantumToAny(
1517             GetPixelRed(image,p),range),range) ? MagickTrue : MagickFalse;
1518           if (status != MagickFalse)
1519             {
1520               channel_statistics[RedChannel].depth++;
1521               continue;
1522             }
1523         }
1524       if (channel_statistics[GreenChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
1525         {
1526           depth=channel_statistics[GreenChannel].depth;
1527           range=GetQuantumRange(depth);
1528           status=GetPixelGreen(image,p) != ScaleAnyToQuantum(ScaleQuantumToAny(
1529             GetPixelGreen(image,p),range),range) ? MagickTrue : MagickFalse;
1530           if (status != MagickFalse)
1531             {
1532               channel_statistics[GreenChannel].depth++;
1533               continue;
1534             }
1535         }
1536       if (channel_statistics[BlueChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
1537         {
1538           depth=channel_statistics[BlueChannel].depth;
1539           range=GetQuantumRange(depth);
1540           status=GetPixelBlue(image,p) != ScaleAnyToQuantum(ScaleQuantumToAny(
1541             GetPixelBlue(image,p),range),range) ? MagickTrue : MagickFalse;
1542           if (status != MagickFalse)
1543             {
1544               channel_statistics[BlueChannel].depth++;
1545               continue;
1546             }
1547         }
1548       if (image->matte != MagickFalse)
1549         {
1550           if (channel_statistics[AlphaChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
1551             {
1552               depth=channel_statistics[AlphaChannel].depth;
1553               range=GetQuantumRange(depth);
1554               status=GetPixelAlpha(image,p) != ScaleAnyToQuantum(
1555                 ScaleQuantumToAny(GetPixelAlpha(image,p),range),range) ?
1556                 MagickTrue : MagickFalse;
1557               if (status != MagickFalse)
1558                 {
1559                   channel_statistics[AlphaChannel].depth++;
1560                   continue;
1561                 }
1562             }
1563           }
1564       if (image->colorspace == CMYKColorspace)
1565         {
1566           if (channel_statistics[BlackChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
1567             {
1568               depth=channel_statistics[BlackChannel].depth;
1569               range=GetQuantumRange(depth);
1570               status=GetPixelBlack(image,p) != ScaleAnyToQuantum(
1571                 ScaleQuantumToAny(GetPixelBlack(image,p),range),range) ?
1572                 MagickTrue : MagickFalse;
1573               if (status != MagickFalse)
1574                 {
1575                   channel_statistics[BlackChannel].depth++;
1576                   continue;
1577                 }
1578             }
1579         }
1580       if ((double) GetPixelRed(image,p) < channel_statistics[RedChannel].minima)
1581         channel_statistics[RedChannel].minima=(double) GetPixelRed(image,p);
1582       if ((double) GetPixelRed(image,p) > channel_statistics[RedChannel].maxima)
1583         channel_statistics[RedChannel].maxima=(double) GetPixelRed(image,p);
1584       channel_statistics[RedChannel].sum+=GetPixelRed(image,p);
1585       channel_statistics[RedChannel].sum_squared+=(double)
1586         GetPixelRed(image,p)*GetPixelRed(image,p);
1587       channel_statistics[RedChannel].sum_cubed+=(double)
1588         GetPixelRed(image,p)*GetPixelRed(image,p)*GetPixelRed(image,p);
1589       channel_statistics[RedChannel].sum_fourth_power+=(double)
1590         GetPixelRed(image,p)*GetPixelRed(image,p)*GetPixelRed(image,p)*
1591         GetPixelRed(image,p);
1592       if ((double) GetPixelGreen(image,p) < channel_statistics[GreenChannel].minima)
1593         channel_statistics[GreenChannel].minima=(double)
1594           GetPixelGreen(image,p);
1595       if ((double) GetPixelGreen(image,p) > channel_statistics[GreenChannel].maxima)
1596         channel_statistics[GreenChannel].maxima=(double) GetPixelGreen(image,p);
1597       channel_statistics[GreenChannel].sum+=GetPixelGreen(image,p);
1598       channel_statistics[GreenChannel].sum_squared+=(double)
1599         GetPixelGreen(image,p)*GetPixelGreen(image,p);
1600       channel_statistics[GreenChannel].sum_cubed+=(double)
1601         GetPixelGreen(image,p)*GetPixelGreen(image,p)*GetPixelGreen(image,p);
1602       channel_statistics[GreenChannel].sum_fourth_power+=(double)
1603         GetPixelGreen(image,p)*GetPixelGreen(image,p)*GetPixelGreen(image,p)*
1604         GetPixelGreen(image,p);
1605       if ((double) GetPixelBlue(image,p) < channel_statistics[BlueChannel].minima)
1606         channel_statistics[BlueChannel].minima=(double) GetPixelBlue(image,p);
1607       if ((double) GetPixelBlue(image,p) > channel_statistics[BlueChannel].maxima)
1608         channel_statistics[BlueChannel].maxima=(double) GetPixelBlue(image,p);
1609       channel_statistics[BlueChannel].sum+=GetPixelBlue(image,p);
1610       channel_statistics[BlueChannel].sum_squared+=(double)
1611         GetPixelBlue(image,p)*GetPixelBlue(image,p);
1612       channel_statistics[BlueChannel].sum_cubed+=(double)
1613         GetPixelBlue(image,p)*GetPixelBlue(image,p)*GetPixelBlue(image,p);
1614       channel_statistics[BlueChannel].sum_fourth_power+=(double)
1615         GetPixelBlue(image,p)*GetPixelBlue(image,p)*GetPixelBlue(image,p)*
1616         GetPixelBlue(image,p);
1617       if (image->colorspace == CMYKColorspace)
1618         {
1619           if ((double) GetPixelBlack(image,p) < channel_statistics[BlackChannel].minima)
1620             channel_statistics[BlackChannel].minima=(double)
1621               GetPixelBlack(image,p);
1622           if ((double) GetPixelBlack(image,p) > channel_statistics[BlackChannel].maxima)
1623             channel_statistics[BlackChannel].maxima=(double)
1624               GetPixelBlack(image,p);
1625           channel_statistics[BlackChannel].sum+=GetPixelBlack(image,p);
1626           channel_statistics[BlackChannel].sum_squared+=(double)
1627             GetPixelBlack(image,p)*GetPixelBlack(image,p);
1628           channel_statistics[BlackChannel].sum_cubed+=(double)
1629             GetPixelBlack(image,p)*GetPixelBlack(image,p)*
1630             GetPixelBlack(image,p);
1631           channel_statistics[BlackChannel].sum_fourth_power+=(double)
1632             GetPixelBlack(image,p)*GetPixelBlack(image,p)*
1633             GetPixelBlack(image,p)*GetPixelBlack(image,p);
1634         }
1635       if (image->matte != MagickFalse)
1636         {
1637           if ((double) GetPixelAlpha(image,p) < channel_statistics[AlphaChannel].minima)
1638             channel_statistics[AlphaChannel].minima=(double)
1639               GetPixelAlpha(image,p);
1640           if ((double) GetPixelAlpha(image,p) > channel_statistics[AlphaChannel].maxima)
1641             channel_statistics[AlphaChannel].maxima=(double)
1642               GetPixelAlpha(image,p);
1643           channel_statistics[AlphaChannel].sum+=GetPixelAlpha(image,p);
1644           channel_statistics[AlphaChannel].sum_squared+=(double)
1645             GetPixelAlpha(image,p)*GetPixelAlpha(image,p);
1646           channel_statistics[AlphaChannel].sum_cubed+=(double)
1647             GetPixelAlpha(image,p)*GetPixelAlpha(image,p)*
1648             GetPixelAlpha(image,p);
1649           channel_statistics[AlphaChannel].sum_fourth_power+=(double)
1650             GetPixelAlpha(image,p)*GetPixelAlpha(image,p)*
1651             GetPixelAlpha(image,p)*GetPixelAlpha(image,p);
1652         }
1653       x++;
1654       p+=GetPixelChannels(image);
1655     }
1656   }
1657   area=(double) image->columns*image->rows;
1658   for (i=0; i < (ssize_t) CompositeChannels; i++)
1659   {
1660     channel_statistics[i].sum/=area;
1661     channel_statistics[i].sum_squared/=area;
1662     channel_statistics[i].sum_cubed/=area;
1663     channel_statistics[i].sum_fourth_power/=area;
1664     channel_statistics[i].mean=channel_statistics[i].sum;
1665     channel_statistics[i].variance=channel_statistics[i].sum_squared;
1666     channel_statistics[i].standard_deviation=sqrt(
1667       channel_statistics[i].variance-(channel_statistics[i].mean*
1668       channel_statistics[i].mean));
1669   }
1670   for (i=0; i < (ssize_t) CompositeChannels; i++)
1671   {
1672     channel_statistics[CompositeChannels].depth=(size_t) MagickMax((double)
1673       channel_statistics[CompositeChannels].depth,(double)
1674       channel_statistics[i].depth);
1675     channel_statistics[CompositeChannels].minima=MagickMin(
1676       channel_statistics[CompositeChannels].minima,
1677       channel_statistics[i].minima);
1678     channel_statistics[CompositeChannels].maxima=MagickMax(
1679       channel_statistics[CompositeChannels].maxima,
1680       channel_statistics[i].maxima);
1681     channel_statistics[CompositeChannels].sum+=channel_statistics[i].sum;
1682     channel_statistics[CompositeChannels].sum_squared+=
1683       channel_statistics[i].sum_squared;
1684     channel_statistics[CompositeChannels].sum_cubed+=
1685       channel_statistics[i].sum_cubed;
1686     channel_statistics[CompositeChannels].sum_fourth_power+=
1687       channel_statistics[i].sum_fourth_power;
1688     channel_statistics[CompositeChannels].mean+=channel_statistics[i].mean;
1689     channel_statistics[CompositeChannels].variance+=
1690       channel_statistics[i].variance-channel_statistics[i].mean*
1691       channel_statistics[i].mean;
1692     channel_statistics[CompositeChannels].standard_deviation+=
1693       channel_statistics[i].variance-channel_statistics[i].mean*
1694       channel_statistics[i].mean;
1695   }
1696   channels=3;
1697   if (image->matte != MagickFalse)
1698     channels++;
1699   if (image->colorspace == CMYKColorspace)
1700     channels++;
1701   channel_statistics[CompositeChannels].sum/=channels;
1702   channel_statistics[CompositeChannels].sum_squared/=channels;
1703   channel_statistics[CompositeChannels].sum_cubed/=channels;
1704   channel_statistics[CompositeChannels].sum_fourth_power/=channels;
1705   channel_statistics[CompositeChannels].mean/=channels;
1706   channel_statistics[CompositeChannels].variance/=channels;
1707   channel_statistics[CompositeChannels].standard_deviation=
1708     sqrt(channel_statistics[CompositeChannels].standard_deviation/channels);
1709   channel_statistics[CompositeChannels].kurtosis/=channels;
1710   channel_statistics[CompositeChannels].skewness/=channels;
1711   for (i=0; i <= (ssize_t) CompositeChannels; i++)
1712   {
1713     if (channel_statistics[i].standard_deviation == 0.0)
1714       continue;
1715     channel_statistics[i].skewness=(channel_statistics[i].sum_cubed-
1716       3.0*channel_statistics[i].mean*channel_statistics[i].sum_squared+
1717       2.0*channel_statistics[i].mean*channel_statistics[i].mean*
1718       channel_statistics[i].mean)/(channel_statistics[i].standard_deviation*
1719       channel_statistics[i].standard_deviation*
1720       channel_statistics[i].standard_deviation);
1721     channel_statistics[i].kurtosis=(channel_statistics[i].sum_fourth_power-
1722       4.0*channel_statistics[i].mean*channel_statistics[i].sum_cubed+
1723       6.0*channel_statistics[i].mean*channel_statistics[i].mean*
1724       channel_statistics[i].sum_squared-3.0*channel_statistics[i].mean*
1725       channel_statistics[i].mean*1.0*channel_statistics[i].mean*
1726       channel_statistics[i].mean)/(channel_statistics[i].standard_deviation*
1727       channel_statistics[i].standard_deviation*
1728       channel_statistics[i].standard_deviation*
1729       channel_statistics[i].standard_deviation)-3.0;
1730   }
1731   return(channel_statistics);
1732 }