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,sizeof(*pixels));
 165   if (pixels == (PixelChannels **) NULL)
 166     return((PixelChannels **) 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]=(PixelChannels *) AcquireQuantumMemory(length,sizeof(**pixels));
 177     if (pixels[i] == (PixelChannels *) 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 PixelChannels
 205     *color_1,
 206     *color_2;
 207
 208   MagickRealType
 209     distance;
 210
 211   register ssize_t
 212     i;
 213
 214   color_1=(const PixelChannels *) x;
 215   color_2=(const PixelChannels *) 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   PixelChannels
 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 == (PixelChannels **) 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     {
 496 #if   defined(MAGICKCORE_OPENMP_SUPPORT)
 497       #pragma omp parallel for schedule(dynamic) shared(progress,status)
 498 #endif
 499       for (y=0; y < (ssize_t) evaluate_image->rows; y++)
 500       {
 501         CacheView
 502           *image_view;
 503
 504         const Image
 505           *next;
 506
 507         const int
 508           id = GetOpenMPThreadId();
 509
 510         register PixelChannels
 511           *evaluate_pixel;
 512
 513         register Quantum
 514           *restrict q;
 515
 516         register ssize_t
 517           x;
 518
 519         if (status == MagickFalse)
 520           continue;
 521         q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,
 522           evaluate_image->columns,1,exception);
 523         if (q == (Quantum *) NULL)
 524           {
 525             status=MagickFalse;
 526             continue;
 527           }
 528         evaluate_pixel=evaluate_pixels[id];
 529         for (x=0; x < (ssize_t) evaluate_image->columns; x++)
 530         {
 531           register ssize_t
 532             j,
 533             k;
 534
 535           for (j=0; j < (ssize_t) number_images; j++)
 536             for (k=0; k < MaxPixelChannels; k++)
 537               evaluate_pixel[j].channel[k]=0.0;
 538           next=images;
 539           for (j=0; j < (ssize_t) number_images; j++)
 540           {
 541             register const Quantum
 542               *p;
 543
 544             register ssize_t
 545               i;
 546
 547             image_view=AcquireCacheView(next);
 548             p=GetCacheViewVirtualPixels(image_view,x,y,1,1,exception);
 549             if (p == (const Quantum *) NULL)
 550               {
 551                 image_view=DestroyCacheView(image_view);
 552                 break;
 553               }
 554             for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++)
 555             {
 556               PixelChannel
 557                 channel;
 558
 559               PixelTrait
 560                 evaluate_traits,
 561                 traits;
 562
 563               evaluate_traits=GetPixelChannelMapTraits(evaluate_image,
 564                 (PixelChannel) i);
 565               channel=GetPixelChannelMapChannel(evaluate_image,
 566                 (PixelChannel) i);
 567               traits=GetPixelChannelMapTraits(next,channel);
 568               if ((traits == UndefinedPixelTrait) ||
 569                   (evaluate_traits == UndefinedPixelTrait))
 570                 continue;
 571               if ((evaluate_traits & UpdatePixelTrait) == 0)
 572                 continue;
 573               evaluate_pixel[j].channel[i]=ApplyEvaluateOperator(
 574                 random_info[id],GetPixelChannel(evaluate_image,channel,p),op,
 575                 evaluate_pixel[j].channel[i]);
 576             }
 577             image_view=DestroyCacheView(image_view);
 578             next=GetNextImageInList(next);
 579           }
 580           qsort((void *) evaluate_pixel,number_images,sizeof(*evaluate_pixel),
 581             IntensityCompare);
 582           for (k=0; k < (ssize_t) GetPixelChannels(evaluate_image); k++)
 583             q[k]=ClampToQuantum(evaluate_pixel[j/2].channel[k]);
 584           q+=GetPixelChannels(evaluate_image);
 585         }
 586         if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
 587           status=MagickFalse;
 588         if (images->progress_monitor != (MagickProgressMonitor) NULL)
 589           {
 590             MagickBooleanType
 591               proceed;
 592
 593 #if   defined(MAGICKCORE_OPENMP_SUPPORT)
 594             #pragma omp critical (MagickCore_EvaluateImages)
 595 #endif
 596             proceed=SetImageProgress(images,EvaluateImageTag,progress++,
 597               evaluate_image->rows);
 598             if (proceed == MagickFalse)
 599               status=MagickFalse;
 600           }
 601       }
 602     }
 603   else
 604     {
 605 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 606       #pragma omp parallel for schedule(dynamic) shared(progress,status)
 607 #endif
 608       for (y=0; y < (ssize_t) evaluate_image->rows; y++)
 609       {
 610         CacheView
 611           *image_view;
 612
 613         const Image
 614           *next;
 615
 616         const int
 617           id = GetOpenMPThreadId();
 618
 619         register ssize_t
 620           i,
 621           x;
 622
 623         register PixelChannels
 624           *evaluate_pixel;
 625
 626         register Quantum
 627           *restrict q;
 628
 629         ssize_t
 630           j;
 631
 632         if (status == MagickFalse)
 633           continue;
 634         q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,
 635           evaluate_image->columns,1,exception);
 636         if (q == (Quantum *) NULL)
 637           {
 638             status=MagickFalse;
 639             continue;
 640           }
 641         evaluate_pixel=evaluate_pixels[id];
 642         for (j=0; j < (ssize_t) evaluate_image->columns; j++)
 643           for (i=0; i < MaxPixelChannels; i++)
 644             evaluate_pixel[j].channel[i]=0.0;
 645         next=images;
 646         for (j=0; j < (ssize_t) number_images; j++)
 647         {
 648           register const Quantum
 649             *p;
 650
 651           image_view=AcquireCacheView(next);
 652           p=GetCacheViewVirtualPixels(image_view,0,y,next->columns,1,exception);
 653           if (p == (const Quantum *) NULL)
 654             {
 655               image_view=DestroyCacheView(image_view);
 656               break;
 657             }
 658           for (x=0; x < (ssize_t) next->columns; x++)
 659           {
 660             register ssize_t
 661               i;
 662
 663             for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++)
 664             {
 665               PixelChannel
 666                 channel;
 667
 668               PixelTrait
 669                 evaluate_traits,
 670                 traits;
 671
 672               evaluate_traits=GetPixelChannelMapTraits(evaluate_image,
 673                 (PixelChannel) i);
 674               channel=GetPixelChannelMapChannel(evaluate_image,(PixelChannel)
 675                 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),
 684                 j == 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,(PixelChannel) 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,(PixelChannel) 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
 906         (e.g. c0*x^3 + 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):
 943         width, center, 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,(PixelChannel) 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,(PixelChannel) 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=channel_statistics[CompositePixelChannel].standard_deviation;
1195   channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1196     channel_statistics);
1197   return(MagickTrue);
1198 }
1199 \f
1200 /*
1201 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1202 %                                                                             %
1203 %                                                                             %
1204 %                                                                             %
1205 %   G e t I m a g e K u r t o s i s                                           %
1206 %                                                                             %
1207 %                                                                             %
1208 %                                                                             %
1209 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1210 %
1211 %  GetImageKurtosis() returns the kurtosis and skewness of one or more
1212 %  image channels.
1213 %
1214 %  The format of the GetImageKurtosis method is:
1215 %
1216 %      MagickBooleanType GetImageKurtosis(const Image *image,double *kurtosis,
1217 %        double *skewness,ExceptionInfo *exception)
1218 %
1219 %  A description of each parameter follows:
1220 %
1221 %    o image: the image.
1222 %
1223 %    o kurtosis: the kurtosis of the channel.
1224 %
1225 %    o skewness: the skewness of the channel.
1226 %
1227 %    o exception: return any errors or warnings in this structure.
1228 %
1229 */
1230 MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
1231   double *kurtosis,double *skewness,ExceptionInfo *exception)
1232 {
1233   CacheView
1234     *image_view;
1235
1236   double
1237     area,
1238     mean,
1239     standard_deviation,
1240     sum_squares,
1241     sum_cubes,
1242     sum_fourth_power;
1243
1244   MagickBooleanType
1245     status;
1246
1247   ssize_t
1248     y;
1249
1250   assert(image != (Image *) NULL);
1251   assert(image->signature == MagickSignature);
1252   if (image->debug != MagickFalse)
1253     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1254   status=MagickTrue;
1255   *kurtosis=0.0;
1256   *skewness=0.0;
1257   area=0.0;
1258   mean=0.0;
1259   standard_deviation=0.0;
1260   sum_squares=0.0;
1261   sum_cubes=0.0;
1262   sum_fourth_power=0.0;
1263   image_view=AcquireCacheView(image);
1264 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1265   #pragma omp parallel for schedule(dynamic) shared(status) omp_throttle(1)
1266 #endif
1267   for (y=0; y < (ssize_t) image->rows; y++)
1268   {
1269     register const Quantum
1270       *restrict p;
1271
1272     register ssize_t
1273       x;
1274
1275     if (status == MagickFalse)
1276       continue;
1277     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1278     if (p == (const Quantum *) NULL)
1279       {
1280         status=MagickFalse;
1281         continue;
1282       }
1283     for (x=0; x < (ssize_t) image->columns; x++)
1284     {
1285       register ssize_t
1286         i;
1287
1288       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1289       {
1290         PixelTrait
1291           traits;
1292
1293         traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
1294         if (traits == UndefinedPixelTrait)
1295           continue;
1296         if ((traits & UpdatePixelTrait) == 0)
1297           continue;
1298 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1299         #pragma omp critical (MagickCore_GetImageKurtosis)
1300 #endif
1301         {
1302           mean+=p[i];
1303           sum_squares+=(double) p[i]*p[i];
1304           sum_cubes+=(double) p[i]*p[i]*p[i];
1305           sum_fourth_power+=(double) p[i]*p[i]*p[i]*p[i];
1306           area++;
1307         }
1308       }
1309       p+=GetPixelChannels(image);
1310     }
1311   }
1312   image_view=DestroyCacheView(image_view);
1313   if (area != 0.0)
1314     {
1315       mean/=area;
1316       sum_squares/=area;
1317       sum_cubes/=area;
1318       sum_fourth_power/=area;
1319     }
1320   standard_deviation=sqrt(sum_squares-(mean*mean));
1321   if (standard_deviation != 0.0)
1322     {
1323       *kurtosis=sum_fourth_power-4.0*mean*sum_cubes+6.0*mean*mean*sum_squares-
1324         3.0*mean*mean*mean*mean;
1325       *kurtosis/=standard_deviation*standard_deviation*standard_deviation*
1326         standard_deviation;
1327       *kurtosis-=3.0;
1328       *skewness=sum_cubes-3.0*mean*sum_squares+2.0*mean*mean*mean;
1329       *skewness/=standard_deviation*standard_deviation*standard_deviation;
1330     }
1331   return(status);
1332 }
1333 \f
1334 /*
1335 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1336 %                                                                             %
1337 %                                                                             %
1338 %                                                                             %
1339 %   G e t I m a g e R a n g e                                                 %
1340 %                                                                             %
1341 %                                                                             %
1342 %                                                                             %
1343 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1344 %
1345 %  GetImageRange() returns the range of one or more image channels.
1346 %
1347 %  The format of the GetImageRange method is:
1348 %
1349 %      MagickBooleanType GetImageRange(const Image *image,double *minima,
1350 %        double *maxima,ExceptionInfo *exception)
1351 %
1352 %  A description of each parameter follows:
1353 %
1354 %    o image: the image.
1355 %
1356 %    o minima: the minimum value in the channel.
1357 %
1358 %    o maxima: the maximum value in the channel.
1359 %
1360 %    o exception: return any errors or warnings in this structure.
1361 %
1362 */
1363 MagickExport MagickBooleanType GetImageRange(const Image *image,double *minima,
1364   double *maxima,ExceptionInfo *exception)
1365 {
1366   CacheView
1367     *image_view;
1368
1369   MagickBooleanType
1370     status;
1371
1372   ssize_t
1373     y;
1374
1375   assert(image != (Image *) NULL);
1376   assert(image->signature == MagickSignature);
1377   if (image->debug != MagickFalse)
1378     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1379   status=MagickTrue;
1380   *maxima=(-MagickHuge);
1381   *minima=MagickHuge;
1382   image_view=AcquireCacheView(image);
1383 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1384   #pragma omp parallel for schedule(dynamic) shared(status) omp_throttle(1)
1385 #endif
1386   for (y=0; y < (ssize_t) image->rows; y++)
1387   {
1388     register const Quantum
1389       *restrict p;
1390
1391     register ssize_t
1392       x;
1393
1394     if (status == MagickFalse)
1395       continue;
1396     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1397     if (p == (const Quantum *) NULL)
1398       {
1399         status=MagickFalse;
1400         continue;
1401       }
1402     for (x=0; x < (ssize_t) image->columns; x++)
1403     {
1404       register ssize_t
1405         i;
1406
1407       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1408       {
1409         PixelTrait
1410           traits;
1411
1412         traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
1413         if (traits == UndefinedPixelTrait)
1414           continue;
1415         if ((traits & UpdatePixelTrait) == 0)
1416           continue;
1417 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1418         #pragma omp critical (MagickCore_GetImageRange)
1419 #endif
1420         {
1421           if (p[i] < *minima)
1422             *minima=(double) p[i];
1423           if (p[i] > *maxima)
1424             *maxima=(double) p[i];
1425         }
1426       }
1427       p+=GetPixelChannels(image);
1428     }
1429   }
1430   image_view=DestroyCacheView(image_view);
1431   return(status);
1432 }
1433 \f
1434 /*
1435 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1436 %                                                                             %
1437 %                                                                             %
1438 %                                                                             %
1439 %   G e t I m a g e S t a t i s t i c s                                       %
1440 %                                                                             %
1441 %                                                                             %
1442 %                                                                             %
1443 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1444 %
1445 %  GetImageStatistics() returns statistics for each channel in the
1446 %  image.  The statistics include the channel depth, its minima, maxima, mean,
1447 %  standard deviation, kurtosis and skewness.  You can access the red channel
1448 %  mean, for example, like this:
1449 %
1450 %      channel_statistics=GetImageStatistics(image,exception);
1451 %      red_mean=channel_statistics[RedPixelChannel].mean;
1452 %
1453 %  Use MagickRelinquishMemory() to free the statistics buffer.
1454 %
1455 %  The format of the GetImageStatistics method is:
1456 %
1457 %      ChannelStatistics *GetImageStatistics(const Image *image,
1458 %        ExceptionInfo *exception)
1459 %
1460 %  A description of each parameter follows:
1461 %
1462 %    o image: the image.
1463 %
1464 %    o exception: return any errors or warnings in this structure.
1465 %
1466 */
1467
1468 static size_t GetImageChannels(const Image *image)
1469 {
1470   register ssize_t
1471     i;
1472
1473   size_t
1474     channels;
1475
1476   channels=0;
1477   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1478   {
1479     PixelTrait
1480       traits;
1481
1482     traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
1483     if ((traits & UpdatePixelTrait) != 0)
1484       channels++;
1485   }
1486   return(channels);
1487 }
1488
1489 MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
1490   ExceptionInfo *exception)
1491 {
1492   ChannelStatistics
1493     *channel_statistics;
1494
1495   double
1496     area;
1497
1498   MagickStatusType
1499     status;
1500
1501   QuantumAny
1502     range;
1503
1504   register ssize_t
1505     i;
1506
1507   size_t
1508     channels,
1509     depth;
1510
1511   ssize_t
1512     y;
1513
1514   assert(image != (Image *) NULL);
1515   assert(image->signature == MagickSignature);
1516   if (image->debug != MagickFalse)
1517     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1518   channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(
1519     MaxPixelChannels+1,sizeof(*channel_statistics));
1520   if (channel_statistics == (ChannelStatistics *) NULL)
1521     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1522   (void) ResetMagickMemory(channel_statistics,0,(MaxPixelChannels+1)*
1523     sizeof(*channel_statistics));
1524   for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
1525   {
1526     channel_statistics[i].depth=1;
1527     channel_statistics[i].maxima=(-MagickHuge);
1528     channel_statistics[i].minima=MagickHuge;
1529   }
1530   for (y=0; y < (ssize_t) image->rows; y++)
1531   {
1532     register const Quantum
1533       *restrict p;
1534
1535     register ssize_t
1536       x;
1537
1538     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1539     if (p == (const Quantum *) NULL)
1540       break;
1541     for (x=0; x < (ssize_t) image->columns; x++)
1542     {
1543       register ssize_t
1544         i;
1545
1546       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1547       {
1548         PixelTrait
1549           traits;
1550
1551         traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
1552         if (traits == UndefinedPixelTrait)
1553           continue;
1554         if (channel_statistics[i].depth != MAGICKCORE_QUANTUM_DEPTH)
1555           {
1556             depth=channel_statistics[i].depth;
1557             range=GetQuantumRange(depth);
1558             status=p[i] != ScaleAnyToQuantum(ScaleQuantumToAny(p[i],range),
1559               range) ? MagickTrue : MagickFalse;
1560             if (status != MagickFalse)
1561               {
1562                 channel_statistics[i].depth++;
1563                 continue;
1564               }
1565           }
1566         if ((double) p[i] < channel_statistics[i].minima)
1567           channel_statistics[i].minima=(double) p[i];
1568         if ((double) p[i] > channel_statistics[i].maxima)
1569           channel_statistics[i].maxima=(double) p[i];
1570         channel_statistics[i].sum+=p[i];
1571         channel_statistics[i].sum_squared+=(double) p[i]*p[i];
1572         channel_statistics[i].sum_cubed+=(double) p[i]*p[i]*p[i];
1573         channel_statistics[i].sum_fourth_power+=(double) p[i]*p[i]*p[i]*p[i];
1574       }
1575       p+=GetPixelChannels(image);
1576     }
1577   }
1578   area=(double) image->columns*image->rows;
1579   for (i=0; i < (ssize_t) MaxPixelChannels; i++)
1580   {
1581     channel_statistics[i].sum/=area;
1582     channel_statistics[i].sum_squared/=area;
1583     channel_statistics[i].sum_cubed/=area;
1584     channel_statistics[i].sum_fourth_power/=area;
1585     channel_statistics[i].mean=channel_statistics[i].sum;
1586     channel_statistics[i].variance=channel_statistics[i].sum_squared;
1587     channel_statistics[i].standard_deviation=sqrt(
1588       channel_statistics[i].variance-(channel_statistics[i].mean*
1589       channel_statistics[i].mean));
1590   }
1591   for (i=0; i < (ssize_t) MaxPixelChannels; i++)
1592   {
1593     channel_statistics[CompositePixelChannel].depth=(size_t) MagickMax((double)
1594       channel_statistics[CompositePixelChannel].depth,(double)
1595       channel_statistics[i].depth);
1596     channel_statistics[CompositePixelChannel].minima=MagickMin(
1597       channel_statistics[CompositePixelChannel].minima,
1598       channel_statistics[i].minima);
1599     channel_statistics[CompositePixelChannel].maxima=MagickMax(
1600       channel_statistics[CompositePixelChannel].maxima,
1601       channel_statistics[i].maxima);
1602     channel_statistics[CompositePixelChannel].sum+=channel_statistics[i].sum;
1603     channel_statistics[CompositePixelChannel].sum_squared+=
1604       channel_statistics[i].sum_squared;
1605     channel_statistics[CompositePixelChannel].sum_cubed+=
1606       channel_statistics[i].sum_cubed;
1607     channel_statistics[CompositePixelChannel].sum_fourth_power+=
1608       channel_statistics[i].sum_fourth_power;
1609     channel_statistics[CompositePixelChannel].mean+=channel_statistics[i].mean;
1610     channel_statistics[CompositePixelChannel].variance+=
1611       channel_statistics[i].variance-channel_statistics[i].mean*
1612       channel_statistics[i].mean;
1613     channel_statistics[CompositePixelChannel].standard_deviation+=
1614       channel_statistics[i].variance-channel_statistics[i].mean*
1615       channel_statistics[i].mean;
1616   }
1617   channels=GetImageChannels(image);
1618   channel_statistics[CompositePixelChannel].sum/=channels;
1619   channel_statistics[CompositePixelChannel].sum_squared/=channels;
1620   channel_statistics[CompositePixelChannel].sum_cubed/=channels;
1621   channel_statistics[CompositePixelChannel].sum_fourth_power/=channels;
1622   channel_statistics[CompositePixelChannel].mean/=channels;
1623   channel_statistics[CompositePixelChannel].variance/=channels;
1624   channel_statistics[CompositePixelChannel].standard_deviation=
1625     sqrt(channel_statistics[CompositePixelChannel].standard_deviation/channels);
1626   channel_statistics[CompositePixelChannel].kurtosis/=channels;
1627   channel_statistics[CompositePixelChannel].skewness/=channels;
1628   for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
1629   {
1630     if (channel_statistics[i].standard_deviation == 0.0)
1631       continue;
1632     channel_statistics[i].skewness=(channel_statistics[i].sum_cubed-
1633       3.0*channel_statistics[i].mean*channel_statistics[i].sum_squared+
1634       2.0*channel_statistics[i].mean*channel_statistics[i].mean*
1635       channel_statistics[i].mean)/(channel_statistics[i].standard_deviation*
1636       channel_statistics[i].standard_deviation*
1637       channel_statistics[i].standard_deviation);
1638     channel_statistics[i].kurtosis=(channel_statistics[i].sum_fourth_power-
1639       4.0*channel_statistics[i].mean*channel_statistics[i].sum_cubed+
1640       6.0*channel_statistics[i].mean*channel_statistics[i].mean*
1641       channel_statistics[i].sum_squared-3.0*channel_statistics[i].mean*
1642       channel_statistics[i].mean*1.0*channel_statistics[i].mean*
1643       channel_statistics[i].mean)/(channel_statistics[i].standard_deviation*
1644       channel_statistics[i].standard_deviation*
1645       channel_statistics[i].standard_deviation*
1646       channel_statistics[i].standard_deviation)-3.0;
1647   }
1648   return(channel_statistics);
1649 }