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