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