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-2009 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(_OPENMP) && (_OPENMP >= 200203)
 231   #pragma omp parallel for schedule(static,1) 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(_OPENMP) && (_OPENMP >= 200203)
 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 B o u n d i n g B o x                                     %
 349 %                                                                             %
 350 %                                                                             %
 351 %                                                                             %
 352 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 353 %
 354 %  GetImageBoundingBox() returns the bounding box of an image canvas.
 355 %
 356 %  The format of the GetImageBoundingBox method is:
 357 %
 358 %      RectangleInfo GetImageBoundingBox(const Image *image,
 359 %        ExceptionInfo *exception)
 360 %
 361 %  A description of each parameter follows:
 362 %
 363 %    o bounds: Method GetImageBoundingBox returns the bounding box of an
 364 %      image canvas.
 365 %
 366 %    o image: the image.
 367 %
 368 %    o exception: return any errors or warnings in this structure.
 369 %
 370 */
 371 MagickExport RectangleInfo GetImageBoundingBox(const Image *image,
 372   ExceptionInfo *exception)
 373 {
 374   long
 375     y;
 376
 377   MagickBooleanType
 378     status;
 379
 380   MagickPixelPacket
 381     target[3],
 382     zero;
 383
 384   RectangleInfo
 385     bounds;
 386
 387   register const PixelPacket
 388     *p;
 389
 390   CacheView
 391     *image_view;
 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   bounds.width=0;
 398   bounds.height=0;
 399   bounds.x=(long) image->columns;
 400   bounds.y=(long) image->rows;
 401   GetMagickPixelPacket(image,&target[0]);
 402   image_view=AcquireCacheView(image);
 403   p=GetCacheViewVirtualPixels(image_view,0,0,1,1,exception);
 404   if (p == (const PixelPacket *) NULL)
 405     {
 406       image_view=DestroyCacheView(image_view);
 407       return(bounds);
 408     }
 409   SetMagickPixelPacket(image,p,GetCacheViewAuthenticIndexQueue(image_view),
 410     &target[0]);
 411   GetMagickPixelPacket(image,&target[1]);
 412   p=GetCacheViewVirtualPixels(image_view,(long) image->columns-1,0,1,1,
 413     exception);
 414   SetMagickPixelPacket(image,p,GetCacheViewAuthenticIndexQueue(image_view),
 415     &target[1]);
 416   GetMagickPixelPacket(image,&target[2]);
 417   p=GetCacheViewVirtualPixels(image_view,0,(long) image->rows-1,1,1,exception);
 418   SetMagickPixelPacket(image,p,GetCacheViewAuthenticIndexQueue(image_view),
 419     &target[2]);
 420   status=MagickTrue;
 421   GetMagickPixelPacket(image,&zero);
 422 #if defined(_OPENMP) && (_OPENMP >= 200203)
 423   #pragma omp parallel for schedule(static,1) shared(status)
 424 #endif
 425   for (y=0; y < (long) image->rows; y++)
 426   {
 427     MagickPixelPacket
 428       pixel;
 429
 430     RectangleInfo
 431       bounding_box;
 432
 433     register const IndexPacket
 434       *__restrict indexes;
 435
 436     register const PixelPacket
 437       *__restrict p;
 438
 439     register long
 440       x;
 441
 442     if (status == MagickFalse)
 443       continue;
 444 #if defined(HAVE_OPENMP)
 445 #  pragma omp critical (MagickCore_GetImageBoundingBox)
 446 #endif
 447     bounding_box=bounds;
 448     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
 449     if (p == (const PixelPacket *) NULL)
 450       {
 451         status=MagickFalse;
 452         continue;
 453       }
 454     indexes=GetCacheViewVirtualIndexQueue(image_view);
 455     pixel=zero;
 456     for (x=0; x < (long) image->columns; x++)
 457     {
 458       SetMagickPixelPacket(image,p,indexes+x,&pixel);
 459       if ((x < bounding_box.x) &&
 460           (IsMagickColorSimilar(&pixel,&target[0]) == MagickFalse))
 461         bounding_box.x=x;
 462       if ((x > (long) bounding_box.width) &&
 463           (IsMagickColorSimilar(&pixel,&target[1]) == MagickFalse))
 464         bounding_box.width=(unsigned long) x;
 465       if ((y < bounding_box.y) &&
 466           (IsMagickColorSimilar(&pixel,&target[0]) == MagickFalse))
 467         bounding_box.y=y;
 468       if ((y > (long) bounding_box.height) &&
 469           (IsMagickColorSimilar(&pixel,&target[2]) == MagickFalse))
 470         bounding_box.height=(unsigned long) y;
 471       p++;
 472     }
 473 #if defined(HAVE_OPENMP)
 474 #  pragma omp critical (MagickCore_GetImageBoundingBox)
 475 #endif
 476     {
 477       if (bounding_box.x < bounds.x)
 478         bounds.x=bounding_box.x;
 479       if (bounding_box.y < bounds.y)
 480         bounds.y=bounding_box.y;
 481       if (bounding_box.width > bounds.width)
 482         bounds.width=bounding_box.width;
 483       if (bounding_box.height > bounds.height)
 484         bounds.height=bounding_box.height;
 485     }
 486   }
 487   image_view=DestroyCacheView(image_view);
 488   if ((bounds.width == 0) || (bounds.height == 0))
 489     (void) ThrowMagickException(exception,GetMagickModule(),OptionWarning,
 490       "GeometryDoesNotContainImage","`%s'",image->filename);
 491   else
 492     {
 493       bounds.width-=(bounds.x-1);
 494       bounds.height-=(bounds.y-1);
 495     }
 496   return(bounds);
 497 }
 498 \f
 499 /*
 500 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 501 %                                                                             %
 502 %                                                                             %
 503 %                                                                             %
 504 %   G e t I m a g e C h a n n e l D e p t h                                   %
 505 %                                                                             %
 506 %                                                                             %
 507 %                                                                             %
 508 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 509 %
 510 %  GetImageChannelDepth() returns the depth of a particular image channel.
 511 %
 512 %  The format of the GetImageChannelDepth method is:
 513 %
 514 %      unsigned long GetImageDepth(const Image *image,ExceptionInfo *exception)
 515 %      unsigned long GetImageChannelDepth(const Image *image,
 516 %        const ChannelType channel,ExceptionInfo *exception)
 517 %
 518 %  A description of each parameter follows:
 519 %
 520 %    o image: the image.
 521 %
 522 %    o channel: the channel.
 523 %
 524 %    o exception: return any errors or warnings in this structure.
 525 %
 526 */
 527
 528 MagickExport unsigned long GetImageDepth(const Image *image,
 529   ExceptionInfo *exception)
 530 {
 531   return(GetImageChannelDepth(image,AllChannels,exception));
 532 }
 533
 534 MagickExport unsigned long GetImageChannelDepth(const Image *image,
 535   const ChannelType channel,ExceptionInfo *exception)
 536 {
 537   long
 538     y;
 539
 540   MagickBooleanType
 541     status;
 542
 543   register long
 544     id;
 545
 546   unsigned long
 547     *current_depth,
 548     depth,
 549     number_threads;
 550
 551   CacheView
 552     *image_view;
 553
 554   /*
 555     Compute image depth.
 556   */
 557   assert(image != (Image *) NULL);
 558   assert(image->signature == MagickSignature);
 559   if (image->debug != MagickFalse)
 560     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 561   number_threads=GetOpenMPMaximumThreads();
 562   current_depth=(unsigned long *) AcquireQuantumMemory(number_threads,
 563     sizeof(*current_depth));
 564   if (current_depth == (unsigned long *) NULL)
 565     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
 566   status=MagickTrue;
 567   for (id=0; id < (long) number_threads; id++)
 568     current_depth[id]=1;
 569   if ((image->storage_class == PseudoClass) && (image->matte == MagickFalse))
 570     {
 571       register const PixelPacket
 572         *__restrict p;
 573
 574       register long
 575         i;
 576
 577       p=image->colormap;
 578 #if defined(_OPENMP) && (_OPENMP >= 200203)
 579   #pragma omp parallel for schedule(static,1) shared(status)
 580 #endif
 581       for (i=0; i < (long) image->colors; i++)
 582       {
 583         if (status == MagickFalse)
 584           continue;
 585         id=GetOpenMPThreadId();
 586         while (current_depth[id] < MAGICKCORE_QUANTUM_DEPTH)
 587         {
 588           MagickStatusType
 589             status;
 590
 591           QuantumAny
 592             range;
 593
 594           status=0;
 595           range=GetQuantumRange(current_depth[id]);
 596           if ((channel & RedChannel) != 0)
 597             status|=p->red != ScaleAnyToQuantum(ScaleQuantumToAny(p->red,
 598               range),range);
 599           if ((channel & GreenChannel) != 0)
 600             status|=p->green != ScaleAnyToQuantum(ScaleQuantumToAny(p->green,
 601               range),range);
 602           if ((channel & BlueChannel) != 0)
 603             status|=p->blue != ScaleAnyToQuantum(ScaleQuantumToAny(p->blue,
 604               range),range);
 605           if (status == 0)
 606             break;
 607           current_depth[id]++;
 608         }
 609         p++;
 610       }
 611       depth=current_depth[0];
 612       for (id=1; id < (long) number_threads; id++)
 613         if (depth < current_depth[id])
 614           depth=current_depth[id];
 615       current_depth=(unsigned long *) RelinquishMagickMemory(current_depth);
 616       return(depth);
 617     }
 618   image_view=AcquireCacheView(image);
 619 #if defined(_OPENMP) && (_OPENMP >= 200203)
 620   #pragma omp parallel for schedule(static,1) shared(status)
 621 #endif
 622   for (y=0; y < (long) image->rows; y++)
 623   {
 624     register const IndexPacket
 625       *__restrict indexes;
 626
 627     register const PixelPacket
 628       *__restrict p;
 629
 630     register long
 631       id,
 632       x;
 633
 634     if (status == MagickFalse)
 635       continue;
 636     id=GetOpenMPThreadId();
 637     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
 638     if (p == (const PixelPacket *) NULL)
 639       continue;
 640     indexes=GetCacheViewVirtualIndexQueue(image_view);
 641     for (x=0; x < (long) image->columns; x++)
 642     {
 643       while (current_depth[id] < MAGICKCORE_QUANTUM_DEPTH)
 644       {
 645         MagickStatusType
 646           status;
 647
 648         QuantumAny
 649           range;
 650
 651         status=0;
 652         range=GetQuantumRange(current_depth[id]);
 653         if ((channel & RedChannel) != 0)
 654           status|=p->red != ScaleAnyToQuantum(ScaleQuantumToAny(p->red,range),
 655             range);
 656         if ((channel & GreenChannel) != 0)
 657           status|=p->green != ScaleAnyToQuantum(ScaleQuantumToAny(p->green,
 658             range),range);
 659         if ((channel & BlueChannel) != 0)
 660           status|=p->blue != ScaleAnyToQuantum(ScaleQuantumToAny(p->blue,range),
 661             range);
 662         if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
 663           status|=p->opacity != ScaleAnyToQuantum(ScaleQuantumToAny(p->opacity,
 664             range),range);
 665         if (((channel & IndexChannel) != 0) &&
 666             (image->colorspace == CMYKColorspace))
 667           status|=indexes[x] != ScaleAnyToQuantum(ScaleQuantumToAny(indexes[x],
 668             range),range);
 669         if (status == 0)
 670           break;
 671         current_depth[id]++;
 672       }
 673       p++;
 674     }
 675     if (current_depth[id] == MAGICKCORE_QUANTUM_DEPTH)
 676       status=MagickFalse;
 677   }
 678   image_view=DestroyCacheView(image_view);
 679   depth=current_depth[0];
 680   for (id=1; id < (long) number_threads; id++)
 681     if (depth < current_depth[id])
 682       depth=current_depth[id];
 683   current_depth=(unsigned long *) RelinquishMagickMemory(current_depth);
 684   return(depth);
 685 }
 686 \f
 687 /*
 688 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 689 %                                                                             %
 690 %                                                                             %
 691 %                                                                             %
 692 +   G e t I m a g e C h a n n e l E x t r e m a                               %
 693 %                                                                             %
 694 %                                                                             %
 695 %                                                                             %
 696 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 697 %
 698 %  GetImageChannelExtrema() returns the extrema of one or more image channels.
 699 %
 700 %  The format of the GetImageChannelExtrema method is:
 701 %
 702 %      MagickBooleanType GetImageChannelExtrema(const Image *image,
 703 %        const ChannelType channel,unsigned long *minima,unsigned long *maxima,
 704 %        ExceptionInfo *exception)
 705 %
 706 %  A description of each parameter follows:
 707 %
 708 %    o image: the image.
 709 %
 710 %    o channel: the channel.
 711 %
 712 %    o minima: the minimum value in the channel.
 713 %
 714 %    o maxima: the maximum value in the channel.
 715 %
 716 %    o exception: return any errors or warnings in this structure.
 717 %
 718 */
 719
 720 MagickExport MagickBooleanType GetImageExtrema(const Image *image,
 721   unsigned long *minima,unsigned long *maxima,ExceptionInfo *exception)
 722 {
 723   return(GetImageChannelExtrema(image,AllChannels,minima,maxima,exception));
 724 }
 725
 726 MagickExport MagickBooleanType GetImageChannelExtrema(const Image *image,
 727   const ChannelType channel,unsigned long *minima,unsigned long *maxima,
 728   ExceptionInfo *exception)
 729 {
 730   double
 731     max,
 732     min;
 733
 734   MagickBooleanType
 735     status;
 736
 737   assert(image != (Image *) NULL);
 738   assert(image->signature == MagickSignature);
 739   if (image->debug != MagickFalse)
 740     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 741   status=GetImageChannelRange(image,channel,&min,&max,exception);
 742   *minima=(unsigned long) (min+0.5);
 743   *maxima=(unsigned long) (max+0.5);
 744   return(status);
 745 }
 746 \f
 747 /*
 748 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 749 %                                                                             %
 750 %                                                                             %
 751 %                                                                             %
 752 %   G e t I m a g e C h a n n e l M e a n                                     %
 753 %                                                                             %
 754 %                                                                             %
 755 %                                                                             %
 756 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 757 %
 758 %  GetImageChannelMean() returns the mean and standard deviation of one or more
 759 %  image channels.
 760 %
 761 %  The format of the GetImageChannelMean method is:
 762 %
 763 %      MagickBooleanType GetImageChannelMean(const Image *image,
 764 %        const ChannelType channel,double *mean,double *standard_deviation,
 765 %        ExceptionInfo *exception)
 766 %
 767 %  A description of each parameter follows:
 768 %
 769 %    o image: the image.
 770 %
 771 %    o channel: the channel.
 772 %
 773 %    o mean: the average value in the channel.
 774 %
 775 %    o standard_deviation: the standard deviation of the channel.
 776 %
 777 %    o exception: return any errors or warnings in this structure.
 778 %
 779 */
 780
 781 MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean,
 782   double *standard_deviation,ExceptionInfo *exception)
 783 {
 784   MagickBooleanType
 785     status;
 786
 787   status=GetImageChannelMean(image,AllChannels,mean,standard_deviation,
 788     exception);
 789   return(status);
 790 }
 791
 792 MagickExport MagickBooleanType GetImageChannelMean(const Image *image,
 793   const ChannelType channel,double *mean,double *standard_deviation,
 794   ExceptionInfo *exception)
 795 {
 796   double
 797     area;
 798
 799   long
 800     y;
 801
 802   assert(image != (Image *) NULL);
 803   assert(image->signature == MagickSignature);
 804   if (image->debug != MagickFalse)
 805     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 806   *mean=0.0;
 807   *standard_deviation=0.0;
 808   area=0.0;
 809   for (y=0; y < (long) image->rows; y++)
 810   {
 811     register const IndexPacket
 812       *__restrict indexes;
 813
 814     register const PixelPacket
 815       *__restrict p;
 816
 817     register long
 818       x;
 819
 820     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
 821     if (p == (const PixelPacket *) NULL)
 822       break;
 823     indexes=GetVirtualIndexQueue(image);
 824     for (x=0; x < (long) image->columns; x++)
 825     {
 826       if ((channel & RedChannel) != 0)
 827         {
 828           *mean+=p->red;
 829           *standard_deviation+=(double) p->red*p->red;
 830           area++;
 831         }
 832       if ((channel & GreenChannel) != 0)
 833         {
 834           *mean+=p->green;
 835           *standard_deviation+=(double) p->green*p->green;
 836           area++;
 837         }
 838       if ((channel & BlueChannel) != 0)
 839         {
 840           *mean+=p->blue;
 841           *standard_deviation+=(double) p->blue*p->blue;
 842           area++;
 843         }
 844       if ((channel & OpacityChannel) != 0)
 845         {
 846           *mean+=p->opacity;
 847           *standard_deviation+=(double) p->opacity*p->opacity;
 848           area++;
 849         }
 850       if (((channel & IndexChannel) != 0) &&
 851           (image->colorspace == CMYKColorspace))
 852         {
 853           *mean+=indexes[x];
 854           *standard_deviation+=(double) indexes[x]*indexes[x];
 855           area++;
 856         }
 857       p++;
 858     }
 859   }
 860   if (y < (long) image->rows)
 861     return(MagickFalse);
 862   if (area != 0)
 863     {
 864       *mean/=area;
 865       *standard_deviation/=area;
 866     }
 867   *standard_deviation=sqrt(*standard_deviation-(*mean*(*mean)));
 868   return(y == (long) image->rows ? MagickTrue : MagickFalse);
 869 }
 870 \f
 871 /*
 872 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 873 %                                                                             %
 874 %                                                                             %
 875 %                                                                             %
 876 %   G e t I m a g e C h a n n e l K u r t o s i s                             %
 877 %                                                                             %
 878 %                                                                             %
 879 %                                                                             %
 880 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 881 %
 882 %  GetImageChannelKurtosis() returns the kurtosis and skewness of one or more
 883 %  image channels.
 884 %
 885 %  The format of the GetImageChannelKurtosis method is:
 886 %
 887 %      MagickBooleanType GetImageChannelKurtosis(const Image *image,
 888 %        const ChannelType channel,double *kurtosis,double *skewness,
 889 %        ExceptionInfo *exception)
 890 %
 891 %  A description of each parameter follows:
 892 %
 893 %    o image: the image.
 894 %
 895 %    o channel: the channel.
 896 %
 897 %    o kurtosis: the kurtosis of the channel.
 898 %
 899 %    o skewness: the skewness of the channel.
 900 %
 901 %    o exception: return any errors or warnings in this structure.
 902 %
 903 */
 904
 905 MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
 906   double *kurtosis,double *skewness,ExceptionInfo *exception)
 907 {
 908   MagickBooleanType
 909     status;
 910
 911   status=GetImageChannelKurtosis(image,AllChannels,kurtosis,skewness,
 912     exception);
 913   return(status);
 914 }
 915
 916 MagickExport MagickBooleanType GetImageChannelKurtosis(const Image *image,
 917   const ChannelType channel,double *kurtosis,double *skewness,
 918   ExceptionInfo *exception)
 919 {
 920   double
 921     area,
 922     mean,
 923     standard_deviation,
 924     sum_squares,
 925     sum_cubes,
 926     sum_fourth_power;
 927
 928   long
 929     y;
 930
 931   assert(image != (Image *) NULL);
 932   assert(image->signature == MagickSignature);
 933   if (image->debug != MagickFalse)
 934     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 935   *kurtosis=0.0;
 936   *skewness=0.0;
 937   area=0.0;
 938   mean=0.0;
 939   standard_deviation=0.0;
 940   sum_squares=0.0;
 941   sum_cubes=0.0;
 942   sum_fourth_power=0.0;
 943   for (y=0; y < (long) image->rows; y++)
 944   {
 945     register const IndexPacket
 946       *__restrict indexes;
 947
 948     register const PixelPacket
 949       *__restrict p;
 950
 951     register long
 952       x;
 953
 954     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
 955     if (p == (const PixelPacket *) NULL)
 956       break;
 957     indexes=GetVirtualIndexQueue(image);
 958     for (x=0; x < (long) image->columns; x++)
 959     {
 960       if ((channel & RedChannel) != 0)
 961         {
 962           mean+=p->red;
 963           sum_squares+=(double) p->red*p->red;
 964           sum_cubes+=(double) p->red*p->red*p->red;
 965           sum_fourth_power+=(double) p->red*p->red*p->red*p->red;
 966           area++;
 967         }
 968       if ((channel & GreenChannel) != 0)
 969         {
 970           mean+=p->green;
 971           sum_squares+=(double) p->green*p->green;
 972           sum_cubes+=(double) p->green*p->green*p->green;
 973           sum_fourth_power+=(double) p->green*p->green*p->green*p->green;
 974           area++;
 975         }
 976       if ((channel & BlueChannel) != 0)
 977         {
 978           mean+=p->blue;
 979           sum_squares+=(double) p->blue*p->blue;
 980           sum_cubes+=(double) p->blue*p->blue*p->blue;
 981           sum_fourth_power+=(double) p->blue*p->blue*p->blue*p->blue;
 982           area++;
 983         }
 984       if ((channel & OpacityChannel) != 0)
 985         {
 986           mean+=p->opacity;
 987           sum_squares+=(double) p->opacity*p->opacity;
 988           sum_cubes+=(double) p->opacity*p->opacity*p->opacity;
 989           sum_fourth_power+=(double) p->opacity*p->opacity*p->opacity*
 990             p->opacity;
 991           area++;
 992         }
 993       if (((channel & IndexChannel) != 0) &&
 994           (image->colorspace == CMYKColorspace))
 995         {
 996           mean+=indexes[x];
 997           sum_squares+=(double) indexes[x]*indexes[x];
 998           sum_cubes+=(double) indexes[x]*indexes[x]*indexes[x];
 999           sum_fourth_power+=(double) indexes[x]*indexes[x]*indexes[x]*
1000             indexes[x];
1001           area++;
1002         }
1003       p++;
1004     }
1005   }
1006   if (y < (long) image->rows)
1007     return(MagickFalse);
1008   if (area != 0.0)
1009     {
1010       mean/=area;
1011       sum_squares/=area;
1012       sum_cubes/=area;
1013       sum_fourth_power/=area;
1014     }
1015   standard_deviation=sqrt(sum_squares-(mean*mean));
1016   if (standard_deviation != 0.0)
1017     {
1018       *kurtosis=sum_fourth_power-4.0*mean*sum_cubes+6.0*mean*mean*sum_squares-
1019         3.0*mean*mean*mean*mean;
1020       *kurtosis/=standard_deviation*standard_deviation*standard_deviation*
1021         standard_deviation;
1022       *kurtosis-=3.0;
1023       *skewness=sum_cubes-3.0*mean*sum_squares+2.0*mean*mean*mean;
1024       *skewness/=standard_deviation*standard_deviation*standard_deviation;
1025     }
1026   return(y == (long) image->rows ? MagickTrue : MagickFalse);
1027 }
1028
1029 /*
1030 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1031 %                                                                             %
1032 %                                                                             %
1033 %                                                                             %
1034 %   G e t I m a g e C h a n n e l R a n g e                                   %
1035 %                                                                             %
1036 %                                                                             %
1037 %                                                                             %
1038 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1039 %
1040 %  GetImageChannelRange() returns the range of one or more image channels.
1041 %
1042 %  The format of the GetImageChannelRange method is:
1043 %
1044 %      MagickBooleanType GetImageChannelRange(const Image *image,
1045 %        const ChannelType channel,double *minima,double *maxima,
1046 %        ExceptionInfo *exception)
1047 %
1048 %  A description of each parameter follows:
1049 %
1050 %    o image: the image.
1051 %
1052 %    o channel: the channel.
1053 %
1054 %    o minima: the minimum value in the channel.
1055 %
1056 %    o maxima: the maximum value in the channel.
1057 %
1058 %    o exception: return any errors or warnings in this structure.
1059 %
1060 */
1061
1062 MagickExport MagickBooleanType GetImageRange(const Image *image,
1063   double *minima,double *maxima,ExceptionInfo *exception)
1064 {
1065   return(GetImageChannelRange(image,AllChannels,minima,maxima,exception));
1066 }
1067
1068 MagickExport MagickBooleanType GetImageChannelRange(const Image *image,
1069   const ChannelType channel,double *minima,double *maxima,
1070   ExceptionInfo *exception)
1071 {
1072   long
1073     y;
1074
1075   MagickPixelPacket
1076     pixel;
1077
1078   assert(image != (Image *) NULL);
1079   assert(image->signature == MagickSignature);
1080   if (image->debug != MagickFalse)
1081     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1082   *maxima=(-1.0E-37);
1083   *minima=1.0E+37;
1084   GetMagickPixelPacket(image,&pixel);
1085   for (y=0; y < (long) image->rows; y++)
1086   {
1087     register const IndexPacket
1088       *__restrict indexes;
1089
1090     register const PixelPacket
1091       *__restrict p;
1092
1093     register long
1094       x;
1095
1096     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1097     if (p == (const PixelPacket *) NULL)
1098       break;
1099     indexes=GetVirtualIndexQueue(image);
1100     for (x=0; x < (long) image->columns; x++)
1101     {
1102       SetMagickPixelPacket(image,p,indexes+x,&pixel);
1103       if ((channel & RedChannel) != 0)
1104         {
1105           if (pixel.red < *minima)
1106             *minima=(double) pixel.red;
1107           if (pixel.red > *maxima)
1108             *maxima=(double) pixel.red;
1109         }
1110       if ((channel & GreenChannel) != 0)
1111         {
1112           if (pixel.green < *minima)
1113             *minima=(double) pixel.green;
1114           if (pixel.green > *maxima)
1115             *maxima=(double) pixel.green;
1116         }
1117       if ((channel & BlueChannel) != 0)
1118         {
1119           if (pixel.blue < *minima)
1120             *minima=(double) pixel.blue;
1121           if (pixel.blue > *maxima)
1122             *maxima=(double) pixel.blue;
1123         }
1124       if ((channel & OpacityChannel) != 0)
1125         {
1126           if (pixel.opacity < *minima)
1127             *minima=(double) pixel.opacity;
1128           if (pixel.opacity > *maxima)
1129             *maxima=(double) pixel.opacity;
1130         }
1131       if (((channel & IndexChannel) != 0) &&
1132           (image->colorspace == CMYKColorspace))
1133         {
1134           if ((double) indexes[x] < *minima)
1135             *minima=(double) indexes[x];
1136           if ((double) indexes[x] > *maxima)
1137             *maxima=(double) indexes[x];
1138         }
1139       p++;
1140     }
1141   }
1142   return(y == (long) image->rows ? MagickTrue : MagickFalse);
1143 }
1144 \f
1145 /*
1146 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1147 %                                                                             %
1148 %                                                                             %
1149 %                                                                             %
1150 %   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                         %
1151 %                                                                             %
1152 %                                                                             %
1153 %                                                                             %
1154 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1155 %
1156 %  GetImageChannelStatistics() returns statistics for each channel in the
1157 %  image.  The statistics include the channel depth, its minima, maxima, mean,
1158 %  standard deviation, kurtosis and skewness.  You can access the red channel
1159 %  mean, for example, like this:
1160 %
1161 %      channel_statistics=GetImageChannelStatistics(image,excepton);
1162 %      red_mean=channel_statistics[RedChannel].mean;
1163 %
1164 %  Use MagickRelinquishMemory() to free the statistics buffer.
1165 %
1166 %  The format of the GetImageChannelStatistics method is:
1167 %
1168 %      ChannelStatistics *GetImageChannelStatistics(const Image *image,
1169 %        ExceptionInfo *exception)
1170 %
1171 %  A description of each parameter follows:
1172 %
1173 %    o image: the image.
1174 %
1175 %    o exception: return any errors or warnings in this structure.
1176 %
1177 */
1178
1179 static inline double MagickMax(const double x,const double y)
1180 {
1181   if (x > y)
1182     return(x);
1183   return(y);
1184 }
1185
1186 static inline double MagickMin(const double x,const double y)
1187 {
1188   if (x < y)
1189     return(x);
1190   return(y);
1191 }
1192
1193 MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
1194   ExceptionInfo *exception)
1195 {
1196   ChannelStatistics
1197     *channel_statistics;
1198
1199   double
1200     area,
1201     sum_squares,
1202     sum_cubes;
1203
1204   long
1205     y;
1206
1207   MagickStatusType
1208     status;
1209
1210   QuantumAny
1211     range;
1212
1213   register long
1214     i;
1215
1216   size_t
1217     length;
1218
1219   unsigned long
1220     channels,
1221     depth;
1222
1223   assert(image != (Image *) NULL);
1224   assert(image->signature == MagickSignature);
1225   if (image->debug != MagickFalse)
1226     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1227   length=AllChannels+1UL;
1228   channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(length,
1229     sizeof(*channel_statistics));
1230   if (channel_statistics == (ChannelStatistics *) NULL)
1231     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1232   (void) ResetMagickMemory(channel_statistics,0,length*
1233     sizeof(*channel_statistics));
1234   for (i=0; i <= AllChannels; i++)
1235   {
1236     channel_statistics[i].depth=1;
1237     channel_statistics[i].maxima=(-1.0E-37);
1238     channel_statistics[i].minima=1.0E+37;
1239     channel_statistics[i].mean=0.0;
1240     channel_statistics[i].standard_deviation=0.0;
1241     channel_statistics[i].kurtosis=0.0;
1242     channel_statistics[i].skewness=0.0;
1243   }
1244   for (y=0; y < (long) image->rows; y++)
1245   {
1246     register const IndexPacket
1247       *__restrict indexes;
1248
1249     register const PixelPacket
1250       *__restrict p;
1251
1252     register long
1253       x;
1254
1255     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1256     if (p == (const PixelPacket *) NULL)
1257       break;
1258     indexes=GetVirtualIndexQueue(image);
1259     for (x=0; x < (long) image->columns; )
1260     {
1261       if (channel_statistics[RedChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
1262         {
1263           depth=channel_statistics[RedChannel].depth;
1264           range=GetQuantumRange(depth);
1265           status=p->red != ScaleAnyToQuantum(ScaleQuantumToAny(p->red,range),
1266             range) ? MagickTrue : MagickFalse;
1267           if (status != MagickFalse)
1268             {
1269               channel_statistics[RedChannel].depth++;
1270               continue;
1271             }
1272         }
1273       if (channel_statistics[GreenChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
1274         {
1275           depth=channel_statistics[GreenChannel].depth;
1276           range=GetQuantumRange(depth);
1277           status=p->green != ScaleAnyToQuantum(ScaleQuantumToAny(p->green,
1278             range),range) ? MagickTrue : MagickFalse;
1279           if (status != MagickFalse)
1280             {
1281               channel_statistics[GreenChannel].depth++;
1282               continue;
1283             }
1284         }
1285       if (channel_statistics[BlueChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
1286         {
1287           depth=channel_statistics[BlueChannel].depth;
1288           range=GetQuantumRange(depth);
1289           status=p->blue != ScaleAnyToQuantum(ScaleQuantumToAny(p->blue,
1290             range),range) ? MagickTrue : MagickFalse;
1291           if (status != MagickFalse)
1292             {
1293               channel_statistics[BlueChannel].depth++;
1294               continue;
1295             }
1296         }
1297       if (image->matte != MagickFalse)
1298         {
1299           if (channel_statistics[OpacityChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
1300             {
1301               depth=channel_statistics[OpacityChannel].depth;
1302               range=GetQuantumRange(depth);
1303               status=p->opacity != ScaleAnyToQuantum(ScaleQuantumToAny(
1304                 p->opacity,range),range) ? MagickTrue : MagickFalse;
1305               if (status != MagickFalse)
1306                 {
1307                   channel_statistics[OpacityChannel].depth++;
1308                   continue;
1309                 }
1310             }
1311           }
1312       if (image->colorspace == CMYKColorspace)
1313         {
1314           if (channel_statistics[BlackChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
1315             {
1316               depth=channel_statistics[BlackChannel].depth;
1317               range=GetQuantumRange(depth);
1318               status=indexes[x] != ScaleAnyToQuantum(ScaleQuantumToAny(
1319                 indexes[x],range),range) ? MagickTrue : MagickFalse;
1320               if (status != MagickFalse)
1321                 {
1322                   channel_statistics[BlackChannel].depth++;
1323                   continue;
1324                 }
1325             }
1326         }
1327       if ((double) p->red < channel_statistics[RedChannel].minima)
1328         channel_statistics[RedChannel].minima=(double) p->red;
1329       if ((double) p->red > channel_statistics[RedChannel].maxima)
1330         channel_statistics[RedChannel].maxima=(double) p->red;
1331       channel_statistics[RedChannel].mean+=p->red;
1332       channel_statistics[RedChannel].standard_deviation+=(double) p->red*p->red;
1333       channel_statistics[RedChannel].kurtosis+=(double) p->red*p->red*
1334         p->red*p->red;
1335       channel_statistics[RedChannel].skewness+=(double) p->red*p->red*p->red;
1336       if ((double) p->green < channel_statistics[GreenChannel].minima)
1337         channel_statistics[GreenChannel].minima=(double) p->green;
1338       if ((double) p->green > channel_statistics[GreenChannel].maxima)
1339         channel_statistics[GreenChannel].maxima=(double) p->green;
1340       channel_statistics[GreenChannel].mean+=p->green;
1341       channel_statistics[GreenChannel].standard_deviation+=(double) p->green*
1342         p->green;
1343       channel_statistics[GreenChannel].kurtosis+=(double) p->green*p->green*
1344         p->green*p->green;
1345       channel_statistics[GreenChannel].skewness+=(double) p->green*p->green*
1346         p->green;
1347       if ((double) p->blue < channel_statistics[BlueChannel].minima)
1348         channel_statistics[BlueChannel].minima=(double) p->blue;
1349       if ((double) p->blue > channel_statistics[BlueChannel].maxima)
1350         channel_statistics[BlueChannel].maxima=(double) p->blue;
1351       channel_statistics[BlueChannel].mean+=p->blue;
1352       channel_statistics[BlueChannel].standard_deviation+=(double) p->blue*
1353         p->blue;
1354       channel_statistics[BlueChannel].kurtosis+=(double) p->blue*p->blue*
1355         p->blue*p->blue;
1356       channel_statistics[BlueChannel].skewness+=(double) p->blue*p->blue*
1357         p->blue;
1358       if (image->matte != MagickFalse)
1359         {
1360           if ((double) p->opacity < channel_statistics[OpacityChannel].minima)
1361             channel_statistics[OpacityChannel].minima=(double) p->opacity;
1362           if ((double) p->opacity > channel_statistics[OpacityChannel].maxima)
1363             channel_statistics[OpacityChannel].maxima=(double) p->opacity;
1364           channel_statistics[OpacityChannel].mean+=p->opacity;
1365           channel_statistics[OpacityChannel].standard_deviation+=(double)
1366             p->opacity*p->opacity;
1367           channel_statistics[OpacityChannel].kurtosis+=(double) p->opacity*
1368             p->opacity*p->opacity*p->opacity;
1369           channel_statistics[OpacityChannel].skewness+=(double) p->opacity*
1370             p->opacity*p->opacity;
1371         }
1372       if (image->colorspace == CMYKColorspace)
1373         {
1374           if ((double) indexes[x] < channel_statistics[BlackChannel].minima)
1375             channel_statistics[BlackChannel].minima=(double) indexes[x];
1376           if ((double) indexes[x] > channel_statistics[BlackChannel].maxima)
1377             channel_statistics[BlackChannel].maxima=(double) indexes[x];
1378           channel_statistics[BlackChannel].mean+=indexes[x];
1379           channel_statistics[BlackChannel].standard_deviation+=(double)
1380             indexes[x]*indexes[x];
1381           channel_statistics[BlackChannel].kurtosis+=(double) indexes[x]*
1382             indexes[x]*indexes[x]*indexes[x];
1383           channel_statistics[BlackChannel].skewness+=(double) indexes[x]*
1384             indexes[x]*indexes[x];
1385         }
1386       x++;
1387       p++;
1388     }
1389   }
1390   area=(double) image->columns*image->rows;
1391   for (i=0; i < AllChannels; i++)
1392   {
1393     channel_statistics[i].mean/=area;
1394     channel_statistics[i].standard_deviation/=area;
1395     channel_statistics[i].kurtosis/=area;
1396     channel_statistics[i].skewness/=area;
1397   }
1398   for (i=0; i < AllChannels; i++)
1399   {
1400     channel_statistics[AllChannels].depth=(unsigned long) MagickMax((double)
1401       channel_statistics[AllChannels].depth,(double)
1402       channel_statistics[i].depth);
1403     channel_statistics[AllChannels].minima=MagickMin(
1404       channel_statistics[AllChannels].minima,channel_statistics[i].minima);
1405     channel_statistics[AllChannels].maxima=MagickMax(
1406       channel_statistics[AllChannels].maxima,channel_statistics[i].maxima);
1407     channel_statistics[AllChannels].mean+=channel_statistics[i].mean;
1408     channel_statistics[AllChannels].standard_deviation+=
1409       channel_statistics[i].standard_deviation;
1410     channel_statistics[AllChannels].kurtosis+=channel_statistics[i].kurtosis;
1411     channel_statistics[AllChannels].skewness+=channel_statistics[i].skewness;
1412   }
1413   channels=4;
1414   if (image->colorspace == CMYKColorspace)
1415     channels++;
1416   channel_statistics[AllChannels].mean/=channels;
1417   channel_statistics[AllChannels].standard_deviation/=channels;
1418   channel_statistics[AllChannels].kurtosis/=channels;
1419   channel_statistics[AllChannels].skewness/=channels;
1420   for (i=0; i <= AllChannels; i++)
1421   {
1422     sum_squares=0.0;
1423     sum_squares=channel_statistics[i].standard_deviation;
1424     sum_cubes=0.0;
1425     sum_cubes=channel_statistics[i].skewness;
1426     channel_statistics[i].standard_deviation=sqrt(
1427       channel_statistics[i].standard_deviation-
1428        (channel_statistics[i].mean*channel_statistics[i].mean));
1429     if (channel_statistics[i].standard_deviation == 0.0)
1430       {
1431         channel_statistics[i].kurtosis=0.0;
1432         channel_statistics[i].skewness=0.0;
1433       }
1434     else
1435       {
1436         channel_statistics[i].skewness=(channel_statistics[i].skewness-
1437           3.0*channel_statistics[i].mean*sum_squares+
1438           2.0*channel_statistics[i].mean*channel_statistics[i].mean*
1439           channel_statistics[i].mean)/
1440           (channel_statistics[i].standard_deviation*
1441            channel_statistics[i].standard_deviation*
1442            channel_statistics[i].standard_deviation);
1443         channel_statistics[i].kurtosis=(channel_statistics[i].kurtosis-
1444           4.0*channel_statistics[i].mean*sum_cubes+
1445           6.0*channel_statistics[i].mean*channel_statistics[i].mean*sum_squares-
1446           3.0*channel_statistics[i].mean*channel_statistics[i].mean*
1447           1.0*channel_statistics[i].mean*channel_statistics[i].mean)/
1448           (channel_statistics[i].standard_deviation*
1449            channel_statistics[i].standard_deviation*
1450            channel_statistics[i].standard_deviation*
1451            channel_statistics[i].standard_deviation)-3.0;
1452       }
1453   }
1454   return(channel_statistics);
1455 }
1456 \f
1457 /*
1458 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1459 %                                                                             %
1460 %                                                                             %
1461 %                                                                             %
1462 %   G e t I m a g e Q u a n t u m D e p t h                                   %
1463 %                                                                             %
1464 %                                                                             %
1465 %                                                                             %
1466 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1467 %
1468 %  GetImageQuantumDepth() returns the depth of the image rounded to a legal
1469 %  quantum depth: 8, 16, or 32.
1470 %
1471 %  The format of the GetImageQuantumDepth method is:
1472 %
1473 %      unsigned long GetImageQuantumDepth(const Image *image,
1474 %        const MagickBooleanType constrain)
1475 %
1476 %  A description of each parameter follows:
1477 %
1478 %    o image: the image.
1479 %
1480 %    o constrain: A value other than MagickFalse, constrains the depth to
1481 %      a maximum of MAGICKCORE_QUANTUM_DEPTH.
1482 %
1483 */
1484 MagickExport unsigned long GetImageQuantumDepth(const Image *image,
1485   const MagickBooleanType constrain)
1486 {
1487   unsigned long
1488     depth;
1489
1490   depth=image->depth;
1491   if (depth <= 8)
1492     depth=8;
1493   else
1494     if (depth <= 16)
1495       depth=16;
1496     else
1497       if (depth <= 32)
1498         depth=32;
1499       else
1500         if (depth <= 64)
1501           depth=64;
1502   if (constrain != MagickFalse)
1503     depth=(unsigned long) MagickMin((double) depth,(double)
1504       MAGICKCORE_QUANTUM_DEPTH);
1505   return(depth);
1506 }
1507 \f
1508 /*
1509 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1510 %                                                                             %
1511 %                                                                             %
1512 %                                                                             %
1513 %   S e t I m a g e C h a n n e l D e p t h                                   %
1514 %                                                                             %
1515 %                                                                             %
1516 %                                                                             %
1517 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1518 %
1519 %  SetImageChannelDepth() sets the depth of the image.
1520 %
1521 %  The format of the SetImageChannelDepth method is:
1522 %
1523 %      MagickBooleanType SetImageDepth(Image *image,const unsigned long depth)
1524 %      MagickBooleanType SetImageChannelDepth(Image *image,
1525 %        const ChannelType channel,const unsigned long depth)
1526 %
1527 %  A description of each parameter follows:
1528 %
1529 %    o image: the image.
1530 %
1531 %    o channel: the channel.
1532 %
1533 %    o depth: the image depth.
1534 %
1535 */
1536
1537 MagickExport MagickBooleanType SetImageDepth(Image *image,
1538   const unsigned long depth)
1539 {
1540   return(SetImageChannelDepth(image,AllChannels,depth));
1541 }
1542
1543 MagickExport MagickBooleanType SetImageChannelDepth(Image *image,
1544   const ChannelType channel,const unsigned long depth)
1545 {
1546   ExceptionInfo
1547     *exception;
1548
1549   long
1550     y;
1551
1552   MagickBooleanType
1553     status;
1554
1555   QuantumAny
1556     range;
1557
1558   CacheView
1559     *image_view;
1560
1561   assert(image != (Image *) NULL);
1562   if (image->debug != MagickFalse)
1563     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
1564   assert(image->signature == MagickSignature);
1565   if (GetImageDepth(image,&image->exception) <= (unsigned long)
1566       MagickMin((double) depth,(double) MAGICKCORE_QUANTUM_DEPTH))
1567     {
1568       image->depth=depth;
1569       return(MagickTrue);
1570     }
1571   /*
1572     Scale pixels to desired depth.
1573   */
1574   status=MagickTrue;
1575   range=GetQuantumRange(depth);
1576   exception=(&image->exception);
1577   image_view=AcquireCacheView(image);
1578 #if defined(_OPENMP) && (_OPENMP >= 200203)
1579   #pragma omp parallel for schedule(static,1) shared(status)
1580 #endif
1581   for (y=0; y < (long) image->rows; y++)
1582   {
1583     register IndexPacket
1584       *__restrict indexes;
1585
1586     register long
1587       x;
1588
1589     register PixelPacket
1590       *__restrict q;
1591
1592     if (status == MagickFalse)
1593       continue;
1594     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,
1595       exception);
1596     if (q == (PixelPacket *) NULL)
1597       {
1598         status=MagickFalse;
1599         continue;
1600       }
1601     indexes=GetCacheViewAuthenticIndexQueue(image_view);
1602     for (x=0; x < (long) image->columns; x++)
1603     {
1604       if ((channel & RedChannel) != 0)
1605         q->red=ScaleAnyToQuantum(ScaleQuantumToAny(q->red,range),range);
1606       if ((channel & GreenChannel) != 0)
1607         q->green=ScaleAnyToQuantum(ScaleQuantumToAny(q->green,range),range);
1608       if ((channel & BlueChannel) != 0)
1609         q->blue=ScaleAnyToQuantum(ScaleQuantumToAny(q->blue,range),range);
1610       if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
1611         q->opacity=ScaleAnyToQuantum(ScaleQuantumToAny(q->opacity,range),range);
1612       if (((channel & IndexChannel) != 0) &&
1613           (image->colorspace == CMYKColorspace))
1614         indexes[x]=ScaleAnyToQuantum(ScaleQuantumToAny(indexes[x],range),range);
1615       q++;
1616     }
1617     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1618       {
1619         status=MagickFalse;
1620         continue;
1621       }
1622   }
1623   image_view=DestroyCacheView(image_view);
1624   if (image->storage_class == PseudoClass)
1625     {
1626       QuantumAny
1627         range;
1628
1629       register long
1630         i;
1631
1632       register PixelPacket
1633         *__restrict p;
1634
1635       p=image->colormap;
1636       range=GetQuantumRange(depth);
1637 #if defined(_OPENMP) && (_OPENMP >= 200203)
1638   #pragma omp parallel for schedule(static,1) shared(status)
1639 #endif
1640       for (i=0; i < (long) image->colors; i++)
1641       {
1642         if ((channel & RedChannel) != 0)
1643           p->red=ScaleAnyToQuantum(ScaleQuantumToAny(p->red,range),range);
1644         if ((channel & GreenChannel) != 0)
1645           p->green=ScaleAnyToQuantum(ScaleQuantumToAny(p->green,range),range);
1646         if ((channel & BlueChannel) != 0)
1647           p->blue=ScaleAnyToQuantum(ScaleQuantumToAny(p->blue,range),range);
1648         if ((channel & OpacityChannel) != 0)
1649           p->opacity=ScaleAnyToQuantum(ScaleQuantumToAny(p->opacity,range),
1650             range);
1651         p++;
1652       }
1653     }
1654   image->depth=depth;
1655   return(status);
1656 }