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