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           op,evaluate_pixel[x].red);
 517         evaluate_pixel[x].green=ApplyEvaluateOperator(random_info[id],p->green,
 518           op,evaluate_pixel[x].green);
 519         evaluate_pixel[x].blue=ApplyEvaluateOperator(random_info[id],p->blue,
 520           op,evaluate_pixel[x].blue);
 521         evaluate_pixel[x].opacity=ApplyEvaluateOperator(random_info[id],
 522           p->opacity,op,evaluate_pixel[x].opacity);
 523         if (evaluate_image->colorspace == CMYKColorspace)
 524           evaluate_pixel[x].index=ApplyEvaluateOperator(random_info[id],
 525             indexes[x],op,evaluate_pixel[x].index);
 526         p++;
 527       }
 528       image_view=DestroyCacheView(image_view);
 529       next=GetNextImageInList(next);
 530     }
 531     for (x=0; x < (long) evaluate_image->columns; x++)
 532     {
 533       q->red=ClampToQuantum(evaluate_pixel[x].red);
 534       q->green=ClampToQuantum(evaluate_pixel[x].green);
 535       q->blue=ClampToQuantum(evaluate_pixel[x].blue);
 536       q->opacity=ClampToQuantum(evaluate_pixel[x].opacity);
 537       if (evaluate_image->colorspace == CMYKColorspace)
 538         evaluate_indexes[x]=ClampToQuantum(evaluate_pixel[x].index);
 539       q++;
 540     }
 541     if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
 542       status=MagickFalse;
 543     if (images->progress_monitor != (MagickProgressMonitor) NULL)
 544       {
 545         MagickBooleanType
 546           proceed;
 547
 548 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 549         #pragma omp critical (MagickCore_EvaluateImages)
 550 #endif
 551         proceed=SetImageProgress(images,EvaluateImageTag,progress++,
 552           evaluate_image->rows);
 553         if (proceed == MagickFalse)
 554           status=MagickFalse;
 555       }
 556   }
 557   evaluate_view=DestroyCacheView(evaluate_view);
 558   evaluate_pixels=DestroyPixelThreadSet(evaluate_pixels);
 559   random_info=DestroyRandomInfoThreadSet(random_info);
 560   if (status == MagickFalse)
 561     evaluate_image=DestroyImage(evaluate_image);
 562   return(evaluate_image);
 563 }
 564
 565 MagickExport MagickBooleanType EvaluateImageChannel(Image *image,
 566   const ChannelType channel,const MagickEvaluateOperator op,const double value,
 567   ExceptionInfo *exception)
 568 {
 569   CacheView
 570     *image_view;
 571
 572   long
 573     progress,
 574     y;
 575
 576   MagickBooleanType
 577     status;
 578
 579   RandomInfo
 580     **restrict random_info;
 581
 582   assert(image != (Image *) NULL);
 583   assert(image->signature == MagickSignature);
 584   if (image->debug != MagickFalse)
 585     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 586   assert(exception != (ExceptionInfo *) NULL);
 587   assert(exception->signature == MagickSignature);
 588   if (SetImageStorageClass(image,DirectClass) == MagickFalse)
 589     {
 590       InheritException(exception,&image->exception);
 591       return(MagickFalse);
 592     }
 593   status=MagickTrue;
 594   progress=0;
 595   random_info=AcquireRandomInfoThreadSet();
 596   image_view=AcquireCacheView(image);
 597 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 598   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
 599 #endif
 600   for (y=0; y < (long) image->rows; y++)
 601   {
 602     register IndexPacket
 603       *restrict indexes;
 604
 605     register long
 606       id,
 607       x;
 608
 609     register PixelPacket
 610       *restrict q;
 611
 612     if (status == MagickFalse)
 613       continue;
 614     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
 615     if (q == (PixelPacket *) NULL)
 616       {
 617         status=MagickFalse;
 618         continue;
 619       }
 620     indexes=GetCacheViewAuthenticIndexQueue(image_view);
 621     id=GetOpenMPThreadId();
 622     for (x=0; x < (long) image->columns; x++)
 623     {
 624       if ((channel & RedChannel) != 0)
 625         q->red=ClampToQuantum(ApplyEvaluateOperator(random_info[id],q->red,op,
 626           value));
 627       if ((channel & GreenChannel) != 0)
 628         q->green=ClampToQuantum(ApplyEvaluateOperator(random_info[id],q->green,
 629           op,value));
 630       if ((channel & BlueChannel) != 0)
 631         q->blue=ClampToQuantum(ApplyEvaluateOperator(random_info[id],q->blue,op,
 632           value));
 633       if ((channel & OpacityChannel) != 0)
 634         {
 635           if (image->matte == MagickFalse)
 636             q->opacity=ClampToQuantum(ApplyEvaluateOperator(random_info[id],
 637               q->opacity,op,value));
 638           else
 639             q->opacity=ClampToQuantum(QuantumRange-ApplyEvaluateOperator(
 640               random_info[id],(Quantum) GetAlphaPixelComponent(q),op,value));
 641         }
 642       if (((channel & IndexChannel) != 0) && (indexes != (IndexPacket *) NULL))
 643         indexes[x]=(IndexPacket) ClampToQuantum(ApplyEvaluateOperator(
 644           random_info[id],indexes[x],op,value));
 645       q++;
 646     }
 647     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
 648       status=MagickFalse;
 649     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 650       {
 651         MagickBooleanType
 652           proceed;
 653
 654 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 655   #pragma omp critical (MagickCore_EvaluateImageChannel)
 656 #endif
 657         proceed=SetImageProgress(image,EvaluateImageTag,progress++,image->rows);
 658         if (proceed == MagickFalse)
 659           status=MagickFalse;
 660       }
 661   }
 662   image_view=DestroyCacheView(image_view);
 663   random_info=DestroyRandomInfoThreadSet(random_info);
 664   return(status);
 665 }
 666 \f
 667 /*
 668 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 669 %                                                                             %
 670 %                                                                             %
 671 %                                                                             %
 672 %     F u n c t i o n I m a g e                                               %
 673 %                                                                             %
 674 %                                                                             %
 675 %                                                                             %
 676 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 677 %
 678 %  FunctionImage() applies a value to the image with an arithmetic, relational,
 679 %  or logical operator to an image. Use these operations to lighten or darken
 680 %  an image, to increase or decrease contrast in an image, or to produce the
 681 %  "negative" of an image.
 682 %
 683 %  The format of the FunctionImageChannel method is:
 684 %
 685 %      MagickBooleanType FunctionImage(Image *image,
 686 %        const MagickFunction function,const long number_parameters,
 687 %        const double *parameters,ExceptionInfo *exception)
 688 %      MagickBooleanType FunctionImageChannel(Image *image,
 689 %        const ChannelType channel,const MagickFunction function,
 690 %        const long number_parameters,const double *argument,
 691 %        ExceptionInfo *exception)
 692 %
 693 %  A description of each parameter follows:
 694 %
 695 %    o image: the image.
 696 %
 697 %    o channel: the channel.
 698 %
 699 %    o function: A channel function.
 700 %
 701 %    o parameters: one or more parameters.
 702 %
 703 %    o exception: return any errors or warnings in this structure.
 704 %
 705 */
 706
 707 static Quantum ApplyFunction(Quantum pixel,const MagickFunction function,
 708   const unsigned long number_parameters,const double *parameters,
 709   ExceptionInfo *exception)
 710 {
 711   MagickRealType
 712     result;
 713
 714   register long
 715     i;
 716
 717   (void) exception;
 718   result=0.0;
 719   switch (function)
 720   {
 721     case PolynomialFunction:
 722     {
 723       /*
 724        * Polynomial
 725        * Parameters:   polynomial constants,  highest to lowest order
 726        *   For example:      c0*x^3 + c1*x^2 + c2*x  + c3
 727        */
 728       result=0.0;
 729       for (i=0; i < (long) number_parameters; i++)
 730         result = result*QuantumScale*pixel + parameters[i];
 731       result *= QuantumRange;
 732       break;
 733     }
 734     case SinusoidFunction:
 735     {
 736       /* Sinusoid Function
 737        * Parameters:   Freq, Phase, Ampl, bias
 738        */
 739       double  freq,phase,ampl,bias;
 740       freq  = ( number_parameters >= 1 ) ? parameters[0] : 1.0;
 741       phase = ( number_parameters >= 2 ) ? parameters[1] : 0.0;
 742       ampl  = ( number_parameters >= 3 ) ? parameters[2] : 0.5;
 743       bias  = ( number_parameters >= 4 ) ? parameters[3] : 0.5;
 744       result=(MagickRealType) (QuantumRange*(ampl*sin((double) (2.0*MagickPI*
 745         (freq*QuantumScale*pixel + phase/360.0) )) + bias ) );
 746       break;
 747     }
 748     case ArcsinFunction:
 749     {
 750       /* Arcsin Function  (peged at range limits for invalid results)
 751        * Parameters:   Width, Center, Range, Bias
 752        */
 753       double  width,range,center,bias;
 754       width  = ( number_parameters >= 1 ) ? parameters[0] : 1.0;
 755       center = ( number_parameters >= 2 ) ? parameters[1] : 0.5;
 756       range  = ( number_parameters >= 3 ) ? parameters[2] : 1.0;
 757       bias   = ( number_parameters >= 4 ) ? parameters[3] : 0.5;
 758       result = 2.0/width*(QuantumScale*pixel - center);
 759       if ( result <= -1.0 )
 760         result = bias - range/2.0;
 761       else if ( result >= 1.0 )
 762         result = bias + range/2.0;
 763       else
 764         result=range/MagickPI*asin((double)result) + bias;
 765       result *= QuantumRange;
 766       break;
 767     }
 768     case ArctanFunction:
 769     {
 770       /* Arctan Function
 771        * Parameters:   Slope, Center, Range, Bias
 772        */
 773       double  slope,range,center,bias;
 774       slope  = ( number_parameters >= 1 ) ? parameters[0] : 1.0;
 775       center = ( number_parameters >= 2 ) ? parameters[1] : 0.5;
 776       range  = ( number_parameters >= 3 ) ? parameters[2] : 1.0;
 777       bias   = ( number_parameters >= 4 ) ? parameters[3] : 0.5;
 778       result = MagickPI*slope*(QuantumScale*pixel - center);
 779       result=(MagickRealType) (QuantumRange*(range/MagickPI*atan((double)
 780                   result) + bias ) );
 781       break;
 782     }
 783     case UndefinedFunction:
 784       break;
 785   }
 786   return(ClampToQuantum(result));
 787 }
 788
 789 MagickExport MagickBooleanType FunctionImage(Image *image,
 790   const MagickFunction function,const unsigned long number_parameters,
 791   const double *parameters,ExceptionInfo *exception)
 792 {
 793   MagickBooleanType
 794     status;
 795
 796   status=FunctionImageChannel(image,AllChannels,function,number_parameters,
 797     parameters,exception);
 798   return(status);
 799 }
 800
 801 MagickExport MagickBooleanType FunctionImageChannel(Image *image,
 802   const ChannelType channel,const MagickFunction function,
 803   const unsigned long number_parameters,const double *parameters,
 804   ExceptionInfo *exception)
 805 {
 806 #define FunctionImageTag  "Function/Image "
 807
 808   CacheView
 809     *image_view;
 810
 811   long
 812     progress,
 813     y;
 814
 815   MagickBooleanType
 816     status;
 817
 818   assert(image != (Image *) NULL);
 819   assert(image->signature == MagickSignature);
 820   if (image->debug != MagickFalse)
 821     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 822   assert(exception != (ExceptionInfo *) NULL);
 823   assert(exception->signature == MagickSignature);
 824   if (SetImageStorageClass(image,DirectClass) == MagickFalse)
 825     {
 826       InheritException(exception,&image->exception);
 827       return(MagickFalse);
 828     }
 829   status=MagickTrue;
 830   progress=0;
 831   image_view=AcquireCacheView(image);
 832 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 833   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
 834 #endif
 835   for (y=0; y < (long) image->rows; y++)
 836   {
 837     register IndexPacket
 838       *restrict indexes;
 839
 840     register long
 841       x;
 842
 843     register PixelPacket
 844       *restrict q;
 845
 846     if (status == MagickFalse)
 847       continue;
 848     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
 849     if (q == (PixelPacket *) NULL)
 850       {
 851         status=MagickFalse;
 852         continue;
 853       }
 854     indexes=GetCacheViewAuthenticIndexQueue(image_view);
 855     for (x=0; x < (long) image->columns; x++)
 856     {
 857       if ((channel & RedChannel) != 0)
 858         q->red=ApplyFunction(q->red,function,number_parameters,parameters,
 859           exception);
 860       if ((channel & GreenChannel) != 0)
 861         q->green=ApplyFunction(q->green,function,number_parameters,parameters,
 862           exception);
 863       if ((channel & BlueChannel) != 0)
 864         q->blue=ApplyFunction(q->blue,function,number_parameters,parameters,
 865           exception);
 866       if ((channel & OpacityChannel) != 0)
 867         {
 868           if (image->matte == MagickFalse)
 869             q->opacity=ApplyFunction(q->opacity,function,number_parameters,
 870               parameters,exception);
 871           else
 872             q->opacity=(Quantum) QuantumRange-ApplyFunction((Quantum)
 873               GetAlphaPixelComponent(q),function,number_parameters,parameters,
 874               exception);
 875         }
 876       if (((channel & IndexChannel) != 0) && (indexes != (IndexPacket *) NULL))
 877         indexes[x]=(IndexPacket) ApplyFunction(indexes[x],function,
 878           number_parameters,parameters,exception);
 879       q++;
 880     }
 881     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
 882       status=MagickFalse;
 883     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 884       {
 885         MagickBooleanType
 886           proceed;
 887
 888 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 889   #pragma omp critical (MagickCore_FunctionImageChannel)
 890 #endif
 891         proceed=SetImageProgress(image,FunctionImageTag,progress++,image->rows);
 892         if (proceed == MagickFalse)
 893           status=MagickFalse;
 894       }
 895   }
 896   image_view=DestroyCacheView(image_view);
 897   return(status);
 898 }
 899 \f
 900 /*
 901 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 902 %                                                                             %
 903 %                                                                             %
 904 %                                                                             %
 905 +   G e t I m a g e C h a n n e l E x t r e m a                               %
 906 %                                                                             %
 907 %                                                                             %
 908 %                                                                             %
 909 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 910 %
 911 %  GetImageChannelExtrema() returns the extrema of one or more image channels.
 912 %
 913 %  The format of the GetImageChannelExtrema method is:
 914 %
 915 %      MagickBooleanType GetImageChannelExtrema(const Image *image,
 916 %        const ChannelType channel,unsigned long *minima,unsigned long *maxima,
 917 %        ExceptionInfo *exception)
 918 %
 919 %  A description of each parameter follows:
 920 %
 921 %    o image: the image.
 922 %
 923 %    o channel: the channel.
 924 %
 925 %    o minima: the minimum value in the channel.
 926 %
 927 %    o maxima: the maximum value in the channel.
 928 %
 929 %    o exception: return any errors or warnings in this structure.
 930 %
 931 */
 932
 933 MagickExport MagickBooleanType GetImageExtrema(const Image *image,
 934   unsigned long *minima,unsigned long *maxima,ExceptionInfo *exception)
 935 {
 936   return(GetImageChannelExtrema(image,AllChannels,minima,maxima,exception));
 937 }
 938
 939 MagickExport MagickBooleanType GetImageChannelExtrema(const Image *image,
 940   const ChannelType channel,unsigned long *minima,unsigned long *maxima,
 941   ExceptionInfo *exception)
 942 {
 943   double
 944     max,
 945     min;
 946
 947   MagickBooleanType
 948     status;
 949
 950   assert(image != (Image *) NULL);
 951   assert(image->signature == MagickSignature);
 952   if (image->debug != MagickFalse)
 953     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 954   status=GetImageChannelRange(image,channel,&min,&max,exception);
 955   *minima=(unsigned long) (min+0.5);
 956   *maxima=(unsigned long) (max+0.5);
 957   return(status);
 958 }
 959 \f
 960 /*
 961 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 962 %                                                                             %
 963 %                                                                             %
 964 %                                                                             %
 965 %   G e t I m a g e C h a n n e l M e a n                                     %
 966 %                                                                             %
 967 %                                                                             %
 968 %                                                                             %
 969 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 970 %
 971 %  GetImageChannelMean() returns the mean and standard deviation of one or more
 972 %  image channels.
 973 %
 974 %  The format of the GetImageChannelMean method is:
 975 %
 976 %      MagickBooleanType GetImageChannelMean(const Image *image,
 977 %        const ChannelType channel,double *mean,double *standard_deviation,
 978 %        ExceptionInfo *exception)
 979 %
 980 %  A description of each parameter follows:
 981 %
 982 %    o image: the image.
 983 %
 984 %    o channel: the channel.
 985 %
 986 %    o mean: the average value in the channel.
 987 %
 988 %    o standard_deviation: the standard deviation of the channel.
 989 %
 990 %    o exception: return any errors or warnings in this structure.
 991 %
 992 */
 993
 994 MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean,
 995   double *standard_deviation,ExceptionInfo *exception)
 996 {
 997   MagickBooleanType
 998     status;
 999
1000   status=GetImageChannelMean(image,AllChannels,mean,standard_deviation,
1001     exception);
1002   return(status);
1003 }
1004
1005 MagickExport MagickBooleanType GetImageChannelMean(const Image *image,
1006   const ChannelType channel,double *mean,double *standard_deviation,
1007   ExceptionInfo *exception)
1008 {
1009   double
1010     area;
1011
1012   long
1013     y;
1014
1015   assert(image != (Image *) NULL);
1016   assert(image->signature == MagickSignature);
1017   if (image->debug != MagickFalse)
1018     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1019   *mean=0.0;
1020   *standard_deviation=0.0;
1021   area=0.0;
1022   for (y=0; y < (long) image->rows; y++)
1023   {
1024     register const IndexPacket
1025       *restrict indexes;
1026
1027     register const PixelPacket
1028       *restrict p;
1029
1030     register long
1031       x;
1032
1033     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1034     if (p == (const PixelPacket *) NULL)
1035       break;
1036     indexes=GetVirtualIndexQueue(image);
1037     for (x=0; x < (long) image->columns; x++)
1038     {
1039       if ((channel & RedChannel) != 0)
1040         {
1041           *mean+=GetRedPixelComponent(p);
1042           *standard_deviation+=(double) p->red*GetRedPixelComponent(p);
1043           area++;
1044         }
1045       if ((channel & GreenChannel) != 0)
1046         {
1047           *mean+=GetGreenPixelComponent(p);
1048           *standard_deviation+=(double) p->green*GetGreenPixelComponent(p);
1049           area++;
1050         }
1051       if ((channel & BlueChannel) != 0)
1052         {
1053           *mean+=GetBluePixelComponent(p);
1054           *standard_deviation+=(double) p->blue*GetBluePixelComponent(p);
1055           area++;
1056         }
1057       if ((channel & OpacityChannel) != 0)
1058         {
1059           *mean+=GetOpacityPixelComponent(p);
1060           *standard_deviation+=(double) p->opacity*GetOpacityPixelComponent(p);
1061           area++;
1062         }
1063       if (((channel & IndexChannel) != 0) &&
1064           (image->colorspace == CMYKColorspace))
1065         {
1066           *mean+=indexes[x];
1067           *standard_deviation+=(double) indexes[x]*indexes[x];
1068           area++;
1069         }
1070       p++;
1071     }
1072   }
1073   if (y < (long) image->rows)
1074     return(MagickFalse);
1075   if (area != 0)
1076     {
1077       *mean/=area;
1078       *standard_deviation/=area;
1079     }
1080   *standard_deviation=sqrt(*standard_deviation-(*mean*(*mean)));
1081   return(y == (long) image->rows ? MagickTrue : MagickFalse);
1082 }
1083 \f
1084 /*
1085 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1086 %                                                                             %
1087 %                                                                             %
1088 %                                                                             %
1089 %   G e t I m a g e C h a n n e l K u r t o s i s                             %
1090 %                                                                             %
1091 %                                                                             %
1092 %                                                                             %
1093 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1094 %
1095 %  GetImageChannelKurtosis() returns the kurtosis and skewness of one or more
1096 %  image channels.
1097 %
1098 %  The format of the GetImageChannelKurtosis method is:
1099 %
1100 %      MagickBooleanType GetImageChannelKurtosis(const Image *image,
1101 %        const ChannelType channel,double *kurtosis,double *skewness,
1102 %        ExceptionInfo *exception)
1103 %
1104 %  A description of each parameter follows:
1105 %
1106 %    o image: the image.
1107 %
1108 %    o channel: the channel.
1109 %
1110 %    o kurtosis: the kurtosis of the channel.
1111 %
1112 %    o skewness: the skewness of the channel.
1113 %
1114 %    o exception: return any errors or warnings in this structure.
1115 %
1116 */
1117
1118 MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
1119   double *kurtosis,double *skewness,ExceptionInfo *exception)
1120 {
1121   MagickBooleanType
1122     status;
1123
1124   status=GetImageChannelKurtosis(image,AllChannels,kurtosis,skewness,
1125     exception);
1126   return(status);
1127 }
1128
1129 MagickExport MagickBooleanType GetImageChannelKurtosis(const Image *image,
1130   const ChannelType channel,double *kurtosis,double *skewness,
1131   ExceptionInfo *exception)
1132 {
1133   double
1134     area,
1135     mean,
1136     standard_deviation,
1137     sum_squares,
1138     sum_cubes,
1139     sum_fourth_power;
1140
1141   long
1142     y;
1143
1144   assert(image != (Image *) NULL);
1145   assert(image->signature == MagickSignature);
1146   if (image->debug != MagickFalse)
1147     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1148   *kurtosis=0.0;
1149   *skewness=0.0;
1150   area=0.0;
1151   mean=0.0;
1152   standard_deviation=0.0;
1153   sum_squares=0.0;
1154   sum_cubes=0.0;
1155   sum_fourth_power=0.0;
1156   for (y=0; y < (long) image->rows; y++)
1157   {
1158     register const IndexPacket
1159       *restrict indexes;
1160
1161     register const PixelPacket
1162       *restrict p;
1163
1164     register long
1165       x;
1166
1167     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1168     if (p == (const PixelPacket *) NULL)
1169       break;
1170     indexes=GetVirtualIndexQueue(image);
1171     for (x=0; x < (long) image->columns; x++)
1172     {
1173       if ((channel & RedChannel) != 0)
1174         {
1175           mean+=GetRedPixelComponent(p);
1176           sum_squares+=(double) p->red*GetRedPixelComponent(p);
1177           sum_cubes+=(double) p->red*p->red*GetRedPixelComponent(p);
1178           sum_fourth_power+=(double) p->red*p->red*p->red*
1179             GetRedPixelComponent(p);
1180           area++;
1181         }
1182       if ((channel & GreenChannel) != 0)
1183         {
1184           mean+=GetGreenPixelComponent(p);
1185           sum_squares+=(double) p->green*GetGreenPixelComponent(p);
1186           sum_cubes+=(double) p->green*p->green*GetGreenPixelComponent(p);
1187           sum_fourth_power+=(double) p->green*p->green*p->green*
1188             GetGreenPixelComponent(p);
1189           area++;
1190         }
1191       if ((channel & BlueChannel) != 0)
1192         {
1193           mean+=GetBluePixelComponent(p);
1194           sum_squares+=(double) p->blue*GetBluePixelComponent(p);
1195           sum_cubes+=(double) p->blue*p->blue*GetBluePixelComponent(p);
1196           sum_fourth_power+=(double) p->blue*p->blue*p->blue*
1197             GetBluePixelComponent(p);
1198           area++;
1199         }
1200       if ((channel & OpacityChannel) != 0)
1201         {
1202           mean+=GetOpacityPixelComponent(p);
1203           sum_squares+=(double) p->opacity*GetOpacityPixelComponent(p);
1204           sum_cubes+=(double) p->opacity*p->opacity*GetOpacityPixelComponent(p);
1205           sum_fourth_power+=(double) p->opacity*p->opacity*p->opacity*
1206             GetOpacityPixelComponent(p);
1207           area++;
1208         }
1209       if (((channel & IndexChannel) != 0) &&
1210           (image->colorspace == CMYKColorspace))
1211         {
1212           mean+=indexes[x];
1213           sum_squares+=(double) indexes[x]*indexes[x];
1214           sum_cubes+=(double) indexes[x]*indexes[x]*indexes[x];
1215           sum_fourth_power+=(double) indexes[x]*indexes[x]*indexes[x]*
1216             indexes[x];
1217           area++;
1218         }
1219       p++;
1220     }
1221   }
1222   if (y < (long) image->rows)
1223     return(MagickFalse);
1224   if (area != 0.0)
1225     {
1226       mean/=area;
1227       sum_squares/=area;
1228       sum_cubes/=area;
1229       sum_fourth_power/=area;
1230     }
1231   standard_deviation=sqrt(sum_squares-(mean*mean));
1232   if (standard_deviation != 0.0)
1233     {
1234       *kurtosis=sum_fourth_power-4.0*mean*sum_cubes+6.0*mean*mean*sum_squares-
1235         3.0*mean*mean*mean*mean;
1236       *kurtosis/=standard_deviation*standard_deviation*standard_deviation*
1237         standard_deviation;
1238       *kurtosis-=3.0;
1239       *skewness=sum_cubes-3.0*mean*sum_squares+2.0*mean*mean*mean;
1240       *skewness/=standard_deviation*standard_deviation*standard_deviation;
1241     }
1242   return(y == (long) image->rows ? MagickTrue : MagickFalse);
1243 }
1244 \f
1245 /*
1246 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1247 %                                                                             %
1248 %                                                                             %
1249 %                                                                             %
1250 %   G e t I m a g e C h a n n e l R a n g e                                   %
1251 %                                                                             %
1252 %                                                                             %
1253 %                                                                             %
1254 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1255 %
1256 %  GetImageChannelRange() returns the range of one or more image channels.
1257 %
1258 %  The format of the GetImageChannelRange method is:
1259 %
1260 %      MagickBooleanType GetImageChannelRange(const Image *image,
1261 %        const ChannelType channel,double *minima,double *maxima,
1262 %        ExceptionInfo *exception)
1263 %
1264 %  A description of each parameter follows:
1265 %
1266 %    o image: the image.
1267 %
1268 %    o channel: the channel.
1269 %
1270 %    o minima: the minimum value in the channel.
1271 %
1272 %    o maxima: the maximum value in the channel.
1273 %
1274 %    o exception: return any errors or warnings in this structure.
1275 %
1276 */
1277
1278 MagickExport MagickBooleanType GetImageRange(const Image *image,
1279   double *minima,double *maxima,ExceptionInfo *exception)
1280 {
1281   return(GetImageChannelRange(image,AllChannels,minima,maxima,exception));
1282 }
1283
1284 MagickExport MagickBooleanType GetImageChannelRange(const Image *image,
1285   const ChannelType channel,double *minima,double *maxima,
1286   ExceptionInfo *exception)
1287 {
1288   long
1289     y;
1290
1291   MagickPixelPacket
1292     pixel;
1293
1294   assert(image != (Image *) NULL);
1295   assert(image->signature == MagickSignature);
1296   if (image->debug != MagickFalse)
1297     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1298   *maxima=(-1.0E-37);
1299   *minima=1.0E+37;
1300   GetMagickPixelPacket(image,&pixel);
1301   for (y=0; y < (long) image->rows; y++)
1302   {
1303     register const IndexPacket
1304       *restrict indexes;
1305
1306     register const PixelPacket
1307       *restrict p;
1308
1309     register long
1310       x;
1311
1312     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1313     if (p == (const PixelPacket *) NULL)
1314       break;
1315     indexes=GetVirtualIndexQueue(image);
1316     for (x=0; x < (long) image->columns; x++)
1317     {
1318       SetMagickPixelPacket(image,p,indexes+x,&pixel);
1319       if ((channel & RedChannel) != 0)
1320         {
1321           if (pixel.red < *minima)
1322             *minima=(double) pixel.red;
1323           if (pixel.red > *maxima)
1324             *maxima=(double) pixel.red;
1325         }
1326       if ((channel & GreenChannel) != 0)
1327         {
1328           if (pixel.green < *minima)
1329             *minima=(double) pixel.green;
1330           if (pixel.green > *maxima)
1331             *maxima=(double) pixel.green;
1332         }
1333       if ((channel & BlueChannel) != 0)
1334         {
1335           if (pixel.blue < *minima)
1336             *minima=(double) pixel.blue;
1337           if (pixel.blue > *maxima)
1338             *maxima=(double) pixel.blue;
1339         }
1340       if ((channel & OpacityChannel) != 0)
1341         {
1342           if (pixel.opacity < *minima)
1343             *minima=(double) pixel.opacity;
1344           if (pixel.opacity > *maxima)
1345             *maxima=(double) pixel.opacity;
1346         }
1347       if (((channel & IndexChannel) != 0) &&
1348           (image->colorspace == CMYKColorspace))
1349         {
1350           if ((double) indexes[x] < *minima)
1351             *minima=(double) indexes[x];
1352           if ((double) indexes[x] > *maxima)
1353             *maxima=(double) indexes[x];
1354         }
1355       p++;
1356     }
1357   }
1358   return(y == (long) image->rows ? MagickTrue : MagickFalse);
1359 }
1360 \f
1361 /*
1362 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1363 %                                                                             %
1364 %                                                                             %
1365 %                                                                             %
1366 %   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                         %
1367 %                                                                             %
1368 %                                                                             %
1369 %                                                                             %
1370 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1371 %
1372 %  GetImageChannelStatistics() returns statistics for each channel in the
1373 %  image.  The statistics include the channel depth, its minima, maxima, mean,
1374 %  standard deviation, kurtosis and skewness.  You can access the red channel
1375 %  mean, for example, like this:
1376 %
1377 %      channel_statistics=GetImageChannelStatistics(image,excepton);
1378 %      red_mean=channel_statistics[RedChannel].mean;
1379 %
1380 %  Use MagickRelinquishMemory() to free the statistics buffer.
1381 %
1382 %  The format of the GetImageChannelStatistics method is:
1383 %
1384 %      ChannelStatistics *GetImageChannelStatistics(const Image *image,
1385 %        ExceptionInfo *exception)
1386 %
1387 %  A description of each parameter follows:
1388 %
1389 %    o image: the image.
1390 %
1391 %    o exception: return any errors or warnings in this structure.
1392 %
1393 */
1394 MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
1395   ExceptionInfo *exception)
1396 {
1397   ChannelStatistics
1398     *channel_statistics;
1399
1400   double
1401     area,
1402     sum_squares,
1403     sum_cubes;
1404
1405   long
1406     y;
1407
1408   MagickStatusType
1409     status;
1410
1411   QuantumAny
1412     range;
1413
1414   register long
1415     i;
1416
1417   size_t
1418     length;
1419
1420   unsigned long
1421     channels,
1422     depth;
1423
1424   assert(image != (Image *) NULL);
1425   assert(image->signature == MagickSignature);
1426   if (image->debug != MagickFalse)
1427     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1428   length=AllChannels+1UL;
1429   channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(length,
1430     sizeof(*channel_statistics));
1431   if (channel_statistics == (ChannelStatistics *) NULL)
1432     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1433   (void) ResetMagickMemory(channel_statistics,0,length*
1434     sizeof(*channel_statistics));
1435   for (i=0; i <= AllChannels; i++)
1436   {
1437     channel_statistics[i].depth=1;
1438     channel_statistics[i].maxima=(-1.0E-37);
1439     channel_statistics[i].minima=1.0E+37;
1440     channel_statistics[i].mean=0.0;
1441     channel_statistics[i].standard_deviation=0.0;
1442     channel_statistics[i].kurtosis=0.0;
1443     channel_statistics[i].skewness=0.0;
1444   }
1445   for (y=0; y < (long) image->rows; y++)
1446   {
1447     register const IndexPacket
1448       *restrict indexes;
1449
1450     register const PixelPacket
1451       *restrict p;
1452
1453     register long
1454       x;
1455
1456     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1457     if (p == (const PixelPacket *) NULL)
1458       break;
1459     indexes=GetVirtualIndexQueue(image);
1460     for (x=0; x < (long) image->columns; )
1461     {
1462       if (channel_statistics[RedChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
1463         {
1464           depth=channel_statistics[RedChannel].depth;
1465           range=GetQuantumRange(depth);
1466           status=p->red != ScaleAnyToQuantum(ScaleQuantumToAny(p->red,range),
1467             range) ? MagickTrue : MagickFalse;
1468           if (status != MagickFalse)
1469             {
1470               channel_statistics[RedChannel].depth++;
1471               continue;
1472             }
1473         }
1474       if (channel_statistics[GreenChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
1475         {
1476           depth=channel_statistics[GreenChannel].depth;
1477           range=GetQuantumRange(depth);
1478           status=p->green != ScaleAnyToQuantum(ScaleQuantumToAny(p->green,
1479             range),range) ? MagickTrue : MagickFalse;
1480           if (status != MagickFalse)
1481             {
1482               channel_statistics[GreenChannel].depth++;
1483               continue;
1484             }
1485         }
1486       if (channel_statistics[BlueChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
1487         {
1488           depth=channel_statistics[BlueChannel].depth;
1489           range=GetQuantumRange(depth);
1490           status=p->blue != ScaleAnyToQuantum(ScaleQuantumToAny(p->blue,
1491             range),range) ? MagickTrue : MagickFalse;
1492           if (status != MagickFalse)
1493             {
1494               channel_statistics[BlueChannel].depth++;
1495               continue;
1496             }
1497         }
1498       if (image->matte != MagickFalse)
1499         {
1500           if (channel_statistics[OpacityChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
1501             {
1502               depth=channel_statistics[OpacityChannel].depth;
1503               range=GetQuantumRange(depth);
1504               status=p->opacity != ScaleAnyToQuantum(ScaleQuantumToAny(
1505                 p->opacity,range),range) ? MagickTrue : MagickFalse;
1506               if (status != MagickFalse)
1507                 {
1508                   channel_statistics[OpacityChannel].depth++;
1509                   continue;
1510                 }
1511             }
1512           }
1513       if (image->colorspace == CMYKColorspace)
1514         {
1515           if (channel_statistics[BlackChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
1516             {
1517               depth=channel_statistics[BlackChannel].depth;
1518               range=GetQuantumRange(depth);
1519               status=indexes[x] != ScaleAnyToQuantum(ScaleQuantumToAny(
1520                 indexes[x],range),range) ? MagickTrue : MagickFalse;
1521               if (status != MagickFalse)
1522                 {
1523                   channel_statistics[BlackChannel].depth++;
1524                   continue;
1525                 }
1526             }
1527         }
1528       if ((double) p->red < channel_statistics[RedChannel].minima)
1529         channel_statistics[RedChannel].minima=(double) GetRedPixelComponent(p);
1530       if ((double) p->red > channel_statistics[RedChannel].maxima)
1531         channel_statistics[RedChannel].maxima=(double) GetRedPixelComponent(p);
1532       channel_statistics[RedChannel].mean+=GetRedPixelComponent(p);
1533       channel_statistics[RedChannel].standard_deviation+=(double) p->red*
1534         GetRedPixelComponent(p);
1535       channel_statistics[RedChannel].kurtosis+=(double) p->red*p->red*
1536         p->red*GetRedPixelComponent(p);
1537       channel_statistics[RedChannel].skewness+=(double) p->red*p->red*
1538         GetRedPixelComponent(p);
1539       if ((double) p->green < channel_statistics[GreenChannel].minima)
1540         channel_statistics[GreenChannel].minima=(double)
1541           GetGreenPixelComponent(p);
1542       if ((double) p->green > channel_statistics[GreenChannel].maxima)
1543         channel_statistics[GreenChannel].maxima=(double)
1544           GetGreenPixelComponent(p);
1545       channel_statistics[GreenChannel].mean+=GetGreenPixelComponent(p);
1546       channel_statistics[GreenChannel].standard_deviation+=(double) p->green*
1547         GetGreenPixelComponent(p);
1548       channel_statistics[GreenChannel].kurtosis+=(double) p->green*p->green*
1549         p->green*GetGreenPixelComponent(p);
1550       channel_statistics[GreenChannel].skewness+=(double) p->green*p->green*
1551         GetGreenPixelComponent(p);
1552       if ((double) p->blue < channel_statistics[BlueChannel].minima)
1553         channel_statistics[BlueChannel].minima=(double)
1554           GetBluePixelComponent(p);
1555       if ((double) p->blue > channel_statistics[BlueChannel].maxima)
1556         channel_statistics[BlueChannel].maxima=(double)
1557           GetBluePixelComponent(p);
1558       channel_statistics[BlueChannel].mean+=GetBluePixelComponent(p);
1559       channel_statistics[BlueChannel].standard_deviation+=(double) p->blue*
1560         GetBluePixelComponent(p);
1561       channel_statistics[BlueChannel].kurtosis+=(double) p->blue*p->blue*
1562         p->blue*GetBluePixelComponent(p);
1563       channel_statistics[BlueChannel].skewness+=(double) p->blue*p->blue*
1564         GetBluePixelComponent(p);
1565       if (image->matte != MagickFalse)
1566         {
1567           if ((double) p->opacity < channel_statistics[OpacityChannel].minima)
1568             channel_statistics[OpacityChannel].minima=(double)
1569               GetOpacityPixelComponent(p);
1570           if ((double) p->opacity > channel_statistics[OpacityChannel].maxima)
1571             channel_statistics[OpacityChannel].maxima=(double)
1572               GetOpacityPixelComponent(p);
1573           channel_statistics[OpacityChannel].mean+=GetOpacityPixelComponent(p);
1574           channel_statistics[OpacityChannel].standard_deviation+=(double)
1575             p->opacity*GetOpacityPixelComponent(p);
1576           channel_statistics[OpacityChannel].kurtosis+=(double) p->opacity*
1577             p->opacity*p->opacity*GetOpacityPixelComponent(p);
1578           channel_statistics[OpacityChannel].skewness+=(double) p->opacity*
1579             p->opacity*GetOpacityPixelComponent(p);
1580         }
1581       if (image->colorspace == CMYKColorspace)
1582         {
1583           if ((double) indexes[x] < channel_statistics[BlackChannel].minima)
1584             channel_statistics[BlackChannel].minima=(double) indexes[x];
1585           if ((double) indexes[x] > channel_statistics[BlackChannel].maxima)
1586             channel_statistics[BlackChannel].maxima=(double) indexes[x];
1587           channel_statistics[BlackChannel].mean+=indexes[x];
1588           channel_statistics[BlackChannel].standard_deviation+=(double)
1589             indexes[x]*indexes[x];
1590           channel_statistics[BlackChannel].kurtosis+=(double) indexes[x]*
1591             indexes[x]*indexes[x]*indexes[x];
1592           channel_statistics[BlackChannel].skewness+=(double) indexes[x]*
1593             indexes[x]*indexes[x];
1594         }
1595       x++;
1596       p++;
1597     }
1598   }
1599   area=(double) image->columns*image->rows;
1600   for (i=0; i < AllChannels; i++)
1601   {
1602     channel_statistics[i].mean/=area;
1603     channel_statistics[i].standard_deviation/=area;
1604     channel_statistics[i].kurtosis/=area;
1605     channel_statistics[i].skewness/=area;
1606   }
1607   for (i=0; i < AllChannels; i++)
1608   {
1609     channel_statistics[AllChannels].depth=(unsigned long) MagickMax((double)
1610       channel_statistics[AllChannels].depth,(double)
1611       channel_statistics[i].depth);
1612     channel_statistics[AllChannels].minima=MagickMin(
1613       channel_statistics[AllChannels].minima,channel_statistics[i].minima);
1614     channel_statistics[AllChannels].maxima=MagickMax(
1615       channel_statistics[AllChannels].maxima,channel_statistics[i].maxima);
1616     channel_statistics[AllChannels].mean+=channel_statistics[i].mean;
1617     channel_statistics[AllChannels].standard_deviation+=
1618       channel_statistics[i].standard_deviation;
1619     channel_statistics[AllChannels].kurtosis+=channel_statistics[i].kurtosis;
1620     channel_statistics[AllChannels].skewness+=channel_statistics[i].skewness;
1621   }
1622   channels=4;
1623   if (image->colorspace == CMYKColorspace)
1624     channels++;
1625   channel_statistics[AllChannels].mean/=channels;
1626   channel_statistics[AllChannels].standard_deviation/=channels;
1627   channel_statistics[AllChannels].kurtosis/=channels;
1628   channel_statistics[AllChannels].skewness/=channels;
1629   for (i=0; i <= AllChannels; i++)
1630   {
1631     sum_squares=0.0;
1632     sum_squares=channel_statistics[i].standard_deviation;
1633     sum_cubes=0.0;
1634     sum_cubes=channel_statistics[i].skewness;
1635     channel_statistics[i].standard_deviation=sqrt(
1636       channel_statistics[i].standard_deviation-
1637        (channel_statistics[i].mean*channel_statistics[i].mean));
1638     if (channel_statistics[i].standard_deviation == 0.0)
1639       {
1640         channel_statistics[i].kurtosis=0.0;
1641         channel_statistics[i].skewness=0.0;
1642       }
1643     else
1644       {
1645         channel_statistics[i].skewness=(channel_statistics[i].skewness-
1646           3.0*channel_statistics[i].mean*sum_squares+
1647           2.0*channel_statistics[i].mean*channel_statistics[i].mean*
1648           channel_statistics[i].mean)/
1649           (channel_statistics[i].standard_deviation*
1650            channel_statistics[i].standard_deviation*
1651            channel_statistics[i].standard_deviation);
1652         channel_statistics[i].kurtosis=(channel_statistics[i].kurtosis-
1653           4.0*channel_statistics[i].mean*sum_cubes+
1654           6.0*channel_statistics[i].mean*channel_statistics[i].mean*sum_squares-
1655           3.0*channel_statistics[i].mean*channel_statistics[i].mean*
1656           1.0*channel_statistics[i].mean*channel_statistics[i].mean)/
1657           (channel_statistics[i].standard_deviation*
1658            channel_statistics[i].standard_deviation*
1659            channel_statistics[i].standard_deviation*
1660            channel_statistics[i].standard_deviation)-3.0;
1661       }
1662   }
1663   return(channel_statistics);
1664 }