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