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 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/segment.h"
  83 #include "magick/semaphore.h"
  84 #include "magick/signature-private.h"
  85 #include "magick/statistic.h"
  86 #include "magick/string_.h"
  87 #include "magick/thread-private.h"
  88 #include "magick/timer.h"
  89 #include "magick/utility.h"
  90 #include "magick/version.h"
  91 \f
  92 /*
  93 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  94 %                                                                             %
  95 %                                                                             %
  96 %                                                                             %
  97 %     A v e r a g e I m a g e s                                               %
  98 %                                                                             %
  99 %                                                                             %
 100 %                                                                             %
 101 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 102 %
 103 %  AverageImages() takes a set of images and averages them together.  Each
 104 %  image in the set must have the same width and height.  AverageImages()
 105 %  returns a single image with each corresponding pixel component of each
 106 %  image averaged.   On failure, a NULL image is returned and exception
 107 %  describes the reason for the failure.
 108 %
 109 %  The format of the AverageImages method is:
 110 %
 111 %      Image *AverageImages(Image *image,ExceptionInfo *exception)
 112 %
 113 %  A description of each parameter follows:
 114 %
 115 %    o image: the image sequence.
 116 %
 117 %    o exception: return any errors or warnings in this structure.
 118 %
 119 */
 120
 121 static MagickPixelPacket **DestroyPixelThreadSet(MagickPixelPacket **pixels)
 122 {
 123   register long
 124     i;
 125
 126   assert(pixels != (MagickPixelPacket **) NULL);
 127   for (i=0; i < (long) GetOpenMPMaximumThreads(); i++)
 128     if (pixels[i] != (MagickPixelPacket *) NULL)
 129       pixels[i]=(MagickPixelPacket *) RelinquishMagickMemory(pixels[i]);
 130   pixels=(MagickPixelPacket **) RelinquishAlignedMemory(pixels);
 131   return(pixels);
 132 }
 133
 134 static MagickPixelPacket **AcquirePixelThreadSet(const Image *image)
 135 {
 136   register long
 137     i,
 138     j;
 139
 140   MagickPixelPacket
 141     **pixels;
 142
 143   unsigned long
 144     number_threads;
 145
 146   number_threads=GetOpenMPMaximumThreads();
 147   pixels=(MagickPixelPacket **) AcquireAlignedMemory(number_threads,
 148     sizeof(*pixels));
 149   if (pixels == (MagickPixelPacket **) NULL)
 150     return((MagickPixelPacket **) NULL);
 151   (void) ResetMagickMemory(pixels,0,number_threads*sizeof(*pixels));
 152   for (i=0; i < (long) number_threads; i++)
 153   {
 154     pixels[i]=(MagickPixelPacket *) AcquireQuantumMemory(image->columns,
 155       sizeof(**pixels));
 156     if (pixels[i] == (MagickPixelPacket *) NULL)
 157       return(DestroyPixelThreadSet(pixels));
 158     for (j=0; j < (long) image->columns; j++)
 159       GetMagickPixelPacket(image,&pixels[i][j]);
 160   }
 161   return(pixels);
 162 }
 163
 164 MagickExport Image *AverageImages(const Image *image,ExceptionInfo *exception)
 165 {
 166 #define AverageImageTag  "Average/Image"
 167
 168   CacheView
 169     *average_view;
 170
 171   const Image
 172     *next;
 173
 174   Image
 175     *average_image;
 176
 177   long
 178     progress,
 179     y;
 180
 181   MagickBooleanType
 182     status;
 183
 184   MagickPixelPacket
 185     **restrict average_pixels,
 186     zero;
 187
 188   unsigned long
 189     number_images;
 190
 191   /*
 192     Ensure the image are the same size.
 193   */
 194   assert(image != (Image *) NULL);
 195   assert(image->signature == MagickSignature);
 196   if (image->debug != MagickFalse)
 197     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 198   assert(exception != (ExceptionInfo *) NULL);
 199   assert(exception->signature == MagickSignature);
 200   for (next=image; next != (Image *) NULL; next=GetNextImageInList(next))
 201     if ((next->columns != image->columns) || (next->rows != image->rows))
 202       ThrowImageException(OptionError,"ImageWidthsOrHeightsDiffer");
 203   /*
 204     Initialize average next attributes.
 205   */
 206   average_image=CloneImage(image,image->columns,image->rows,MagickTrue,
 207     exception);
 208   if (average_image == (Image *) NULL)
 209     return((Image *) NULL);
 210   if (SetImageStorageClass(average_image,DirectClass) == MagickFalse)
 211     {
 212       InheritException(exception,&average_image->exception);
 213       average_image=DestroyImage(average_image);
 214       return((Image *) NULL);
 215     }
 216   average_pixels=AcquirePixelThreadSet(image);
 217   if (average_pixels == (MagickPixelPacket **) NULL)
 218     {
 219       average_image=DestroyImage(average_image);
 220       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 221     }
 222   /*
 223     Average image pixels.
 224   */
 225   status=MagickTrue;
 226   progress=0;
 227   GetMagickPixelPacket(image,&zero);
 228   number_images=GetImageListLength(image);
 229   average_view=AcquireCacheView(average_image);
 230 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 231   #pragma omp parallel for schedule(dynamic) shared(progress,status)
 232 #endif
 233   for (y=0; y < (long) average_image->rows; y++)
 234   {
 235     CacheView
 236       *image_view;
 237
 238     const Image
 239       *next;
 240
 241     MagickPixelPacket
 242       pixel;
 243
 244     register IndexPacket
 245       *restrict average_indexes;
 246
 247     register long
 248       i,
 249       id,
 250       x;
 251
 252     register MagickPixelPacket
 253       *average_pixel;
 254
 255     register PixelPacket
 256       *restrict q;
 257
 258     if (status == MagickFalse)
 259       continue;
 260     q=QueueCacheViewAuthenticPixels(average_view,0,y,average_image->columns,1,
 261       exception);
 262     if (q == (PixelPacket *) NULL)
 263       {
 264         status=MagickFalse;
 265         continue;
 266       }
 267     average_indexes=GetCacheViewAuthenticIndexQueue(average_view);
 268     pixel=zero;
 269     id=GetOpenMPThreadId();
 270     average_pixel=average_pixels[id];
 271     for (x=0; x < (long) average_image->columns; x++)
 272       average_pixel[x]=zero;
 273     next=image;
 274     for (i=0; i < (long) number_images; i++)
 275     {
 276       register const IndexPacket
 277         *indexes;
 278
 279       register const PixelPacket
 280         *p;
 281
 282       image_view=AcquireCacheView(next);
 283       p=GetCacheViewVirtualPixels(image_view,0,y,next->columns,1,exception);
 284       if (p == (const PixelPacket *) NULL)
 285         {
 286           image_view=DestroyCacheView(image_view);
 287           break;
 288         }
 289       indexes=GetCacheViewVirtualIndexQueue(image_view);
 290       for (x=0; x < (long) next->columns; x++)
 291       {
 292         SetMagickPixelPacket(next,p,indexes+x,&pixel);
 293         average_pixel[x].red+=QuantumScale*pixel.red;
 294         average_pixel[x].green+=QuantumScale*pixel.green;
 295         average_pixel[x].blue+=QuantumScale*pixel.blue;
 296         average_pixel[x].opacity+=QuantumScale*pixel.opacity;
 297         if (average_image->colorspace == CMYKColorspace)
 298           average_pixel[x].index+=QuantumScale*pixel.index;
 299         p++;
 300       }
 301       image_view=DestroyCacheView(image_view);
 302       next=GetNextImageInList(next);
 303     }
 304     for (x=0; x < (long) average_image->columns; x++)
 305     {
 306       average_pixel[x].red=(MagickRealType) (QuantumRange*
 307         average_pixel[x].red/number_images);
 308       average_pixel[x].green=(MagickRealType) (QuantumRange*
 309         average_pixel[x].green/number_images);
 310       average_pixel[x].blue=(MagickRealType) (QuantumRange*
 311         average_pixel[x].blue/number_images);
 312       average_pixel[x].opacity=(MagickRealType) (QuantumRange*
 313         average_pixel[x].opacity/number_images);
 314       if (average_image->colorspace == CMYKColorspace)
 315         average_pixel[x].index=(MagickRealType) (QuantumRange*
 316           average_pixel[x].index/number_images);
 317       SetPixelPacket(average_image,&average_pixel[x],q,average_indexes+x);
 318       q++;
 319     }
 320     if (SyncCacheViewAuthenticPixels(average_view,exception) == MagickFalse)
 321       status=MagickFalse;
 322     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 323       {
 324         MagickBooleanType
 325           proceed;
 326
 327 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 328         #pragma omp critical (MagickCore_AverageImages)
 329 #endif
 330         proceed=SetImageProgress(image,AverageImageTag,progress++,
 331           average_image->rows);
 332         if (proceed == MagickFalse)
 333           status=MagickFalse;
 334       }
 335   }
 336   average_view=DestroyCacheView(average_view);
 337   average_pixels=DestroyPixelThreadSet(average_pixels);
 338   if (status == MagickFalse)
 339     average_image=DestroyImage(average_image);
 340   return(average_image);
 341 }
 342 \f
 343 /*
 344 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 345 %                                                                             %
 346 %                                                                             %
 347 %                                                                             %
 348 +   G e t I m a g e C h a n n e l E x t r e m a                               %
 349 %                                                                             %
 350 %                                                                             %
 351 %                                                                             %
 352 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 353 %
 354 %  GetImageChannelExtrema() returns the extrema of one or more image channels.
 355 %
 356 %  The format of the GetImageChannelExtrema method is:
 357 %
 358 %      MagickBooleanType GetImageChannelExtrema(const Image *image,
 359 %        const ChannelType channel,unsigned long *minima,unsigned long *maxima,
 360 %        ExceptionInfo *exception)
 361 %
 362 %  A description of each parameter follows:
 363 %
 364 %    o image: the image.
 365 %
 366 %    o channel: the channel.
 367 %
 368 %    o minima: the minimum value in the channel.
 369 %
 370 %    o maxima: the maximum value in the channel.
 371 %
 372 %    o exception: return any errors or warnings in this structure.
 373 %
 374 */
 375
 376 MagickExport MagickBooleanType GetImageExtrema(const Image *image,
 377   unsigned long *minima,unsigned long *maxima,ExceptionInfo *exception)
 378 {
 379   return(GetImageChannelExtrema(image,AllChannels,minima,maxima,exception));
 380 }
 381
 382 MagickExport MagickBooleanType GetImageChannelExtrema(const Image *image,
 383   const ChannelType channel,unsigned long *minima,unsigned long *maxima,
 384   ExceptionInfo *exception)
 385 {
 386   double
 387     max,
 388     min;
 389
 390   MagickBooleanType
 391     status;
 392
 393   assert(image != (Image *) NULL);
 394   assert(image->signature == MagickSignature);
 395   if (image->debug != MagickFalse)
 396     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 397   status=GetImageChannelRange(image,channel,&min,&max,exception);
 398   *minima=(unsigned long) (min+0.5);
 399   *maxima=(unsigned long) (max+0.5);
 400   return(status);
 401 }
 402 \f
 403 /*
 404 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 405 %                                                                             %
 406 %                                                                             %
 407 %                                                                             %
 408 %   G e t I m a g e C h a n n e l M e a n                                     %
 409 %                                                                             %
 410 %                                                                             %
 411 %                                                                             %
 412 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 413 %
 414 %  GetImageChannelMean() returns the mean and standard deviation of one or more
 415 %  image channels.
 416 %
 417 %  The format of the GetImageChannelMean method is:
 418 %
 419 %      MagickBooleanType GetImageChannelMean(const Image *image,
 420 %        const ChannelType channel,double *mean,double *standard_deviation,
 421 %        ExceptionInfo *exception)
 422 %
 423 %  A description of each parameter follows:
 424 %
 425 %    o image: the image.
 426 %
 427 %    o channel: the channel.
 428 %
 429 %    o mean: the average value in the channel.
 430 %
 431 %    o standard_deviation: the standard deviation of the channel.
 432 %
 433 %    o exception: return any errors or warnings in this structure.
 434 %
 435 */
 436
 437 MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean,
 438   double *standard_deviation,ExceptionInfo *exception)
 439 {
 440   MagickBooleanType
 441     status;
 442
 443   status=GetImageChannelMean(image,AllChannels,mean,standard_deviation,
 444     exception);
 445   return(status);
 446 }
 447
 448 MagickExport MagickBooleanType GetImageChannelMean(const Image *image,
 449   const ChannelType channel,double *mean,double *standard_deviation,
 450   ExceptionInfo *exception)
 451 {
 452   double
 453     area;
 454
 455   long
 456     y;
 457
 458   assert(image != (Image *) NULL);
 459   assert(image->signature == MagickSignature);
 460   if (image->debug != MagickFalse)
 461     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 462   *mean=0.0;
 463   *standard_deviation=0.0;
 464   area=0.0;
 465   for (y=0; y < (long) image->rows; y++)
 466   {
 467     register const IndexPacket
 468       *restrict indexes;
 469
 470     register const PixelPacket
 471       *restrict p;
 472
 473     register long
 474       x;
 475
 476     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
 477     if (p == (const PixelPacket *) NULL)
 478       break;
 479     indexes=GetVirtualIndexQueue(image);
 480     for (x=0; x < (long) image->columns; x++)
 481     {
 482       if ((channel & RedChannel) != 0)
 483         {
 484           *mean+=GetRedPixelComponent(p);
 485           *standard_deviation+=(double) p->red*GetRedPixelComponent(p);
 486           area++;
 487         }
 488       if ((channel & GreenChannel) != 0)
 489         {
 490           *mean+=GetGreenPixelComponent(p);
 491           *standard_deviation+=(double) p->green*GetGreenPixelComponent(p);
 492           area++;
 493         }
 494       if ((channel & BlueChannel) != 0)
 495         {
 496           *mean+=GetBluePixelComponent(p);
 497           *standard_deviation+=(double) p->blue*GetBluePixelComponent(p);
 498           area++;
 499         }
 500       if ((channel & OpacityChannel) != 0)
 501         {
 502           *mean+=GetOpacityPixelComponent(p);
 503           *standard_deviation+=(double) p->opacity*GetOpacityPixelComponent(p);
 504           area++;
 505         }
 506       if (((channel & IndexChannel) != 0) &&
 507           (image->colorspace == CMYKColorspace))
 508         {
 509           *mean+=indexes[x];
 510           *standard_deviation+=(double) indexes[x]*indexes[x];
 511           area++;
 512         }
 513       p++;
 514     }
 515   }
 516   if (y < (long) image->rows)
 517     return(MagickFalse);
 518   if (area != 0)
 519     {
 520       *mean/=area;
 521       *standard_deviation/=area;
 522     }
 523   *standard_deviation=sqrt(*standard_deviation-(*mean*(*mean)));
 524   return(y == (long) image->rows ? MagickTrue : MagickFalse);
 525 }
 526 \f
 527 /*
 528 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 529 %                                                                             %
 530 %                                                                             %
 531 %                                                                             %
 532 %   G e t I m a g e C h a n n e l K u r t o s i s                             %
 533 %                                                                             %
 534 %                                                                             %
 535 %                                                                             %
 536 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 537 %
 538 %  GetImageChannelKurtosis() returns the kurtosis and skewness of one or more
 539 %  image channels.
 540 %
 541 %  The format of the GetImageChannelKurtosis method is:
 542 %
 543 %      MagickBooleanType GetImageChannelKurtosis(const Image *image,
 544 %        const ChannelType channel,double *kurtosis,double *skewness,
 545 %        ExceptionInfo *exception)
 546 %
 547 %  A description of each parameter follows:
 548 %
 549 %    o image: the image.
 550 %
 551 %    o channel: the channel.
 552 %
 553 %    o kurtosis: the kurtosis of the channel.
 554 %
 555 %    o skewness: the skewness of the channel.
 556 %
 557 %    o exception: return any errors or warnings in this structure.
 558 %
 559 */
 560
 561 MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
 562   double *kurtosis,double *skewness,ExceptionInfo *exception)
 563 {
 564   MagickBooleanType
 565     status;
 566
 567   status=GetImageChannelKurtosis(image,AllChannels,kurtosis,skewness,
 568     exception);
 569   return(status);
 570 }
 571
 572 MagickExport MagickBooleanType GetImageChannelKurtosis(const Image *image,
 573   const ChannelType channel,double *kurtosis,double *skewness,
 574   ExceptionInfo *exception)
 575 {
 576   double
 577     area,
 578     mean,
 579     standard_deviation,
 580     sum_squares,
 581     sum_cubes,
 582     sum_fourth_power;
 583
 584   long
 585     y;
 586
 587   assert(image != (Image *) NULL);
 588   assert(image->signature == MagickSignature);
 589   if (image->debug != MagickFalse)
 590     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 591   *kurtosis=0.0;
 592   *skewness=0.0;
 593   area=0.0;
 594   mean=0.0;
 595   standard_deviation=0.0;
 596   sum_squares=0.0;
 597   sum_cubes=0.0;
 598   sum_fourth_power=0.0;
 599   for (y=0; y < (long) image->rows; y++)
 600   {
 601     register const IndexPacket
 602       *restrict indexes;
 603
 604     register const PixelPacket
 605       *restrict p;
 606
 607     register long
 608       x;
 609
 610     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
 611     if (p == (const PixelPacket *) NULL)
 612       break;
 613     indexes=GetVirtualIndexQueue(image);
 614     for (x=0; x < (long) image->columns; x++)
 615     {
 616       if ((channel & RedChannel) != 0)
 617         {
 618           mean+=GetRedPixelComponent(p);
 619           sum_squares+=(double) p->red*GetRedPixelComponent(p);
 620           sum_cubes+=(double) p->red*p->red*GetRedPixelComponent(p);
 621           sum_fourth_power+=(double) p->red*p->red*p->red*
 622             GetRedPixelComponent(p);
 623           area++;
 624         }
 625       if ((channel & GreenChannel) != 0)
 626         {
 627           mean+=GetGreenPixelComponent(p);
 628           sum_squares+=(double) p->green*GetGreenPixelComponent(p);
 629           sum_cubes+=(double) p->green*p->green*GetGreenPixelComponent(p);
 630           sum_fourth_power+=(double) p->green*p->green*p->green*
 631             GetGreenPixelComponent(p);
 632           area++;
 633         }
 634       if ((channel & BlueChannel) != 0)
 635         {
 636           mean+=GetBluePixelComponent(p);
 637           sum_squares+=(double) p->blue*GetBluePixelComponent(p);
 638           sum_cubes+=(double) p->blue*p->blue*GetBluePixelComponent(p);
 639           sum_fourth_power+=(double) p->blue*p->blue*p->blue*
 640             GetBluePixelComponent(p);
 641           area++;
 642         }
 643       if ((channel & OpacityChannel) != 0)
 644         {
 645           mean+=GetOpacityPixelComponent(p);
 646           sum_squares+=(double) p->opacity*GetOpacityPixelComponent(p);
 647           sum_cubes+=(double) p->opacity*p->opacity*GetOpacityPixelComponent(p);
 648           sum_fourth_power+=(double) p->opacity*p->opacity*p->opacity*
 649             GetOpacityPixelComponent(p);
 650           area++;
 651         }
 652       if (((channel & IndexChannel) != 0) &&
 653           (image->colorspace == CMYKColorspace))
 654         {
 655           mean+=indexes[x];
 656           sum_squares+=(double) indexes[x]*indexes[x];
 657           sum_cubes+=(double) indexes[x]*indexes[x]*indexes[x];
 658           sum_fourth_power+=(double) indexes[x]*indexes[x]*indexes[x]*
 659             indexes[x];
 660           area++;
 661         }
 662       p++;
 663     }
 664   }
 665   if (y < (long) image->rows)
 666     return(MagickFalse);
 667   if (area != 0.0)
 668     {
 669       mean/=area;
 670       sum_squares/=area;
 671       sum_cubes/=area;
 672       sum_fourth_power/=area;
 673     }
 674   standard_deviation=sqrt(sum_squares-(mean*mean));
 675   if (standard_deviation != 0.0)
 676     {
 677       *kurtosis=sum_fourth_power-4.0*mean*sum_cubes+6.0*mean*mean*sum_squares-
 678         3.0*mean*mean*mean*mean;
 679       *kurtosis/=standard_deviation*standard_deviation*standard_deviation*
 680         standard_deviation;
 681       *kurtosis-=3.0;
 682       *skewness=sum_cubes-3.0*mean*sum_squares+2.0*mean*mean*mean;
 683       *skewness/=standard_deviation*standard_deviation*standard_deviation;
 684     }
 685   return(y == (long) image->rows ? MagickTrue : MagickFalse);
 686 }
 687 \f
 688 /*
 689 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 690 %                                                                             %
 691 %                                                                             %
 692 %                                                                             %
 693 %   G e t I m a g e C h a n n e l R a n g e                                   %
 694 %                                                                             %
 695 %                                                                             %
 696 %                                                                             %
 697 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 698 %
 699 %  GetImageChannelRange() returns the range of one or more image channels.
 700 %
 701 %  The format of the GetImageChannelRange method is:
 702 %
 703 %      MagickBooleanType GetImageChannelRange(const Image *image,
 704 %        const ChannelType channel,double *minima,double *maxima,
 705 %        ExceptionInfo *exception)
 706 %
 707 %  A description of each parameter follows:
 708 %
 709 %    o image: the image.
 710 %
 711 %    o channel: the channel.
 712 %
 713 %    o minima: the minimum value in the channel.
 714 %
 715 %    o maxima: the maximum value in the channel.
 716 %
 717 %    o exception: return any errors or warnings in this structure.
 718 %
 719 */
 720
 721 MagickExport MagickBooleanType GetImageRange(const Image *image,
 722   double *minima,double *maxima,ExceptionInfo *exception)
 723 {
 724   return(GetImageChannelRange(image,AllChannels,minima,maxima,exception));
 725 }
 726
 727 MagickExport MagickBooleanType GetImageChannelRange(const Image *image,
 728   const ChannelType channel,double *minima,double *maxima,
 729   ExceptionInfo *exception)
 730 {
 731   long
 732     y;
 733
 734   MagickPixelPacket
 735     pixel;
 736
 737   assert(image != (Image *) NULL);
 738   assert(image->signature == MagickSignature);
 739   if (image->debug != MagickFalse)
 740     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 741   *maxima=(-1.0E-37);
 742   *minima=1.0E+37;
 743   GetMagickPixelPacket(image,&pixel);
 744   for (y=0; y < (long) image->rows; y++)
 745   {
 746     register const IndexPacket
 747       *restrict indexes;
 748
 749     register const PixelPacket
 750       *restrict p;
 751
 752     register long
 753       x;
 754
 755     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
 756     if (p == (const PixelPacket *) NULL)
 757       break;
 758     indexes=GetVirtualIndexQueue(image);
 759     for (x=0; x < (long) image->columns; x++)
 760     {
 761       SetMagickPixelPacket(image,p,indexes+x,&pixel);
 762       if ((channel & RedChannel) != 0)
 763         {
 764           if (pixel.red < *minima)
 765             *minima=(double) pixel.red;
 766           if (pixel.red > *maxima)
 767             *maxima=(double) pixel.red;
 768         }
 769       if ((channel & GreenChannel) != 0)
 770         {
 771           if (pixel.green < *minima)
 772             *minima=(double) pixel.green;
 773           if (pixel.green > *maxima)
 774             *maxima=(double) pixel.green;
 775         }
 776       if ((channel & BlueChannel) != 0)
 777         {
 778           if (pixel.blue < *minima)
 779             *minima=(double) pixel.blue;
 780           if (pixel.blue > *maxima)
 781             *maxima=(double) pixel.blue;
 782         }
 783       if ((channel & OpacityChannel) != 0)
 784         {
 785           if (pixel.opacity < *minima)
 786             *minima=(double) pixel.opacity;
 787           if (pixel.opacity > *maxima)
 788             *maxima=(double) pixel.opacity;
 789         }
 790       if (((channel & IndexChannel) != 0) &&
 791           (image->colorspace == CMYKColorspace))
 792         {
 793           if ((double) indexes[x] < *minima)
 794             *minima=(double) indexes[x];
 795           if ((double) indexes[x] > *maxima)
 796             *maxima=(double) indexes[x];
 797         }
 798       p++;
 799     }
 800   }
 801   return(y == (long) image->rows ? MagickTrue : MagickFalse);
 802 }
 803 \f
 804 /*
 805 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 806 %                                                                             %
 807 %                                                                             %
 808 %                                                                             %
 809 %   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                         %
 810 %                                                                             %
 811 %                                                                             %
 812 %                                                                             %
 813 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 814 %
 815 %  GetImageChannelStatistics() returns statistics for each channel in the
 816 %  image.  The statistics include the channel depth, its minima, maxima, mean,
 817 %  standard deviation, kurtosis and skewness.  You can access the red channel
 818 %  mean, for example, like this:
 819 %
 820 %      channel_statistics=GetImageChannelStatistics(image,excepton);
 821 %      red_mean=channel_statistics[RedChannel].mean;
 822 %
 823 %  Use MagickRelinquishMemory() to free the statistics buffer.
 824 %
 825 %  The format of the GetImageChannelStatistics method is:
 826 %
 827 %      ChannelStatistics *GetImageChannelStatistics(const Image *image,
 828 %        ExceptionInfo *exception)
 829 %
 830 %  A description of each parameter follows:
 831 %
 832 %    o image: the image.
 833 %
 834 %    o exception: return any errors or warnings in this structure.
 835 %
 836 */
 837
 838 static inline double MagickMax(const double x,const double y)
 839 {
 840   if (x > y)
 841     return(x);
 842   return(y);
 843 }
 844
 845 static inline double MagickMin(const double x,const double y)
 846 {
 847   if (x < y)
 848     return(x);
 849   return(y);
 850 }
 851
 852 MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
 853   ExceptionInfo *exception)
 854 {
 855   ChannelStatistics
 856     *channel_statistics;
 857
 858   double
 859     area,
 860     sum_squares,
 861     sum_cubes;
 862
 863   long
 864     y;
 865
 866   MagickStatusType
 867     status;
 868
 869   QuantumAny
 870     range;
 871
 872   register long
 873     i;
 874
 875   size_t
 876     length;
 877
 878   unsigned long
 879     channels,
 880     depth;
 881
 882   assert(image != (Image *) NULL);
 883   assert(image->signature == MagickSignature);
 884   if (image->debug != MagickFalse)
 885     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 886   length=AllChannels+1UL;
 887   channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(length,
 888     sizeof(*channel_statistics));
 889   if (channel_statistics == (ChannelStatistics *) NULL)
 890     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
 891   (void) ResetMagickMemory(channel_statistics,0,length*
 892     sizeof(*channel_statistics));
 893   for (i=0; i <= AllChannels; i++)
 894   {
 895     channel_statistics[i].depth=1;
 896     channel_statistics[i].maxima=(-1.0E-37);
 897     channel_statistics[i].minima=1.0E+37;
 898     channel_statistics[i].mean=0.0;
 899     channel_statistics[i].standard_deviation=0.0;
 900     channel_statistics[i].kurtosis=0.0;
 901     channel_statistics[i].skewness=0.0;
 902   }
 903   for (y=0; y < (long) image->rows; y++)
 904   {
 905     register const IndexPacket
 906       *restrict indexes;
 907
 908     register const PixelPacket
 909       *restrict p;
 910
 911     register long
 912       x;
 913
 914     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
 915     if (p == (const PixelPacket *) NULL)
 916       break;
 917     indexes=GetVirtualIndexQueue(image);
 918     for (x=0; x < (long) image->columns; )
 919     {
 920       if (channel_statistics[RedChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
 921         {
 922           depth=channel_statistics[RedChannel].depth;
 923           range=GetQuantumRange(depth);
 924           status=p->red != ScaleAnyToQuantum(ScaleQuantumToAny(p->red,range),
 925             range) ? MagickTrue : MagickFalse;
 926           if (status != MagickFalse)
 927             {
 928               channel_statistics[RedChannel].depth++;
 929               continue;
 930             }
 931         }
 932       if (channel_statistics[GreenChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
 933         {
 934           depth=channel_statistics[GreenChannel].depth;
 935           range=GetQuantumRange(depth);
 936           status=p->green != ScaleAnyToQuantum(ScaleQuantumToAny(p->green,
 937             range),range) ? MagickTrue : MagickFalse;
 938           if (status != MagickFalse)
 939             {
 940               channel_statistics[GreenChannel].depth++;
 941               continue;
 942             }
 943         }
 944       if (channel_statistics[BlueChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
 945         {
 946           depth=channel_statistics[BlueChannel].depth;
 947           range=GetQuantumRange(depth);
 948           status=p->blue != ScaleAnyToQuantum(ScaleQuantumToAny(p->blue,
 949             range),range) ? MagickTrue : MagickFalse;
 950           if (status != MagickFalse)
 951             {
 952               channel_statistics[BlueChannel].depth++;
 953               continue;
 954             }
 955         }
 956       if (image->matte != MagickFalse)
 957         {
 958           if (channel_statistics[OpacityChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
 959             {
 960               depth=channel_statistics[OpacityChannel].depth;
 961               range=GetQuantumRange(depth);
 962               status=p->opacity != ScaleAnyToQuantum(ScaleQuantumToAny(
 963                 p->opacity,range),range) ? MagickTrue : MagickFalse;
 964               if (status != MagickFalse)
 965                 {
 966                   channel_statistics[OpacityChannel].depth++;
 967                   continue;
 968                 }
 969             }
 970           }
 971       if (image->colorspace == CMYKColorspace)
 972         {
 973           if (channel_statistics[BlackChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
 974             {
 975               depth=channel_statistics[BlackChannel].depth;
 976               range=GetQuantumRange(depth);
 977               status=indexes[x] != ScaleAnyToQuantum(ScaleQuantumToAny(
 978                 indexes[x],range),range) ? MagickTrue : MagickFalse;
 979               if (status != MagickFalse)
 980                 {
 981                   channel_statistics[BlackChannel].depth++;
 982                   continue;
 983                 }
 984             }
 985         }
 986       if ((double) p->red < channel_statistics[RedChannel].minima)
 987         channel_statistics[RedChannel].minima=(double) GetRedPixelComponent(p);
 988       if ((double) p->red > channel_statistics[RedChannel].maxima)
 989         channel_statistics[RedChannel].maxima=(double) GetRedPixelComponent(p);
 990       channel_statistics[RedChannel].mean+=GetRedPixelComponent(p);
 991       channel_statistics[RedChannel].standard_deviation+=(double) p->red*
 992         GetRedPixelComponent(p);
 993       channel_statistics[RedChannel].kurtosis+=(double) p->red*p->red*
 994         p->red*GetRedPixelComponent(p);
 995       channel_statistics[RedChannel].skewness+=(double) p->red*p->red*
 996         GetRedPixelComponent(p);
 997       if ((double) p->green < channel_statistics[GreenChannel].minima)
 998         channel_statistics[GreenChannel].minima=(double)
 999           GetGreenPixelComponent(p);
1000       if ((double) p->green > channel_statistics[GreenChannel].maxima)
1001         channel_statistics[GreenChannel].maxima=(double)
1002           GetGreenPixelComponent(p);
1003       channel_statistics[GreenChannel].mean+=GetGreenPixelComponent(p);
1004       channel_statistics[GreenChannel].standard_deviation+=(double) p->green*
1005         GetGreenPixelComponent(p);
1006       channel_statistics[GreenChannel].kurtosis+=(double) p->green*p->green*
1007         p->green*GetGreenPixelComponent(p);
1008       channel_statistics[GreenChannel].skewness+=(double) p->green*p->green*
1009         GetGreenPixelComponent(p);
1010       if ((double) p->blue < channel_statistics[BlueChannel].minima)
1011         channel_statistics[BlueChannel].minima=(double)
1012           GetBluePixelComponent(p);
1013       if ((double) p->blue > channel_statistics[BlueChannel].maxima)
1014         channel_statistics[BlueChannel].maxima=(double)
1015           GetBluePixelComponent(p);
1016       channel_statistics[BlueChannel].mean+=GetBluePixelComponent(p);
1017       channel_statistics[BlueChannel].standard_deviation+=(double) p->blue*
1018         GetBluePixelComponent(p);
1019       channel_statistics[BlueChannel].kurtosis+=(double) p->blue*p->blue*
1020         p->blue*GetBluePixelComponent(p);
1021       channel_statistics[BlueChannel].skewness+=(double) p->blue*p->blue*
1022         GetBluePixelComponent(p);
1023       if (image->matte != MagickFalse)
1024         {
1025           if ((double) p->opacity < channel_statistics[OpacityChannel].minima)
1026             channel_statistics[OpacityChannel].minima=(double)
1027               GetOpacityPixelComponent(p);
1028           if ((double) p->opacity > channel_statistics[OpacityChannel].maxima)
1029             channel_statistics[OpacityChannel].maxima=(double)
1030               GetOpacityPixelComponent(p);
1031           channel_statistics[OpacityChannel].mean+=GetOpacityPixelComponent(p);
1032           channel_statistics[OpacityChannel].standard_deviation+=(double)
1033             p->opacity*GetOpacityPixelComponent(p);
1034           channel_statistics[OpacityChannel].kurtosis+=(double) p->opacity*
1035             p->opacity*p->opacity*GetOpacityPixelComponent(p);
1036           channel_statistics[OpacityChannel].skewness+=(double) p->opacity*
1037             p->opacity*GetOpacityPixelComponent(p);
1038         }
1039       if (image->colorspace == CMYKColorspace)
1040         {
1041           if ((double) indexes[x] < channel_statistics[BlackChannel].minima)
1042             channel_statistics[BlackChannel].minima=(double) indexes[x];
1043           if ((double) indexes[x] > channel_statistics[BlackChannel].maxima)
1044             channel_statistics[BlackChannel].maxima=(double) indexes[x];
1045           channel_statistics[BlackChannel].mean+=indexes[x];
1046           channel_statistics[BlackChannel].standard_deviation+=(double)
1047             indexes[x]*indexes[x];
1048           channel_statistics[BlackChannel].kurtosis+=(double) indexes[x]*
1049             indexes[x]*indexes[x]*indexes[x];
1050           channel_statistics[BlackChannel].skewness+=(double) indexes[x]*
1051             indexes[x]*indexes[x];
1052         }
1053       x++;
1054       p++;
1055     }
1056   }
1057   area=(double) image->columns*image->rows;
1058   for (i=0; i < AllChannels; i++)
1059   {
1060     channel_statistics[i].mean/=area;
1061     channel_statistics[i].standard_deviation/=area;
1062     channel_statistics[i].kurtosis/=area;
1063     channel_statistics[i].skewness/=area;
1064   }
1065   for (i=0; i < AllChannels; i++)
1066   {
1067     channel_statistics[AllChannels].depth=(unsigned long) MagickMax((double)
1068       channel_statistics[AllChannels].depth,(double)
1069       channel_statistics[i].depth);
1070     channel_statistics[AllChannels].minima=MagickMin(
1071       channel_statistics[AllChannels].minima,channel_statistics[i].minima);
1072     channel_statistics[AllChannels].maxima=MagickMax(
1073       channel_statistics[AllChannels].maxima,channel_statistics[i].maxima);
1074     channel_statistics[AllChannels].mean+=channel_statistics[i].mean;
1075     channel_statistics[AllChannels].standard_deviation+=
1076       channel_statistics[i].standard_deviation;
1077     channel_statistics[AllChannels].kurtosis+=channel_statistics[i].kurtosis;
1078     channel_statistics[AllChannels].skewness+=channel_statistics[i].skewness;
1079   }
1080   channels=4;
1081   if (image->colorspace == CMYKColorspace)
1082     channels++;
1083   channel_statistics[AllChannels].mean/=channels;
1084   channel_statistics[AllChannels].standard_deviation/=channels;
1085   channel_statistics[AllChannels].kurtosis/=channels;
1086   channel_statistics[AllChannels].skewness/=channels;
1087   for (i=0; i <= AllChannels; i++)
1088   {
1089     sum_squares=0.0;
1090     sum_squares=channel_statistics[i].standard_deviation;
1091     sum_cubes=0.0;
1092     sum_cubes=channel_statistics[i].skewness;
1093     channel_statistics[i].standard_deviation=sqrt(
1094       channel_statistics[i].standard_deviation-
1095        (channel_statistics[i].mean*channel_statistics[i].mean));
1096     if (channel_statistics[i].standard_deviation == 0.0)
1097       {
1098         channel_statistics[i].kurtosis=0.0;
1099         channel_statistics[i].skewness=0.0;
1100       }
1101     else
1102       {
1103         channel_statistics[i].skewness=(channel_statistics[i].skewness-
1104           3.0*channel_statistics[i].mean*sum_squares+
1105           2.0*channel_statistics[i].mean*channel_statistics[i].mean*
1106           channel_statistics[i].mean)/
1107           (channel_statistics[i].standard_deviation*
1108            channel_statistics[i].standard_deviation*
1109            channel_statistics[i].standard_deviation);
1110         channel_statistics[i].kurtosis=(channel_statistics[i].kurtosis-
1111           4.0*channel_statistics[i].mean*sum_cubes+
1112           6.0*channel_statistics[i].mean*channel_statistics[i].mean*sum_squares-
1113           3.0*channel_statistics[i].mean*channel_statistics[i].mean*
1114           1.0*channel_statistics[i].mean*channel_statistics[i].mean)/
1115           (channel_statistics[i].standard_deviation*
1116            channel_statistics[i].standard_deviation*
1117            channel_statistics[i].standard_deviation*
1118            channel_statistics[i].standard_deviation)-3.0;
1119       }
1120   }
1121   return(channel_statistics);
1122 }