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