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