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

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