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