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