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