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