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

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