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,ClampToQuantum(evaluate_pixel[x].black),
 687             q);
 688         if (evaluate_image->matte == MagickFalse)
 689           SetPixelAlpha(evaluate_image,ClampToQuantum(evaluate_pixel[x].alpha),
 690             q);
 691         else
 692           SetPixelAlpha(evaluate_image,ClampToQuantum(evaluate_pixel[x].alpha),
 693             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(ApplyEvaluateOperator(
 780           random_info[id],GetPixelRed(image,q),op,value)),q);
 781       if ((channel & GreenChannel) != 0)
 782         SetPixelGreen(image,ClampToQuantum(ApplyEvaluateOperator(
 783           random_info[id],GetPixelGreen(image,q),op,value)),q);
 784       if ((channel & BlueChannel) != 0)
 785         SetPixelBlue(image,ClampToQuantum(ApplyEvaluateOperator(
 786           random_info[id],GetPixelBlue(image,q),op,value)),q);
 787       if (((channel & BlackChannel) != 0) &&
 788           (image->colorspace == CMYKColorspace))
 789         SetPixelBlack(image,ClampToQuantum(ApplyEvaluateOperator(
 790           random_info[id],GetPixelBlack(image,q),op,value)),q);
 791       if ((channel & AlphaChannel) != 0)
 792         {
 793           if (image->matte == MagickFalse)
 794             SetPixelAlpha(image,ClampToQuantum(ApplyEvaluateOperator(
 795               random_info[id],GetPixelAlpha(image,q),op,value)),q);
 796           else
 797             SetPixelAlpha(image,ClampToQuantum(ApplyEvaluateOperator(
 798               random_info[id],GetPixelAlpha(image,q),op,value)),q);
 799         }
 800       q+=GetPixelChannels(image);
 801     }
 802     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
 803       status=MagickFalse;
 804     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 805       {
 806         MagickBooleanType
 807           proceed;
 808
 809 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 810   #pragma omp critical (MagickCore_EvaluateImageChannel)
 811 #endif
 812         proceed=SetImageProgress(image,EvaluateImageTag,progress++,image->rows);
 813         if (proceed == MagickFalse)
 814           status=MagickFalse;
 815       }
 816   }
 817   image_view=DestroyCacheView(image_view);
 818   random_info=DestroyRandomInfoThreadSet(random_info);
 819   return(status);
 820 }
 821 \f
 822 /*
 823 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 824 %                                                                             %
 825 %                                                                             %
 826 %                                                                             %
 827 %     F u n c t i o n I m a g e                                               %
 828 %                                                                             %
 829 %                                                                             %
 830 %                                                                             %
 831 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 832 %
 833 %  FunctionImage() applies a value to the image with an arithmetic, relational,
 834 %  or logical operator to an image. Use these operations to lighten or darken
 835 %  an image, to increase or decrease contrast in an image, or to produce the
 836 %  "negative" of an image.
 837 %
 838 %  The format of the FunctionImageChannel method is:
 839 %
 840 %      MagickBooleanType FunctionImage(Image *image,
 841 %        const MagickFunction function,const ssize_t number_parameters,
 842 %        const double *parameters,ExceptionInfo *exception)
 843 %      MagickBooleanType FunctionImageChannel(Image *image,
 844 %        const ChannelType channel,const MagickFunction function,
 845 %        const ssize_t number_parameters,const double *argument,
 846 %        ExceptionInfo *exception)
 847 %
 848 %  A description of each parameter follows:
 849 %
 850 %    o image: the image.
 851 %
 852 %    o channel: the channel.
 853 %
 854 %    o function: A channel function.
 855 %
 856 %    o parameters: one or more parameters.
 857 %
 858 %    o exception: return any errors or warnings in this structure.
 859 %
 860 */
 861
 862 static Quantum ApplyFunction(Quantum pixel,const MagickFunction function,
 863   const size_t number_parameters,const double *parameters,
 864   ExceptionInfo *exception)
 865 {
 866   MagickRealType
 867     result;
 868
 869   register ssize_t
 870     i;
 871
 872   (void) exception;
 873   result=0.0;
 874   switch (function)
 875   {
 876     case PolynomialFunction:
 877     {
 878       /*
 879        * Polynomial
 880        * Parameters:   polynomial constants,  highest to lowest order
 881        *   For example:      c0*x^3 + c1*x^2 + c2*x  + c3
 882        */
 883       result=0.0;
 884       for (i=0; i < (ssize_t) number_parameters; i++)
 885         result = result*QuantumScale*pixel + parameters[i];
 886       result *= QuantumRange;
 887       break;
 888     }
 889     case SinusoidFunction:
 890     {
 891       /* Sinusoid Function
 892        * Parameters:   Freq, Phase, Ampl, bias
 893        */
 894       double  freq,phase,ampl,bias;
 895       freq  = ( number_parameters >= 1 ) ? parameters[0] : 1.0;
 896       phase = ( number_parameters >= 2 ) ? parameters[1] : 0.0;
 897       ampl  = ( number_parameters >= 3 ) ? parameters[2] : 0.5;
 898       bias  = ( number_parameters >= 4 ) ? parameters[3] : 0.5;
 899       result=(MagickRealType) (QuantumRange*(ampl*sin((double) (2.0*MagickPI*
 900         (freq*QuantumScale*pixel + phase/360.0) )) + bias ) );
 901       break;
 902     }
 903     case ArcsinFunction:
 904     {
 905       /* Arcsin Function  (peged at range limits for invalid results)
 906        * Parameters:   Width, Center, Range, Bias
 907        */
 908       double  width,range,center,bias;
 909       width  = ( number_parameters >= 1 ) ? parameters[0] : 1.0;
 910       center = ( number_parameters >= 2 ) ? parameters[1] : 0.5;
 911       range  = ( number_parameters >= 3 ) ? parameters[2] : 1.0;
 912       bias   = ( number_parameters >= 4 ) ? parameters[3] : 0.5;
 913       result = 2.0/width*(QuantumScale*pixel - center);
 914       if ( result <= -1.0 )
 915         result = bias - range/2.0;
 916       else if ( result >= 1.0 )
 917         result = bias + range/2.0;
 918       else
 919         result=(MagickRealType) (range/MagickPI*asin((double) result)+bias);
 920       result *= QuantumRange;
 921       break;
 922     }
 923     case ArctanFunction:
 924     {
 925       /* Arctan Function
 926        * Parameters:   Slope, Center, Range, Bias
 927        */
 928       double  slope,range,center,bias;
 929       slope  = ( number_parameters >= 1 ) ? parameters[0] : 1.0;
 930       center = ( number_parameters >= 2 ) ? parameters[1] : 0.5;
 931       range  = ( number_parameters >= 3 ) ? parameters[2] : 1.0;
 932       bias   = ( number_parameters >= 4 ) ? parameters[3] : 0.5;
 933       result=(MagickRealType) (MagickPI*slope*(QuantumScale*pixel-center));
 934       result=(MagickRealType) (QuantumRange*(range/MagickPI*atan((double)
 935                   result) + bias ) );
 936       break;
 937     }
 938     case UndefinedFunction:
 939       break;
 940   }
 941   return(ClampToQuantum(result));
 942 }
 943
 944 MagickExport MagickBooleanType FunctionImage(Image *image,
 945   const MagickFunction function,const size_t number_parameters,
 946   const double *parameters,ExceptionInfo *exception)
 947 {
 948   MagickBooleanType
 949     status;
 950
 951   status=FunctionImageChannel(image,CompositeChannels,function,number_parameters,
 952     parameters,exception);
 953   return(status);
 954 }
 955
 956 MagickExport MagickBooleanType FunctionImageChannel(Image *image,
 957   const ChannelType channel,const MagickFunction function,
 958   const size_t number_parameters,const double *parameters,
 959   ExceptionInfo *exception)
 960 {
 961 #define FunctionImageTag  "Function/Image "
 962
 963   CacheView
 964     *image_view;
 965
 966   MagickBooleanType
 967     status;
 968
 969   MagickOffsetType
 970     progress;
 971
 972   ssize_t
 973     y;
 974
 975   assert(image != (Image *) NULL);
 976   assert(image->signature == MagickSignature);
 977   if (image->debug != MagickFalse)
 978     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 979   assert(exception != (ExceptionInfo *) NULL);
 980   assert(exception->signature == MagickSignature);
 981   if (SetImageStorageClass(image,DirectClass) == MagickFalse)
 982     {
 983       InheritException(exception,&image->exception);
 984       return(MagickFalse);
 985     }
 986   status=MagickTrue;
 987   progress=0;
 988   image_view=AcquireCacheView(image);
 989 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 990   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
 991 #endif
 992   for (y=0; y < (ssize_t) image->rows; y++)
 993   {
 994     register ssize_t
 995       x;
 996
 997     register Quantum
 998       *restrict q;
 999
1000     if (status == MagickFalse)
1001       continue;
1002     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
1003     if (q == (const Quantum *) NULL)
1004       {
1005         status=MagickFalse;
1006         continue;
1007       }
1008     for (x=0; x < (ssize_t) image->columns; x++)
1009     {
1010       if ((channel & RedChannel) != 0)
1011         SetPixelRed(image,ApplyFunction(GetPixelRed(image,q),function,
1012           number_parameters,parameters,exception),q);
1013       if ((channel & GreenChannel) != 0)
1014         SetPixelGreen(image,ApplyFunction(GetPixelGreen(image,q),function,
1015           number_parameters,parameters,exception),q);
1016       if ((channel & BlueChannel) != 0)
1017         SetPixelBlue(image,ApplyFunction(GetPixelBlue(image,q),function,
1018           number_parameters,parameters,exception),q);
1019       if (((channel & BlackChannel) != 0) &&
1020           (image->colorspace == CMYKColorspace))
1021         SetPixelBlack(image,ApplyFunction(GetPixelBlack(image,q),function,
1022           number_parameters,parameters,exception),q);
1023       if ((channel & AlphaChannel) != 0)
1024         {
1025           if (image->matte == MagickFalse)
1026             SetPixelAlpha(image,ApplyFunction(GetPixelAlpha(image,q),function,
1027               number_parameters,parameters,exception),q);
1028           else
1029             SetPixelAlpha(image,ApplyFunction(GetPixelAlpha(image,q),function,
1030               number_parameters,parameters,exception),q);
1031         }
1032       q+=GetPixelChannels(image);
1033     }
1034     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1035       status=MagickFalse;
1036     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1037       {
1038         MagickBooleanType
1039           proceed;
1040
1041 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1042   #pragma omp critical (MagickCore_FunctionImageChannel)
1043 #endif
1044         proceed=SetImageProgress(image,FunctionImageTag,progress++,image->rows);
1045         if (proceed == MagickFalse)
1046           status=MagickFalse;
1047       }
1048   }
1049   image_view=DestroyCacheView(image_view);
1050   return(status);
1051 }
1052 \f
1053 /*
1054 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1055 %                                                                             %
1056 %                                                                             %
1057 %                                                                             %
1058 +   G e t I m a g e C h a n n e l E x t r e m a                               %
1059 %                                                                             %
1060 %                                                                             %
1061 %                                                                             %
1062 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1063 %
1064 %  GetImageChannelExtrema() returns the extrema of one or more image channels.
1065 %
1066 %  The format of the GetImageChannelExtrema method is:
1067 %
1068 %      MagickBooleanType GetImageChannelExtrema(const Image *image,
1069 %        const ChannelType channel,size_t *minima,size_t *maxima,
1070 %        ExceptionInfo *exception)
1071 %
1072 %  A description of each parameter follows:
1073 %
1074 %    o image: the image.
1075 %
1076 %    o channel: the channel.
1077 %
1078 %    o minima: the minimum value in the channel.
1079 %
1080 %    o maxima: the maximum value in the channel.
1081 %
1082 %    o exception: return any errors or warnings in this structure.
1083 %
1084 */
1085
1086 MagickExport MagickBooleanType GetImageExtrema(const Image *image,
1087   size_t *minima,size_t *maxima,ExceptionInfo *exception)
1088 {
1089   return(GetImageChannelExtrema(image,CompositeChannels,minima,maxima,exception));
1090 }
1091
1092 MagickExport MagickBooleanType GetImageChannelExtrema(const Image *image,
1093   const ChannelType channel,size_t *minima,size_t *maxima,
1094   ExceptionInfo *exception)
1095 {
1096   double
1097     max,
1098     min;
1099
1100   MagickBooleanType
1101     status;
1102
1103   assert(image != (Image *) NULL);
1104   assert(image->signature == MagickSignature);
1105   if (image->debug != MagickFalse)
1106     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1107   status=GetImageChannelRange(image,channel,&min,&max,exception);
1108   *minima=(size_t) ceil(min-0.5);
1109   *maxima=(size_t) floor(max+0.5);
1110   return(status);
1111 }
1112 \f
1113 /*
1114 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1115 %                                                                             %
1116 %                                                                             %
1117 %                                                                             %
1118 %   G e t I m a g e C h a n n e l M e a n                                     %
1119 %                                                                             %
1120 %                                                                             %
1121 %                                                                             %
1122 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1123 %
1124 %  GetImageChannelMean() returns the mean and standard deviation of one or more
1125 %  image channels.
1126 %
1127 %  The format of the GetImageChannelMean method is:
1128 %
1129 %      MagickBooleanType GetImageChannelMean(const Image *image,
1130 %        const ChannelType channel,double *mean,double *standard_deviation,
1131 %        ExceptionInfo *exception)
1132 %
1133 %  A description of each parameter follows:
1134 %
1135 %    o image: the image.
1136 %
1137 %    o channel: the channel.
1138 %
1139 %    o mean: the average value in the channel.
1140 %
1141 %    o standard_deviation: the standard deviation of the channel.
1142 %
1143 %    o exception: return any errors or warnings in this structure.
1144 %
1145 */
1146
1147 MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean,
1148   double *standard_deviation,ExceptionInfo *exception)
1149 {
1150   MagickBooleanType
1151     status;
1152
1153   status=GetImageChannelMean(image,CompositeChannels,mean,standard_deviation,
1154     exception);
1155   return(status);
1156 }
1157
1158 MagickExport MagickBooleanType GetImageChannelMean(const Image *image,
1159   const ChannelType channel,double *mean,double *standard_deviation,
1160   ExceptionInfo *exception)
1161 {
1162   ChannelStatistics
1163     *channel_statistics;
1164
1165   size_t
1166     channels;
1167
1168   assert(image != (Image *) NULL);
1169   assert(image->signature == MagickSignature);
1170   if (image->debug != MagickFalse)
1171     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1172   channel_statistics=GetImageChannelStatistics(image,exception);
1173   if (channel_statistics == (ChannelStatistics *) NULL)
1174     return(MagickFalse);
1175   channels=0;
1176   channel_statistics[CompositeChannels].mean=0.0;
1177   channel_statistics[CompositeChannels].standard_deviation=0.0;
1178   if ((channel & RedChannel) != 0)
1179     {
1180       channel_statistics[CompositeChannels].mean+=
1181         channel_statistics[RedChannel].mean;
1182       channel_statistics[CompositeChannels].standard_deviation+=
1183         channel_statistics[RedChannel].variance-
1184         channel_statistics[RedChannel].mean*
1185         channel_statistics[RedChannel].mean;
1186       channels++;
1187     }
1188   if ((channel & GreenChannel) != 0)
1189     {
1190       channel_statistics[CompositeChannels].mean+=
1191         channel_statistics[GreenChannel].mean;
1192       channel_statistics[CompositeChannels].standard_deviation+=
1193         channel_statistics[GreenChannel].variance-
1194         channel_statistics[GreenChannel].mean*
1195         channel_statistics[GreenChannel].mean;
1196       channels++;
1197     }
1198   if ((channel & BlueChannel) != 0)
1199     {
1200       channel_statistics[CompositeChannels].mean+=
1201         channel_statistics[BlueChannel].mean;
1202       channel_statistics[CompositeChannels].standard_deviation+=
1203         channel_statistics[BlueChannel].variance-
1204         channel_statistics[BlueChannel].mean*
1205         channel_statistics[BlueChannel].mean;
1206       channels++;
1207     }
1208   if (((channel & BlackChannel) != 0) &&
1209       (image->colorspace == CMYKColorspace))
1210     {
1211       channel_statistics[CompositeChannels].mean+=
1212         channel_statistics[BlackChannel].mean;
1213       channel_statistics[CompositeChannels].standard_deviation+=
1214         channel_statistics[BlackChannel].variance-
1215         channel_statistics[BlackChannel].mean*
1216         channel_statistics[BlackChannel].mean;
1217       channels++;
1218     }
1219   if (((channel & AlphaChannel) != 0) &&
1220       (image->matte != MagickFalse))
1221     {
1222       channel_statistics[CompositeChannels].mean+=
1223         channel_statistics[AlphaChannel].mean;
1224       channel_statistics[CompositeChannels].standard_deviation+=
1225         channel_statistics[AlphaChannel].variance-
1226         channel_statistics[AlphaChannel].mean*
1227         channel_statistics[AlphaChannel].mean;
1228       channels++;
1229     }
1230   channel_statistics[CompositeChannels].mean/=channels;
1231   channel_statistics[CompositeChannels].standard_deviation=
1232     sqrt(channel_statistics[CompositeChannels].standard_deviation/channels);
1233   *mean=channel_statistics[CompositeChannels].mean;
1234   *standard_deviation=channel_statistics[CompositeChannels].standard_deviation;
1235   channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1236     channel_statistics);
1237   return(MagickTrue);
1238 }
1239 \f
1240 /*
1241 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1242 %                                                                             %
1243 %                                                                             %
1244 %                                                                             %
1245 %   G e t I m a g e C h a n n e l K u r t o s i s                             %
1246 %                                                                             %
1247 %                                                                             %
1248 %                                                                             %
1249 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1250 %
1251 %  GetImageChannelKurtosis() returns the kurtosis and skewness of one or more
1252 %  image channels.
1253 %
1254 %  The format of the GetImageChannelKurtosis method is:
1255 %
1256 %      MagickBooleanType GetImageChannelKurtosis(const Image *image,
1257 %        const ChannelType channel,double *kurtosis,double *skewness,
1258 %        ExceptionInfo *exception)
1259 %
1260 %  A description of each parameter follows:
1261 %
1262 %    o image: the image.
1263 %
1264 %    o channel: the channel.
1265 %
1266 %    o kurtosis: the kurtosis of the channel.
1267 %
1268 %    o skewness: the skewness of the channel.
1269 %
1270 %    o exception: return any errors or warnings in this structure.
1271 %
1272 */
1273
1274 MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
1275   double *kurtosis,double *skewness,ExceptionInfo *exception)
1276 {
1277   MagickBooleanType
1278     status;
1279
1280   status=GetImageChannelKurtosis(image,CompositeChannels,kurtosis,skewness,
1281     exception);
1282   return(status);
1283 }
1284
1285 MagickExport MagickBooleanType GetImageChannelKurtosis(const Image *image,
1286   const ChannelType channel,double *kurtosis,double *skewness,
1287   ExceptionInfo *exception)
1288 {
1289   double
1290     area,
1291     mean,
1292     standard_deviation,
1293     sum_squares,
1294     sum_cubes,
1295     sum_fourth_power;
1296
1297   ssize_t
1298     y;
1299
1300   assert(image != (Image *) NULL);
1301   assert(image->signature == MagickSignature);
1302   if (image->debug != MagickFalse)
1303     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1304   *kurtosis=0.0;
1305   *skewness=0.0;
1306   area=0.0;
1307   mean=0.0;
1308   standard_deviation=0.0;
1309   sum_squares=0.0;
1310   sum_cubes=0.0;
1311   sum_fourth_power=0.0;
1312   for (y=0; y < (ssize_t) image->rows; y++)
1313   {
1314     register const Quantum
1315       *restrict p;
1316
1317     register ssize_t
1318       x;
1319
1320     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1321     if (p == (const Quantum *) NULL)
1322       break;
1323     for (x=0; x < (ssize_t) image->columns; x++)
1324     {
1325       if ((channel & RedChannel) != 0)
1326         {
1327           mean+=GetPixelRed(image,p);
1328           sum_squares+=(double) GetPixelRed(image,p)*GetPixelRed(image,p);
1329           sum_cubes+=(double) GetPixelRed(image,p)*GetPixelRed(image,p)*
1330             GetPixelRed(image,p);
1331           sum_fourth_power+=(double) GetPixelRed(image,p)*
1332             GetPixelRed(image,p)*GetPixelRed(image,p)*GetPixelRed(image,p);
1333           area++;
1334         }
1335       if ((channel & GreenChannel) != 0)
1336         {
1337           mean+=GetPixelGreen(image,p);
1338           sum_squares+=(double) GetPixelGreen(image,p)*GetPixelGreen(image,p);
1339           sum_cubes+=(double) GetPixelGreen(image,p)*GetPixelGreen(image,p)*
1340             GetPixelGreen(image,p);
1341           sum_fourth_power+=(double) GetPixelGreen(image,p)*
1342             GetPixelGreen(image,p)*GetPixelGreen(image,p)*
1343             GetPixelGreen(image,p);
1344           area++;
1345         }
1346       if ((channel & BlueChannel) != 0)
1347         {
1348           mean+=GetPixelBlue(image,p);
1349           sum_squares+=(double) GetPixelBlue(image,p)*GetPixelBlue(image,p);
1350           sum_cubes+=(double) GetPixelBlue(image,p)*GetPixelBlue(image,p)*
1351             GetPixelBlue(image,p);
1352           sum_fourth_power+=(double) GetPixelBlue(image,p)*
1353             GetPixelBlue(image,p)*GetPixelBlue(image,p)*
1354             GetPixelBlue(image,p);
1355           area++;
1356         }
1357       if (((channel & BlackChannel) != 0) &&
1358           (image->colorspace == CMYKColorspace))
1359         {
1360           mean+=GetPixelBlack(image,p);
1361           sum_squares+=(double) GetPixelBlack(image,p)*GetPixelBlack(image,p);
1362           sum_cubes+=(double) GetPixelBlack(image,p)*GetPixelBlack(image,p)*
1363             GetPixelBlack(image,p);
1364           sum_fourth_power+=(double) GetPixelBlack(image,p)*
1365             GetPixelBlack(image,p)*GetPixelBlack(image,p)*
1366             GetPixelBlack(image,p);
1367           area++;
1368         }
1369       if ((channel & AlphaChannel) != 0)
1370         {
1371           mean+=GetPixelAlpha(image,p);
1372           sum_squares+=(double) GetPixelAlpha(image,p)*GetPixelAlpha(image,p);
1373           sum_cubes+=(double) GetPixelAlpha(image,p)*GetPixelAlpha(image,p)*
1374             GetPixelAlpha(image,p);
1375           sum_fourth_power+=(double) GetPixelAlpha(image,p)*
1376             GetPixelAlpha(image,p)*GetPixelAlpha(image,p)*
1377             GetPixelAlpha(image,p);
1378           area++;
1379         }
1380       p+=GetPixelChannels(image);
1381     }
1382   }
1383   if (y < (ssize_t) image->rows)
1384     return(MagickFalse);
1385   if (area != 0.0)
1386     {
1387       mean/=area;
1388       sum_squares/=area;
1389       sum_cubes/=area;
1390       sum_fourth_power/=area;
1391     }
1392   standard_deviation=sqrt(sum_squares-(mean*mean));
1393   if (standard_deviation != 0.0)
1394     {
1395       *kurtosis=sum_fourth_power-4.0*mean*sum_cubes+6.0*mean*mean*sum_squares-
1396         3.0*mean*mean*mean*mean;
1397       *kurtosis/=standard_deviation*standard_deviation*standard_deviation*
1398         standard_deviation;
1399       *kurtosis-=3.0;
1400       *skewness=sum_cubes-3.0*mean*sum_squares+2.0*mean*mean*mean;
1401       *skewness/=standard_deviation*standard_deviation*standard_deviation;
1402     }
1403   return(y == (ssize_t) image->rows ? MagickTrue : MagickFalse);
1404 }
1405 \f
1406 /*
1407 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1408 %                                                                             %
1409 %                                                                             %
1410 %                                                                             %
1411 %   G e t I m a g e C h a n n e l R a n g e                                   %
1412 %                                                                             %
1413 %                                                                             %
1414 %                                                                             %
1415 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1416 %
1417 %  GetImageChannelRange() returns the range of one or more image channels.
1418 %
1419 %  The format of the GetImageChannelRange method is:
1420 %
1421 %      MagickBooleanType GetImageChannelRange(const Image *image,
1422 %        const ChannelType channel,double *minima,double *maxima,
1423 %        ExceptionInfo *exception)
1424 %
1425 %  A description of each parameter follows:
1426 %
1427 %    o image: the image.
1428 %
1429 %    o channel: the channel.
1430 %
1431 %    o minima: the minimum value in the channel.
1432 %
1433 %    o maxima: the maximum value in the channel.
1434 %
1435 %    o exception: return any errors or warnings in this structure.
1436 %
1437 */
1438
1439 MagickExport MagickBooleanType GetImageRange(const Image *image,
1440   double *minima,double *maxima,ExceptionInfo *exception)
1441 {
1442   return(GetImageChannelRange(image,CompositeChannels,minima,maxima,exception));
1443 }
1444
1445 MagickExport MagickBooleanType GetImageChannelRange(const Image *image,
1446   const ChannelType channel,double *minima,double *maxima,
1447   ExceptionInfo *exception)
1448 {
1449   PixelInfo
1450     pixel;
1451
1452   ssize_t
1453     y;
1454
1455   assert(image != (Image *) NULL);
1456   assert(image->signature == MagickSignature);
1457   if (image->debug != MagickFalse)
1458     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1459   *maxima=(-1.0E-37);
1460   *minima=1.0E+37;
1461   GetPixelInfo(image,&pixel);
1462   for (y=0; y < (ssize_t) image->rows; y++)
1463   {
1464     register const Quantum
1465       *restrict p;
1466
1467     register ssize_t
1468       x;
1469
1470     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1471     if (p == (const Quantum *) NULL)
1472       break;
1473     for (x=0; x < (ssize_t) image->columns; x++)
1474     {
1475       SetPixelInfo(image,p,&pixel);
1476       if ((channel & RedChannel) != 0)
1477         {
1478           if (pixel.red < *minima)
1479             *minima=(double) pixel.red;
1480           if (pixel.red > *maxima)
1481             *maxima=(double) pixel.red;
1482         }
1483       if ((channel & GreenChannel) != 0)
1484         {
1485           if (pixel.green < *minima)
1486             *minima=(double) pixel.green;
1487           if (pixel.green > *maxima)
1488             *maxima=(double) pixel.green;
1489         }
1490       if ((channel & BlueChannel) != 0)
1491         {
1492           if (pixel.blue < *minima)
1493             *minima=(double) pixel.blue;
1494           if (pixel.blue > *maxima)
1495             *maxima=(double) pixel.blue;
1496         }
1497       if (((channel & BlackChannel) != 0) &&
1498           (image->colorspace == CMYKColorspace))
1499         {
1500           if (pixel.black < *minima)
1501             *minima=(double) pixel.black;
1502           if (pixel.black > *maxima)
1503             *maxima=(double) pixel.black;
1504         }
1505       if ((channel & AlphaChannel) != 0)
1506         {
1507           if (pixel.alpha < *minima)
1508             *minima=(double) pixel.alpha;
1509           if (pixel.alpha > *maxima)
1510             *maxima=(double) pixel.alpha;
1511         }
1512       p+=GetPixelChannels(image);
1513     }
1514   }
1515   return(y == (ssize_t) image->rows ? MagickTrue : MagickFalse);
1516 }
1517 \f
1518 /*
1519 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1520 %                                                                             %
1521 %                                                                             %
1522 %                                                                             %
1523 %   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                         %
1524 %                                                                             %
1525 %                                                                             %
1526 %                                                                             %
1527 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1528 %
1529 %  GetImageChannelStatistics() returns statistics for each channel in the
1530 %  image.  The statistics include the channel depth, its minima, maxima, mean,
1531 %  standard deviation, kurtosis and skewness.  You can access the red channel
1532 %  mean, for example, like this:
1533 %
1534 %      channel_statistics=GetImageChannelStatistics(image,exception);
1535 %      red_mean=channel_statistics[RedChannel].mean;
1536 %
1537 %  Use MagickRelinquishMemory() to free the statistics buffer.
1538 %
1539 %  The format of the GetImageChannelStatistics method is:
1540 %
1541 %      ChannelStatistics *GetImageChannelStatistics(const Image *image,
1542 %        ExceptionInfo *exception)
1543 %
1544 %  A description of each parameter follows:
1545 %
1546 %    o image: the image.
1547 %
1548 %    o exception: return any errors or warnings in this structure.
1549 %
1550 */
1551 MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
1552   ExceptionInfo *exception)
1553 {
1554   ChannelStatistics
1555     *channel_statistics;
1556
1557   double
1558     area;
1559
1560   MagickStatusType
1561     status;
1562
1563   QuantumAny
1564     range;
1565
1566   register ssize_t
1567     i;
1568
1569   size_t
1570     channels,
1571     depth,
1572     length;
1573
1574   ssize_t
1575     y;
1576
1577   assert(image != (Image *) NULL);
1578   assert(image->signature == MagickSignature);
1579   if (image->debug != MagickFalse)
1580     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1581   length=CompositeChannels+1UL;
1582   channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(length,
1583     sizeof(*channel_statistics));
1584   if (channel_statistics == (ChannelStatistics *) NULL)
1585     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1586   (void) ResetMagickMemory(channel_statistics,0,length*
1587     sizeof(*channel_statistics));
1588   for (i=0; i <= (ssize_t) CompositeChannels; i++)
1589   {
1590     channel_statistics[i].depth=1;
1591     channel_statistics[i].maxima=(-1.0E-37);
1592     channel_statistics[i].minima=1.0E+37;
1593   }
1594   for (y=0; y < (ssize_t) image->rows; y++)
1595   {
1596     register const Quantum
1597       *restrict p;
1598
1599     register ssize_t
1600       x;
1601
1602     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1603     if (p == (const Quantum *) NULL)
1604       break;
1605     for (x=0; x < (ssize_t) image->columns; )
1606     {
1607       if (channel_statistics[RedChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
1608         {
1609           depth=channel_statistics[RedChannel].depth;
1610           range=GetQuantumRange(depth);
1611           status=GetPixelRed(image,p) != ScaleAnyToQuantum(ScaleQuantumToAny(
1612             GetPixelRed(image,p),range),range) ? MagickTrue : MagickFalse;
1613           if (status != MagickFalse)
1614             {
1615               channel_statistics[RedChannel].depth++;
1616               continue;
1617             }
1618         }
1619       if (channel_statistics[GreenChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
1620         {
1621           depth=channel_statistics[GreenChannel].depth;
1622           range=GetQuantumRange(depth);
1623           status=GetPixelGreen(image,p) != ScaleAnyToQuantum(ScaleQuantumToAny(
1624             GetPixelGreen(image,p),range),range) ? MagickTrue : MagickFalse;
1625           if (status != MagickFalse)
1626             {
1627               channel_statistics[GreenChannel].depth++;
1628               continue;
1629             }
1630         }
1631       if (channel_statistics[BlueChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
1632         {
1633           depth=channel_statistics[BlueChannel].depth;
1634           range=GetQuantumRange(depth);
1635           status=GetPixelBlue(image,p) != ScaleAnyToQuantum(ScaleQuantumToAny(
1636             GetPixelBlue(image,p),range),range) ? MagickTrue : MagickFalse;
1637           if (status != MagickFalse)
1638             {
1639               channel_statistics[BlueChannel].depth++;
1640               continue;
1641             }
1642         }
1643       if (image->matte != MagickFalse)
1644         {
1645           if (channel_statistics[AlphaChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
1646             {
1647               depth=channel_statistics[AlphaChannel].depth;
1648               range=GetQuantumRange(depth);
1649               status=GetPixelAlpha(image,p) != ScaleAnyToQuantum(
1650                 ScaleQuantumToAny(GetPixelAlpha(image,p),range),range) ?
1651                 MagickTrue : MagickFalse;
1652               if (status != MagickFalse)
1653                 {
1654                   channel_statistics[AlphaChannel].depth++;
1655                   continue;
1656                 }
1657             }
1658           }
1659       if (image->colorspace == CMYKColorspace)
1660         {
1661           if (channel_statistics[BlackChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
1662             {
1663               depth=channel_statistics[BlackChannel].depth;
1664               range=GetQuantumRange(depth);
1665               status=GetPixelBlack(image,p) != ScaleAnyToQuantum(
1666                 ScaleQuantumToAny(GetPixelBlack(image,p),range),range) ?
1667                 MagickTrue : MagickFalse;
1668               if (status != MagickFalse)
1669                 {
1670                   channel_statistics[BlackChannel].depth++;
1671                   continue;
1672                 }
1673             }
1674         }
1675       if ((double) GetPixelRed(image,p) < channel_statistics[RedChannel].minima)
1676         channel_statistics[RedChannel].minima=(double) GetPixelRed(image,p);
1677       if ((double) GetPixelRed(image,p) > channel_statistics[RedChannel].maxima)
1678         channel_statistics[RedChannel].maxima=(double) GetPixelRed(image,p);
1679       channel_statistics[RedChannel].sum+=GetPixelRed(image,p);
1680       channel_statistics[RedChannel].sum_squared+=(double)
1681         GetPixelRed(image,p)*GetPixelRed(image,p);
1682       channel_statistics[RedChannel].sum_cubed+=(double)
1683         GetPixelRed(image,p)*GetPixelRed(image,p)*GetPixelRed(image,p);
1684       channel_statistics[RedChannel].sum_fourth_power+=(double)
1685         GetPixelRed(image,p)*GetPixelRed(image,p)*GetPixelRed(image,p)*
1686         GetPixelRed(image,p);
1687       if ((double) GetPixelGreen(image,p) < channel_statistics[GreenChannel].minima)
1688         channel_statistics[GreenChannel].minima=(double)
1689           GetPixelGreen(image,p);
1690       if ((double) GetPixelGreen(image,p) > channel_statistics[GreenChannel].maxima)
1691         channel_statistics[GreenChannel].maxima=(double) GetPixelGreen(image,p);
1692       channel_statistics[GreenChannel].sum+=GetPixelGreen(image,p);
1693       channel_statistics[GreenChannel].sum_squared+=(double)
1694         GetPixelGreen(image,p)*GetPixelGreen(image,p);
1695       channel_statistics[GreenChannel].sum_cubed+=(double)
1696         GetPixelGreen(image,p)*GetPixelGreen(image,p)*GetPixelGreen(image,p);
1697       channel_statistics[GreenChannel].sum_fourth_power+=(double)
1698         GetPixelGreen(image,p)*GetPixelGreen(image,p)*GetPixelGreen(image,p)*
1699         GetPixelGreen(image,p);
1700       if ((double) GetPixelBlue(image,p) < channel_statistics[BlueChannel].minima)
1701         channel_statistics[BlueChannel].minima=(double) GetPixelBlue(image,p);
1702       if ((double) GetPixelBlue(image,p) > channel_statistics[BlueChannel].maxima)
1703         channel_statistics[BlueChannel].maxima=(double) GetPixelBlue(image,p);
1704       channel_statistics[BlueChannel].sum+=GetPixelBlue(image,p);
1705       channel_statistics[BlueChannel].sum_squared+=(double)
1706         GetPixelBlue(image,p)*GetPixelBlue(image,p);
1707       channel_statistics[BlueChannel].sum_cubed+=(double)
1708         GetPixelBlue(image,p)*GetPixelBlue(image,p)*GetPixelBlue(image,p);
1709       channel_statistics[BlueChannel].sum_fourth_power+=(double)
1710         GetPixelBlue(image,p)*GetPixelBlue(image,p)*GetPixelBlue(image,p)*
1711         GetPixelBlue(image,p);
1712       if (image->colorspace == CMYKColorspace)
1713         {
1714           if ((double) GetPixelBlack(image,p) < channel_statistics[BlackChannel].minima)
1715             channel_statistics[BlackChannel].minima=(double)
1716               GetPixelBlack(image,p);
1717           if ((double) GetPixelBlack(image,p) > channel_statistics[BlackChannel].maxima)
1718             channel_statistics[BlackChannel].maxima=(double)
1719               GetPixelBlack(image,p);
1720           channel_statistics[BlackChannel].sum+=GetPixelBlack(image,p);
1721           channel_statistics[BlackChannel].sum_squared+=(double)
1722             GetPixelBlack(image,p)*GetPixelBlack(image,p);
1723           channel_statistics[BlackChannel].sum_cubed+=(double)
1724             GetPixelBlack(image,p)*GetPixelBlack(image,p)*
1725             GetPixelBlack(image,p);
1726           channel_statistics[BlackChannel].sum_fourth_power+=(double)
1727             GetPixelBlack(image,p)*GetPixelBlack(image,p)*
1728             GetPixelBlack(image,p)*GetPixelBlack(image,p);
1729         }
1730       if (image->matte != MagickFalse)
1731         {
1732           if ((double) GetPixelAlpha(image,p) < channel_statistics[AlphaChannel].minima)
1733             channel_statistics[AlphaChannel].minima=(double)
1734               GetPixelAlpha(image,p);
1735           if ((double) GetPixelAlpha(image,p) > channel_statistics[AlphaChannel].maxima)
1736             channel_statistics[AlphaChannel].maxima=(double)
1737               GetPixelAlpha(image,p);
1738           channel_statistics[AlphaChannel].sum+=GetPixelAlpha(image,p);
1739           channel_statistics[AlphaChannel].sum_squared+=(double)
1740             GetPixelAlpha(image,p)*GetPixelAlpha(image,p);
1741           channel_statistics[AlphaChannel].sum_cubed+=(double)
1742             GetPixelAlpha(image,p)*GetPixelAlpha(image,p)*
1743             GetPixelAlpha(image,p);
1744           channel_statistics[AlphaChannel].sum_fourth_power+=(double)
1745             GetPixelAlpha(image,p)*GetPixelAlpha(image,p)*
1746             GetPixelAlpha(image,p)*GetPixelAlpha(image,p);
1747         }
1748       x++;
1749       p+=GetPixelChannels(image);
1750     }
1751   }
1752   area=(double) image->columns*image->rows;
1753   for (i=0; i < (ssize_t) CompositeChannels; i++)
1754   {
1755     channel_statistics[i].sum/=area;
1756     channel_statistics[i].sum_squared/=area;
1757     channel_statistics[i].sum_cubed/=area;
1758     channel_statistics[i].sum_fourth_power/=area;
1759     channel_statistics[i].mean=channel_statistics[i].sum;
1760     channel_statistics[i].variance=channel_statistics[i].sum_squared;
1761     channel_statistics[i].standard_deviation=sqrt(
1762       channel_statistics[i].variance-(channel_statistics[i].mean*
1763       channel_statistics[i].mean));
1764   }
1765   for (i=0; i < (ssize_t) CompositeChannels; i++)
1766   {
1767     channel_statistics[CompositeChannels].depth=(size_t) MagickMax((double)
1768       channel_statistics[CompositeChannels].depth,(double)
1769       channel_statistics[i].depth);
1770     channel_statistics[CompositeChannels].minima=MagickMin(
1771       channel_statistics[CompositeChannels].minima,
1772       channel_statistics[i].minima);
1773     channel_statistics[CompositeChannels].maxima=MagickMax(
1774       channel_statistics[CompositeChannels].maxima,
1775       channel_statistics[i].maxima);
1776     channel_statistics[CompositeChannels].sum+=channel_statistics[i].sum;
1777     channel_statistics[CompositeChannels].sum_squared+=
1778       channel_statistics[i].sum_squared;
1779     channel_statistics[CompositeChannels].sum_cubed+=
1780       channel_statistics[i].sum_cubed;
1781     channel_statistics[CompositeChannels].sum_fourth_power+=
1782       channel_statistics[i].sum_fourth_power;
1783     channel_statistics[CompositeChannels].mean+=channel_statistics[i].mean;
1784     channel_statistics[CompositeChannels].variance+=
1785       channel_statistics[i].variance-channel_statistics[i].mean*
1786       channel_statistics[i].mean;
1787     channel_statistics[CompositeChannels].standard_deviation+=
1788       channel_statistics[i].variance-channel_statistics[i].mean*
1789       channel_statistics[i].mean;
1790   }
1791   channels=3;
1792   if (image->matte != MagickFalse)
1793     channels++;
1794   if (image->colorspace == CMYKColorspace)
1795     channels++;
1796   channel_statistics[CompositeChannels].sum/=channels;
1797   channel_statistics[CompositeChannels].sum_squared/=channels;
1798   channel_statistics[CompositeChannels].sum_cubed/=channels;
1799   channel_statistics[CompositeChannels].sum_fourth_power/=channels;
1800   channel_statistics[CompositeChannels].mean/=channels;
1801   channel_statistics[CompositeChannels].variance/=channels;
1802   channel_statistics[CompositeChannels].standard_deviation=
1803     sqrt(channel_statistics[CompositeChannels].standard_deviation/channels);
1804   channel_statistics[CompositeChannels].kurtosis/=channels;
1805   channel_statistics[CompositeChannels].skewness/=channels;
1806   for (i=0; i <= (ssize_t) CompositeChannels; i++)
1807   {
1808     if (channel_statistics[i].standard_deviation == 0.0)
1809       continue;
1810     channel_statistics[i].skewness=(channel_statistics[i].sum_cubed-
1811       3.0*channel_statistics[i].mean*channel_statistics[i].sum_squared+
1812       2.0*channel_statistics[i].mean*channel_statistics[i].mean*
1813       channel_statistics[i].mean)/(channel_statistics[i].standard_deviation*
1814       channel_statistics[i].standard_deviation*
1815       channel_statistics[i].standard_deviation);
1816     channel_statistics[i].kurtosis=(channel_statistics[i].sum_fourth_power-
1817       4.0*channel_statistics[i].mean*channel_statistics[i].sum_cubed+
1818       6.0*channel_statistics[i].mean*channel_statistics[i].mean*
1819       channel_statistics[i].sum_squared-3.0*channel_statistics[i].mean*
1820       channel_statistics[i].mean*1.0*channel_statistics[i].mean*
1821       channel_statistics[i].mean)/(channel_statistics[i].standard_deviation*
1822       channel_statistics[i].standard_deviation*
1823       channel_statistics[i].standard_deviation*
1824       channel_statistics[i].standard_deviation)-3.0;
1825   }
1826   return(channel_statistics);
1827 }