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