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     **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+=p->red;
 485           *standard_deviation+=(double) p->red*p->red;
 486           area++;
 487         }
 488       if ((channel & GreenChannel) != 0)
 489         {
 490           *mean+=p->green;
 491           *standard_deviation+=(double) p->green*p->green;
 492           area++;
 493         }
 494       if ((channel & BlueChannel) != 0)
 495         {
 496           *mean+=p->blue;
 497           *standard_deviation+=(double) p->blue*p->blue;
 498           area++;
 499         }
 500       if ((channel & OpacityChannel) != 0)
 501         {
 502           *mean+=p->opacity;
 503           *standard_deviation+=(double) p->opacity*p->opacity;
 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+=p->red;
 619           sum_squares+=(double) p->red*p->red;
 620           sum_cubes+=(double) p->red*p->red*p->red;
 621           sum_fourth_power+=(double) p->red*p->red*p->red*p->red;
 622           area++;
 623         }
 624       if ((channel & GreenChannel) != 0)
 625         {
 626           mean+=p->green;
 627           sum_squares+=(double) p->green*p->green;
 628           sum_cubes+=(double) p->green*p->green*p->green;
 629           sum_fourth_power+=(double) p->green*p->green*p->green*p->green;
 630           area++;
 631         }
 632       if ((channel & BlueChannel) != 0)
 633         {
 634           mean+=p->blue;
 635           sum_squares+=(double) p->blue*p->blue;
 636           sum_cubes+=(double) p->blue*p->blue*p->blue;
 637           sum_fourth_power+=(double) p->blue*p->blue*p->blue*p->blue;
 638           area++;
 639         }
 640       if ((channel & OpacityChannel) != 0)
 641         {
 642           mean+=p->opacity;
 643           sum_squares+=(double) p->opacity*p->opacity;
 644           sum_cubes+=(double) p->opacity*p->opacity*p->opacity;
 645           sum_fourth_power+=(double) p->opacity*p->opacity*p->opacity*
 646             p->opacity;
 647           area++;
 648         }
 649       if (((channel & IndexChannel) != 0) &&
 650           (image->colorspace == CMYKColorspace))
 651         {
 652           mean+=indexes[x];
 653           sum_squares+=(double) indexes[x]*indexes[x];
 654           sum_cubes+=(double) indexes[x]*indexes[x]*indexes[x];
 655           sum_fourth_power+=(double) indexes[x]*indexes[x]*indexes[x]*
 656             indexes[x];
 657           area++;
 658         }
 659       p++;
 660     }
 661   }
 662   if (y < (long) image->rows)
 663     return(MagickFalse);
 664   if (area != 0.0)
 665     {
 666       mean/=area;
 667       sum_squares/=area;
 668       sum_cubes/=area;
 669       sum_fourth_power/=area;
 670     }
 671   standard_deviation=sqrt(sum_squares-(mean*mean));
 672   if (standard_deviation != 0.0)
 673     {
 674       *kurtosis=sum_fourth_power-4.0*mean*sum_cubes+6.0*mean*mean*sum_squares-
 675         3.0*mean*mean*mean*mean;
 676       *kurtosis/=standard_deviation*standard_deviation*standard_deviation*
 677         standard_deviation;
 678       *kurtosis-=3.0;
 679       *skewness=sum_cubes-3.0*mean*sum_squares+2.0*mean*mean*mean;
 680       *skewness/=standard_deviation*standard_deviation*standard_deviation;
 681     }
 682   return(y == (long) image->rows ? MagickTrue : MagickFalse);
 683 }
 684
 685 /*
 686 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 687 %                                                                             %
 688 %                                                                             %
 689 %                                                                             %
 690 %   G e t I m a g e C h a n n e l R a n g e                                   %
 691 %                                                                             %
 692 %                                                                             %
 693 %                                                                             %
 694 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 695 %
 696 %  GetImageChannelRange() returns the range of one or more image channels.
 697 %
 698 %  The format of the GetImageChannelRange method is:
 699 %
 700 %      MagickBooleanType GetImageChannelRange(const Image *image,
 701 %        const ChannelType channel,double *minima,double *maxima,
 702 %        ExceptionInfo *exception)
 703 %
 704 %  A description of each parameter follows:
 705 %
 706 %    o image: the image.
 707 %
 708 %    o channel: the channel.
 709 %
 710 %    o minima: the minimum value in the channel.
 711 %
 712 %    o maxima: the maximum value in the channel.
 713 %
 714 %    o exception: return any errors or warnings in this structure.
 715 %
 716 */
 717
 718 MagickExport MagickBooleanType GetImageRange(const Image *image,
 719   double *minima,double *maxima,ExceptionInfo *exception)
 720 {
 721   return(GetImageChannelRange(image,AllChannels,minima,maxima,exception));
 722 }
 723
 724 MagickExport MagickBooleanType GetImageChannelRange(const Image *image,
 725   const ChannelType channel,double *minima,double *maxima,
 726   ExceptionInfo *exception)
 727 {
 728   long
 729     y;
 730
 731   MagickPixelPacket
 732     pixel;
 733
 734   assert(image != (Image *) NULL);
 735   assert(image->signature == MagickSignature);
 736   if (image->debug != MagickFalse)
 737     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 738   *maxima=(-1.0E-37);
 739   *minima=1.0E+37;
 740   GetMagickPixelPacket(image,&pixel);
 741   for (y=0; y < (long) image->rows; y++)
 742   {
 743     register const IndexPacket
 744       *restrict indexes;
 745
 746     register const PixelPacket
 747       *restrict p;
 748
 749     register long
 750       x;
 751
 752     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
 753     if (p == (const PixelPacket *) NULL)
 754       break;
 755     indexes=GetVirtualIndexQueue(image);
 756     for (x=0; x < (long) image->columns; x++)
 757     {
 758       SetMagickPixelPacket(image,p,indexes+x,&pixel);
 759       if ((channel & RedChannel) != 0)
 760         {
 761           if (pixel.red < *minima)
 762             *minima=(double) pixel.red;
 763           if (pixel.red > *maxima)
 764             *maxima=(double) pixel.red;
 765         }
 766       if ((channel & GreenChannel) != 0)
 767         {
 768           if (pixel.green < *minima)
 769             *minima=(double) pixel.green;
 770           if (pixel.green > *maxima)
 771             *maxima=(double) pixel.green;
 772         }
 773       if ((channel & BlueChannel) != 0)
 774         {
 775           if (pixel.blue < *minima)
 776             *minima=(double) pixel.blue;
 777           if (pixel.blue > *maxima)
 778             *maxima=(double) pixel.blue;
 779         }
 780       if ((channel & OpacityChannel) != 0)
 781         {
 782           if (pixel.opacity < *minima)
 783             *minima=(double) pixel.opacity;
 784           if (pixel.opacity > *maxima)
 785             *maxima=(double) pixel.opacity;
 786         }
 787       if (((channel & IndexChannel) != 0) &&
 788           (image->colorspace == CMYKColorspace))
 789         {
 790           if ((double) indexes[x] < *minima)
 791             *minima=(double) indexes[x];
 792           if ((double) indexes[x] > *maxima)
 793             *maxima=(double) indexes[x];
 794         }
 795       p++;
 796     }
 797   }
 798   return(y == (long) image->rows ? MagickTrue : MagickFalse);
 799 }
 800 \f
 801 /*
 802 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 803 %                                                                             %
 804 %                                                                             %
 805 %                                                                             %
 806 %   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                         %
 807 %                                                                             %
 808 %                                                                             %
 809 %                                                                             %
 810 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 811 %
 812 %  GetImageChannelStatistics() returns statistics for each channel in the
 813 %  image.  The statistics include the channel depth, its minima, maxima, mean,
 814 %  standard deviation, kurtosis and skewness.  You can access the red channel
 815 %  mean, for example, like this:
 816 %
 817 %      channel_statistics=GetImageChannelStatistics(image,excepton);
 818 %      red_mean=channel_statistics[RedChannel].mean;
 819 %
 820 %  Use MagickRelinquishMemory() to free the statistics buffer.
 821 %
 822 %  The format of the GetImageChannelStatistics method is:
 823 %
 824 %      ChannelStatistics *GetImageChannelStatistics(const Image *image,
 825 %        ExceptionInfo *exception)
 826 %
 827 %  A description of each parameter follows:
 828 %
 829 %    o image: the image.
 830 %
 831 %    o exception: return any errors or warnings in this structure.
 832 %
 833 */
 834
 835 static inline double MagickMax(const double x,const double y)
 836 {
 837   if (x > y)
 838     return(x);
 839   return(y);
 840 }
 841
 842 static inline double MagickMin(const double x,const double y)
 843 {
 844   if (x < y)
 845     return(x);
 846   return(y);
 847 }
 848
 849 MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
 850   ExceptionInfo *exception)
 851 {
 852   ChannelStatistics
 853     *channel_statistics;
 854
 855   double
 856     area,
 857     sum_squares,
 858     sum_cubes;
 859
 860   long
 861     y;
 862
 863   MagickStatusType
 864     status;
 865
 866   QuantumAny
 867     range;
 868
 869   register long
 870     i;
 871
 872   size_t
 873     length;
 874
 875   unsigned long
 876     channels,
 877     depth;
 878
 879   assert(image != (Image *) NULL);
 880   assert(image->signature == MagickSignature);
 881   if (image->debug != MagickFalse)
 882     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 883   length=AllChannels+1UL;
 884   channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(length,
 885     sizeof(*channel_statistics));
 886   if (channel_statistics == (ChannelStatistics *) NULL)
 887     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
 888   (void) ResetMagickMemory(channel_statistics,0,length*
 889     sizeof(*channel_statistics));
 890   for (i=0; i <= AllChannels; i++)
 891   {
 892     channel_statistics[i].depth=1;
 893     channel_statistics[i].maxima=(-1.0E-37);
 894     channel_statistics[i].minima=1.0E+37;
 895     channel_statistics[i].mean=0.0;
 896     channel_statistics[i].standard_deviation=0.0;
 897     channel_statistics[i].kurtosis=0.0;
 898     channel_statistics[i].skewness=0.0;
 899   }
 900   for (y=0; y < (long) image->rows; y++)
 901   {
 902     register const IndexPacket
 903       *restrict indexes;
 904
 905     register const PixelPacket
 906       *restrict p;
 907
 908     register long
 909       x;
 910
 911     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
 912     if (p == (const PixelPacket *) NULL)
 913       break;
 914     indexes=GetVirtualIndexQueue(image);
 915     for (x=0; x < (long) image->columns; )
 916     {
 917       if (channel_statistics[RedChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
 918         {
 919           depth=channel_statistics[RedChannel].depth;
 920           range=GetQuantumRange(depth);
 921           status=p->red != ScaleAnyToQuantum(ScaleQuantumToAny(p->red,range),
 922             range) ? MagickTrue : MagickFalse;
 923           if (status != MagickFalse)
 924             {
 925               channel_statistics[RedChannel].depth++;
 926               continue;
 927             }
 928         }
 929       if (channel_statistics[GreenChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
 930         {
 931           depth=channel_statistics[GreenChannel].depth;
 932           range=GetQuantumRange(depth);
 933           status=p->green != ScaleAnyToQuantum(ScaleQuantumToAny(p->green,
 934             range),range) ? MagickTrue : MagickFalse;
 935           if (status != MagickFalse)
 936             {
 937               channel_statistics[GreenChannel].depth++;
 938               continue;
 939             }
 940         }
 941       if (channel_statistics[BlueChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
 942         {
 943           depth=channel_statistics[BlueChannel].depth;
 944           range=GetQuantumRange(depth);
 945           status=p->blue != ScaleAnyToQuantum(ScaleQuantumToAny(p->blue,
 946             range),range) ? MagickTrue : MagickFalse;
 947           if (status != MagickFalse)
 948             {
 949               channel_statistics[BlueChannel].depth++;
 950               continue;
 951             }
 952         }
 953       if (image->matte != MagickFalse)
 954         {
 955           if (channel_statistics[OpacityChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
 956             {
 957               depth=channel_statistics[OpacityChannel].depth;
 958               range=GetQuantumRange(depth);
 959               status=p->opacity != ScaleAnyToQuantum(ScaleQuantumToAny(
 960                 p->opacity,range),range) ? MagickTrue : MagickFalse;
 961               if (status != MagickFalse)
 962                 {
 963                   channel_statistics[OpacityChannel].depth++;
 964                   continue;
 965                 }
 966             }
 967           }
 968       if (image->colorspace == CMYKColorspace)
 969         {
 970           if (channel_statistics[BlackChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
 971             {
 972               depth=channel_statistics[BlackChannel].depth;
 973               range=GetQuantumRange(depth);
 974               status=indexes[x] != ScaleAnyToQuantum(ScaleQuantumToAny(
 975                 indexes[x],range),range) ? MagickTrue : MagickFalse;
 976               if (status != MagickFalse)
 977                 {
 978                   channel_statistics[BlackChannel].depth++;
 979                   continue;
 980                 }
 981             }
 982         }
 983       if ((double) p->red < channel_statistics[RedChannel].minima)
 984         channel_statistics[RedChannel].minima=(double) p->red;
 985       if ((double) p->red > channel_statistics[RedChannel].maxima)
 986         channel_statistics[RedChannel].maxima=(double) p->red;
 987       channel_statistics[RedChannel].mean+=p->red;
 988       channel_statistics[RedChannel].standard_deviation+=(double) p->red*p->red;
 989       channel_statistics[RedChannel].kurtosis+=(double) p->red*p->red*
 990         p->red*p->red;
 991       channel_statistics[RedChannel].skewness+=(double) p->red*p->red*p->red;
 992       if ((double) p->green < channel_statistics[GreenChannel].minima)
 993         channel_statistics[GreenChannel].minima=(double) p->green;
 994       if ((double) p->green > channel_statistics[GreenChannel].maxima)
 995         channel_statistics[GreenChannel].maxima=(double) p->green;
 996       channel_statistics[GreenChannel].mean+=p->green;
 997       channel_statistics[GreenChannel].standard_deviation+=(double) p->green*
 998         p->green;
 999       channel_statistics[GreenChannel].kurtosis+=(double) p->green*p->green*
1000         p->green*p->green;
1001       channel_statistics[GreenChannel].skewness+=(double) p->green*p->green*
1002         p->green;
1003       if ((double) p->blue < channel_statistics[BlueChannel].minima)
1004         channel_statistics[BlueChannel].minima=(double) p->blue;
1005       if ((double) p->blue > channel_statistics[BlueChannel].maxima)
1006         channel_statistics[BlueChannel].maxima=(double) p->blue;
1007       channel_statistics[BlueChannel].mean+=p->blue;
1008       channel_statistics[BlueChannel].standard_deviation+=(double) p->blue*
1009         p->blue;
1010       channel_statistics[BlueChannel].kurtosis+=(double) p->blue*p->blue*
1011         p->blue*p->blue;
1012       channel_statistics[BlueChannel].skewness+=(double) p->blue*p->blue*
1013         p->blue;
1014       if (image->matte != MagickFalse)
1015         {
1016           if ((double) p->opacity < channel_statistics[OpacityChannel].minima)
1017             channel_statistics[OpacityChannel].minima=(double) p->opacity;
1018           if ((double) p->opacity > channel_statistics[OpacityChannel].maxima)
1019             channel_statistics[OpacityChannel].maxima=(double) p->opacity;
1020           channel_statistics[OpacityChannel].mean+=p->opacity;
1021           channel_statistics[OpacityChannel].standard_deviation+=(double)
1022             p->opacity*p->opacity;
1023           channel_statistics[OpacityChannel].kurtosis+=(double) p->opacity*
1024             p->opacity*p->opacity*p->opacity;
1025           channel_statistics[OpacityChannel].skewness+=(double) p->opacity*
1026             p->opacity*p->opacity;
1027         }
1028       if (image->colorspace == CMYKColorspace)
1029         {
1030           if ((double) indexes[x] < channel_statistics[BlackChannel].minima)
1031             channel_statistics[BlackChannel].minima=(double) indexes[x];
1032           if ((double) indexes[x] > channel_statistics[BlackChannel].maxima)
1033             channel_statistics[BlackChannel].maxima=(double) indexes[x];
1034           channel_statistics[BlackChannel].mean+=indexes[x];
1035           channel_statistics[BlackChannel].standard_deviation+=(double)
1036             indexes[x]*indexes[x];
1037           channel_statistics[BlackChannel].kurtosis+=(double) indexes[x]*
1038             indexes[x]*indexes[x]*indexes[x];
1039           channel_statistics[BlackChannel].skewness+=(double) indexes[x]*
1040             indexes[x]*indexes[x];
1041         }
1042       x++;
1043       p++;
1044     }
1045   }
1046   area=(double) image->columns*image->rows;
1047   for (i=0; i < AllChannels; i++)
1048   {
1049     channel_statistics[i].mean/=area;
1050     channel_statistics[i].standard_deviation/=area;
1051     channel_statistics[i].kurtosis/=area;
1052     channel_statistics[i].skewness/=area;
1053   }
1054   for (i=0; i < AllChannels; i++)
1055   {
1056     channel_statistics[AllChannels].depth=(unsigned long) MagickMax((double)
1057       channel_statistics[AllChannels].depth,(double)
1058       channel_statistics[i].depth);
1059     channel_statistics[AllChannels].minima=MagickMin(
1060       channel_statistics[AllChannels].minima,channel_statistics[i].minima);
1061     channel_statistics[AllChannels].maxima=MagickMax(
1062       channel_statistics[AllChannels].maxima,channel_statistics[i].maxima);
1063     channel_statistics[AllChannels].mean+=channel_statistics[i].mean;
1064     channel_statistics[AllChannels].standard_deviation+=
1065       channel_statistics[i].standard_deviation;
1066     channel_statistics[AllChannels].kurtosis+=channel_statistics[i].kurtosis;
1067     channel_statistics[AllChannels].skewness+=channel_statistics[i].skewness;
1068   }
1069   channels=4;
1070   if (image->colorspace == CMYKColorspace)
1071     channels++;
1072   channel_statistics[AllChannels].mean/=channels;
1073   channel_statistics[AllChannels].standard_deviation/=channels;
1074   channel_statistics[AllChannels].kurtosis/=channels;
1075   channel_statistics[AllChannels].skewness/=channels;
1076   for (i=0; i <= AllChannels; i++)
1077   {
1078     sum_squares=0.0;
1079     sum_squares=channel_statistics[i].standard_deviation;
1080     sum_cubes=0.0;
1081     sum_cubes=channel_statistics[i].skewness;
1082     channel_statistics[i].standard_deviation=sqrt(
1083       channel_statistics[i].standard_deviation-
1084        (channel_statistics[i].mean*channel_statistics[i].mean));
1085     if (channel_statistics[i].standard_deviation == 0.0)
1086       {
1087         channel_statistics[i].kurtosis=0.0;
1088         channel_statistics[i].skewness=0.0;
1089       }
1090     else
1091       {
1092         channel_statistics[i].skewness=(channel_statistics[i].skewness-
1093           3.0*channel_statistics[i].mean*sum_squares+
1094           2.0*channel_statistics[i].mean*channel_statistics[i].mean*
1095           channel_statistics[i].mean)/
1096           (channel_statistics[i].standard_deviation*
1097            channel_statistics[i].standard_deviation*
1098            channel_statistics[i].standard_deviation);
1099         channel_statistics[i].kurtosis=(channel_statistics[i].kurtosis-
1100           4.0*channel_statistics[i].mean*sum_cubes+
1101           6.0*channel_statistics[i].mean*channel_statistics[i].mean*sum_squares-
1102           3.0*channel_statistics[i].mean*channel_statistics[i].mean*
1103           1.0*channel_statistics[i].mean*channel_statistics[i].mean)/
1104           (channel_statistics[i].standard_deviation*
1105            channel_statistics[i].standard_deviation*
1106            channel_statistics[i].standard_deviation*
1107            channel_statistics[i].standard_deviation)-3.0;
1108       }
1109   }
1110   return(channel_statistics);
1111 }