granicus.if.org Git - imagemagick/blob - magick/statistic.c

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