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-2012 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/gem-private.h"
  67 #include "MagickCore/geometry.h"
  68 #include "MagickCore/list.h"
  69 #include "MagickCore/image-private.h"
  70 #include "MagickCore/magic.h"
  71 #include "MagickCore/magick.h"
  72 #include "MagickCore/memory_.h"
  73 #include "MagickCore/module.h"
  74 #include "MagickCore/monitor.h"
  75 #include "MagickCore/monitor-private.h"
  76 #include "MagickCore/option.h"
  77 #include "MagickCore/paint.h"
  78 #include "MagickCore/pixel-accessor.h"
  79 #include "MagickCore/profile.h"
  80 #include "MagickCore/quantize.h"
  81 #include "MagickCore/quantum-private.h"
  82 #include "MagickCore/random_.h"
  83 #include "MagickCore/random-private.h"
  84 #include "MagickCore/resource_.h"
  85 #include "MagickCore/segment.h"
  86 #include "MagickCore/semaphore.h"
  87 #include "MagickCore/signature-private.h"
  88 #include "MagickCore/statistic.h"
  89 #include "MagickCore/string_.h"
  90 #include "MagickCore/thread-private.h"
  91 #include "MagickCore/timer.h"
  92 #include "MagickCore/utility.h"
  93 #include "MagickCore/version.h"
  94 \f
  95 /*
  96 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  97 %                                                                             %
  98 %                                                                             %
  99 %                                                                             %
 100 %     E v a l u a t e I m a g e                                               %
 101 %                                                                             %
 102 %                                                                             %
 103 %                                                                             %
 104 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 105 %
 106 %  EvaluateImage() applies a value to the image with an arithmetic, relational,
 107 %  or logical operator to an image. Use these operations to lighten or darken
 108 %  an image, to increase or decrease contrast in an image, or to produce the
 109 %  "negative" of an image.
 110 %
 111 %  The format of the EvaluateImage method is:
 112 %
 113 %      MagickBooleanType EvaluateImage(Image *image,
 114 %        const MagickEvaluateOperator op,const double value,
 115 %        ExceptionInfo *exception)
 116 %      MagickBooleanType EvaluateImages(Image *images,
 117 %        const MagickEvaluateOperator op,const double value,
 118 %        ExceptionInfo *exception)
 119 %
 120 %  A description of each parameter follows:
 121 %
 122 %    o image: the image.
 123 %
 124 %    o op: A channel op.
 125 %
 126 %    o value: A value value.
 127 %
 128 %    o exception: return any errors or warnings in this structure.
 129 %
 130 */
 131
 132 typedef struct _PixelChannels
 133 {
 134   double
 135     channel[CompositePixelChannel];
 136 } PixelChannels;
 137
 138 static PixelChannels **DestroyPixelThreadSet(PixelChannels **pixels)
 139 {
 140   register ssize_t
 141     i;
 142
 143   assert(pixels != (PixelChannels **) NULL);
 144   for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
 145     if (pixels[i] != (PixelChannels *) NULL)
 146       pixels[i]=(PixelChannels *) RelinquishMagickMemory(pixels[i]);
 147   pixels=(PixelChannels **) RelinquishMagickMemory(pixels);
 148   return(pixels);
 149 }
 150
 151 static PixelChannels **AcquirePixelThreadSet(const Image *image,
 152   const size_t number_images)
 153 {
 154   register ssize_t
 155     i;
 156
 157   PixelChannels
 158     **pixels;
 159
 160   size_t
 161     length,
 162     number_threads;
 163
 164   number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
 165   pixels=(PixelChannels **) AcquireQuantumMemory(number_threads,
 166     sizeof(*pixels));
 167   if (pixels == (PixelChannels **) NULL)
 168     return((PixelChannels **) NULL);
 169   (void) ResetMagickMemory(pixels,0,number_threads*sizeof(*pixels));
 170   for (i=0; i < (ssize_t) number_threads; i++)
 171   {
 172     register ssize_t
 173       j;
 174
 175     length=image->columns;
 176     if (length < number_images)
 177       length=number_images;
 178     pixels[i]=(PixelChannels *) AcquireQuantumMemory(length,sizeof(**pixels));
 179     if (pixels[i] == (PixelChannels *) NULL)
 180       return(DestroyPixelThreadSet(pixels));
 181     for (j=0; j < (ssize_t) length; j++)
 182     {
 183       register ssize_t
 184         k;
 185
 186       for (k=0; k < MaxPixelChannels; k++)
 187         pixels[i][j].channel[k]=0.0;
 188     }
 189   }
 190   return(pixels);
 191 }
 192
 193 static inline double EvaluateMax(const double x,const double y)
 194 {
 195   if (x > y)
 196     return(x);
 197   return(y);
 198 }
 199
 200 #if defined(__cplusplus) || defined(c_plusplus)
 201 extern "C" {
 202 #endif
 203
 204 static int IntensityCompare(const void *x,const void *y)
 205 {
 206   const PixelChannels
 207     *color_1,
 208     *color_2;
 209
 210   double
 211     distance;
 212
 213   register ssize_t
 214     i;
 215
 216   color_1=(const PixelChannels *) x;
 217   color_2=(const PixelChannels *) y;
 218   distance=0.0;
 219   for (i=0; i < MaxPixelChannels; i++)
 220     distance+=color_1->channel[i]-(double) color_2->channel[i];
 221   return(distance < 0 ? -1 : distance > 0 ? 1 : 0);
 222 }
 223
 224 #if defined(__cplusplus) || defined(c_plusplus)
 225 }
 226 #endif
 227
 228 static inline double MagickMin(const double x,const double y)
 229 {
 230   if (x < y)
 231     return(x);
 232   return(y);
 233 }
 234
 235 static double ApplyEvaluateOperator(RandomInfo *random_info,
 236   Quantum pixel,const MagickEvaluateOperator op,const double value)
 237 {
 238   double
 239     result;
 240
 241   result=0.0;
 242   switch (op)
 243   {
 244     case UndefinedEvaluateOperator:
 245       break;
 246     case AbsEvaluateOperator:
 247     {
 248       result=(double) fabs((double) (pixel+value));
 249       break;
 250     }
 251     case AddEvaluateOperator:
 252     {
 253       result=(double) (pixel+value);
 254       break;
 255     }
 256     case AddModulusEvaluateOperator:
 257     {
 258       /*
 259         This returns a 'floored modulus' of the addition which is a positive
 260         result.  It differs from % or fmod() that returns a 'truncated modulus'
 261         result, where floor() is replaced by trunc() and could return a
 262         negative result (which is clipped).
 263       */
 264       result=pixel+value;
 265       result-=(QuantumRange+1.0)*floor((double) result/(QuantumRange+1.0));
 266       break;
 267     }
 268     case AndEvaluateOperator:
 269     {
 270       result=(double) ((size_t) pixel & (size_t) (value+0.5));
 271       break;
 272     }
 273     case CosineEvaluateOperator:
 274     {
 275       result=(double) (QuantumRange*(0.5*cos((double) (2.0*MagickPI*
 276         QuantumScale*pixel*value))+0.5));
 277       break;
 278     }
 279     case DivideEvaluateOperator:
 280     {
 281       result=pixel/(value == 0.0 ? 1.0 : value);
 282       break;
 283     }
 284     case ExponentialEvaluateOperator:
 285     {
 286       result=(double) (QuantumRange*exp((double) (value*QuantumScale*
 287         pixel)));
 288       break;
 289     }
 290     case GaussianNoiseEvaluateOperator:
 291     {
 292       result=(double) GenerateDifferentialNoise(random_info,pixel,
 293         GaussianNoise,value);
 294       break;
 295     }
 296     case ImpulseNoiseEvaluateOperator:
 297     {
 298       result=(double) GenerateDifferentialNoise(random_info,pixel,
 299         ImpulseNoise,value);
 300       break;
 301     }
 302     case LaplacianNoiseEvaluateOperator:
 303     {
 304       result=(double) GenerateDifferentialNoise(random_info,pixel,
 305         LaplacianNoise,value);
 306       break;
 307     }
 308     case LeftShiftEvaluateOperator:
 309     {
 310       result=(double) ((size_t) pixel << (size_t) (value+0.5));
 311       break;
 312     }
 313     case LogEvaluateOperator:
 314     {
 315       if ((QuantumScale*pixel) >= MagickEpsilon)
 316         result=(double) (QuantumRange*log((double) (QuantumScale*value*
 317           pixel+1.0))/log((double) (value+1.0)));
 318       break;
 319     }
 320     case MaxEvaluateOperator:
 321     {
 322       result=(double) EvaluateMax((double) pixel,value);
 323       break;
 324     }
 325     case MeanEvaluateOperator:
 326     {
 327       result=(double) (pixel+value);
 328       break;
 329     }
 330     case MedianEvaluateOperator:
 331     {
 332       result=(double) (pixel+value);
 333       break;
 334     }
 335     case MinEvaluateOperator:
 336     {
 337       result=(double) MagickMin((double) pixel,value);
 338       break;
 339     }
 340     case MultiplicativeNoiseEvaluateOperator:
 341     {
 342       result=(double) GenerateDifferentialNoise(random_info,pixel,
 343         MultiplicativeGaussianNoise,value);
 344       break;
 345     }
 346     case MultiplyEvaluateOperator:
 347     {
 348       result=(double) (value*pixel);
 349       break;
 350     }
 351     case OrEvaluateOperator:
 352     {
 353       result=(double) ((size_t) pixel | (size_t) (value+0.5));
 354       break;
 355     }
 356     case PoissonNoiseEvaluateOperator:
 357     {
 358       result=(double) GenerateDifferentialNoise(random_info,pixel,
 359         PoissonNoise,value);
 360       break;
 361     }
 362     case PowEvaluateOperator:
 363     {
 364       result=(double) (QuantumRange*pow((double) (QuantumScale*pixel),
 365         (double) value));
 366       break;
 367     }
 368     case RightShiftEvaluateOperator:
 369     {
 370       result=(double) ((size_t) pixel >> (size_t) (value+0.5));
 371       break;
 372     }
 373     case SetEvaluateOperator:
 374     {
 375       result=value;
 376       break;
 377     }
 378     case SineEvaluateOperator:
 379     {
 380       result=(double) (QuantumRange*(0.5*sin((double) (2.0*MagickPI*
 381         QuantumScale*pixel*value))+0.5));
 382       break;
 383     }
 384     case SubtractEvaluateOperator:
 385     {
 386       result=(double) (pixel-value);
 387       break;
 388     }
 389     case SumEvaluateOperator:
 390     {
 391       result=(double) (pixel+value);
 392       break;
 393     }
 394     case ThresholdEvaluateOperator:
 395     {
 396       result=(double) (((double) pixel <= value) ? 0 :
 397         QuantumRange);
 398       break;
 399     }
 400     case ThresholdBlackEvaluateOperator:
 401     {
 402       result=(double) (((double) pixel <= value) ? 0 : pixel);
 403       break;
 404     }
 405     case ThresholdWhiteEvaluateOperator:
 406     {
 407       result=(double) (((double) pixel > value) ? QuantumRange :
 408         pixel);
 409       break;
 410     }
 411     case UniformNoiseEvaluateOperator:
 412     {
 413       result=(double) GenerateDifferentialNoise(random_info,pixel,
 414         UniformNoise,value);
 415       break;
 416     }
 417     case XorEvaluateOperator:
 418     {
 419       result=(double) ((size_t) pixel ^ (size_t) (value+0.5));
 420       break;
 421     }
 422   }
 423   return(result);
 424 }
 425
 426 MagickExport Image *EvaluateImages(const Image *images,
 427   const MagickEvaluateOperator op,ExceptionInfo *exception)
 428 {
 429 #define EvaluateImageTag  "Evaluate/Image"
 430
 431   CacheView
 432     *evaluate_view;
 433
 434   const Image
 435     *next;
 436
 437   Image
 438     *image;
 439
 440   MagickBooleanType
 441     status;
 442
 443   MagickOffsetType
 444     progress;
 445
 446   PixelChannels
 447     **restrict evaluate_pixels;
 448
 449   RandomInfo
 450     **restrict random_info;
 451
 452   size_t
 453     number_images;
 454
 455   ssize_t
 456     y;
 457
 458   unsigned long
 459     key;
 460
 461   /*
 462     Ensure the image are the same size.
 463   */
 464   assert(images != (Image *) NULL);
 465   assert(images->signature == MagickSignature);
 466   if (images->debug != MagickFalse)
 467     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
 468   assert(exception != (ExceptionInfo *) NULL);
 469   assert(exception->signature == MagickSignature);
 470   for (next=images; next != (Image *) NULL; next=GetNextImageInList(next))
 471     if ((next->columns != images->columns) || (next->rows != images->rows))
 472       {
 473         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
 474           "ImageWidthsOrHeightsDiffer","'%s'",images->filename);
 475         return((Image *) NULL);
 476       }
 477   /*
 478     Initialize evaluate next attributes.
 479   */
 480   image=CloneImage(images,images->columns,images->rows,MagickTrue,
 481     exception);
 482   if (image == (Image *) NULL)
 483     return((Image *) NULL);
 484   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
 485     {
 486       image=DestroyImage(image);
 487       return((Image *) NULL);
 488     }
 489   number_images=GetImageListLength(images);
 490   evaluate_pixels=AcquirePixelThreadSet(images,number_images);
 491   if (evaluate_pixels == (PixelChannels **) NULL)
 492     {
 493       image=DestroyImage(image);
 494       (void) ThrowMagickException(exception,GetMagickModule(),
 495         ResourceLimitError,"MemoryAllocationFailed","'%s'",images->filename);
 496       return((Image *) NULL);
 497     }
 498   /*
 499     Evaluate image pixels.
 500   */
 501   status=MagickTrue;
 502   progress=0;
 503   random_info=AcquireRandomInfoThreadSet();
 504   key=GetRandomSecretKey(random_info[0]);
 505   evaluate_view=AcquireAuthenticCacheView(image,exception);
 506   if (op == MedianEvaluateOperator)
 507     {
 508 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 509       #pragma omp parallel for schedule(static) shared(progress,status) \
 510         dynamic_number_threads(image,image->columns,image->rows,key == ~0UL)
 511 #endif
 512       for (y=0; y < (ssize_t) image->rows; y++)
 513       {
 514         CacheView
 515           *image_view;
 516
 517         const Image
 518           *next;
 519
 520         const int
 521           id = GetOpenMPThreadId();
 522
 523         register PixelChannels
 524           *evaluate_pixel;
 525
 526         register Quantum
 527           *restrict q;
 528
 529         register ssize_t
 530           x;
 531
 532         if (status == MagickFalse)
 533           continue;
 534         q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,
 535           image->columns,1,exception);
 536         if (q == (Quantum *) NULL)
 537           {
 538             status=MagickFalse;
 539             continue;
 540           }
 541         evaluate_pixel=evaluate_pixels[id];
 542         for (x=0; x < (ssize_t) image->columns; x++)
 543         {
 544           register ssize_t
 545             j,
 546             k;
 547
 548           for (j=0; j < (ssize_t) number_images; j++)
 549             for (k=0; k < MaxPixelChannels; k++)
 550               evaluate_pixel[j].channel[k]=0.0;
 551           next=images;
 552           for (j=0; j < (ssize_t) number_images; j++)
 553           {
 554             register const Quantum
 555               *p;
 556
 557             register ssize_t
 558               i;
 559
 560             image_view=AcquireVirtualCacheView(next,exception);
 561             p=GetCacheViewVirtualPixels(image_view,x,y,1,1,exception);
 562             if (p == (const Quantum *) NULL)
 563               {
 564                 image_view=DestroyCacheView(image_view);
 565                 break;
 566               }
 567             for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 568             {
 569               PixelChannel
 570                 channel;
 571
 572               PixelTrait
 573                 evaluate_traits,
 574                 traits;
 575
 576               channel=GetPixelChannelMapChannel(image,i);
 577               evaluate_traits=GetPixelChannelMapTraits(image,channel);
 578               traits=GetPixelChannelMapTraits(next,channel);
 579               if ((traits == UndefinedPixelTrait) ||
 580                   (evaluate_traits == UndefinedPixelTrait))
 581                 continue;
 582               if ((evaluate_traits & UpdatePixelTrait) == 0)
 583                 continue;
 584               evaluate_pixel[j].channel[i]=ApplyEvaluateOperator(
 585                 random_info[id],GetPixelChannel(image,channel,p),op,
 586                 evaluate_pixel[j].channel[i]);
 587             }
 588             image_view=DestroyCacheView(image_view);
 589             next=GetNextImageInList(next);
 590           }
 591           qsort((void *) evaluate_pixel,number_images,sizeof(*evaluate_pixel),
 592             IntensityCompare);
 593           for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
 594             q[k]=ClampToQuantum(evaluate_pixel[j/2].channel[k]);
 595           q+=GetPixelChannels(image);
 596         }
 597         if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
 598           status=MagickFalse;
 599         if (images->progress_monitor != (MagickProgressMonitor) NULL)
 600           {
 601             MagickBooleanType
 602               proceed;
 603
 604 #if   defined(MAGICKCORE_OPENMP_SUPPORT)
 605             #pragma omp critical (MagickCore_EvaluateImages)
 606 #endif
 607             proceed=SetImageProgress(images,EvaluateImageTag,progress++,
 608               image->rows);
 609             if (proceed == MagickFalse)
 610               status=MagickFalse;
 611           }
 612       }
 613     }
 614   else
 615     {
 616 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 617       #pragma omp parallel for schedule(static) shared(progress,status) \
 618         dynamic_number_threads(image,image->columns,image->rows,key == ~0UL)
 619 #endif
 620       for (y=0; y < (ssize_t) image->rows; y++)
 621       {
 622         CacheView
 623           *image_view;
 624
 625         const Image
 626           *next;
 627
 628         const int
 629           id = GetOpenMPThreadId();
 630
 631         register ssize_t
 632           i,
 633           x;
 634
 635         register PixelChannels
 636           *evaluate_pixel;
 637
 638         register Quantum
 639           *restrict q;
 640
 641         ssize_t
 642           j;
 643
 644         if (status == MagickFalse)
 645           continue;
 646         q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,
 647           image->columns,1,exception);
 648         if (q == (Quantum *) NULL)
 649           {
 650             status=MagickFalse;
 651             continue;
 652           }
 653         evaluate_pixel=evaluate_pixels[id];
 654         for (j=0; j < (ssize_t) image->columns; j++)
 655           for (i=0; i < MaxPixelChannels; i++)
 656             evaluate_pixel[j].channel[i]=0.0;
 657         next=images;
 658         for (j=0; j < (ssize_t) number_images; j++)
 659         {
 660           register const Quantum
 661             *p;
 662
 663           image_view=AcquireVirtualCacheView(next,exception);
 664           p=GetCacheViewVirtualPixels(image_view,0,y,next->columns,1,exception);
 665           if (p == (const Quantum *) NULL)
 666             {
 667               image_view=DestroyCacheView(image_view);
 668               break;
 669             }
 670           for (x=0; x < (ssize_t) next->columns; x++)
 671           {
 672             register ssize_t
 673               i;
 674
 675             if (GetPixelMask(next,p) != 0)
 676               {
 677                 p+=GetPixelChannels(next);
 678                 continue;
 679               }
 680             for (i=0; i < (ssize_t) GetPixelChannels(next); i++)
 681             {
 682               PixelChannel
 683                 channel;
 684
 685               PixelTrait
 686                 evaluate_traits,
 687                 traits;
 688
 689               channel=GetPixelChannelMapChannel(image,i);
 690               traits=GetPixelChannelMapTraits(next,channel);
 691               evaluate_traits=GetPixelChannelMapTraits(image,channel);
 692               if ((traits == UndefinedPixelTrait) ||
 693                   (evaluate_traits == UndefinedPixelTrait))
 694                 continue;
 695               if ((traits & UpdatePixelTrait) == 0)
 696                 continue;
 697               evaluate_pixel[x].channel[i]=ApplyEvaluateOperator(
 698                 random_info[id],GetPixelChannel(image,channel,p),j ==
 699                 0 ? AddEvaluateOperator : op,evaluate_pixel[x].channel[i]);
 700             }
 701             p+=GetPixelChannels(next);
 702           }
 703           image_view=DestroyCacheView(image_view);
 704           next=GetNextImageInList(next);
 705         }
 706         for (x=0; x < (ssize_t) image->columns; x++)
 707         {
 708           register ssize_t
 709              i;
 710
 711           switch (op)
 712           {
 713             case MeanEvaluateOperator:
 714             {
 715               for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 716                 evaluate_pixel[x].channel[i]/=(double) number_images;
 717               break;
 718             }
 719             case MultiplyEvaluateOperator:
 720             {
 721               for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 722               {
 723                 register ssize_t
 724                   j;
 725
 726                 for (j=0; j < (ssize_t) (number_images-1); j++)
 727                   evaluate_pixel[x].channel[i]*=QuantumScale;
 728               }
 729               break;
 730             }
 731             default:
 732               break;
 733           }
 734         }
 735         for (x=0; x < (ssize_t) image->columns; x++)
 736         {
 737           register ssize_t
 738             i;
 739
 740           if (GetPixelMask(image,q) != 0)
 741             {
 742               q+=GetPixelChannels(image);
 743               continue;
 744             }
 745           for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 746           {
 747             PixelChannel
 748               channel;
 749
 750             PixelTrait
 751               traits;
 752
 753             channel=GetPixelChannelMapChannel(image,i);
 754             traits=GetPixelChannelMapTraits(image,channel);
 755             if (traits == UndefinedPixelTrait)
 756               continue;
 757             if ((traits & UpdatePixelTrait) == 0)
 758               continue;
 759             q[i]=ClampToQuantum(evaluate_pixel[x].channel[i]);
 760           }
 761           q+=GetPixelChannels(image);
 762         }
 763         if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
 764           status=MagickFalse;
 765         if (images->progress_monitor != (MagickProgressMonitor) NULL)
 766           {
 767             MagickBooleanType
 768               proceed;
 769
 770 #if   defined(MAGICKCORE_OPENMP_SUPPORT)
 771             #pragma omp critical (MagickCore_EvaluateImages)
 772 #endif
 773             proceed=SetImageProgress(images,EvaluateImageTag,progress++,
 774               image->rows);
 775             if (proceed == MagickFalse)
 776               status=MagickFalse;
 777           }
 778       }
 779     }
 780   evaluate_view=DestroyCacheView(evaluate_view);
 781   evaluate_pixels=DestroyPixelThreadSet(evaluate_pixels);
 782   random_info=DestroyRandomInfoThreadSet(random_info);
 783   if (status == MagickFalse)
 784     image=DestroyImage(image);
 785   return(image);
 786 }
 787
 788 MagickExport MagickBooleanType EvaluateImage(Image *image,
 789   const MagickEvaluateOperator op,const double value,ExceptionInfo *exception)
 790 {
 791   CacheView
 792     *image_view;
 793
 794   MagickBooleanType
 795     status;
 796
 797   MagickOffsetType
 798     progress;
 799
 800   RandomInfo
 801     **restrict random_info;
 802
 803   ssize_t
 804     y;
 805
 806   unsigned long
 807     key;
 808
 809   assert(image != (Image *) NULL);
 810   assert(image->signature == MagickSignature);
 811   if (image->debug != MagickFalse)
 812     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 813   assert(exception != (ExceptionInfo *) NULL);
 814   assert(exception->signature == MagickSignature);
 815   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
 816     return(MagickFalse);
 817   status=MagickTrue;
 818   progress=0;
 819   random_info=AcquireRandomInfoThreadSet();
 820   key=GetRandomSecretKey(random_info[0]);
 821   image_view=AcquireAuthenticCacheView(image,exception);
 822 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 823   #pragma omp parallel for schedule(static,4) shared(progress,status) \
 824     dynamic_number_threads(image,image->columns,image->rows,key == ~0UL)
 825 #endif
 826   for (y=0; y < (ssize_t) image->rows; y++)
 827   {
 828     const int
 829       id = GetOpenMPThreadId();
 830
 831     register Quantum
 832       *restrict q;
 833
 834     register ssize_t
 835       x;
 836
 837     if (status == MagickFalse)
 838       continue;
 839     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
 840     if (q == (Quantum *) NULL)
 841       {
 842         status=MagickFalse;
 843         continue;
 844       }
 845     for (x=0; x < (ssize_t) image->columns; x++)
 846     {
 847       register ssize_t
 848         i;
 849
 850       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 851       {
 852         PixelChannel
 853           channel;
 854
 855         PixelTrait
 856           traits;
 857
 858         channel=GetPixelChannelMapChannel(image,i);
 859         traits=GetPixelChannelMapTraits(image,channel);
 860         if (traits == UndefinedPixelTrait)
 861           continue;
 862         if (((traits & CopyPixelTrait) != 0) ||
 863             (GetPixelMask(image,q) != 0))
 864           continue;
 865         q[i]=ClampToQuantum(ApplyEvaluateOperator(random_info[id],q[i],op,
 866           value));
 867       }
 868       q+=GetPixelChannels(image);
 869     }
 870     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
 871       status=MagickFalse;
 872     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 873       {
 874         MagickBooleanType
 875           proceed;
 876
 877 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 878         #pragma omp critical (MagickCore_EvaluateImage)
 879 #endif
 880         proceed=SetImageProgress(image,EvaluateImageTag,progress++,image->rows);
 881         if (proceed == MagickFalse)
 882           status=MagickFalse;
 883       }
 884   }
 885   image_view=DestroyCacheView(image_view);
 886   random_info=DestroyRandomInfoThreadSet(random_info);
 887   return(status);
 888 }
 889 \f
 890 /*
 891 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 892 %                                                                             %
 893 %                                                                             %
 894 %                                                                             %
 895 %     F u n c t i o n I m a g e                                               %
 896 %                                                                             %
 897 %                                                                             %
 898 %                                                                             %
 899 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 900 %
 901 %  FunctionImage() applies a value to the image with an arithmetic, relational,
 902 %  or logical operator to an image. Use these operations to lighten or darken
 903 %  an image, to increase or decrease contrast in an image, or to produce the
 904 %  "negative" of an image.
 905 %
 906 %  The format of the FunctionImage method is:
 907 %
 908 %      MagickBooleanType FunctionImage(Image *image,
 909 %        const MagickFunction function,const ssize_t number_parameters,
 910 %        const double *parameters,ExceptionInfo *exception)
 911 %
 912 %  A description of each parameter follows:
 913 %
 914 %    o image: the image.
 915 %
 916 %    o function: A channel function.
 917 %
 918 %    o parameters: one or more parameters.
 919 %
 920 %    o exception: return any errors or warnings in this structure.
 921 %
 922 */
 923
 924 static Quantum ApplyFunction(Quantum pixel,const MagickFunction function,
 925   const size_t number_parameters,const double *parameters,
 926   ExceptionInfo *exception)
 927 {
 928   double
 929     result;
 930
 931   register ssize_t
 932     i;
 933
 934   (void) exception;
 935   result=0.0;
 936   switch (function)
 937   {
 938     case PolynomialFunction:
 939     {
 940       /*
 941         Polynomial: polynomial constants, highest to lowest order (e.g. c0*x^3+
 942         c1*x^2+c2*x+c3).
 943       */
 944       result=0.0;
 945       for (i=0; i < (ssize_t) number_parameters; i++)
 946         result=result*QuantumScale*pixel+parameters[i];
 947       result*=QuantumRange;
 948       break;
 949     }
 950     case SinusoidFunction:
 951     {
 952       double
 953         amplitude,
 954         bias,
 955         frequency,
 956         phase;
 957
 958       /*
 959         Sinusoid: frequency, phase, amplitude, bias.
 960       */
 961       frequency=(number_parameters >= 1) ? parameters[0] : 1.0;
 962       phase=(number_parameters >= 2) ? parameters[1] : 0.0;
 963       amplitude=(number_parameters >= 3) ? parameters[2] : 0.5;
 964       bias=(number_parameters >= 4) ? parameters[3] : 0.5;
 965       result=(double) (QuantumRange*(amplitude*sin((double) (2.0*
 966         MagickPI*(frequency*QuantumScale*pixel+phase/360.0)))+bias));
 967       break;
 968     }
 969     case ArcsinFunction:
 970     {
 971       double
 972         bias,
 973         center,
 974         range,
 975         width;
 976
 977       /*
 978         Arcsin (peged at range limits for invalid results): width, center,
 979         range, and bias.
 980       */
 981       width=(number_parameters >= 1) ? parameters[0] : 1.0;
 982       center=(number_parameters >= 2) ? parameters[1] : 0.5;
 983       range=(number_parameters >= 3) ? parameters[2] : 1.0;
 984       bias=(number_parameters >= 4) ? parameters[3] : 0.5;
 985       result=2.0/width*(QuantumScale*pixel-center);
 986       if ( result <= -1.0 )
 987         result=bias-range/2.0;
 988       else
 989         if (result >= 1.0)
 990           result=bias+range/2.0;
 991         else
 992           result=(double) (range/MagickPI*asin((double) result)+bias);
 993       result*=QuantumRange;
 994       break;
 995     }
 996     case ArctanFunction:
 997     {
 998       double
 999         center,
1000         bias,
1001         range,
1002         slope;
1003
1004       /*
1005         Arctan: slope, center, range, and bias.
1006       */
1007       slope=(number_parameters >= 1) ? parameters[0] : 1.0;
1008       center=(number_parameters >= 2) ? parameters[1] : 0.5;
1009       range=(number_parameters >= 3) ? parameters[2] : 1.0;
1010       bias=(number_parameters >= 4) ? parameters[3] : 0.5;
1011       result=(double) (MagickPI*slope*(QuantumScale*pixel-center));
1012       result=(double) (QuantumRange*(range/MagickPI*atan((double)
1013         result)+bias));
1014       break;
1015     }
1016     case UndefinedFunction:
1017       break;
1018   }
1019   return(ClampToQuantum(result));
1020 }
1021
1022 MagickExport MagickBooleanType FunctionImage(Image *image,
1023   const MagickFunction function,const size_t number_parameters,
1024   const double *parameters,ExceptionInfo *exception)
1025 {
1026 #define FunctionImageTag  "Function/Image "
1027
1028   CacheView
1029     *image_view;
1030
1031   MagickBooleanType
1032     status;
1033
1034   MagickOffsetType
1035     progress;
1036
1037   ssize_t
1038     y;
1039
1040   assert(image != (Image *) NULL);
1041   assert(image->signature == MagickSignature);
1042   if (image->debug != MagickFalse)
1043     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1044   assert(exception != (ExceptionInfo *) NULL);
1045   assert(exception->signature == MagickSignature);
1046   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
1047     return(MagickFalse);
1048   status=MagickTrue;
1049   progress=0;
1050   image_view=AcquireAuthenticCacheView(image,exception);
1051 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1052   #pragma omp parallel for schedule(static,4) shared(progress,status) \
1053     dynamic_number_threads(image,image->columns,image->rows,1)
1054 #endif
1055   for (y=0; y < (ssize_t) image->rows; y++)
1056   {
1057     register Quantum
1058       *restrict q;
1059
1060     register ssize_t
1061       x;
1062
1063     if (status == MagickFalse)
1064       continue;
1065     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
1066     if (q == (Quantum *) NULL)
1067       {
1068         status=MagickFalse;
1069         continue;
1070       }
1071     for (x=0; x < (ssize_t) image->columns; x++)
1072     {
1073       register ssize_t
1074         i;
1075
1076       if (GetPixelMask(image,q) != 0)
1077         {
1078           q+=GetPixelChannels(image);
1079           continue;
1080         }
1081       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1082       {
1083         PixelChannel
1084           channel;
1085
1086         PixelTrait
1087           traits;
1088
1089         channel=GetPixelChannelMapChannel(image,i);
1090         traits=GetPixelChannelMapTraits(image,channel);
1091         if (traits == UndefinedPixelTrait)
1092           continue;
1093         if ((traits & UpdatePixelTrait) == 0)
1094           continue;
1095         q[i]=ApplyFunction(q[i],function,number_parameters,parameters,
1096           exception);
1097       }
1098       q+=GetPixelChannels(image);
1099     }
1100     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1101       status=MagickFalse;
1102     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1103       {
1104         MagickBooleanType
1105           proceed;
1106
1107 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1108         #pragma omp critical (MagickCore_FunctionImage)
1109 #endif
1110         proceed=SetImageProgress(image,FunctionImageTag,progress++,image->rows);
1111         if (proceed == MagickFalse)
1112           status=MagickFalse;
1113       }
1114   }
1115   image_view=DestroyCacheView(image_view);
1116   return(status);
1117 }
1118 \f
1119 /*
1120 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1121 %                                                                             %
1122 %                                                                             %
1123 %                                                                             %
1124 %   G e t I m a g e E x t r e m a                                             %
1125 %                                                                             %
1126 %                                                                             %
1127 %                                                                             %
1128 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1129 %
1130 %  GetImageExtrema() returns the extrema of one or more image channels.
1131 %
1132 %  The format of the GetImageExtrema method is:
1133 %
1134 %      MagickBooleanType GetImageExtrema(const Image *image,size_t *minima,
1135 %        size_t *maxima,ExceptionInfo *exception)
1136 %
1137 %  A description of each parameter follows:
1138 %
1139 %    o image: the image.
1140 %
1141 %    o minima: the minimum value in the channel.
1142 %
1143 %    o maxima: the maximum value in the channel.
1144 %
1145 %    o exception: return any errors or warnings in this structure.
1146 %
1147 */
1148 MagickExport MagickBooleanType GetImageExtrema(const Image *image,
1149   size_t *minima,size_t *maxima,ExceptionInfo *exception)
1150 {
1151   double
1152     max,
1153     min;
1154
1155   MagickBooleanType
1156     status;
1157
1158   assert(image != (Image *) NULL);
1159   assert(image->signature == MagickSignature);
1160   if (image->debug != MagickFalse)
1161     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1162   status=GetImageRange(image,&min,&max,exception);
1163   *minima=(size_t) ceil(min-0.5);
1164   *maxima=(size_t) floor(max+0.5);
1165   return(status);
1166 }
1167 \f
1168 /*
1169 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1170 %                                                                             %
1171 %                                                                             %
1172 %                                                                             %
1173 %   G e t I m a g e M e a n                                                   %
1174 %                                                                             %
1175 %                                                                             %
1176 %                                                                             %
1177 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1178 %
1179 %  GetImageMean() returns the mean and standard deviation of one or more image
1180 %  channels.
1181 %
1182 %  The format of the GetImageMean method is:
1183 %
1184 %      MagickBooleanType GetImageMean(const Image *image,double *mean,
1185 %        double *standard_deviation,ExceptionInfo *exception)
1186 %
1187 %  A description of each parameter follows:
1188 %
1189 %    o image: the image.
1190 %
1191 %    o mean: the average value in the channel.
1192 %
1193 %    o standard_deviation: the standard deviation of the channel.
1194 %
1195 %    o exception: return any errors or warnings in this structure.
1196 %
1197 */
1198 MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean,
1199   double *standard_deviation,ExceptionInfo *exception)
1200 {
1201   double
1202     area;
1203
1204   ChannelStatistics
1205     *channel_statistics;
1206
1207   register ssize_t
1208     i;
1209
1210   assert(image != (Image *) NULL);
1211   assert(image->signature == MagickSignature);
1212   if (image->debug != MagickFalse)
1213     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1214   channel_statistics=GetImageStatistics(image,exception);
1215   if (channel_statistics == (ChannelStatistics *) NULL)
1216     return(MagickFalse);
1217   area=0.0;
1218   channel_statistics[CompositePixelChannel].mean=0.0;
1219   channel_statistics[CompositePixelChannel].standard_deviation=0.0;
1220   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1221   {
1222     PixelChannel
1223       channel;
1224
1225     PixelTrait
1226       traits;
1227
1228     channel=GetPixelChannelMapChannel(image,i);
1229     traits=GetPixelChannelMapTraits(image,channel);
1230     if (traits == UndefinedPixelTrait)
1231       continue;
1232     if ((traits & UpdatePixelTrait) == 0)
1233       continue;
1234     channel_statistics[CompositePixelChannel].mean+=channel_statistics[i].mean;
1235     channel_statistics[CompositePixelChannel].standard_deviation+=
1236       channel_statistics[i].variance-channel_statistics[i].mean*
1237       channel_statistics[i].mean;
1238     area++;
1239   }
1240   channel_statistics[CompositePixelChannel].mean/=area;
1241   channel_statistics[CompositePixelChannel].standard_deviation=
1242     sqrt(channel_statistics[CompositePixelChannel].standard_deviation/area);
1243   *mean=channel_statistics[CompositePixelChannel].mean;
1244   *standard_deviation=
1245     channel_statistics[CompositePixelChannel].standard_deviation;
1246   channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1247     channel_statistics);
1248   return(MagickTrue);
1249 }
1250 \f
1251 /*
1252 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1253 %                                                                             %
1254 %                                                                             %
1255 %                                                                             %
1256 %   G e t I m a g e K u r t o s i s                                           %
1257 %                                                                             %
1258 %                                                                             %
1259 %                                                                             %
1260 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1261 %
1262 %  GetImageKurtosis() returns the kurtosis and skewness of one or more
1263 %  image channels.
1264 %
1265 %  The format of the GetImageKurtosis method is:
1266 %
1267 %      MagickBooleanType GetImageKurtosis(const Image *image,double *kurtosis,
1268 %        double *skewness,ExceptionInfo *exception)
1269 %
1270 %  A description of each parameter follows:
1271 %
1272 %    o image: the image.
1273 %
1274 %    o kurtosis: the kurtosis of the channel.
1275 %
1276 %    o skewness: the skewness of the channel.
1277 %
1278 %    o exception: return any errors or warnings in this structure.
1279 %
1280 */
1281 MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
1282   double *kurtosis,double *skewness,ExceptionInfo *exception)
1283 {
1284   CacheView
1285     *image_view;
1286
1287   double
1288     area,
1289     mean,
1290     standard_deviation,
1291     sum_squares,
1292     sum_cubes,
1293     sum_fourth_power;
1294
1295   MagickBooleanType
1296     status;
1297
1298   ssize_t
1299     y;
1300
1301   assert(image != (Image *) NULL);
1302   assert(image->signature == MagickSignature);
1303   if (image->debug != MagickFalse)
1304     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1305   status=MagickTrue;
1306   *kurtosis=0.0;
1307   *skewness=0.0;
1308   area=0.0;
1309   mean=0.0;
1310   standard_deviation=0.0;
1311   sum_squares=0.0;
1312   sum_cubes=0.0;
1313   sum_fourth_power=0.0;
1314   image_view=AcquireVirtualCacheView(image,exception);
1315 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1316   #pragma omp parallel for schedule(static) shared(status) \
1317     dynamic_number_threads(image,image->columns,image->rows,1)
1318 #endif
1319   for (y=0; y < (ssize_t) image->rows; y++)
1320   {
1321     register const Quantum
1322       *restrict p;
1323
1324     register ssize_t
1325       x;
1326
1327     if (status == MagickFalse)
1328       continue;
1329     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1330     if (p == (const Quantum *) NULL)
1331       {
1332         status=MagickFalse;
1333         continue;
1334       }
1335     for (x=0; x < (ssize_t) image->columns; x++)
1336     {
1337       register ssize_t
1338         i;
1339
1340       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1341       {
1342         PixelChannel
1343           channel;
1344
1345         PixelTrait
1346           traits;
1347
1348         channel=GetPixelChannelMapChannel(image,i);
1349         traits=GetPixelChannelMapTraits(image,channel);
1350         if (traits == UndefinedPixelTrait)
1351           continue;
1352         if (((traits & UpdatePixelTrait) == 0) ||
1353             (GetPixelMask(image,p) != 0))
1354           continue;
1355 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1356         #pragma omp critical (MagickCore_GetImageKurtosis)
1357 #endif
1358         {
1359           mean+=p[i];
1360           sum_squares+=(double) p[i]*p[i];
1361           sum_cubes+=(double) p[i]*p[i]*p[i];
1362           sum_fourth_power+=(double) p[i]*p[i]*p[i]*p[i];
1363           area++;
1364         }
1365       }
1366       p+=GetPixelChannels(image);
1367     }
1368   }
1369   image_view=DestroyCacheView(image_view);
1370   if (area != 0.0)
1371     {
1372       mean/=area;
1373       sum_squares/=area;
1374       sum_cubes/=area;
1375       sum_fourth_power/=area;
1376     }
1377   standard_deviation=sqrt(sum_squares-(mean*mean));
1378   if (standard_deviation != 0.0)
1379     {
1380       *kurtosis=sum_fourth_power-4.0*mean*sum_cubes+6.0*mean*mean*sum_squares-
1381         3.0*mean*mean*mean*mean;
1382       *kurtosis/=standard_deviation*standard_deviation*standard_deviation*
1383         standard_deviation;
1384       *kurtosis-=3.0;
1385       *skewness=sum_cubes-3.0*mean*sum_squares+2.0*mean*mean*mean;
1386       *skewness/=standard_deviation*standard_deviation*standard_deviation;
1387     }
1388   return(status);
1389 }
1390 \f
1391 /*
1392 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1393 %                                                                             %
1394 %                                                                             %
1395 %                                                                             %
1396 %   G e t I m a g e R a n g e                                                 %
1397 %                                                                             %
1398 %                                                                             %
1399 %                                                                             %
1400 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1401 %
1402 %  GetImageRange() returns the range of one or more image channels.
1403 %
1404 %  The format of the GetImageRange method is:
1405 %
1406 %      MagickBooleanType GetImageRange(const Image *image,double *minima,
1407 %        double *maxima,ExceptionInfo *exception)
1408 %
1409 %  A description of each parameter follows:
1410 %
1411 %    o image: the image.
1412 %
1413 %    o minima: the minimum value in the channel.
1414 %
1415 %    o maxima: the maximum value in the channel.
1416 %
1417 %    o exception: return any errors or warnings in this structure.
1418 %
1419 */
1420 MagickExport MagickBooleanType GetImageRange(const Image *image,double *minima,
1421   double *maxima,ExceptionInfo *exception)
1422 {
1423   CacheView
1424     *image_view;
1425
1426   MagickBooleanType
1427     initialize,
1428     status;
1429
1430   ssize_t
1431     y;
1432
1433   assert(image != (Image *) NULL);
1434   assert(image->signature == MagickSignature);
1435   if (image->debug != MagickFalse)
1436     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1437   status=MagickTrue;
1438   initialize=MagickTrue;
1439   *maxima=0.0;
1440   *minima=0.0;
1441   image_view=AcquireVirtualCacheView(image,exception);
1442 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1443   #pragma omp parallel for schedule(static) shared(status,initialize) \
1444     dynamic_number_threads(image,image->columns,image->rows,1)
1445 #endif
1446   for (y=0; y < (ssize_t) image->rows; y++)
1447   {
1448     register const Quantum
1449       *restrict p;
1450
1451     register ssize_t
1452       x;
1453
1454     if (status == MagickFalse)
1455       continue;
1456     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1457     if (p == (const Quantum *) NULL)
1458       {
1459         status=MagickFalse;
1460         continue;
1461       }
1462     for (x=0; x < (ssize_t) image->columns; x++)
1463     {
1464       register ssize_t
1465         i;
1466
1467       if (GetPixelMask(image,p) != 0)
1468         {
1469           p+=GetPixelChannels(image);
1470           continue;
1471         }
1472       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1473       {
1474         PixelChannel
1475           channel;
1476
1477         PixelTrait
1478           traits;
1479
1480         channel=GetPixelChannelMapChannel(image,i);
1481         traits=GetPixelChannelMapTraits(image,channel);
1482         if (traits == UndefinedPixelTrait)
1483           continue;
1484         if ((traits & UpdatePixelTrait) == 0)
1485           continue;
1486 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1487         #pragma omp critical (MagickCore_GetImageRange)
1488 #endif
1489         {
1490           if (initialize != MagickFalse)
1491             {
1492               *minima=(double) p[i];
1493               *maxima=(double) p[i];
1494               initialize=MagickFalse;
1495             }
1496           else
1497             {
1498               if ((double) p[i] < *minima)
1499                 *minima=(double) p[i];
1500               if ((double) p[i] > *maxima)
1501                 *maxima=(double) p[i];
1502            }
1503         }
1504       }
1505       p+=GetPixelChannels(image);
1506     }
1507   }
1508   image_view=DestroyCacheView(image_view);
1509   return(status);
1510 }
1511 \f
1512 /*
1513 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1514 %                                                                             %
1515 %                                                                             %
1516 %                                                                             %
1517 %   G e t I m a g e S t a t i s t i c s                                       %
1518 %                                                                             %
1519 %                                                                             %
1520 %                                                                             %
1521 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1522 %
1523 %  GetImageStatistics() returns statistics for each channel in the image.  The
1524 %  statistics include the channel depth, its minima, maxima, mean, standard
1525 %  deviation, kurtosis and skewness.  You can access the red channel mean, for
1526 %  example, like this:
1527 %
1528 %      channel_statistics=GetImageStatistics(image,exception);
1529 %      red_mean=channel_statistics[RedPixelChannel].mean;
1530 %
1531 %  Use MagickRelinquishMemory() to free the statistics buffer.
1532 %
1533 %  The format of the GetImageStatistics method is:
1534 %
1535 %      ChannelStatistics *GetImageStatistics(const Image *image,
1536 %        ExceptionInfo *exception)
1537 %
1538 %  A description of each parameter follows:
1539 %
1540 %    o image: the image.
1541 %
1542 %    o exception: return any errors or warnings in this structure.
1543 %
1544 */
1545
1546 static size_t GetImageChannels(const Image *image)
1547 {
1548   register ssize_t
1549     i;
1550
1551   size_t
1552     channels;
1553
1554   channels=0;
1555   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1556   {
1557     PixelChannel
1558       channel;
1559
1560     PixelTrait
1561       traits;
1562
1563     channel=GetPixelChannelMapChannel(image,i);
1564     traits=GetPixelChannelMapTraits(image,channel);
1565     if ((traits & UpdatePixelTrait) != 0)
1566       channels++;
1567   }
1568   return(channels);
1569 }
1570
1571 MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
1572   ExceptionInfo *exception)
1573 {
1574   ChannelStatistics
1575     *channel_statistics;
1576
1577   MagickStatusType
1578     status;
1579
1580   QuantumAny
1581     range;
1582
1583   register ssize_t
1584     i;
1585
1586   size_t
1587     channels,
1588     depth;
1589
1590   ssize_t
1591     y;
1592
1593   assert(image != (Image *) NULL);
1594   assert(image->signature == MagickSignature);
1595   if (image->debug != MagickFalse)
1596     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1597   channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(
1598     MaxPixelChannels+1,sizeof(*channel_statistics));
1599   if (channel_statistics == (ChannelStatistics *) NULL)
1600     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1601   (void) ResetMagickMemory(channel_statistics,0,(MaxPixelChannels+1)*
1602     sizeof(*channel_statistics));
1603   for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
1604   {
1605     channel_statistics[i].depth=1;
1606     channel_statistics[i].maxima=(-MagickHuge);
1607     channel_statistics[i].minima=MagickHuge;
1608   }
1609   for (y=0; y < (ssize_t) image->rows; y++)
1610   {
1611     register const Quantum
1612       *restrict p;
1613
1614     register ssize_t
1615       x;
1616
1617     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1618     if (p == (const Quantum *) NULL)
1619       break;
1620     for (x=0; x < (ssize_t) image->columns; x++)
1621     {
1622       register ssize_t
1623         i;
1624
1625       if (GetPixelMask(image,p) != 0)
1626         {
1627           p+=GetPixelChannels(image);
1628           continue;
1629         }
1630       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1631       {
1632         PixelChannel
1633           channel;
1634
1635         PixelTrait
1636           traits;
1637
1638         channel=GetPixelChannelMapChannel(image,i);
1639         traits=GetPixelChannelMapTraits(image,channel);
1640         if (traits == UndefinedPixelTrait)
1641           continue;
1642         if (channel_statistics[channel].depth != MAGICKCORE_QUANTUM_DEPTH)
1643           {
1644             depth=channel_statistics[channel].depth;
1645             range=GetQuantumRange(depth);
1646             status=p[i] != ScaleAnyToQuantum(ScaleQuantumToAny(p[i],range),
1647               range) ? MagickTrue : MagickFalse;
1648             if (status != MagickFalse)
1649               {
1650                 channel_statistics[channel].depth++;
1651                 i--;
1652                 continue;
1653               }
1654           }
1655         if ((double) p[i] < channel_statistics[channel].minima)
1656           channel_statistics[channel].minima=(double) p[i];
1657         if ((double) p[i] > channel_statistics[channel].maxima)
1658           channel_statistics[channel].maxima=(double) p[i];
1659         channel_statistics[channel].sum+=p[i];
1660         channel_statistics[channel].sum_squared+=(double) p[i]*p[i];
1661         channel_statistics[channel].sum_cubed+=(double) p[i]*p[i]*p[i];
1662         channel_statistics[channel].sum_fourth_power+=(double) p[i]*p[i]*p[i]*
1663           p[i];
1664         channel_statistics[channel].area++;
1665       }
1666       p+=GetPixelChannels(image);
1667     }
1668   }
1669   for (i=0; i < (ssize_t) MaxPixelChannels; i++)
1670   {
1671     double
1672       area;
1673
1674     area=MagickEpsilonReciprocal(channel_statistics[i].area);
1675     channel_statistics[i].sum*=area;
1676     channel_statistics[i].sum_squared*=area;
1677     channel_statistics[i].sum_cubed*=area;
1678     channel_statistics[i].sum_fourth_power*=area;
1679     channel_statistics[i].mean=channel_statistics[i].sum;
1680     channel_statistics[i].variance=channel_statistics[i].sum_squared;
1681     channel_statistics[i].standard_deviation=sqrt(
1682       channel_statistics[i].variance-(channel_statistics[i].mean*
1683       channel_statistics[i].mean));
1684   }
1685   for (i=0; i < (ssize_t) MaxPixelChannels; i++)
1686   {
1687     channel_statistics[CompositePixelChannel].depth=(size_t) EvaluateMax(
1688       (double) channel_statistics[CompositePixelChannel].depth,(double)
1689       channel_statistics[i].depth);
1690     channel_statistics[CompositePixelChannel].minima=MagickMin(
1691       channel_statistics[CompositePixelChannel].minima,
1692       channel_statistics[i].minima);
1693     channel_statistics[CompositePixelChannel].maxima=EvaluateMax(
1694       channel_statistics[CompositePixelChannel].maxima,
1695       channel_statistics[i].maxima);
1696     channel_statistics[CompositePixelChannel].sum+=channel_statistics[i].sum;
1697     channel_statistics[CompositePixelChannel].sum_squared+=
1698       channel_statistics[i].sum_squared;
1699     channel_statistics[CompositePixelChannel].sum_cubed+=
1700       channel_statistics[i].sum_cubed;
1701     channel_statistics[CompositePixelChannel].sum_fourth_power+=
1702       channel_statistics[i].sum_fourth_power;
1703     channel_statistics[CompositePixelChannel].mean+=channel_statistics[i].mean;
1704     channel_statistics[CompositePixelChannel].variance+=
1705       channel_statistics[i].variance-channel_statistics[i].mean*
1706       channel_statistics[i].mean;
1707     channel_statistics[CompositePixelChannel].standard_deviation+=
1708       channel_statistics[i].variance-channel_statistics[i].mean*
1709       channel_statistics[i].mean;
1710   }
1711   channels=GetImageChannels(image);
1712   channel_statistics[CompositePixelChannel].sum/=channels;
1713   channel_statistics[CompositePixelChannel].sum_squared/=channels;
1714   channel_statistics[CompositePixelChannel].sum_cubed/=channels;
1715   channel_statistics[CompositePixelChannel].sum_fourth_power/=channels;
1716   channel_statistics[CompositePixelChannel].mean/=channels;
1717   channel_statistics[CompositePixelChannel].variance/=channels;
1718   channel_statistics[CompositePixelChannel].standard_deviation=
1719     sqrt(channel_statistics[CompositePixelChannel].standard_deviation/channels);
1720   channel_statistics[CompositePixelChannel].kurtosis/=channels;
1721   channel_statistics[CompositePixelChannel].skewness/=channels;
1722   for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
1723   {
1724     double
1725       standard_deviation;
1726
1727     standard_deviation=MagickEpsilonReciprocal(
1728       channel_statistics[i].standard_deviation);
1729     channel_statistics[i].skewness=(channel_statistics[i].sum_cubed-3.0*
1730       channel_statistics[i].mean*channel_statistics[i].sum_squared+2.0*
1731       channel_statistics[i].mean*channel_statistics[i].mean*
1732       channel_statistics[i].mean)*(standard_deviation*standard_deviation*
1733       standard_deviation);
1734     channel_statistics[i].kurtosis=(channel_statistics[i].sum_fourth_power-4.0*
1735       channel_statistics[i].mean*channel_statistics[i].sum_cubed+6.0*
1736       channel_statistics[i].mean*channel_statistics[i].mean*
1737       channel_statistics[i].sum_squared-3.0*channel_statistics[i].mean*
1738       channel_statistics[i].mean*1.0*channel_statistics[i].mean*
1739       channel_statistics[i].mean)*(standard_deviation*standard_deviation*
1740       standard_deviation*standard_deviation)-3.0;
1741   }
1742   return(channel_statistics);
1743 }
1744 \f
1745 /*
1746 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1747 %                                                                             %
1748 %                                                                             %
1749 %                                                                             %
1750 %     S t a t i s t i c I m a g e                                             %
1751 %                                                                             %
1752 %                                                                             %
1753 %                                                                             %
1754 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1755 %
1756 %  StatisticImage() makes each pixel the min / max / median / mode / etc. of
1757 %  the neighborhood of the specified width and height.
1758 %
1759 %  The format of the StatisticImage method is:
1760 %
1761 %      Image *StatisticImage(const Image *image,const StatisticType type,
1762 %        const size_t width,const size_t height,ExceptionInfo *exception)
1763 %
1764 %  A description of each parameter follows:
1765 %
1766 %    o image: the image.
1767 %
1768 %    o type: the statistic type (median, mode, etc.).
1769 %
1770 %    o width: the width of the pixel neighborhood.
1771 %
1772 %    o height: the height of the pixel neighborhood.
1773 %
1774 %    o exception: return any errors or warnings in this structure.
1775 %
1776 */
1777
1778 typedef struct _SkipNode
1779 {
1780   size_t
1781     next[9],
1782     count,
1783     signature;
1784 } SkipNode;
1785
1786 typedef struct _SkipList
1787 {
1788   ssize_t
1789     level;
1790
1791   SkipNode
1792     *nodes;
1793 } SkipList;
1794
1795 typedef struct _PixelList
1796 {
1797   size_t
1798     length,
1799     seed;
1800
1801   SkipList
1802     skip_list;
1803
1804   size_t
1805     signature;
1806 } PixelList;
1807
1808 static PixelList *DestroyPixelList(PixelList *pixel_list)
1809 {
1810   if (pixel_list == (PixelList *) NULL)
1811     return((PixelList *) NULL);
1812   if (pixel_list->skip_list.nodes != (SkipNode *) NULL)
1813     pixel_list->skip_list.nodes=(SkipNode *) RelinquishMagickMemory(
1814       pixel_list->skip_list.nodes);
1815   pixel_list=(PixelList *) RelinquishMagickMemory(pixel_list);
1816   return(pixel_list);
1817 }
1818
1819 static PixelList **DestroyPixelListThreadSet(PixelList **pixel_list)
1820 {
1821   register ssize_t
1822     i;
1823
1824   assert(pixel_list != (PixelList **) NULL);
1825   for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
1826     if (pixel_list[i] != (PixelList *) NULL)
1827       pixel_list[i]=DestroyPixelList(pixel_list[i]);
1828   pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list);
1829   return(pixel_list);
1830 }
1831
1832 static PixelList *AcquirePixelList(const size_t width,const size_t height)
1833 {
1834   PixelList
1835     *pixel_list;
1836
1837   pixel_list=(PixelList *) AcquireMagickMemory(sizeof(*pixel_list));
1838   if (pixel_list == (PixelList *) NULL)
1839     return(pixel_list);
1840   (void) ResetMagickMemory((void *) pixel_list,0,sizeof(*pixel_list));
1841   pixel_list->length=width*height;
1842   pixel_list->skip_list.nodes=(SkipNode *) AcquireQuantumMemory(65537UL,
1843     sizeof(*pixel_list->skip_list.nodes));
1844   if (pixel_list->skip_list.nodes == (SkipNode *) NULL)
1845     return(DestroyPixelList(pixel_list));
1846   (void) ResetMagickMemory(pixel_list->skip_list.nodes,0,65537UL*
1847     sizeof(*pixel_list->skip_list.nodes));
1848   pixel_list->signature=MagickSignature;
1849   return(pixel_list);
1850 }
1851
1852 static PixelList **AcquirePixelListThreadSet(const size_t width,
1853   const size_t height)
1854 {
1855   PixelList
1856     **pixel_list;
1857
1858   register ssize_t
1859     i;
1860
1861   size_t
1862     number_threads;
1863
1864   number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
1865   pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
1866     sizeof(*pixel_list));
1867   if (pixel_list == (PixelList **) NULL)
1868     return((PixelList **) NULL);
1869   (void) ResetMagickMemory(pixel_list,0,number_threads*sizeof(*pixel_list));
1870   for (i=0; i < (ssize_t) number_threads; i++)
1871   {
1872     pixel_list[i]=AcquirePixelList(width,height);
1873     if (pixel_list[i] == (PixelList *) NULL)
1874       return(DestroyPixelListThreadSet(pixel_list));
1875   }
1876   return(pixel_list);
1877 }
1878
1879 static void AddNodePixelList(PixelList *pixel_list,const size_t color)
1880 {
1881   register SkipList
1882     *p;
1883
1884   register ssize_t
1885     level;
1886
1887   size_t
1888     search,
1889     update[9];
1890
1891   /*
1892     Initialize the node.
1893   */
1894   p=(&pixel_list->skip_list);
1895   p->nodes[color].signature=pixel_list->signature;
1896   p->nodes[color].count=1;
1897   /*
1898     Determine where it belongs in the list.
1899   */
1900   search=65536UL;
1901   for (level=p->level; level >= 0; level--)
1902   {
1903     while (p->nodes[search].next[level] < color)
1904       search=p->nodes[search].next[level];
1905     update[level]=search;
1906   }
1907   /*
1908     Generate a pseudo-random level for this node.
1909   */
1910   for (level=0; ; level++)
1911   {
1912     pixel_list->seed=(pixel_list->seed*42893621L)+1L;
1913     if ((pixel_list->seed & 0x300) != 0x300)
1914       break;
1915   }
1916   if (level > 8)
1917     level=8;
1918   if (level > (p->level+2))
1919     level=p->level+2;
1920   /*
1921     If we're raising the list's level, link back to the root node.
1922   */
1923   while (level > p->level)
1924   {
1925     p->level++;
1926     update[p->level]=65536UL;
1927   }
1928   /*
1929     Link the node into the skip-list.
1930   */
1931   do
1932   {
1933     p->nodes[color].next[level]=p->nodes[update[level]].next[level];
1934     p->nodes[update[level]].next[level]=color;
1935   } while (level-- > 0);
1936 }
1937
1938 static inline void GetMaximumPixelList(PixelList *pixel_list,Quantum *pixel)
1939 {
1940   register SkipList
1941     *p;
1942
1943   size_t
1944     color,
1945     maximum;
1946
1947   ssize_t
1948     count;
1949
1950   /*
1951     Find the maximum value for each of the color.
1952   */
1953   p=(&pixel_list->skip_list);
1954   color=65536L;
1955   count=0;
1956   maximum=p->nodes[color].next[0];
1957   do
1958   {
1959     color=p->nodes[color].next[0];
1960     if (color > maximum)
1961       maximum=color;
1962     count+=p->nodes[color].count;
1963   } while (count < (ssize_t) pixel_list->length);
1964   *pixel=ScaleShortToQuantum((unsigned short) maximum);
1965 }
1966
1967 static inline void GetMeanPixelList(PixelList *pixel_list,Quantum *pixel)
1968 {
1969   double
1970     sum;
1971
1972   register SkipList
1973     *p;
1974
1975   size_t
1976     color;
1977
1978   ssize_t
1979     count;
1980
1981   /*
1982     Find the mean value for each of the color.
1983   */
1984   p=(&pixel_list->skip_list);
1985   color=65536L;
1986   count=0;
1987   sum=0.0;
1988   do
1989   {
1990     color=p->nodes[color].next[0];
1991     sum+=(double) p->nodes[color].count*color;
1992     count+=p->nodes[color].count;
1993   } while (count < (ssize_t) pixel_list->length);
1994   sum/=pixel_list->length;
1995   *pixel=ScaleShortToQuantum((unsigned short) sum);
1996 }
1997
1998 static inline void GetMedianPixelList(PixelList *pixel_list,Quantum *pixel)
1999 {
2000   register SkipList
2001     *p;
2002
2003   size_t
2004     color;
2005
2006   ssize_t
2007     count;
2008
2009   /*
2010     Find the median value for each of the color.
2011   */
2012   p=(&pixel_list->skip_list);
2013   color=65536L;
2014   count=0;
2015   do
2016   {
2017     color=p->nodes[color].next[0];
2018     count+=p->nodes[color].count;
2019   } while (count <= (ssize_t) (pixel_list->length >> 1));
2020   *pixel=ScaleShortToQuantum((unsigned short) color);
2021 }
2022
2023 static inline void GetMinimumPixelList(PixelList *pixel_list,Quantum *pixel)
2024 {
2025   register SkipList
2026     *p;
2027
2028   size_t
2029     color,
2030     minimum;
2031
2032   ssize_t
2033     count;
2034
2035   /*
2036     Find the minimum value for each of the color.
2037   */
2038   p=(&pixel_list->skip_list);
2039   count=0;
2040   color=65536UL;
2041   minimum=p->nodes[color].next[0];
2042   do
2043   {
2044     color=p->nodes[color].next[0];
2045     if (color < minimum)
2046       minimum=color;
2047     count+=p->nodes[color].count;
2048   } while (count < (ssize_t) pixel_list->length);
2049   *pixel=ScaleShortToQuantum((unsigned short) minimum);
2050 }
2051
2052 static inline void GetModePixelList(PixelList *pixel_list,Quantum *pixel)
2053 {
2054   register SkipList
2055     *p;
2056
2057   size_t
2058     color,
2059     max_count,
2060     mode;
2061
2062   ssize_t
2063     count;
2064
2065   /*
2066     Make each pixel the 'predominant color' of the specified neighborhood.
2067   */
2068   p=(&pixel_list->skip_list);
2069   color=65536L;
2070   mode=color;
2071   max_count=p->nodes[mode].count;
2072   count=0;
2073   do
2074   {
2075     color=p->nodes[color].next[0];
2076     if (p->nodes[color].count > max_count)
2077       {
2078         mode=color;
2079         max_count=p->nodes[mode].count;
2080       }
2081     count+=p->nodes[color].count;
2082   } while (count < (ssize_t) pixel_list->length);
2083   *pixel=ScaleShortToQuantum((unsigned short) mode);
2084 }
2085
2086 static inline void GetNonpeakPixelList(PixelList *pixel_list,Quantum *pixel)
2087 {
2088   register SkipList
2089     *p;
2090
2091   size_t
2092     color,
2093     next,
2094     previous;
2095
2096   ssize_t
2097     count;
2098
2099   /*
2100     Finds the non peak value for each of the colors.
2101   */
2102   p=(&pixel_list->skip_list);
2103   color=65536L;
2104   next=p->nodes[color].next[0];
2105   count=0;
2106   do
2107   {
2108     previous=color;
2109     color=next;
2110     next=p->nodes[color].next[0];
2111     count+=p->nodes[color].count;
2112   } while (count <= (ssize_t) (pixel_list->length >> 1));
2113   if ((previous == 65536UL) && (next != 65536UL))
2114     color=next;
2115   else
2116     if ((previous != 65536UL) && (next == 65536UL))
2117       color=previous;
2118   *pixel=ScaleShortToQuantum((unsigned short) color);
2119 }
2120
2121 static inline void GetStandardDeviationPixelList(PixelList *pixel_list,
2122   Quantum *pixel)
2123 {
2124   double
2125     sum,
2126     sum_squared;
2127
2128   register SkipList
2129     *p;
2130
2131   size_t
2132     color;
2133
2134   ssize_t
2135     count;
2136
2137   /*
2138     Find the standard-deviation value for each of the color.
2139   */
2140   p=(&pixel_list->skip_list);
2141   color=65536L;
2142   count=0;
2143   sum=0.0;
2144   sum_squared=0.0;
2145   do
2146   {
2147     register ssize_t
2148       i;
2149
2150     color=p->nodes[color].next[0];
2151     sum+=(double) p->nodes[color].count*color;
2152     for (i=0; i < (ssize_t) p->nodes[color].count; i++)
2153       sum_squared+=((double) color)*((double) color);
2154     count+=p->nodes[color].count;
2155   } while (count < (ssize_t) pixel_list->length);
2156   sum/=pixel_list->length;
2157   sum_squared/=pixel_list->length;
2158   *pixel=ScaleShortToQuantum((unsigned short) sqrt(sum_squared-(sum*sum)));
2159 }
2160
2161 static inline void InsertPixelList(const Image *image,const Quantum pixel,
2162   PixelList *pixel_list)
2163 {
2164   size_t
2165     signature;
2166
2167   unsigned short
2168     index;
2169
2170   index=ScaleQuantumToShort(pixel);
2171   signature=pixel_list->skip_list.nodes[index].signature;
2172   if (signature == pixel_list->signature)
2173     {
2174       pixel_list->skip_list.nodes[index].count++;
2175       return;
2176     }
2177   AddNodePixelList(pixel_list,index);
2178 }
2179
2180 static inline double MagickAbsoluteValue(const double x)
2181 {
2182   if (x < 0)
2183     return(-x);
2184   return(x);
2185 }
2186
2187 static inline size_t MagickMax(const size_t x,const size_t y)
2188 {
2189   if (x > y)
2190     return(x);
2191   return(y);
2192 }
2193
2194 static void ResetPixelList(PixelList *pixel_list)
2195 {
2196   int
2197     level;
2198
2199   register SkipNode
2200     *root;
2201
2202   register SkipList
2203     *p;
2204
2205   /*
2206     Reset the skip-list.
2207   */
2208   p=(&pixel_list->skip_list);
2209   root=p->nodes+65536UL;
2210   p->level=0;
2211   for (level=0; level < 9; level++)
2212     root->next[level]=65536UL;
2213   pixel_list->seed=pixel_list->signature++;
2214 }
2215
2216 MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
2217   const size_t width,const size_t height,ExceptionInfo *exception)
2218 {
2219 #define StatisticImageTag  "Statistic/Image"
2220
2221   CacheView
2222     *image_view,
2223     *statistic_view;
2224
2225   Image
2226     *statistic_image;
2227
2228   MagickBooleanType
2229     status;
2230
2231   MagickOffsetType
2232     progress;
2233
2234   PixelList
2235     **restrict pixel_list;
2236
2237   ssize_t
2238     center,
2239     y;
2240
2241   /*
2242     Initialize statistics image attributes.
2243   */
2244   assert(image != (Image *) NULL);
2245   assert(image->signature == MagickSignature);
2246   if (image->debug != MagickFalse)
2247     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2248   assert(exception != (ExceptionInfo *) NULL);
2249   assert(exception->signature == MagickSignature);
2250   statistic_image=CloneImage(image,image->columns,image->rows,MagickTrue,
2251     exception);
2252   if (statistic_image == (Image *) NULL)
2253     return((Image *) NULL);
2254   status=SetImageStorageClass(statistic_image,DirectClass,exception);
2255   if (status == MagickFalse)
2256     {
2257       statistic_image=DestroyImage(statistic_image);
2258       return((Image *) NULL);
2259     }
2260   pixel_list=AcquirePixelListThreadSet(MagickMax(width,1),MagickMax(height,1));
2261   if (pixel_list == (PixelList **) NULL)
2262     {
2263       statistic_image=DestroyImage(statistic_image);
2264       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2265     }
2266   /*
2267     Make each pixel the min / max / median / mode / etc. of the neighborhood.
2268   */
2269   center=(ssize_t) GetPixelChannels(image)*(image->columns+MagickMax(width,1))*
2270     (MagickMax(height,1)/2L)+GetPixelChannels(image)*(MagickMax(width,1)/2L);
2271   status=MagickTrue;
2272   progress=0;
2273   image_view=AcquireVirtualCacheView(image,exception);
2274   statistic_view=AcquireAuthenticCacheView(statistic_image,exception);
2275 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2276   #pragma omp parallel for schedule(static,4) shared(progress,status) \
2277     dynamic_number_threads(image,image->columns,image->rows,1)
2278 #endif
2279   for (y=0; y < (ssize_t) statistic_image->rows; y++)
2280   {
2281     const int
2282       id = GetOpenMPThreadId();
2283
2284     register const Quantum
2285       *restrict p;
2286
2287     register Quantum
2288       *restrict q;
2289
2290     register ssize_t
2291       x;
2292
2293     if (status == MagickFalse)
2294       continue;
2295     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) MagickMax(width,1)/2L),y-
2296       (ssize_t) (MagickMax(height,1)/2L),image->columns+MagickMax(width,1),
2297       MagickMax(height,1),exception);
2298     q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns,      1,exception);
2299     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2300       {
2301         status=MagickFalse;
2302         continue;
2303       }
2304     for (x=0; x < (ssize_t) statistic_image->columns; x++)
2305     {
2306       register ssize_t
2307         i;
2308
2309       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2310       {
2311         PixelChannel
2312           channel;
2313
2314         PixelTrait
2315           statistic_traits,
2316           traits;
2317
2318         Quantum
2319           pixel;
2320
2321         register const Quantum
2322           *restrict pixels;
2323
2324         register ssize_t
2325           u;
2326
2327         ssize_t
2328           v;
2329
2330         channel=GetPixelChannelMapChannel(image,i);
2331         traits=GetPixelChannelMapTraits(image,channel);
2332         statistic_traits=GetPixelChannelMapTraits(statistic_image,channel);
2333         if ((traits == UndefinedPixelTrait) ||
2334             (statistic_traits == UndefinedPixelTrait))
2335           continue;
2336         if (((statistic_traits & CopyPixelTrait) != 0) ||
2337             (GetPixelMask(image,p) != 0))
2338           {
2339             SetPixelChannel(statistic_image,channel,p[center+i],q);
2340             continue;
2341           }
2342         pixels=p;
2343         ResetPixelList(pixel_list[id]);
2344         for (v=0; v < (ssize_t) MagickMax(height,1); v++)
2345         {
2346           for (u=0; u < (ssize_t) MagickMax(width,1); u++)
2347           {
2348             InsertPixelList(image,pixels[i],pixel_list[id]);
2349             pixels+=GetPixelChannels(image);
2350           }
2351           pixels+=image->columns*GetPixelChannels(image);
2352         }
2353         switch (type)
2354         {
2355           case GradientStatistic:
2356           {
2357             double
2358               maximum,
2359               minimum;
2360
2361             GetMinimumPixelList(pixel_list[id],&pixel);
2362             minimum=(double) pixel;
2363             GetMaximumPixelList(pixel_list[id],&pixel);
2364             maximum=(double) pixel;
2365             pixel=ClampToQuantum(MagickAbsoluteValue(maximum-minimum));
2366             break;
2367           }
2368           case MaximumStatistic:
2369           {
2370             GetMaximumPixelList(pixel_list[id],&pixel);
2371             break;
2372           }
2373           case MeanStatistic:
2374           {
2375             GetMeanPixelList(pixel_list[id],&pixel);
2376             break;
2377           }
2378           case MedianStatistic:
2379           default:
2380           {
2381             GetMedianPixelList(pixel_list[id],&pixel);
2382             break;
2383           }
2384           case MinimumStatistic:
2385           {
2386             GetMinimumPixelList(pixel_list[id],&pixel);
2387             break;
2388           }
2389           case ModeStatistic:
2390           {
2391             GetModePixelList(pixel_list[id],&pixel);
2392             break;
2393           }
2394           case NonpeakStatistic:
2395           {
2396             GetNonpeakPixelList(pixel_list[id],&pixel);
2397             break;
2398           }
2399           case StandardDeviationStatistic:
2400           {
2401             GetStandardDeviationPixelList(pixel_list[id],&pixel);
2402             break;
2403           }
2404         }
2405         SetPixelChannel(statistic_image,channel,pixel,q);
2406       }
2407       p+=GetPixelChannels(image);
2408       q+=GetPixelChannels(statistic_image);
2409     }
2410     if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
2411       status=MagickFalse;
2412     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2413       {
2414         MagickBooleanType
2415           proceed;
2416
2417 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2418         #pragma omp critical (MagickCore_StatisticImage)
2419 #endif
2420         proceed=SetImageProgress(image,StatisticImageTag,progress++,
2421           image->rows);
2422         if (proceed == MagickFalse)
2423           status=MagickFalse;
2424       }
2425   }
2426   statistic_view=DestroyCacheView(statistic_view);
2427   image_view=DestroyCacheView(image_view);
2428   pixel_list=DestroyPixelListThreadSet(pixel_list);
2429   return(statistic_image);
2430 }