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 %                                   Cristy                                    %
  17 %                                 July 1992                                   %
  18 %                                                                             %
  19 %                                                                             %
  20 %  Copyright 1999-2019 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 %    https://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/accelerate-private.h"
  45 #include "MagickCore/animate.h"
  46 #include "MagickCore/artifact.h"
  47 #include "MagickCore/blob.h"
  48 #include "MagickCore/blob-private.h"
  49 #include "MagickCore/cache.h"
  50 #include "MagickCore/cache-private.h"
  51 #include "MagickCore/cache-view.h"
  52 #include "MagickCore/client.h"
  53 #include "MagickCore/color.h"
  54 #include "MagickCore/color-private.h"
  55 #include "MagickCore/colorspace.h"
  56 #include "MagickCore/colorspace-private.h"
  57 #include "MagickCore/composite.h"
  58 #include "MagickCore/composite-private.h"
  59 #include "MagickCore/compress.h"
  60 #include "MagickCore/constitute.h"
  61 #include "MagickCore/display.h"
  62 #include "MagickCore/draw.h"
  63 #include "MagickCore/enhance.h"
  64 #include "MagickCore/exception.h"
  65 #include "MagickCore/exception-private.h"
  66 #include "MagickCore/gem.h"
  67 #include "MagickCore/gem-private.h"
  68 #include "MagickCore/geometry.h"
  69 #include "MagickCore/list.h"
  70 #include "MagickCore/image-private.h"
  71 #include "MagickCore/magic.h"
  72 #include "MagickCore/magick.h"
  73 #include "MagickCore/memory_.h"
  74 #include "MagickCore/module.h"
  75 #include "MagickCore/monitor.h"
  76 #include "MagickCore/monitor-private.h"
  77 #include "MagickCore/option.h"
  78 #include "MagickCore/paint.h"
  79 #include "MagickCore/pixel-accessor.h"
  80 #include "MagickCore/profile.h"
  81 #include "MagickCore/property.h"
  82 #include "MagickCore/quantize.h"
  83 #include "MagickCore/quantum-private.h"
  84 #include "MagickCore/random_.h"
  85 #include "MagickCore/random-private.h"
  86 #include "MagickCore/resource_.h"
  87 #include "MagickCore/segment.h"
  88 #include "MagickCore/semaphore.h"
  89 #include "MagickCore/signature-private.h"
  90 #include "MagickCore/statistic.h"
  91 #include "MagickCore/string_.h"
  92 #include "MagickCore/thread-private.h"
  93 #include "MagickCore/timer.h"
  94 #include "MagickCore/utility.h"
  95 #include "MagickCore/version.h"
  96 \f
  97 /*
  98 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  99 %                                                                             %
 100 %                                                                             %
 101 %                                                                             %
 102 %     E v a l u a t e I m a g e                                               %
 103 %                                                                             %
 104 %                                                                             %
 105 %                                                                             %
 106 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 107 %
 108 %  EvaluateImage() applies a value to the image with an arithmetic, relational,
 109 %  or logical operator to an image. Use these operations to lighten or darken
 110 %  an image, to increase or decrease contrast in an image, or to produce the
 111 %  "negative" of an image.
 112 %
 113 %  The format of the EvaluateImage method is:
 114 %
 115 %      MagickBooleanType EvaluateImage(Image *image,
 116 %        const MagickEvaluateOperator op,const double value,
 117 %        ExceptionInfo *exception)
 118 %      MagickBooleanType EvaluateImages(Image *images,
 119 %        const MagickEvaluateOperator op,const double value,
 120 %        ExceptionInfo *exception)
 121 %
 122 %  A description of each parameter follows:
 123 %
 124 %    o image: the image.
 125 %
 126 %    o op: A channel op.
 127 %
 128 %    o value: A value value.
 129 %
 130 %    o exception: return any errors or warnings in this structure.
 131 %
 132 */
 133
 134 typedef struct _PixelChannels
 135 {
 136   double
 137     channel[CompositePixelChannel];
 138 } PixelChannels;
 139
 140 static PixelChannels **DestroyPixelThreadSet(const Image *images,
 141   PixelChannels **pixels)
 142 {
 143   register ssize_t
 144     i;
 145
 146   size_t
 147     rows;
 148
 149   assert(pixels != (PixelChannels **) NULL);
 150   rows=MagickMax(GetImageListLength(images),
 151     (size_t) GetMagickResourceLimit(ThreadResource));
 152   for (i=0; i < (ssize_t) rows; i++)
 153     if (pixels[i] != (PixelChannels *) NULL)
 154       pixels[i]=(PixelChannels *) RelinquishMagickMemory(pixels[i]);
 155   pixels=(PixelChannels **) RelinquishMagickMemory(pixels);
 156   return(pixels);
 157 }
 158
 159 static PixelChannels **AcquirePixelThreadSet(const Image *images)
 160 {
 161   const Image
 162     *next;
 163
 164   PixelChannels
 165     **pixels;
 166
 167   register ssize_t
 168     i;
 169
 170   size_t
 171     columns,
 172     rows;
 173
 174   rows=MagickMax(GetImageListLength(images),
 175     (size_t) GetMagickResourceLimit(ThreadResource));
 176   pixels=(PixelChannels **) AcquireQuantumMemory(rows,sizeof(*pixels));
 177   if (pixels == (PixelChannels **) NULL)
 178     return((PixelChannels **) NULL);
 179   (void) memset(pixels,0,rows*sizeof(*pixels));
 180   columns=MagickMax(GetImageListLength(images),MaxPixelChannels);
 181   for (next=images; next != (Image *) NULL; next=next->next)
 182     columns=MagickMax(next->columns,columns);
 183   for (i=0; i < (ssize_t) rows; i++)
 184   {
 185     register ssize_t
 186       j;
 187
 188     pixels[i]=(PixelChannels *) AcquireQuantumMemory(columns,sizeof(**pixels));
 189     if (pixels[i] == (PixelChannels *) NULL)
 190       return(DestroyPixelThreadSet(images,pixels));
 191     for (j=0; j < (ssize_t) columns; j++)
 192     {
 193       register ssize_t
 194         k;
 195
 196       for (k=0; k < MaxPixelChannels; k++)
 197         pixels[i][j].channel[k]=0.0;
 198     }
 199   }
 200   return(pixels);
 201 }
 202
 203 static inline double EvaluateMax(const double x,const double y)
 204 {
 205   if (x > y)
 206     return(x);
 207   return(y);
 208 }
 209
 210 #if defined(__cplusplus) || defined(c_plusplus)
 211 extern "C" {
 212 #endif
 213
 214 static int IntensityCompare(const void *x,const void *y)
 215 {
 216   const PixelChannels
 217     *color_1,
 218     *color_2;
 219
 220   double
 221     distance;
 222
 223   register ssize_t
 224     i;
 225
 226   color_1=(const PixelChannels *) x;
 227   color_2=(const PixelChannels *) y;
 228   distance=0.0;
 229   for (i=0; i < MaxPixelChannels; i++)
 230     distance+=color_1->channel[i]-(double) color_2->channel[i];
 231   return(distance < 0 ? -1 : distance > 0 ? 1 : 0);
 232 }
 233
 234 #if defined(__cplusplus) || defined(c_plusplus)
 235 }
 236 #endif
 237
 238 static double ApplyEvaluateOperator(RandomInfo *random_info,const Quantum pixel,
 239   const MagickEvaluateOperator op,const double value)
 240 {
 241   double
 242     result;
 243
 244   register ssize_t
 245     i;
 246
 247   result=0.0;
 248   switch (op)
 249   {
 250     case UndefinedEvaluateOperator:
 251       break;
 252     case AbsEvaluateOperator:
 253     {
 254       result=(double) fabs((double) (pixel+value));
 255       break;
 256     }
 257     case AddEvaluateOperator:
 258     {
 259       result=(double) (pixel+value);
 260       break;
 261     }
 262     case AddModulusEvaluateOperator:
 263     {
 264       /*
 265         This returns a 'floored modulus' of the addition which is a positive
 266         result.  It differs from % or fmod() that returns a 'truncated modulus'
 267         result, where floor() is replaced by trunc() and could return a
 268         negative result (which is clipped).
 269       */
 270       result=pixel+value;
 271       result-=(QuantumRange+1.0)*floor((double) result/(QuantumRange+1.0));
 272       break;
 273     }
 274     case AndEvaluateOperator:
 275     {
 276       result=(double) ((ssize_t) pixel & (ssize_t) (value+0.5));
 277       break;
 278     }
 279     case CosineEvaluateOperator:
 280     {
 281       result=(double) (QuantumRange*(0.5*cos((double) (2.0*MagickPI*
 282         QuantumScale*pixel*value))+0.5));
 283       break;
 284     }
 285     case DivideEvaluateOperator:
 286     {
 287       result=pixel/(value == 0.0 ? 1.0 : value);
 288       break;
 289     }
 290     case ExponentialEvaluateOperator:
 291     {
 292       result=(double) (QuantumRange*exp((double) (value*QuantumScale*pixel)));
 293       break;
 294     }
 295     case GaussianNoiseEvaluateOperator:
 296     {
 297       result=(double) GenerateDifferentialNoise(random_info,pixel,
 298         GaussianNoise,value);
 299       break;
 300     }
 301     case ImpulseNoiseEvaluateOperator:
 302     {
 303       result=(double) GenerateDifferentialNoise(random_info,pixel,ImpulseNoise,
 304         value);
 305       break;
 306     }
 307     case LaplacianNoiseEvaluateOperator:
 308     {
 309       result=(double) GenerateDifferentialNoise(random_info,pixel,
 310         LaplacianNoise,value);
 311       break;
 312     }
 313     case LeftShiftEvaluateOperator:
 314     {
 315       result=(double) pixel;
 316       for (i=0; i < (ssize_t) value; i++)
 317         result*=2.0;
 318       break;
 319     }
 320     case LogEvaluateOperator:
 321     {
 322       if ((QuantumScale*pixel) >= MagickEpsilon)
 323         result=(double) (QuantumRange*log((double) (QuantumScale*value*pixel+
 324           1.0))/log((double) (value+1.0)));
 325       break;
 326     }
 327     case MaxEvaluateOperator:
 328     {
 329       result=(double) EvaluateMax((double) pixel,value);
 330       break;
 331     }
 332     case MeanEvaluateOperator:
 333     {
 334       result=(double) (pixel+value);
 335       break;
 336     }
 337     case MedianEvaluateOperator:
 338     {
 339       result=(double) (pixel+value);
 340       break;
 341     }
 342     case MinEvaluateOperator:
 343     {
 344       result=(double) MagickMin((double) pixel,value);
 345       break;
 346     }
 347     case MultiplicativeNoiseEvaluateOperator:
 348     {
 349       result=(double) GenerateDifferentialNoise(random_info,pixel,
 350         MultiplicativeGaussianNoise,value);
 351       break;
 352     }
 353     case MultiplyEvaluateOperator:
 354     {
 355       result=(double) (value*pixel);
 356       break;
 357     }
 358     case OrEvaluateOperator:
 359     {
 360       result=(double) ((ssize_t) pixel | (ssize_t) (value+0.5));
 361       break;
 362     }
 363     case PoissonNoiseEvaluateOperator:
 364     {
 365       result=(double) GenerateDifferentialNoise(random_info,pixel,PoissonNoise,
 366         value);
 367       break;
 368     }
 369     case PowEvaluateOperator:
 370     {
 371       if (pixel < 0)
 372         result=(double) -(QuantumRange*pow((double) -(QuantumScale*pixel),
 373           (double) value));
 374       else
 375         result=(double) (QuantumRange*pow((double) (QuantumScale*pixel),
 376           (double) value));
 377       break;
 378     }
 379     case RightShiftEvaluateOperator:
 380     {
 381       result=(double) pixel;
 382       for (i=0; i < (ssize_t) value; i++)
 383         result/=2.0;
 384       break;
 385     }
 386     case RootMeanSquareEvaluateOperator:
 387     {
 388       result=(double) (pixel*pixel+value);
 389       break;
 390     }
 391     case SetEvaluateOperator:
 392     {
 393       result=value;
 394       break;
 395     }
 396     case SineEvaluateOperator:
 397     {
 398       result=(double) (QuantumRange*(0.5*sin((double) (2.0*MagickPI*
 399         QuantumScale*pixel*value))+0.5));
 400       break;
 401     }
 402     case SubtractEvaluateOperator:
 403     {
 404       result=(double) (pixel-value);
 405       break;
 406     }
 407     case SumEvaluateOperator:
 408     {
 409       result=(double) (pixel+value);
 410       break;
 411     }
 412     case ThresholdEvaluateOperator:
 413     {
 414       result=(double) (((double) pixel <= value) ? 0 : QuantumRange);
 415       break;
 416     }
 417     case ThresholdBlackEvaluateOperator:
 418     {
 419       result=(double) (((double) pixel <= value) ? 0 : pixel);
 420       break;
 421     }
 422     case ThresholdWhiteEvaluateOperator:
 423     {
 424       result=(double) (((double) pixel > value) ? QuantumRange : pixel);
 425       break;
 426     }
 427     case UniformNoiseEvaluateOperator:
 428     {
 429       result=(double) GenerateDifferentialNoise(random_info,pixel,UniformNoise,
 430         value);
 431       break;
 432     }
 433     case XorEvaluateOperator:
 434     {
 435       result=(double) ((ssize_t) pixel ^ (ssize_t) (value+0.5));
 436       break;
 437     }
 438   }
 439   return(result);
 440 }
 441
 442 static Image *AcquireImageCanvas(const Image *images,ExceptionInfo *exception)
 443 {
 444   const Image
 445     *p,
 446     *q;
 447
 448   size_t
 449     columns,
 450     rows;
 451
 452   q=images;
 453   columns=images->columns;
 454   rows=images->rows;
 455   for (p=images; p != (Image *) NULL; p=p->next)
 456   {
 457     if (p->number_channels > q->number_channels)
 458       q=p;
 459     if (p->columns > columns)
 460       columns=p->columns;
 461     if (p->rows > rows)
 462       rows=p->rows;
 463   }
 464   return(CloneImage(q,columns,rows,MagickTrue,exception));
 465 }
 466
 467 MagickExport Image *EvaluateImages(const Image *images,
 468   const MagickEvaluateOperator op,ExceptionInfo *exception)
 469 {
 470 #define EvaluateImageTag  "Evaluate/Image"
 471
 472   CacheView
 473     *evaluate_view;
 474
 475   Image
 476     *image;
 477
 478   MagickBooleanType
 479     status;
 480
 481   MagickOffsetType
 482     progress;
 483
 484   PixelChannels
 485     **magick_restrict evaluate_pixels;
 486
 487   RandomInfo
 488     **magick_restrict random_info;
 489
 490   size_t
 491     number_images;
 492
 493   ssize_t
 494     y;
 495
 496 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 497   unsigned long
 498     key;
 499 #endif
 500
 501   assert(images != (Image *) NULL);
 502   assert(images->signature == MagickCoreSignature);
 503   if (images->debug != MagickFalse)
 504     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
 505   assert(exception != (ExceptionInfo *) NULL);
 506   assert(exception->signature == MagickCoreSignature);
 507   image=AcquireImageCanvas(images,exception);
 508   if (image == (Image *) NULL)
 509     return((Image *) NULL);
 510   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
 511     {
 512       image=DestroyImage(image);
 513       return((Image *) NULL);
 514     }
 515   number_images=GetImageListLength(images);
 516   evaluate_pixels=AcquirePixelThreadSet(images);
 517   if (evaluate_pixels == (PixelChannels **) NULL)
 518     {
 519       image=DestroyImage(image);
 520       (void) ThrowMagickException(exception,GetMagickModule(),
 521         ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
 522       return((Image *) NULL);
 523     }
 524   /*
 525     Evaluate image pixels.
 526   */
 527   status=MagickTrue;
 528   progress=0;
 529   random_info=AcquireRandomInfoThreadSet();
 530   evaluate_view=AcquireAuthenticCacheView(image,exception);
 531   if (op == MedianEvaluateOperator)
 532     {
 533 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 534       key=GetRandomSecretKey(random_info[0]);
 535       #pragma omp parallel for schedule(static) shared(progress,status) \
 536         magick_number_threads(image,images,image->rows,key == ~0UL)
 537 #endif
 538       for (y=0; y < (ssize_t) image->rows; y++)
 539       {
 540         CacheView
 541           *image_view;
 542
 543         const Image
 544           *next;
 545
 546         const int
 547           id = GetOpenMPThreadId();
 548
 549         register PixelChannels
 550           *evaluate_pixel;
 551
 552         register Quantum
 553           *magick_restrict q;
 554
 555         register ssize_t
 556           x;
 557
 558         if (status == MagickFalse)
 559           continue;
 560         q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,image->columns,1,
 561           exception);
 562         if (q == (Quantum *) NULL)
 563           {
 564             status=MagickFalse;
 565             continue;
 566           }
 567         evaluate_pixel=evaluate_pixels[id];
 568         for (x=0; x < (ssize_t) image->columns; x++)
 569         {
 570           register ssize_t
 571             j,
 572             k;
 573
 574           for (j=0; j < (ssize_t) number_images; j++)
 575             for (k=0; k < MaxPixelChannels; k++)
 576               evaluate_pixel[j].channel[k]=0.0;
 577           next=images;
 578           for (j=0; j < (ssize_t) number_images; j++)
 579           {
 580             register const Quantum
 581               *p;
 582
 583             register ssize_t
 584               i;
 585
 586             image_view=AcquireVirtualCacheView(next,exception);
 587             p=GetCacheViewVirtualPixels(image_view,x,y,1,1,exception);
 588             if (p == (const Quantum *) NULL)
 589               {
 590                 image_view=DestroyCacheView(image_view);
 591                 break;
 592               }
 593             for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 594             {
 595               PixelChannel channel = GetPixelChannelChannel(image,i);
 596               PixelTrait traits = GetPixelChannelTraits(next,channel);
 597               PixelTrait evaluate_traits = GetPixelChannelTraits(image,channel);
 598               if ((traits == UndefinedPixelTrait) ||
 599                   (evaluate_traits == UndefinedPixelTrait))
 600                 continue;
 601               if ((traits & UpdatePixelTrait) == 0)
 602                 continue;
 603               evaluate_pixel[j].channel[i]=ApplyEvaluateOperator(
 604                 random_info[id],GetPixelChannel(next,channel,p),op,
 605                 evaluate_pixel[j].channel[i]);
 606             }
 607             image_view=DestroyCacheView(image_view);
 608             next=GetNextImageInList(next);
 609           }
 610           qsort((void *) evaluate_pixel,number_images,sizeof(*evaluate_pixel),
 611             IntensityCompare);
 612           for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
 613             q[k]=ClampToQuantum(evaluate_pixel[j/2].channel[k]);
 614           q+=GetPixelChannels(image);
 615         }
 616         if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
 617           status=MagickFalse;
 618         if (images->progress_monitor != (MagickProgressMonitor) NULL)
 619           {
 620             MagickBooleanType
 621               proceed;
 622
 623 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 624             #pragma omp atomic
 625 #endif
 626             progress++;
 627             proceed=SetImageProgress(images,EvaluateImageTag,progress,
 628               image->rows);
 629             if (proceed == MagickFalse)
 630               status=MagickFalse;
 631           }
 632       }
 633     }
 634   else
 635     {
 636 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 637       key=GetRandomSecretKey(random_info[0]);
 638       #pragma omp parallel for schedule(static) shared(progress,status) \
 639         magick_number_threads(image,images,image->rows,key == ~0UL)
 640 #endif
 641       for (y=0; y < (ssize_t) image->rows; y++)
 642       {
 643         CacheView
 644           *image_view;
 645
 646         const Image
 647           *next;
 648
 649         const int
 650           id = GetOpenMPThreadId();
 651
 652         register ssize_t
 653           i,
 654           x;
 655
 656         register PixelChannels
 657           *evaluate_pixel;
 658
 659         register Quantum
 660           *magick_restrict q;
 661
 662         ssize_t
 663           j;
 664
 665         if (status == MagickFalse)
 666           continue;
 667         q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,image->columns,1,
 668           exception);
 669         if (q == (Quantum *) NULL)
 670           {
 671             status=MagickFalse;
 672             continue;
 673           }
 674         evaluate_pixel=evaluate_pixels[id];
 675         for (j=0; j < (ssize_t) image->columns; j++)
 676           for (i=0; i < MaxPixelChannels; i++)
 677             evaluate_pixel[j].channel[i]=0.0;
 678         next=images;
 679         for (j=0; j < (ssize_t) number_images; j++)
 680         {
 681           register const Quantum
 682             *p;
 683
 684           image_view=AcquireVirtualCacheView(next,exception);
 685           p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,
 686             exception);
 687           if (p == (const Quantum *) NULL)
 688             {
 689               image_view=DestroyCacheView(image_view);
 690               break;
 691             }
 692           for (x=0; x < (ssize_t) image->columns; x++)
 693           {
 694             register ssize_t
 695               i;
 696
 697             for (i=0; i < (ssize_t) GetPixelChannels(next); i++)
 698             {
 699               PixelChannel channel = GetPixelChannelChannel(image,i);
 700               PixelTrait traits = GetPixelChannelTraits(next,channel);
 701               PixelTrait evaluate_traits = GetPixelChannelTraits(image,channel);
 702               if ((traits == UndefinedPixelTrait) ||
 703                   (evaluate_traits == UndefinedPixelTrait))
 704                 continue;
 705               if ((traits & UpdatePixelTrait) == 0)
 706                 continue;
 707               evaluate_pixel[x].channel[i]=ApplyEvaluateOperator(
 708                 random_info[id],GetPixelChannel(next,channel,p),j == 0 ?
 709                 AddEvaluateOperator : op,evaluate_pixel[x].channel[i]);
 710             }
 711             p+=GetPixelChannels(next);
 712           }
 713           image_view=DestroyCacheView(image_view);
 714           next=GetNextImageInList(next);
 715         }
 716         for (x=0; x < (ssize_t) image->columns; x++)
 717         {
 718           register ssize_t
 719              i;
 720
 721           switch (op)
 722           {
 723             case MeanEvaluateOperator:
 724             {
 725               for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 726                 evaluate_pixel[x].channel[i]/=(double) number_images;
 727               break;
 728             }
 729             case MultiplyEvaluateOperator:
 730             {
 731               for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 732               {
 733                 register ssize_t
 734                   j;
 735
 736                 for (j=0; j < (ssize_t) (number_images-1); j++)
 737                   evaluate_pixel[x].channel[i]*=QuantumScale;
 738               }
 739               break;
 740             }
 741             case RootMeanSquareEvaluateOperator:
 742             {
 743               for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 744                 evaluate_pixel[x].channel[i]=sqrt(evaluate_pixel[x].channel[i]/
 745                   number_images);
 746               break;
 747             }
 748             default:
 749               break;
 750           }
 751         }
 752         for (x=0; x < (ssize_t) image->columns; x++)
 753         {
 754           register ssize_t
 755             i;
 756
 757           for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 758           {
 759             PixelChannel channel = GetPixelChannelChannel(image,i);
 760             PixelTrait traits = GetPixelChannelTraits(image,channel);
 761             if (traits == UndefinedPixelTrait)
 762               continue;
 763             if ((traits & UpdatePixelTrait) == 0)
 764               continue;
 765             q[i]=ClampToQuantum(evaluate_pixel[x].channel[i]);
 766           }
 767           q+=GetPixelChannels(image);
 768         }
 769         if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
 770           status=MagickFalse;
 771         if (images->progress_monitor != (MagickProgressMonitor) NULL)
 772           {
 773             MagickBooleanType
 774               proceed;
 775
 776 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 777             #pragma omp atomic
 778 #endif
 779             progress++;
 780             proceed=SetImageProgress(images,EvaluateImageTag,progress,
 781               image->rows);
 782             if (proceed == MagickFalse)
 783               status=MagickFalse;
 784           }
 785       }
 786     }
 787   evaluate_view=DestroyCacheView(evaluate_view);
 788   evaluate_pixels=DestroyPixelThreadSet(images,evaluate_pixels);
 789   random_info=DestroyRandomInfoThreadSet(random_info);
 790   if (status == MagickFalse)
 791     image=DestroyImage(image);
 792   return(image);
 793 }
 794
 795 MagickExport MagickBooleanType EvaluateImage(Image *image,
 796   const MagickEvaluateOperator op,const double value,ExceptionInfo *exception)
 797 {
 798   CacheView
 799     *image_view;
 800
 801   MagickBooleanType
 802     status;
 803
 804   MagickOffsetType
 805     progress;
 806
 807   RandomInfo
 808     **magick_restrict random_info;
 809
 810   ssize_t
 811     y;
 812
 813 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 814   unsigned long
 815     key;
 816 #endif
 817
 818   assert(image != (Image *) NULL);
 819   assert(image->signature == MagickCoreSignature);
 820   if (image->debug != MagickFalse)
 821     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 822   assert(exception != (ExceptionInfo *) NULL);
 823   assert(exception->signature == MagickCoreSignature);
 824   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
 825     return(MagickFalse);
 826   status=MagickTrue;
 827   progress=0;
 828   random_info=AcquireRandomInfoThreadSet();
 829   image_view=AcquireAuthenticCacheView(image,exception);
 830 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 831   key=GetRandomSecretKey(random_info[0]);
 832   #pragma omp parallel for schedule(static) shared(progress,status) \
 833     magick_number_threads(image,image,image->rows,key == ~0UL)
 834 #endif
 835   for (y=0; y < (ssize_t) image->rows; y++)
 836   {
 837     const int
 838       id = GetOpenMPThreadId();
 839
 840     register Quantum
 841       *magick_restrict q;
 842
 843     register ssize_t
 844       x;
 845
 846     if (status == MagickFalse)
 847       continue;
 848     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
 849     if (q == (Quantum *) NULL)
 850       {
 851         status=MagickFalse;
 852         continue;
 853       }
 854     for (x=0; x < (ssize_t) image->columns; x++)
 855     {
 856       double
 857         result;
 858
 859       register ssize_t
 860         i;
 861
 862       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 863       {
 864         PixelChannel channel = GetPixelChannelChannel(image,i);
 865         PixelTrait traits = GetPixelChannelTraits(image,channel);
 866         if (traits == UndefinedPixelTrait)
 867           continue;
 868         if ((traits & CopyPixelTrait) != 0)
 869           continue;
 870         if ((traits & UpdatePixelTrait) == 0)
 871           continue;
 872         result=ApplyEvaluateOperator(random_info[id],q[i],op,value);
 873         if (op == MeanEvaluateOperator)
 874           result/=2.0;
 875         q[i]=ClampToQuantum(result);
 876       }
 877       q+=GetPixelChannels(image);
 878     }
 879     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
 880       status=MagickFalse;
 881     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 882       {
 883         MagickBooleanType
 884           proceed;
 885
 886 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 887         #pragma omp atomic
 888 #endif
 889         progress++;
 890         proceed=SetImageProgress(image,EvaluateImageTag,progress,image->rows);
 891         if (proceed == MagickFalse)
 892           status=MagickFalse;
 893       }
 894   }
 895   image_view=DestroyCacheView(image_view);
 896   random_info=DestroyRandomInfoThreadSet(random_info);
 897   return(status);
 898 }
 899 \f
 900 /*
 901 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 902 %                                                                             %
 903 %                                                                             %
 904 %                                                                             %
 905 %     F u n c t i o n I m a g e                                               %
 906 %                                                                             %
 907 %                                                                             %
 908 %                                                                             %
 909 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 910 %
 911 %  FunctionImage() applies a value to the image with an arithmetic, relational,
 912 %  or logical operator to an image. Use these operations to lighten or darken
 913 %  an image, to increase or decrease contrast in an image, or to produce the
 914 %  "negative" of an image.
 915 %
 916 %  The format of the FunctionImage method is:
 917 %
 918 %      MagickBooleanType FunctionImage(Image *image,
 919 %        const MagickFunction function,const ssize_t number_parameters,
 920 %        const double *parameters,ExceptionInfo *exception)
 921 %
 922 %  A description of each parameter follows:
 923 %
 924 %    o image: the image.
 925 %
 926 %    o function: A channel function.
 927 %
 928 %    o parameters: one or more parameters.
 929 %
 930 %    o exception: return any errors or warnings in this structure.
 931 %
 932 */
 933
 934 static Quantum ApplyFunction(Quantum pixel,const MagickFunction function,
 935   const size_t number_parameters,const double *parameters,
 936   ExceptionInfo *exception)
 937 {
 938   double
 939     result;
 940
 941   register ssize_t
 942     i;
 943
 944   (void) exception;
 945   result=0.0;
 946   switch (function)
 947   {
 948     case PolynomialFunction:
 949     {
 950       /*
 951         Polynomial: polynomial constants, highest to lowest order (e.g. c0*x^3+
 952         c1*x^2+c2*x+c3).
 953       */
 954       result=0.0;
 955       for (i=0; i < (ssize_t) number_parameters; i++)
 956         result=result*QuantumScale*pixel+parameters[i];
 957       result*=QuantumRange;
 958       break;
 959     }
 960     case SinusoidFunction:
 961     {
 962       double
 963         amplitude,
 964         bias,
 965         frequency,
 966         phase;
 967
 968       /*
 969         Sinusoid: frequency, phase, amplitude, bias.
 970       */
 971       frequency=(number_parameters >= 1) ? parameters[0] : 1.0;
 972       phase=(number_parameters >= 2) ? parameters[1] : 0.0;
 973       amplitude=(number_parameters >= 3) ? parameters[2] : 0.5;
 974       bias=(number_parameters >= 4) ? parameters[3] : 0.5;
 975       result=(double) (QuantumRange*(amplitude*sin((double) (2.0*
 976         MagickPI*(frequency*QuantumScale*pixel+phase/360.0)))+bias));
 977       break;
 978     }
 979     case ArcsinFunction:
 980     {
 981       double
 982         bias,
 983         center,
 984         range,
 985         width;
 986
 987       /*
 988         Arcsin (peged at range limits for invalid results): width, center,
 989         range, and bias.
 990       */
 991       width=(number_parameters >= 1) ? parameters[0] : 1.0;
 992       center=(number_parameters >= 2) ? parameters[1] : 0.5;
 993       range=(number_parameters >= 3) ? parameters[2] : 1.0;
 994       bias=(number_parameters >= 4) ? parameters[3] : 0.5;
 995       result=2.0/width*(QuantumScale*pixel-center);
 996       if ( result <= -1.0 )
 997         result=bias-range/2.0;
 998       else
 999         if (result >= 1.0)
1000           result=bias+range/2.0;
1001         else
1002           result=(double) (range/MagickPI*asin((double) result)+bias);
1003       result*=QuantumRange;
1004       break;
1005     }
1006     case ArctanFunction:
1007     {
1008       double
1009         center,
1010         bias,
1011         range,
1012         slope;
1013
1014       /*
1015         Arctan: slope, center, range, and bias.
1016       */
1017       slope=(number_parameters >= 1) ? parameters[0] : 1.0;
1018       center=(number_parameters >= 2) ? parameters[1] : 0.5;
1019       range=(number_parameters >= 3) ? parameters[2] : 1.0;
1020       bias=(number_parameters >= 4) ? parameters[3] : 0.5;
1021       result=(double) (MagickPI*slope*(QuantumScale*pixel-center));
1022       result=(double) (QuantumRange*(range/MagickPI*atan((double)
1023         result)+bias));
1024       break;
1025     }
1026     case UndefinedFunction:
1027       break;
1028   }
1029   return(ClampToQuantum(result));
1030 }
1031
1032 MagickExport MagickBooleanType FunctionImage(Image *image,
1033   const MagickFunction function,const size_t number_parameters,
1034   const double *parameters,ExceptionInfo *exception)
1035 {
1036 #define FunctionImageTag  "Function/Image "
1037
1038   CacheView
1039     *image_view;
1040
1041   MagickBooleanType
1042     status;
1043
1044   MagickOffsetType
1045     progress;
1046
1047   ssize_t
1048     y;
1049
1050   assert(image != (Image *) NULL);
1051   assert(image->signature == MagickCoreSignature);
1052   if (image->debug != MagickFalse)
1053     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1054   assert(exception != (ExceptionInfo *) NULL);
1055   assert(exception->signature == MagickCoreSignature);
1056 #if defined(MAGICKCORE_OPENCL_SUPPORT)
1057   if (AccelerateFunctionImage(image,function,number_parameters,parameters,
1058         exception) != MagickFalse)
1059     return(MagickTrue);
1060 #endif
1061   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
1062     return(MagickFalse);
1063   status=MagickTrue;
1064   progress=0;
1065   image_view=AcquireAuthenticCacheView(image,exception);
1066 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1067   #pragma omp parallel for schedule(static) shared(progress,status) \
1068     magick_number_threads(image,image,image->rows,1)
1069 #endif
1070   for (y=0; y < (ssize_t) image->rows; y++)
1071   {
1072     register Quantum
1073       *magick_restrict q;
1074
1075     register ssize_t
1076       x;
1077
1078     if (status == MagickFalse)
1079       continue;
1080     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
1081     if (q == (Quantum *) NULL)
1082       {
1083         status=MagickFalse;
1084         continue;
1085       }
1086     for (x=0; x < (ssize_t) image->columns; x++)
1087     {
1088       register ssize_t
1089         i;
1090
1091       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1092       {
1093         PixelChannel channel = GetPixelChannelChannel(image,i);
1094         PixelTrait traits = GetPixelChannelTraits(image,channel);
1095         if (traits == UndefinedPixelTrait)
1096           continue;
1097         if ((traits & UpdatePixelTrait) == 0)
1098           continue;
1099         q[i]=ApplyFunction(q[i],function,number_parameters,parameters,
1100           exception);
1101       }
1102       q+=GetPixelChannels(image);
1103     }
1104     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1105       status=MagickFalse;
1106     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1107       {
1108         MagickBooleanType
1109           proceed;
1110
1111 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1112         #pragma omp atomic
1113 #endif
1114         progress++;
1115         proceed=SetImageProgress(image,FunctionImageTag,progress,image->rows);
1116         if (proceed == MagickFalse)
1117           status=MagickFalse;
1118       }
1119   }
1120   image_view=DestroyCacheView(image_view);
1121   return(status);
1122 }
1123 \f
1124 /*
1125 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1126 %                                                                             %
1127 %                                                                             %
1128 %                                                                             %
1129 %   G e t I m a g e E n t r o p y                                             %
1130 %                                                                             %
1131 %                                                                             %
1132 %                                                                             %
1133 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1134 %
1135 %  GetImageEntropy() returns the entropy of one or more image channels.
1136 %
1137 %  The format of the GetImageEntropy method is:
1138 %
1139 %      MagickBooleanType GetImageEntropy(const Image *image,double *entropy,
1140 %        ExceptionInfo *exception)
1141 %
1142 %  A description of each parameter follows:
1143 %
1144 %    o image: the image.
1145 %
1146 %    o entropy: the average entropy of the selected channels.
1147 %
1148 %    o exception: return any errors or warnings in this structure.
1149 %
1150 */
1151 MagickExport MagickBooleanType GetImageEntropy(const Image *image,
1152   double *entropy,ExceptionInfo *exception)
1153 {
1154   ChannelStatistics
1155     *channel_statistics;
1156
1157   assert(image != (Image *) NULL);
1158   assert(image->signature == MagickCoreSignature);
1159   if (image->debug != MagickFalse)
1160     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1161   channel_statistics=GetImageStatistics(image,exception);
1162   if (channel_statistics == (ChannelStatistics *) NULL)
1163     return(MagickFalse);
1164   *entropy=channel_statistics[CompositePixelChannel].entropy;
1165   channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1166     channel_statistics);
1167   return(MagickTrue);
1168 }
1169 \f
1170 /*
1171 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1172 %                                                                             %
1173 %                                                                             %
1174 %                                                                             %
1175 %   G e t I m a g e E x t r e m a                                             %
1176 %                                                                             %
1177 %                                                                             %
1178 %                                                                             %
1179 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1180 %
1181 %  GetImageExtrema() returns the extrema of one or more image channels.
1182 %
1183 %  The format of the GetImageExtrema method is:
1184 %
1185 %      MagickBooleanType GetImageExtrema(const Image *image,size_t *minima,
1186 %        size_t *maxima,ExceptionInfo *exception)
1187 %
1188 %  A description of each parameter follows:
1189 %
1190 %    o image: the image.
1191 %
1192 %    o minima: the minimum value in the channel.
1193 %
1194 %    o maxima: the maximum value in the channel.
1195 %
1196 %    o exception: return any errors or warnings in this structure.
1197 %
1198 */
1199 MagickExport MagickBooleanType GetImageExtrema(const Image *image,
1200   size_t *minima,size_t *maxima,ExceptionInfo *exception)
1201 {
1202   double
1203     max,
1204     min;
1205
1206   MagickBooleanType
1207     status;
1208
1209   assert(image != (Image *) NULL);
1210   assert(image->signature == MagickCoreSignature);
1211   if (image->debug != MagickFalse)
1212     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1213   status=GetImageRange(image,&min,&max,exception);
1214   *minima=(size_t) ceil(min-0.5);
1215   *maxima=(size_t) floor(max+0.5);
1216   return(status);
1217 }
1218 \f
1219 /*
1220 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1221 %                                                                             %
1222 %                                                                             %
1223 %                                                                             %
1224 %   G e t I m a g e K u r t o s i s                                           %
1225 %                                                                             %
1226 %                                                                             %
1227 %                                                                             %
1228 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1229 %
1230 %  GetImageKurtosis() returns the kurtosis and skewness of one or more image
1231 %  channels.
1232 %
1233 %  The format of the GetImageKurtosis method is:
1234 %
1235 %      MagickBooleanType GetImageKurtosis(const Image *image,double *kurtosis,
1236 %        double *skewness,ExceptionInfo *exception)
1237 %
1238 %  A description of each parameter follows:
1239 %
1240 %    o image: the image.
1241 %
1242 %    o kurtosis: the kurtosis of the channel.
1243 %
1244 %    o skewness: the skewness of the channel.
1245 %
1246 %    o exception: return any errors or warnings in this structure.
1247 %
1248 */
1249 MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
1250   double *kurtosis,double *skewness,ExceptionInfo *exception)
1251 {
1252   ChannelStatistics
1253     *channel_statistics;
1254
1255   assert(image != (Image *) NULL);
1256   assert(image->signature == MagickCoreSignature);
1257   if (image->debug != MagickFalse)
1258     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1259   channel_statistics=GetImageStatistics(image,exception);
1260   if (channel_statistics == (ChannelStatistics *) NULL)
1261     return(MagickFalse);
1262   *kurtosis=channel_statistics[CompositePixelChannel].kurtosis;
1263   *skewness=channel_statistics[CompositePixelChannel].skewness;
1264   channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1265     channel_statistics);
1266   return(MagickTrue);
1267 }
1268 \f
1269 /*
1270 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1271 %                                                                             %
1272 %                                                                             %
1273 %                                                                             %
1274 %   G e t I m a g e M e a n                                                   %
1275 %                                                                             %
1276 %                                                                             %
1277 %                                                                             %
1278 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1279 %
1280 %  GetImageMean() returns the mean and standard deviation of one or more image
1281 %  channels.
1282 %
1283 %  The format of the GetImageMean method is:
1284 %
1285 %      MagickBooleanType GetImageMean(const Image *image,double *mean,
1286 %        double *standard_deviation,ExceptionInfo *exception)
1287 %
1288 %  A description of each parameter follows:
1289 %
1290 %    o image: the image.
1291 %
1292 %    o mean: the average value in the channel.
1293 %
1294 %    o standard_deviation: the standard deviation of the channel.
1295 %
1296 %    o exception: return any errors or warnings in this structure.
1297 %
1298 */
1299 MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean,
1300   double *standard_deviation,ExceptionInfo *exception)
1301 {
1302   ChannelStatistics
1303     *channel_statistics;
1304
1305   assert(image != (Image *) NULL);
1306   assert(image->signature == MagickCoreSignature);
1307   if (image->debug != MagickFalse)
1308     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1309   channel_statistics=GetImageStatistics(image,exception);
1310   if (channel_statistics == (ChannelStatistics *) NULL)
1311     return(MagickFalse);
1312   *mean=channel_statistics[CompositePixelChannel].mean;
1313   *standard_deviation=
1314     channel_statistics[CompositePixelChannel].standard_deviation;
1315   channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1316     channel_statistics);
1317   return(MagickTrue);
1318 }
1319 \f
1320 /*
1321 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1322 %                                                                             %
1323 %                                                                             %
1324 %                                                                             %
1325 %   G e t I m a g e M o m e n t s                                             %
1326 %                                                                             %
1327 %                                                                             %
1328 %                                                                             %
1329 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1330 %
1331 %  GetImageMoments() returns the normalized moments of one or more image
1332 %  channels.
1333 %
1334 %  The format of the GetImageMoments method is:
1335 %
1336 %      ChannelMoments *GetImageMoments(const Image *image,
1337 %        ExceptionInfo *exception)
1338 %
1339 %  A description of each parameter follows:
1340 %
1341 %    o image: the image.
1342 %
1343 %    o exception: return any errors or warnings in this structure.
1344 %
1345 */
1346
1347 static size_t GetImageChannels(const Image *image)
1348 {
1349   register ssize_t
1350     i;
1351
1352   size_t
1353     channels;
1354
1355   channels=0;
1356   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1357   {
1358     PixelChannel channel = GetPixelChannelChannel(image,i);
1359     PixelTrait traits = GetPixelChannelTraits(image,channel);
1360     if (traits == UndefinedPixelTrait)
1361       continue;
1362     if ((traits & UpdatePixelTrait) == 0)
1363       continue;
1364     channels++;
1365   }
1366   return((size_t) (channels == 0 ? 1 : channels));
1367 }
1368
1369 MagickExport ChannelMoments *GetImageMoments(const Image *image,
1370   ExceptionInfo *exception)
1371 {
1372 #define MaxNumberImageMoments  8
1373
1374   CacheView
1375     *image_view;
1376
1377   ChannelMoments
1378     *channel_moments;
1379
1380   double
1381     M00[MaxPixelChannels+1],
1382     M01[MaxPixelChannels+1],
1383     M02[MaxPixelChannels+1],
1384     M03[MaxPixelChannels+1],
1385     M10[MaxPixelChannels+1],
1386     M11[MaxPixelChannels+1],
1387     M12[MaxPixelChannels+1],
1388     M20[MaxPixelChannels+1],
1389     M21[MaxPixelChannels+1],
1390     M22[MaxPixelChannels+1],
1391     M30[MaxPixelChannels+1];
1392
1393   PointInfo
1394     centroid[MaxPixelChannels+1];
1395
1396   ssize_t
1397     channel,
1398     y;
1399
1400   assert(image != (Image *) NULL);
1401   assert(image->signature == MagickCoreSignature);
1402   if (image->debug != MagickFalse)
1403     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1404   channel_moments=(ChannelMoments *) AcquireQuantumMemory(MaxPixelChannels+1,
1405     sizeof(*channel_moments));
1406   if (channel_moments == (ChannelMoments *) NULL)
1407     return(channel_moments);
1408   (void) memset(channel_moments,0,(MaxPixelChannels+1)*
1409     sizeof(*channel_moments));
1410   (void) memset(centroid,0,sizeof(centroid));
1411   (void) memset(M00,0,sizeof(M00));
1412   (void) memset(M01,0,sizeof(M01));
1413   (void) memset(M02,0,sizeof(M02));
1414   (void) memset(M03,0,sizeof(M03));
1415   (void) memset(M10,0,sizeof(M10));
1416   (void) memset(M11,0,sizeof(M11));
1417   (void) memset(M12,0,sizeof(M12));
1418   (void) memset(M20,0,sizeof(M20));
1419   (void) memset(M21,0,sizeof(M21));
1420   (void) memset(M22,0,sizeof(M22));
1421   (void) memset(M30,0,sizeof(M30));
1422   image_view=AcquireVirtualCacheView(image,exception);
1423   for (y=0; y < (ssize_t) image->rows; y++)
1424   {
1425     register const Quantum
1426       *magick_restrict p;
1427
1428     register ssize_t
1429       x;
1430
1431     /*
1432       Compute center of mass (centroid).
1433     */
1434     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1435     if (p == (const Quantum *) NULL)
1436       break;
1437     for (x=0; x < (ssize_t) image->columns; x++)
1438     {
1439       register ssize_t
1440         i;
1441
1442       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1443       {
1444         PixelChannel channel = GetPixelChannelChannel(image,i);
1445         PixelTrait traits = GetPixelChannelTraits(image,channel);
1446         if (traits == UndefinedPixelTrait)
1447           continue;
1448         if ((traits & UpdatePixelTrait) == 0)
1449           continue;
1450         M00[channel]+=QuantumScale*p[i];
1451         M00[MaxPixelChannels]+=QuantumScale*p[i];
1452         M10[channel]+=x*QuantumScale*p[i];
1453         M10[MaxPixelChannels]+=x*QuantumScale*p[i];
1454         M01[channel]+=y*QuantumScale*p[i];
1455         M01[MaxPixelChannels]+=y*QuantumScale*p[i];
1456       }
1457       p+=GetPixelChannels(image);
1458     }
1459   }
1460   for (channel=0; channel <= MaxPixelChannels; channel++)
1461   {
1462     /*
1463        Compute center of mass (centroid).
1464     */
1465     if (M00[channel] < MagickEpsilon)
1466       {
1467         M00[channel]+=MagickEpsilon;
1468         centroid[channel].x=(double) image->columns/2.0;
1469         centroid[channel].y=(double) image->rows/2.0;
1470         continue;
1471       }
1472     M00[channel]+=MagickEpsilon;
1473     centroid[channel].x=M10[channel]/M00[channel];
1474     centroid[channel].y=M01[channel]/M00[channel];
1475   }
1476   for (y=0; y < (ssize_t) image->rows; y++)
1477   {
1478     register const Quantum
1479       *magick_restrict p;
1480
1481     register ssize_t
1482       x;
1483
1484     /*
1485       Compute the image moments.
1486     */
1487     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1488     if (p == (const Quantum *) NULL)
1489       break;
1490     for (x=0; x < (ssize_t) image->columns; x++)
1491     {
1492       register ssize_t
1493         i;
1494
1495       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1496       {
1497         PixelChannel channel = GetPixelChannelChannel(image,i);
1498         PixelTrait traits = GetPixelChannelTraits(image,channel);
1499         if (traits == UndefinedPixelTrait)
1500           continue;
1501         if ((traits & UpdatePixelTrait) == 0)
1502           continue;
1503         M11[channel]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
1504           QuantumScale*p[i];
1505         M11[MaxPixelChannels]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
1506           QuantumScale*p[i];
1507         M20[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1508           QuantumScale*p[i];
1509         M20[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1510           QuantumScale*p[i];
1511         M02[channel]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
1512           QuantumScale*p[i];
1513         M02[MaxPixelChannels]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
1514           QuantumScale*p[i];
1515         M21[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1516           (y-centroid[channel].y)*QuantumScale*p[i];
1517         M21[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1518           (y-centroid[channel].y)*QuantumScale*p[i];
1519         M12[channel]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
1520           (y-centroid[channel].y)*QuantumScale*p[i];
1521         M12[MaxPixelChannels]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
1522           (y-centroid[channel].y)*QuantumScale*p[i];
1523         M22[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1524           (y-centroid[channel].y)*(y-centroid[channel].y)*QuantumScale*p[i];
1525         M22[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1526           (y-centroid[channel].y)*(y-centroid[channel].y)*QuantumScale*p[i];
1527         M30[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1528           (x-centroid[channel].x)*QuantumScale*p[i];
1529         M30[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1530           (x-centroid[channel].x)*QuantumScale*p[i];
1531         M03[channel]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
1532           (y-centroid[channel].y)*QuantumScale*p[i];
1533         M03[MaxPixelChannels]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
1534           (y-centroid[channel].y)*QuantumScale*p[i];
1535       }
1536       p+=GetPixelChannels(image);
1537     }
1538   }
1539   M00[MaxPixelChannels]/=GetImageChannels(image);
1540   M01[MaxPixelChannels]/=GetImageChannels(image);
1541   M02[MaxPixelChannels]/=GetImageChannels(image);
1542   M03[MaxPixelChannels]/=GetImageChannels(image);
1543   M10[MaxPixelChannels]/=GetImageChannels(image);
1544   M11[MaxPixelChannels]/=GetImageChannels(image);
1545   M12[MaxPixelChannels]/=GetImageChannels(image);
1546   M20[MaxPixelChannels]/=GetImageChannels(image);
1547   M21[MaxPixelChannels]/=GetImageChannels(image);
1548   M22[MaxPixelChannels]/=GetImageChannels(image);
1549   M30[MaxPixelChannels]/=GetImageChannels(image);
1550   for (channel=0; channel <= MaxPixelChannels; channel++)
1551   {
1552     /*
1553       Compute elliptical angle, major and minor axes, eccentricity, & intensity.
1554     */
1555     channel_moments[channel].centroid=centroid[channel];
1556     channel_moments[channel].ellipse_axis.x=sqrt((2.0/M00[channel])*
1557       ((M20[channel]+M02[channel])+sqrt(4.0*M11[channel]*M11[channel]+
1558       (M20[channel]-M02[channel])*(M20[channel]-M02[channel]))));
1559     channel_moments[channel].ellipse_axis.y=sqrt((2.0/M00[channel])*
1560       ((M20[channel]+M02[channel])-sqrt(4.0*M11[channel]*M11[channel]+
1561       (M20[channel]-M02[channel])*(M20[channel]-M02[channel]))));
1562     channel_moments[channel].ellipse_angle=RadiansToDegrees(0.5*atan(2.0*
1563       M11[channel]/(M20[channel]-M02[channel]+MagickEpsilon)));
1564     if (fabs(M11[channel]) < MagickEpsilon)
1565       {
1566         if (fabs(M20[channel]-M02[channel]) < MagickEpsilon)
1567           channel_moments[channel].ellipse_angle+=0.0;
1568         else
1569           if ((M20[channel]-M02[channel]) < 0.0)
1570             channel_moments[channel].ellipse_angle+=90.0;
1571           else
1572             channel_moments[channel].ellipse_angle+=0.0;
1573       }
1574     else
1575       if (M11[channel] < 0.0)
1576         {
1577           if (fabs(M20[channel]-M02[channel]) < MagickEpsilon)
1578             channel_moments[channel].ellipse_angle+=0.0;
1579           else
1580             if ((M20[channel]-M02[channel]) < 0.0)
1581               channel_moments[channel].ellipse_angle+=90.0;
1582             else
1583               channel_moments[channel].ellipse_angle+=180.0;
1584         }
1585       else
1586         {
1587           if (fabs(M20[channel]-M02[channel]) < MagickEpsilon)
1588             channel_moments[channel].ellipse_angle+=0.0;
1589           else
1590             if ((M20[channel]-M02[channel]) < 0.0)
1591               channel_moments[channel].ellipse_angle+=90.0;
1592             else
1593               channel_moments[channel].ellipse_angle+=0.0;
1594        }
1595     channel_moments[channel].ellipse_eccentricity=sqrt(1.0-(
1596       channel_moments[channel].ellipse_axis.y/
1597       (channel_moments[channel].ellipse_axis.x+MagickEpsilon)));
1598     channel_moments[channel].ellipse_intensity=M00[channel]/
1599       (MagickPI*channel_moments[channel].ellipse_axis.x*
1600       channel_moments[channel].ellipse_axis.y+MagickEpsilon);
1601   }
1602   for (channel=0; channel <= MaxPixelChannels; channel++)
1603   {
1604     /*
1605       Normalize image moments.
1606     */
1607     M10[channel]=0.0;
1608     M01[channel]=0.0;
1609     M11[channel]/=pow(M00[channel],1.0+(1.0+1.0)/2.0);
1610     M20[channel]/=pow(M00[channel],1.0+(2.0+0.0)/2.0);
1611     M02[channel]/=pow(M00[channel],1.0+(0.0+2.0)/2.0);
1612     M21[channel]/=pow(M00[channel],1.0+(2.0+1.0)/2.0);
1613     M12[channel]/=pow(M00[channel],1.0+(1.0+2.0)/2.0);
1614     M22[channel]/=pow(M00[channel],1.0+(2.0+2.0)/2.0);
1615     M30[channel]/=pow(M00[channel],1.0+(3.0+0.0)/2.0);
1616     M03[channel]/=pow(M00[channel],1.0+(0.0+3.0)/2.0);
1617     M00[channel]=1.0;
1618   }
1619   image_view=DestroyCacheView(image_view);
1620   for (channel=0; channel <= MaxPixelChannels; channel++)
1621   {
1622     /*
1623       Compute Hu invariant moments.
1624     */
1625     channel_moments[channel].invariant[0]=M20[channel]+M02[channel];
1626     channel_moments[channel].invariant[1]=(M20[channel]-M02[channel])*
1627       (M20[channel]-M02[channel])+4.0*M11[channel]*M11[channel];
1628     channel_moments[channel].invariant[2]=(M30[channel]-3.0*M12[channel])*
1629       (M30[channel]-3.0*M12[channel])+(3.0*M21[channel]-M03[channel])*
1630       (3.0*M21[channel]-M03[channel]);
1631     channel_moments[channel].invariant[3]=(M30[channel]+M12[channel])*
1632       (M30[channel]+M12[channel])+(M21[channel]+M03[channel])*
1633       (M21[channel]+M03[channel]);
1634     channel_moments[channel].invariant[4]=(M30[channel]-3.0*M12[channel])*
1635       (M30[channel]+M12[channel])*((M30[channel]+M12[channel])*
1636       (M30[channel]+M12[channel])-3.0*(M21[channel]+M03[channel])*
1637       (M21[channel]+M03[channel]))+(3.0*M21[channel]-M03[channel])*
1638       (M21[channel]+M03[channel])*(3.0*(M30[channel]+M12[channel])*
1639       (M30[channel]+M12[channel])-(M21[channel]+M03[channel])*
1640       (M21[channel]+M03[channel]));
1641     channel_moments[channel].invariant[5]=(M20[channel]-M02[channel])*
1642       ((M30[channel]+M12[channel])*(M30[channel]+M12[channel])-
1643       (M21[channel]+M03[channel])*(M21[channel]+M03[channel]))+
1644       4.0*M11[channel]*(M30[channel]+M12[channel])*(M21[channel]+M03[channel]);
1645     channel_moments[channel].invariant[6]=(3.0*M21[channel]-M03[channel])*
1646       (M30[channel]+M12[channel])*((M30[channel]+M12[channel])*
1647       (M30[channel]+M12[channel])-3.0*(M21[channel]+M03[channel])*
1648       (M21[channel]+M03[channel]))-(M30[channel]-3*M12[channel])*
1649       (M21[channel]+M03[channel])*(3.0*(M30[channel]+M12[channel])*
1650       (M30[channel]+M12[channel])-(M21[channel]+M03[channel])*
1651       (M21[channel]+M03[channel]));
1652     channel_moments[channel].invariant[7]=M11[channel]*((M30[channel]+
1653       M12[channel])*(M30[channel]+M12[channel])-(M03[channel]+M21[channel])*
1654       (M03[channel]+M21[channel]))-(M20[channel]-M02[channel])*
1655       (M30[channel]+M12[channel])*(M03[channel]+M21[channel]);
1656   }
1657   if (y < (ssize_t) image->rows)
1658     channel_moments=(ChannelMoments *) RelinquishMagickMemory(channel_moments);
1659   return(channel_moments);
1660 }
1661 \f
1662 /*
1663 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1664 %                                                                             %
1665 %                                                                             %
1666 %                                                                             %
1667 %   G e t I m a g e C h a n n e l P e r c e p t u a l H a s h                 %
1668 %                                                                             %
1669 %                                                                             %
1670 %                                                                             %
1671 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1672 %
1673 %  GetImagePerceptualHash() returns the perceptual hash of one or more
1674 %  image channels.
1675 %
1676 %  The format of the GetImagePerceptualHash method is:
1677 %
1678 %      ChannelPerceptualHash *GetImagePerceptualHash(const Image *image,
1679 %        ExceptionInfo *exception)
1680 %
1681 %  A description of each parameter follows:
1682 %
1683 %    o image: the image.
1684 %
1685 %    o exception: return any errors or warnings in this structure.
1686 %
1687 */
1688
1689 static inline double MagickLog10(const double x)
1690 {
1691 #define Log10Epsilon  (1.0e-11)
1692
1693   if (fabs(x) < Log10Epsilon)
1694     return(log10(Log10Epsilon));
1695   return(log10(fabs(x)));
1696 }
1697
1698 MagickExport ChannelPerceptualHash *GetImagePerceptualHash(const Image *image,
1699   ExceptionInfo *exception)
1700 {
1701   ChannelPerceptualHash
1702     *perceptual_hash;
1703
1704   char
1705     *colorspaces,
1706     *q;
1707
1708   const char
1709     *artifact;
1710
1711   MagickBooleanType
1712     status;
1713
1714   register char
1715     *p;
1716
1717   register ssize_t
1718     i;
1719
1720   perceptual_hash=(ChannelPerceptualHash *) AcquireQuantumMemory(
1721     MaxPixelChannels+1UL,sizeof(*perceptual_hash));
1722   if (perceptual_hash == (ChannelPerceptualHash *) NULL)
1723     return((ChannelPerceptualHash *) NULL);
1724   artifact=GetImageArtifact(image,"phash:colorspaces");
1725   if (artifact != NULL)
1726     colorspaces=AcquireString(artifact);
1727   else
1728     colorspaces=AcquireString("sRGB,HCLp");
1729   perceptual_hash[0].number_colorspaces=0;
1730   perceptual_hash[0].number_channels=0;
1731   q=colorspaces;
1732   for (i=0; (p=StringToken(",",&q)) != (char *) NULL; i++)
1733   {
1734     ChannelMoments
1735       *moments;
1736
1737     Image
1738       *hash_image;
1739
1740     size_t
1741       j;
1742
1743     ssize_t
1744       channel,
1745       colorspace;
1746
1747     if (i >= MaximumNumberOfPerceptualColorspaces)
1748       break;
1749     colorspace=ParseCommandOption(MagickColorspaceOptions,MagickFalse,p);
1750     if (colorspace < 0)
1751       break;
1752     perceptual_hash[0].colorspace[i]=(ColorspaceType) colorspace;
1753     hash_image=BlurImage(image,0.0,1.0,exception);
1754     if (hash_image == (Image *) NULL)
1755       break;
1756     hash_image->depth=8;
1757     status=TransformImageColorspace(hash_image,(ColorspaceType) colorspace,
1758       exception);
1759     if (status == MagickFalse)
1760       break;
1761     moments=GetImageMoments(hash_image,exception);
1762     perceptual_hash[0].number_colorspaces++;
1763     perceptual_hash[0].number_channels+=GetImageChannels(hash_image);
1764     hash_image=DestroyImage(hash_image);
1765     if (moments == (ChannelMoments *) NULL)
1766       break;
1767     for (channel=0; channel <= MaxPixelChannels; channel++)
1768       for (j=0; j < MaximumNumberOfImageMoments; j++)
1769         perceptual_hash[channel].phash[i][j]=
1770           (-MagickLog10(moments[channel].invariant[j]));
1771     moments=(ChannelMoments *) RelinquishMagickMemory(moments);
1772   }
1773   colorspaces=DestroyString(colorspaces);
1774   return(perceptual_hash);
1775 }
1776 \f
1777 /*
1778 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1779 %                                                                             %
1780 %                                                                             %
1781 %                                                                             %
1782 %   G e t I m a g e R a n g e                                                 %
1783 %                                                                             %
1784 %                                                                             %
1785 %                                                                             %
1786 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1787 %
1788 %  GetImageRange() returns the range of one or more image channels.
1789 %
1790 %  The format of the GetImageRange method is:
1791 %
1792 %      MagickBooleanType GetImageRange(const Image *image,double *minima,
1793 %        double *maxima,ExceptionInfo *exception)
1794 %
1795 %  A description of each parameter follows:
1796 %
1797 %    o image: the image.
1798 %
1799 %    o minima: the minimum value in the channel.
1800 %
1801 %    o maxima: the maximum value in the channel.
1802 %
1803 %    o exception: return any errors or warnings in this structure.
1804 %
1805 */
1806 MagickExport MagickBooleanType GetImageRange(const Image *image,double *minima,
1807   double *maxima,ExceptionInfo *exception)
1808 {
1809   CacheView
1810     *image_view;
1811
1812   MagickBooleanType
1813     initialize,
1814     status;
1815
1816   ssize_t
1817     y;
1818
1819   assert(image != (Image *) NULL);
1820   assert(image->signature == MagickCoreSignature);
1821   if (image->debug != MagickFalse)
1822     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1823   status=MagickTrue;
1824   initialize=MagickTrue;
1825   *maxima=0.0;
1826   *minima=0.0;
1827   image_view=AcquireVirtualCacheView(image,exception);
1828 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1829   #pragma omp parallel for schedule(static) shared(status,initialize) \
1830     magick_number_threads(image,image,image->rows,1)
1831 #endif
1832   for (y=0; y < (ssize_t) image->rows; y++)
1833   {
1834     double
1835       row_maxima = 0.0,
1836       row_minima = 0.0;
1837
1838     MagickBooleanType
1839       row_initialize;
1840
1841     register const Quantum
1842       *magick_restrict p;
1843
1844     register ssize_t
1845       x;
1846
1847     if (status == MagickFalse)
1848       continue;
1849     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1850     if (p == (const Quantum *) NULL)
1851       {
1852         status=MagickFalse;
1853         continue;
1854       }
1855     row_initialize=MagickTrue;
1856     for (x=0; x < (ssize_t) image->columns; x++)
1857     {
1858       register ssize_t
1859         i;
1860
1861       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1862       {
1863         PixelChannel channel = GetPixelChannelChannel(image,i);
1864         PixelTrait traits = GetPixelChannelTraits(image,channel);
1865         if (traits == UndefinedPixelTrait)
1866           continue;
1867         if ((traits & UpdatePixelTrait) == 0)
1868           continue;
1869                                 if (row_initialize != MagickFalse)
1870           {
1871             row_minima=(double) p[i];
1872             row_maxima=(double) p[i];
1873             row_initialize=MagickFalse;
1874           }
1875         else
1876           {
1877             if ((double) p[i] < row_minima)
1878               row_minima=(double) p[i];
1879             if ((double) p[i] > row_maxima)
1880               row_maxima=(double) p[i];
1881          }
1882       }
1883       p+=GetPixelChannels(image);
1884     }
1885 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1886 #pragma omp critical (MagickCore_GetImageRange)
1887 #endif
1888     {
1889       if (initialize != MagickFalse)
1890         {
1891           *minima=row_minima;
1892           *maxima=row_maxima;
1893           initialize=MagickFalse;
1894         }
1895       else
1896         {
1897           if (row_minima < *minima)
1898             *minima=row_minima;
1899           if (row_maxima > *maxima)
1900             *maxima=row_maxima;
1901         }
1902     }
1903   }
1904   image_view=DestroyCacheView(image_view);
1905   return(status);
1906 }
1907 \f
1908 /*
1909 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1910 %                                                                             %
1911 %                                                                             %
1912 %                                                                             %
1913 %   G e t I m a g e S t a t i s t i c s                                       %
1914 %                                                                             %
1915 %                                                                             %
1916 %                                                                             %
1917 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1918 %
1919 %  GetImageStatistics() returns statistics for each channel in the image.  The
1920 %  statistics include the channel depth, its minima, maxima, mean, standard
1921 %  deviation, kurtosis and skewness.  You can access the red channel mean, for
1922 %  example, like this:
1923 %
1924 %      channel_statistics=GetImageStatistics(image,exception);
1925 %      red_mean=channel_statistics[RedPixelChannel].mean;
1926 %
1927 %  Use MagickRelinquishMemory() to free the statistics buffer.
1928 %
1929 %  The format of the GetImageStatistics method is:
1930 %
1931 %      ChannelStatistics *GetImageStatistics(const Image *image,
1932 %        ExceptionInfo *exception)
1933 %
1934 %  A description of each parameter follows:
1935 %
1936 %    o image: the image.
1937 %
1938 %    o exception: return any errors or warnings in this structure.
1939 %
1940 */
1941 MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
1942   ExceptionInfo *exception)
1943 {
1944   ChannelStatistics
1945     *channel_statistics;
1946
1947   double
1948     area,
1949     *histogram,
1950     standard_deviation;
1951
1952   MagickStatusType
1953     status;
1954
1955   QuantumAny
1956     range;
1957
1958   register ssize_t
1959     i;
1960
1961   size_t
1962     depth;
1963
1964   ssize_t
1965     y;
1966
1967   assert(image != (Image *) NULL);
1968   assert(image->signature == MagickCoreSignature);
1969   if (image->debug != MagickFalse)
1970     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1971   histogram=(double *) AcquireQuantumMemory(MaxMap+1UL,GetPixelChannels(image)*
1972     sizeof(*histogram));
1973   channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(
1974     MaxPixelChannels+1,sizeof(*channel_statistics));
1975   if ((channel_statistics == (ChannelStatistics *) NULL) ||
1976       (histogram == (double *) NULL))
1977     {
1978       if (histogram != (double *) NULL)
1979         histogram=(double *) RelinquishMagickMemory(histogram);
1980       if (channel_statistics != (ChannelStatistics *) NULL)
1981         channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1982           channel_statistics);
1983       return(channel_statistics);
1984     }
1985   (void) memset(channel_statistics,0,(MaxPixelChannels+1)*
1986     sizeof(*channel_statistics));
1987   for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
1988   {
1989     channel_statistics[i].depth=1;
1990     channel_statistics[i].maxima=(-MagickMaximumValue);
1991     channel_statistics[i].minima=MagickMaximumValue;
1992   }
1993   (void) memset(histogram,0,(MaxMap+1)*GetPixelChannels(image)*
1994     sizeof(*histogram));
1995   for (y=0; y < (ssize_t) image->rows; y++)
1996   {
1997     register const Quantum
1998       *magick_restrict p;
1999
2000     register ssize_t
2001       x;
2002
2003     /*
2004       Compute pixel statistics.
2005     */
2006     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
2007     if (p == (const Quantum *) NULL)
2008       break;
2009     for (x=0; x < (ssize_t) image->columns; x++)
2010     {
2011       register ssize_t
2012         i;
2013
2014       if (GetPixelReadMask(image,p) <= (QuantumRange/2))
2015         {
2016           p+=GetPixelChannels(image);
2017           continue;
2018         }
2019       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2020       {
2021         PixelChannel channel = GetPixelChannelChannel(image,i);
2022         PixelTrait traits = GetPixelChannelTraits(image,channel);
2023         if (traits == UndefinedPixelTrait)
2024           continue;
2025         if ((traits & UpdatePixelTrait) == 0)
2026           continue;
2027         if (channel_statistics[channel].depth != MAGICKCORE_QUANTUM_DEPTH)
2028           {
2029             depth=channel_statistics[channel].depth;
2030             range=GetQuantumRange(depth);
2031             status=p[i] != ScaleAnyToQuantum(ScaleQuantumToAny(p[i],range),
2032               range) ? MagickTrue : MagickFalse;
2033             if (status != MagickFalse)
2034               {
2035                 channel_statistics[channel].depth++;
2036                 i--;
2037                 continue;
2038               }
2039           }
2040         if ((double) p[i] < channel_statistics[channel].minima)
2041           channel_statistics[channel].minima=(double) p[i];
2042         if ((double) p[i] > channel_statistics[channel].maxima)
2043           channel_statistics[channel].maxima=(double) p[i];
2044         channel_statistics[channel].sum+=p[i];
2045         channel_statistics[channel].sum_squared+=(double) p[i]*p[i];
2046         channel_statistics[channel].sum_cubed+=(double) p[i]*p[i]*p[i];
2047         channel_statistics[channel].sum_fourth_power+=(double) p[i]*p[i]*p[i]*
2048           p[i];
2049         channel_statistics[channel].area++;
2050         if ((double) p[i] < channel_statistics[CompositePixelChannel].minima)
2051           channel_statistics[CompositePixelChannel].minima=(double) p[i];
2052         if ((double) p[i] > channel_statistics[CompositePixelChannel].maxima)
2053           channel_statistics[CompositePixelChannel].maxima=(double) p[i];
2054         histogram[GetPixelChannels(image)*ScaleQuantumToMap(
2055           ClampToQuantum((double) p[i]))+i]++;
2056         channel_statistics[CompositePixelChannel].sum+=(double) p[i];
2057         channel_statistics[CompositePixelChannel].sum_squared+=(double)
2058           p[i]*p[i];
2059         channel_statistics[CompositePixelChannel].sum_cubed+=(double)
2060           p[i]*p[i]*p[i];
2061         channel_statistics[CompositePixelChannel].sum_fourth_power+=(double)
2062           p[i]*p[i]*p[i]*p[i];
2063         channel_statistics[CompositePixelChannel].area++;
2064       }
2065       p+=GetPixelChannels(image);
2066     }
2067   }
2068   for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
2069   {
2070     /*
2071       Normalize pixel statistics.
2072     */
2073     area=PerceptibleReciprocal(channel_statistics[i].area);
2074     channel_statistics[i].sum*=area;
2075     channel_statistics[i].sum_squared*=area;
2076     channel_statistics[i].sum_cubed*=area;
2077     channel_statistics[i].sum_fourth_power*=area;
2078     channel_statistics[i].mean=channel_statistics[i].sum;
2079     channel_statistics[i].variance=channel_statistics[i].sum_squared;
2080     standard_deviation=sqrt(channel_statistics[i].variance-
2081       (channel_statistics[i].mean*channel_statistics[i].mean));
2082     standard_deviation=sqrt(PerceptibleReciprocal(channel_statistics[i].area-
2083       1.0)*channel_statistics[i].area*standard_deviation*standard_deviation);
2084     channel_statistics[i].standard_deviation=standard_deviation;
2085   }
2086   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2087   {
2088     double
2089       number_bins;
2090
2091     register ssize_t
2092       j;
2093
2094     /*
2095       Compute pixel entropy.
2096     */
2097     PixelChannel channel = GetPixelChannelChannel(image,i);
2098     number_bins=0.0;
2099     for (j=0; j <= (ssize_t) MaxMap; j++)
2100       if (histogram[GetPixelChannels(image)*j+i] > 0.0)
2101         number_bins++;
2102     area=PerceptibleReciprocal(channel_statistics[channel].area);
2103     for (j=0; j <= (ssize_t) MaxMap; j++)
2104     {
2105       double
2106         count;
2107
2108       count=area*histogram[GetPixelChannels(image)*j+i];
2109       channel_statistics[channel].entropy+=-count*MagickLog10(count)*
2110         PerceptibleReciprocal(MagickLog10(number_bins));
2111       channel_statistics[CompositePixelChannel].entropy+=-count*
2112         MagickLog10(count)*PerceptibleReciprocal(MagickLog10(number_bins))/
2113         GetPixelChannels(image);
2114     }
2115   }
2116   histogram=(double *) RelinquishMagickMemory(histogram);
2117   for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
2118   {
2119     /*
2120       Compute kurtosis & skewness statistics.
2121     */
2122     standard_deviation=PerceptibleReciprocal(
2123       channel_statistics[i].standard_deviation);
2124     channel_statistics[i].skewness=(channel_statistics[i].sum_cubed-3.0*
2125       channel_statistics[i].mean*channel_statistics[i].sum_squared+2.0*
2126       channel_statistics[i].mean*channel_statistics[i].mean*
2127       channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2128       standard_deviation);
2129     channel_statistics[i].kurtosis=(channel_statistics[i].sum_fourth_power-4.0*
2130       channel_statistics[i].mean*channel_statistics[i].sum_cubed+6.0*
2131       channel_statistics[i].mean*channel_statistics[i].mean*
2132       channel_statistics[i].sum_squared-3.0*channel_statistics[i].mean*
2133       channel_statistics[i].mean*1.0*channel_statistics[i].mean*
2134       channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2135       standard_deviation*standard_deviation)-3.0;
2136   }
2137   channel_statistics[CompositePixelChannel].mean=0.0;
2138   channel_statistics[CompositePixelChannel].standard_deviation=0.0;
2139   channel_statistics[CompositePixelChannel].entropy=0.0;
2140   for (i=0; i < (ssize_t) MaxPixelChannels; i++)
2141   {
2142     channel_statistics[CompositePixelChannel].mean+=
2143       channel_statistics[i].mean;
2144     channel_statistics[CompositePixelChannel].standard_deviation+=
2145       channel_statistics[i].standard_deviation;
2146     channel_statistics[CompositePixelChannel].entropy+=
2147       channel_statistics[i].entropy;
2148   }
2149   channel_statistics[CompositePixelChannel].mean/=(double)
2150     GetImageChannels(image);
2151   channel_statistics[CompositePixelChannel].standard_deviation/=(double)
2152     GetImageChannels(image);
2153   channel_statistics[CompositePixelChannel].entropy/=(double)
2154     GetImageChannels(image);
2155   if (y < (ssize_t) image->rows)
2156     channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
2157       channel_statistics);
2158   return(channel_statistics);
2159 }
2160 \f
2161 /*
2162 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2163 %                                                                             %
2164 %                                                                             %
2165 %                                                                             %
2166 %     P o l y n o m i a l I m a g e                                           %
2167 %                                                                             %
2168 %                                                                             %
2169 %                                                                             %
2170 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2171 %
2172 %  PolynomialImage() returns a new image where each pixel is the sum of the
2173 %  pixels in the image sequence after applying its corresponding terms
2174 %  (coefficient and degree pairs).
2175 %
2176 %  The format of the PolynomialImage method is:
2177 %
2178 %      Image *PolynomialImage(const Image *images,const size_t number_terms,
2179 %        const double *terms,ExceptionInfo *exception)
2180 %
2181 %  A description of each parameter follows:
2182 %
2183 %    o images: the image sequence.
2184 %
2185 %    o number_terms: the number of terms in the list.  The actual list length
2186 %      is 2 x number_terms + 1 (the constant).
2187 %
2188 %    o terms: the list of polynomial coefficients and degree pairs and a
2189 %      constant.
2190 %
2191 %    o exception: return any errors or warnings in this structure.
2192 %
2193 */
2194 MagickExport Image *PolynomialImage(const Image *images,
2195   const size_t number_terms,const double *terms,ExceptionInfo *exception)
2196 {
2197 #define PolynomialImageTag  "Polynomial/Image"
2198
2199   CacheView
2200     *polynomial_view;
2201
2202   Image
2203     *image;
2204
2205   MagickBooleanType
2206     status;
2207
2208   MagickOffsetType
2209     progress;
2210
2211   PixelChannels
2212     **magick_restrict polynomial_pixels;
2213
2214   size_t
2215     number_images;
2216
2217   ssize_t
2218     y;
2219
2220   assert(images != (Image *) NULL);
2221   assert(images->signature == MagickCoreSignature);
2222   if (images->debug != MagickFalse)
2223     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
2224   assert(exception != (ExceptionInfo *) NULL);
2225   assert(exception->signature == MagickCoreSignature);
2226   image=AcquireImageCanvas(images,exception);
2227   if (image == (Image *) NULL)
2228     return((Image *) NULL);
2229   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
2230     {
2231       image=DestroyImage(image);
2232       return((Image *) NULL);
2233     }
2234   number_images=GetImageListLength(images);
2235   polynomial_pixels=AcquirePixelThreadSet(images);
2236   if (polynomial_pixels == (PixelChannels **) NULL)
2237     {
2238       image=DestroyImage(image);
2239       (void) ThrowMagickException(exception,GetMagickModule(),
2240         ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
2241       return((Image *) NULL);
2242     }
2243   /*
2244     Polynomial image pixels.
2245   */
2246   status=MagickTrue;
2247   progress=0;
2248   polynomial_view=AcquireAuthenticCacheView(image,exception);
2249 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2250   #pragma omp parallel for schedule(static) shared(progress,status) \
2251     magick_number_threads(image,image,image->rows,1)
2252 #endif
2253   for (y=0; y < (ssize_t) image->rows; y++)
2254   {
2255     CacheView
2256       *image_view;
2257
2258     const Image
2259       *next;
2260
2261     const int
2262       id = GetOpenMPThreadId();
2263
2264     register ssize_t
2265       i,
2266       x;
2267
2268     register PixelChannels
2269       *polynomial_pixel;
2270
2271     register Quantum
2272       *magick_restrict q;
2273
2274     ssize_t
2275       j;
2276
2277     if (status == MagickFalse)
2278       continue;
2279     q=QueueCacheViewAuthenticPixels(polynomial_view,0,y,image->columns,1,
2280       exception);
2281     if (q == (Quantum *) NULL)
2282       {
2283         status=MagickFalse;
2284         continue;
2285       }
2286     polynomial_pixel=polynomial_pixels[id];
2287     for (j=0; j < (ssize_t) image->columns; j++)
2288       for (i=0; i < MaxPixelChannels; i++)
2289         polynomial_pixel[j].channel[i]=0.0;
2290     next=images;
2291     for (j=0; j < (ssize_t) number_images; j++)
2292     {
2293       register const Quantum
2294         *p;
2295
2296       if (j >= (ssize_t) number_terms)
2297         continue;
2298       image_view=AcquireVirtualCacheView(next,exception);
2299       p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2300       if (p == (const Quantum *) NULL)
2301         {
2302           image_view=DestroyCacheView(image_view);
2303           break;
2304         }
2305       for (x=0; x < (ssize_t) image->columns; x++)
2306       {
2307         register ssize_t
2308           i;
2309
2310         for (i=0; i < (ssize_t) GetPixelChannels(next); i++)
2311         {
2312           MagickRealType
2313             coefficient,
2314             degree;
2315
2316           PixelChannel channel = GetPixelChannelChannel(image,i);
2317           PixelTrait traits = GetPixelChannelTraits(next,channel);
2318           PixelTrait polynomial_traits=GetPixelChannelTraits(image,channel);
2319           if ((traits == UndefinedPixelTrait) ||
2320               (polynomial_traits == UndefinedPixelTrait))
2321             continue;
2322           if ((traits & UpdatePixelTrait) == 0)
2323             continue;
2324           coefficient=(MagickRealType) terms[2*j];
2325           degree=(MagickRealType) terms[(j << 1)+1];
2326           polynomial_pixel[x].channel[i]+=coefficient*
2327             pow(QuantumScale*GetPixelChannel(image,channel,p),degree);
2328         }
2329         p+=GetPixelChannels(next);
2330       }
2331       image_view=DestroyCacheView(image_view);
2332       next=GetNextImageInList(next);
2333     }
2334     for (x=0; x < (ssize_t) image->columns; x++)
2335     {
2336       register ssize_t
2337         i;
2338
2339       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2340       {
2341         PixelChannel channel = GetPixelChannelChannel(image,i);
2342         PixelTrait traits = GetPixelChannelTraits(image,channel);
2343         if (traits == UndefinedPixelTrait)
2344           continue;
2345         if ((traits & UpdatePixelTrait) == 0)
2346           continue;
2347         q[i]=ClampToQuantum(QuantumRange*polynomial_pixel[x].channel[i]);
2348       }
2349       q+=GetPixelChannels(image);
2350     }
2351     if (SyncCacheViewAuthenticPixels(polynomial_view,exception) == MagickFalse)
2352       status=MagickFalse;
2353     if (images->progress_monitor != (MagickProgressMonitor) NULL)
2354       {
2355         MagickBooleanType
2356           proceed;
2357
2358 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2359         #pragma omp atomic
2360 #endif
2361         progress++;
2362         proceed=SetImageProgress(images,PolynomialImageTag,progress,
2363           image->rows);
2364         if (proceed == MagickFalse)
2365           status=MagickFalse;
2366       }
2367   }
2368   polynomial_view=DestroyCacheView(polynomial_view);
2369   polynomial_pixels=DestroyPixelThreadSet(images,polynomial_pixels);
2370   if (status == MagickFalse)
2371     image=DestroyImage(image);
2372   return(image);
2373 }
2374 \f
2375 /*
2376 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2377 %                                                                             %
2378 %                                                                             %
2379 %                                                                             %
2380 %     S t a t i s t i c I m a g e                                             %
2381 %                                                                             %
2382 %                                                                             %
2383 %                                                                             %
2384 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2385 %
2386 %  StatisticImage() makes each pixel the min / max / median / mode / etc. of
2387 %  the neighborhood of the specified width and height.
2388 %
2389 %  The format of the StatisticImage method is:
2390 %
2391 %      Image *StatisticImage(const Image *image,const StatisticType type,
2392 %        const size_t width,const size_t height,ExceptionInfo *exception)
2393 %
2394 %  A description of each parameter follows:
2395 %
2396 %    o image: the image.
2397 %
2398 %    o type: the statistic type (median, mode, etc.).
2399 %
2400 %    o width: the width of the pixel neighborhood.
2401 %
2402 %    o height: the height of the pixel neighborhood.
2403 %
2404 %    o exception: return any errors or warnings in this structure.
2405 %
2406 */
2407
2408 typedef struct _SkipNode
2409 {
2410   size_t
2411     next[9],
2412     count,
2413     signature;
2414 } SkipNode;
2415
2416 typedef struct _SkipList
2417 {
2418   ssize_t
2419     level;
2420
2421   SkipNode
2422     *nodes;
2423 } SkipList;
2424
2425 typedef struct _PixelList
2426 {
2427   size_t
2428     length,
2429     seed;
2430
2431   SkipList
2432     skip_list;
2433
2434   size_t
2435     signature;
2436 } PixelList;
2437
2438 static PixelList *DestroyPixelList(PixelList *pixel_list)
2439 {
2440   if (pixel_list == (PixelList *) NULL)
2441     return((PixelList *) NULL);
2442   if (pixel_list->skip_list.nodes != (SkipNode *) NULL)
2443     pixel_list->skip_list.nodes=(SkipNode *) RelinquishAlignedMemory(
2444       pixel_list->skip_list.nodes);
2445   pixel_list=(PixelList *) RelinquishMagickMemory(pixel_list);
2446   return(pixel_list);
2447 }
2448
2449 static PixelList **DestroyPixelListThreadSet(PixelList **pixel_list)
2450 {
2451   register ssize_t
2452     i;
2453
2454   assert(pixel_list != (PixelList **) NULL);
2455   for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
2456     if (pixel_list[i] != (PixelList *) NULL)
2457       pixel_list[i]=DestroyPixelList(pixel_list[i]);
2458   pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list);
2459   return(pixel_list);
2460 }
2461
2462 static PixelList *AcquirePixelList(const size_t width,const size_t height)
2463 {
2464   PixelList
2465     *pixel_list;
2466
2467   pixel_list=(PixelList *) AcquireMagickMemory(sizeof(*pixel_list));
2468   if (pixel_list == (PixelList *) NULL)
2469     return(pixel_list);
2470   (void) memset((void *) pixel_list,0,sizeof(*pixel_list));
2471   pixel_list->length=width*height;
2472   pixel_list->skip_list.nodes=(SkipNode *) AcquireAlignedMemory(65537UL,
2473     sizeof(*pixel_list->skip_list.nodes));
2474   if (pixel_list->skip_list.nodes == (SkipNode *) NULL)
2475     return(DestroyPixelList(pixel_list));
2476   (void) memset(pixel_list->skip_list.nodes,0,65537UL*
2477     sizeof(*pixel_list->skip_list.nodes));
2478   pixel_list->signature=MagickCoreSignature;
2479   return(pixel_list);
2480 }
2481
2482 static PixelList **AcquirePixelListThreadSet(const size_t width,
2483   const size_t height)
2484 {
2485   PixelList
2486     **pixel_list;
2487
2488   register ssize_t
2489     i;
2490
2491   size_t
2492     number_threads;
2493
2494   number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
2495   pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
2496     sizeof(*pixel_list));
2497   if (pixel_list == (PixelList **) NULL)
2498     return((PixelList **) NULL);
2499   (void) memset(pixel_list,0,number_threads*sizeof(*pixel_list));
2500   for (i=0; i < (ssize_t) number_threads; i++)
2501   {
2502     pixel_list[i]=AcquirePixelList(width,height);
2503     if (pixel_list[i] == (PixelList *) NULL)
2504       return(DestroyPixelListThreadSet(pixel_list));
2505   }
2506   return(pixel_list);
2507 }
2508
2509 static void AddNodePixelList(PixelList *pixel_list,const size_t color)
2510 {
2511   register SkipList
2512     *p;
2513
2514   register ssize_t
2515     level;
2516
2517   size_t
2518     search,
2519     update[9];
2520
2521   /*
2522     Initialize the node.
2523   */
2524   p=(&pixel_list->skip_list);
2525   p->nodes[color].signature=pixel_list->signature;
2526   p->nodes[color].count=1;
2527   /*
2528     Determine where it belongs in the list.
2529   */
2530   search=65536UL;
2531   for (level=p->level; level >= 0; level--)
2532   {
2533     while (p->nodes[search].next[level] < color)
2534       search=p->nodes[search].next[level];
2535     update[level]=search;
2536   }
2537   /*
2538     Generate a pseudo-random level for this node.
2539   */
2540   for (level=0; ; level++)
2541   {
2542     pixel_list->seed=(pixel_list->seed*42893621L)+1L;
2543     if ((pixel_list->seed & 0x300) != 0x300)
2544       break;
2545   }
2546   if (level > 8)
2547     level=8;
2548   if (level > (p->level+2))
2549     level=p->level+2;
2550   /*
2551     If we're raising the list's level, link back to the root node.
2552   */
2553   while (level > p->level)
2554   {
2555     p->level++;
2556     update[p->level]=65536UL;
2557   }
2558   /*
2559     Link the node into the skip-list.
2560   */
2561   do
2562   {
2563     p->nodes[color].next[level]=p->nodes[update[level]].next[level];
2564     p->nodes[update[level]].next[level]=color;
2565   } while (level-- > 0);
2566 }
2567
2568 static inline void GetMaximumPixelList(PixelList *pixel_list,Quantum *pixel)
2569 {
2570   register SkipList
2571     *p;
2572
2573   size_t
2574     color,
2575     maximum;
2576
2577   ssize_t
2578     count;
2579
2580   /*
2581     Find the maximum value for each of the color.
2582   */
2583   p=(&pixel_list->skip_list);
2584   color=65536L;
2585   count=0;
2586   maximum=p->nodes[color].next[0];
2587   do
2588   {
2589     color=p->nodes[color].next[0];
2590     if (color > maximum)
2591       maximum=color;
2592     count+=p->nodes[color].count;
2593   } while (count < (ssize_t) pixel_list->length);
2594   *pixel=ScaleShortToQuantum((unsigned short) maximum);
2595 }
2596
2597 static inline void GetMeanPixelList(PixelList *pixel_list,Quantum *pixel)
2598 {
2599   double
2600     sum;
2601
2602   register SkipList
2603     *p;
2604
2605   size_t
2606     color;
2607
2608   ssize_t
2609     count;
2610
2611   /*
2612     Find the mean value for each of the color.
2613   */
2614   p=(&pixel_list->skip_list);
2615   color=65536L;
2616   count=0;
2617   sum=0.0;
2618   do
2619   {
2620     color=p->nodes[color].next[0];
2621     sum+=(double) p->nodes[color].count*color;
2622     count+=p->nodes[color].count;
2623   } while (count < (ssize_t) pixel_list->length);
2624   sum/=pixel_list->length;
2625   *pixel=ScaleShortToQuantum((unsigned short) sum);
2626 }
2627
2628 static inline void GetMedianPixelList(PixelList *pixel_list,Quantum *pixel)
2629 {
2630   register SkipList
2631     *p;
2632
2633   size_t
2634     color;
2635
2636   ssize_t
2637     count;
2638
2639   /*
2640     Find the median value for each of the color.
2641   */
2642   p=(&pixel_list->skip_list);
2643   color=65536L;
2644   count=0;
2645   do
2646   {
2647     color=p->nodes[color].next[0];
2648     count+=p->nodes[color].count;
2649   } while (count <= (ssize_t) (pixel_list->length >> 1));
2650   *pixel=ScaleShortToQuantum((unsigned short) color);
2651 }
2652
2653 static inline void GetMinimumPixelList(PixelList *pixel_list,Quantum *pixel)
2654 {
2655   register SkipList
2656     *p;
2657
2658   size_t
2659     color,
2660     minimum;
2661
2662   ssize_t
2663     count;
2664
2665   /*
2666     Find the minimum value for each of the color.
2667   */
2668   p=(&pixel_list->skip_list);
2669   count=0;
2670   color=65536UL;
2671   minimum=p->nodes[color].next[0];
2672   do
2673   {
2674     color=p->nodes[color].next[0];
2675     if (color < minimum)
2676       minimum=color;
2677     count+=p->nodes[color].count;
2678   } while (count < (ssize_t) pixel_list->length);
2679   *pixel=ScaleShortToQuantum((unsigned short) minimum);
2680 }
2681
2682 static inline void GetModePixelList(PixelList *pixel_list,Quantum *pixel)
2683 {
2684   register SkipList
2685     *p;
2686
2687   size_t
2688     color,
2689     max_count,
2690     mode;
2691
2692   ssize_t
2693     count;
2694
2695   /*
2696     Make each pixel the 'predominant color' of the specified neighborhood.
2697   */
2698   p=(&pixel_list->skip_list);
2699   color=65536L;
2700   mode=color;
2701   max_count=p->nodes[mode].count;
2702   count=0;
2703   do
2704   {
2705     color=p->nodes[color].next[0];
2706     if (p->nodes[color].count > max_count)
2707       {
2708         mode=color;
2709         max_count=p->nodes[mode].count;
2710       }
2711     count+=p->nodes[color].count;
2712   } while (count < (ssize_t) pixel_list->length);
2713   *pixel=ScaleShortToQuantum((unsigned short) mode);
2714 }
2715
2716 static inline void GetNonpeakPixelList(PixelList *pixel_list,Quantum *pixel)
2717 {
2718   register SkipList
2719     *p;
2720
2721   size_t
2722     color,
2723     next,
2724     previous;
2725
2726   ssize_t
2727     count;
2728
2729   /*
2730     Finds the non peak value for each of the colors.
2731   */
2732   p=(&pixel_list->skip_list);
2733   color=65536L;
2734   next=p->nodes[color].next[0];
2735   count=0;
2736   do
2737   {
2738     previous=color;
2739     color=next;
2740     next=p->nodes[color].next[0];
2741     count+=p->nodes[color].count;
2742   } while (count <= (ssize_t) (pixel_list->length >> 1));
2743   if ((previous == 65536UL) && (next != 65536UL))
2744     color=next;
2745   else
2746     if ((previous != 65536UL) && (next == 65536UL))
2747       color=previous;
2748   *pixel=ScaleShortToQuantum((unsigned short) color);
2749 }
2750
2751 static inline void GetRootMeanSquarePixelList(PixelList *pixel_list,
2752   Quantum *pixel)
2753 {
2754   double
2755     sum;
2756
2757   register SkipList
2758     *p;
2759
2760   size_t
2761     color;
2762
2763   ssize_t
2764     count;
2765
2766   /*
2767     Find the root mean square value for each of the color.
2768   */
2769   p=(&pixel_list->skip_list);
2770   color=65536L;
2771   count=0;
2772   sum=0.0;
2773   do
2774   {
2775     color=p->nodes[color].next[0];
2776     sum+=(double) (p->nodes[color].count*color*color);
2777     count+=p->nodes[color].count;
2778   } while (count < (ssize_t) pixel_list->length);
2779   sum/=pixel_list->length;
2780   *pixel=ScaleShortToQuantum((unsigned short) sqrt(sum));
2781 }
2782
2783 static inline void GetStandardDeviationPixelList(PixelList *pixel_list,
2784   Quantum *pixel)
2785 {
2786   double
2787     sum,
2788     sum_squared;
2789
2790   register SkipList
2791     *p;
2792
2793   size_t
2794     color;
2795
2796   ssize_t
2797     count;
2798
2799   /*
2800     Find the standard-deviation value for each of the color.
2801   */
2802   p=(&pixel_list->skip_list);
2803   color=65536L;
2804   count=0;
2805   sum=0.0;
2806   sum_squared=0.0;
2807   do
2808   {
2809     register ssize_t
2810       i;
2811
2812     color=p->nodes[color].next[0];
2813     sum+=(double) p->nodes[color].count*color;
2814     for (i=0; i < (ssize_t) p->nodes[color].count; i++)
2815       sum_squared+=((double) color)*((double) color);
2816     count+=p->nodes[color].count;
2817   } while (count < (ssize_t) pixel_list->length);
2818   sum/=pixel_list->length;
2819   sum_squared/=pixel_list->length;
2820   *pixel=ScaleShortToQuantum((unsigned short) sqrt(sum_squared-(sum*sum)));
2821 }
2822
2823 static inline void InsertPixelList(const Quantum pixel,PixelList *pixel_list)
2824 {
2825   size_t
2826     signature;
2827
2828   unsigned short
2829     index;
2830
2831   index=ScaleQuantumToShort(pixel);
2832   signature=pixel_list->skip_list.nodes[index].signature;
2833   if (signature == pixel_list->signature)
2834     {
2835       pixel_list->skip_list.nodes[index].count++;
2836       return;
2837     }
2838   AddNodePixelList(pixel_list,index);
2839 }
2840
2841 static void ResetPixelList(PixelList *pixel_list)
2842 {
2843   int
2844     level;
2845
2846   register SkipNode
2847     *root;
2848
2849   register SkipList
2850     *p;
2851
2852   /*
2853     Reset the skip-list.
2854   */
2855   p=(&pixel_list->skip_list);
2856   root=p->nodes+65536UL;
2857   p->level=0;
2858   for (level=0; level < 9; level++)
2859     root->next[level]=65536UL;
2860   pixel_list->seed=pixel_list->signature++;
2861 }
2862
2863 MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
2864   const size_t width,const size_t height,ExceptionInfo *exception)
2865 {
2866 #define StatisticImageTag  "Statistic/Image"
2867
2868   CacheView
2869     *image_view,
2870     *statistic_view;
2871
2872   Image
2873     *statistic_image;
2874
2875   MagickBooleanType
2876     status;
2877
2878   MagickOffsetType
2879     progress;
2880
2881   PixelList
2882     **magick_restrict pixel_list;
2883
2884   ssize_t
2885     center,
2886     y;
2887
2888   /*
2889     Initialize statistics image attributes.
2890   */
2891   assert(image != (Image *) NULL);
2892   assert(image->signature == MagickCoreSignature);
2893   if (image->debug != MagickFalse)
2894     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2895   assert(exception != (ExceptionInfo *) NULL);
2896   assert(exception->signature == MagickCoreSignature);
2897   statistic_image=CloneImage(image,0,0,MagickTrue,
2898     exception);
2899   if (statistic_image == (Image *) NULL)
2900     return((Image *) NULL);
2901   status=SetImageStorageClass(statistic_image,DirectClass,exception);
2902   if (status == MagickFalse)
2903     {
2904       statistic_image=DestroyImage(statistic_image);
2905       return((Image *) NULL);
2906     }
2907   pixel_list=AcquirePixelListThreadSet(MagickMax(width,1),MagickMax(height,1));
2908   if (pixel_list == (PixelList **) NULL)
2909     {
2910       statistic_image=DestroyImage(statistic_image);
2911       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2912     }
2913   /*
2914     Make each pixel the min / max / median / mode / etc. of the neighborhood.
2915   */
2916   center=(ssize_t) GetPixelChannels(image)*(image->columns+MagickMax(width,1))*
2917     (MagickMax(height,1)/2L)+GetPixelChannels(image)*(MagickMax(width,1)/2L);
2918   status=MagickTrue;
2919   progress=0;
2920   image_view=AcquireVirtualCacheView(image,exception);
2921   statistic_view=AcquireAuthenticCacheView(statistic_image,exception);
2922 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2923   #pragma omp parallel for schedule(static) shared(progress,status) \
2924     magick_number_threads(image,statistic_image,statistic_image->rows,1)
2925 #endif
2926   for (y=0; y < (ssize_t) statistic_image->rows; y++)
2927   {
2928     const int
2929       id = GetOpenMPThreadId();
2930
2931     register const Quantum
2932       *magick_restrict p;
2933
2934     register Quantum
2935       *magick_restrict q;
2936
2937     register ssize_t
2938       x;
2939
2940     if (status == MagickFalse)
2941       continue;
2942     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) MagickMax(width,1)/2L),y-
2943       (ssize_t) (MagickMax(height,1)/2L),image->columns+MagickMax(width,1),
2944       MagickMax(height,1),exception);
2945     q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns,      1,exception);
2946     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2947       {
2948         status=MagickFalse;
2949         continue;
2950       }
2951     for (x=0; x < (ssize_t) statistic_image->columns; x++)
2952     {
2953       register ssize_t
2954         i;
2955
2956       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2957       {
2958         Quantum
2959           pixel;
2960
2961         register const Quantum
2962           *magick_restrict pixels;
2963
2964         register ssize_t
2965           u;
2966
2967         ssize_t
2968           v;
2969
2970         PixelChannel channel = GetPixelChannelChannel(image,i);
2971         PixelTrait traits = GetPixelChannelTraits(image,channel);
2972         PixelTrait statistic_traits=GetPixelChannelTraits(statistic_image,
2973           channel);
2974         if ((traits == UndefinedPixelTrait) ||
2975             (statistic_traits == UndefinedPixelTrait))
2976           continue;
2977         if (((statistic_traits & CopyPixelTrait) != 0) ||
2978             (GetPixelWriteMask(image,p) <= (QuantumRange/2)))
2979           {
2980             SetPixelChannel(statistic_image,channel,p[center+i],q);
2981             continue;
2982           }
2983         if ((statistic_traits & UpdatePixelTrait) == 0)
2984           continue;
2985         pixels=p;
2986         ResetPixelList(pixel_list[id]);
2987         for (v=0; v < (ssize_t) MagickMax(height,1); v++)
2988         {
2989           for (u=0; u < (ssize_t) MagickMax(width,1); u++)
2990           {
2991             InsertPixelList(pixels[i],pixel_list[id]);
2992             pixels+=GetPixelChannels(image);
2993           }
2994           pixels+=GetPixelChannels(image)*image->columns;
2995         }
2996         switch (type)
2997         {
2998           case GradientStatistic:
2999           {
3000             double
3001               maximum,
3002               minimum;
3003
3004             GetMinimumPixelList(pixel_list[id],&pixel);
3005             minimum=(double) pixel;
3006             GetMaximumPixelList(pixel_list[id],&pixel);
3007             maximum=(double) pixel;
3008             pixel=ClampToQuantum(MagickAbsoluteValue(maximum-minimum));
3009             break;
3010           }
3011           case MaximumStatistic:
3012           {
3013             GetMaximumPixelList(pixel_list[id],&pixel);
3014             break;
3015           }
3016           case MeanStatistic:
3017           {
3018             GetMeanPixelList(pixel_list[id],&pixel);
3019             break;
3020           }
3021           case MedianStatistic:
3022           default:
3023           {
3024             GetMedianPixelList(pixel_list[id],&pixel);
3025             break;
3026           }
3027           case MinimumStatistic:
3028           {
3029             GetMinimumPixelList(pixel_list[id],&pixel);
3030             break;
3031           }
3032           case ModeStatistic:
3033           {
3034             GetModePixelList(pixel_list[id],&pixel);
3035             break;
3036           }
3037           case NonpeakStatistic:
3038           {
3039             GetNonpeakPixelList(pixel_list[id],&pixel);
3040             break;
3041           }
3042           case RootMeanSquareStatistic:
3043           {
3044             GetRootMeanSquarePixelList(pixel_list[id],&pixel);
3045             break;
3046           }
3047           case StandardDeviationStatistic:
3048           {
3049             GetStandardDeviationPixelList(pixel_list[id],&pixel);
3050             break;
3051           }
3052         }
3053         SetPixelChannel(statistic_image,channel,pixel,q);
3054       }
3055       p+=GetPixelChannels(image);
3056       q+=GetPixelChannels(statistic_image);
3057     }
3058     if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
3059       status=MagickFalse;
3060     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3061       {
3062         MagickBooleanType
3063           proceed;
3064
3065 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3066         #pragma omp atomic
3067 #endif
3068         progress++;
3069         proceed=SetImageProgress(image,StatisticImageTag,progress,image->rows);
3070         if (proceed == MagickFalse)
3071           status=MagickFalse;
3072       }
3073   }
3074   statistic_view=DestroyCacheView(statistic_view);
3075   image_view=DestroyCacheView(image_view);
3076   pixel_list=DestroyPixelListThreadSet(pixel_list);
3077   if (status == MagickFalse)
3078     statistic_image=DestroyImage(statistic_image);
3079   return(statistic_image);
3080 }