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