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-2018 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) memset(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 static Image *AcquireImageCanvas(const Image *images,ExceptionInfo *exception)
 420 {
 421   const Image
 422     *p,
 423     *q;
 424
 425   size_t
 426     columns,
 427     rows;
 428
 429   q=images;
 430   columns=images->columns;
 431   rows=images->rows;
 432   for (p=images; p != (Image *) NULL; p=p->next)
 433   {
 434     if (p->number_channels > q->number_channels)
 435       q=p;
 436     if (p->columns > columns)
 437       columns=p->columns;
 438     if (p->rows > rows)
 439       rows=p->rows;
 440   }
 441   return(CloneImage(q,columns,rows,MagickTrue,exception));
 442 }
 443
 444 MagickExport Image *EvaluateImages(const Image *images,
 445   const MagickEvaluateOperator op,ExceptionInfo *exception)
 446 {
 447 #define EvaluateImageTag  "Evaluate/Image"
 448
 449   CacheView
 450     *evaluate_view;
 451
 452   Image
 453     *image;
 454
 455   MagickBooleanType
 456     status;
 457
 458   MagickOffsetType
 459     progress;
 460
 461   PixelChannels
 462     **magick_restrict evaluate_pixels;
 463
 464   RandomInfo
 465     **magick_restrict random_info;
 466
 467   size_t
 468     number_images;
 469
 470   ssize_t
 471     y;
 472
 473 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 474   unsigned long
 475     key;
 476 #endif
 477
 478   assert(images != (Image *) NULL);
 479   assert(images->signature == MagickCoreSignature);
 480   if (images->debug != MagickFalse)
 481     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
 482   assert(exception != (ExceptionInfo *) NULL);
 483   assert(exception->signature == MagickCoreSignature);
 484   image=AcquireImageCanvas(images,exception);
 485   if (image == (Image *) NULL)
 486     return((Image *) NULL);
 487   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
 488     {
 489       image=DestroyImage(image);
 490       return((Image *) NULL);
 491     }
 492   number_images=GetImageListLength(images);
 493   evaluate_pixels=AcquirePixelThreadSet(images);
 494   if (evaluate_pixels == (PixelChannels **) NULL)
 495     {
 496       image=DestroyImage(image);
 497       (void) ThrowMagickException(exception,GetMagickModule(),
 498         ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
 499       return((Image *) NULL);
 500     }
 501   /*
 502     Evaluate image pixels.
 503   */
 504   status=MagickTrue;
 505   progress=0;
 506   random_info=AcquireRandomInfoThreadSet();
 507   evaluate_view=AcquireAuthenticCacheView(image,exception);
 508   if (op == MedianEvaluateOperator)
 509     {
 510 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 511       key=GetRandomSecretKey(random_info[0]);
 512       #pragma omp parallel for schedule(static) shared(progress,status) \
 513         magick_number_threads(image,images,image->rows,key == ~0UL)
 514 #endif
 515       for (y=0; y < (ssize_t) image->rows; y++)
 516       {
 517         CacheView
 518           *image_view;
 519
 520         const Image
 521           *next;
 522
 523         const int
 524           id = GetOpenMPThreadId();
 525
 526         register PixelChannels
 527           *evaluate_pixel;
 528
 529         register Quantum
 530           *magick_restrict q;
 531
 532         register ssize_t
 533           x;
 534
 535         if (status == MagickFalse)
 536           continue;
 537         q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,image->columns,1,
 538           exception);
 539         if (q == (Quantum *) NULL)
 540           {
 541             status=MagickFalse;
 542             continue;
 543           }
 544         evaluate_pixel=evaluate_pixels[id];
 545         for (x=0; x < (ssize_t) image->columns; x++)
 546         {
 547           register ssize_t
 548             j,
 549             k;
 550
 551           for (j=0; j < (ssize_t) number_images; j++)
 552             for (k=0; k < MaxPixelChannels; k++)
 553               evaluate_pixel[j].channel[k]=0.0;
 554           next=images;
 555           for (j=0; j < (ssize_t) number_images; j++)
 556           {
 557             register const Quantum
 558               *p;
 559
 560             register ssize_t
 561               i;
 562
 563             image_view=AcquireVirtualCacheView(next,exception);
 564             p=GetCacheViewVirtualPixels(image_view,x,y,1,1,exception);
 565             if (p == (const Quantum *) NULL)
 566               {
 567                 image_view=DestroyCacheView(image_view);
 568                 break;
 569               }
 570             for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 571             {
 572               PixelChannel channel = GetPixelChannelChannel(image,i);
 573               PixelTrait evaluate_traits=GetPixelChannelTraits(image,channel);
 574               PixelTrait traits = GetPixelChannelTraits(next,channel);
 575               if ((traits == UndefinedPixelTrait) ||
 576                   (evaluate_traits == UndefinedPixelTrait))
 577                 continue;
 578               if ((traits & UpdatePixelTrait) == 0)
 579                 continue;
 580               evaluate_pixel[j].channel[i]=ApplyEvaluateOperator(
 581                 random_info[id],GetPixelChannel(image,channel,p),op,
 582                 evaluate_pixel[j].channel[i]);
 583             }
 584             image_view=DestroyCacheView(image_view);
 585             next=GetNextImageInList(next);
 586           }
 587           qsort((void *) evaluate_pixel,number_images,sizeof(*evaluate_pixel),
 588             IntensityCompare);
 589           for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
 590             q[k]=ClampToQuantum(evaluate_pixel[j/2].channel[k]);
 591           q+=GetPixelChannels(image);
 592         }
 593         if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
 594           status=MagickFalse;
 595         if (images->progress_monitor != (MagickProgressMonitor) NULL)
 596           {
 597             MagickBooleanType
 598               proceed;
 599
 600 #if   defined(MAGICKCORE_OPENMP_SUPPORT)
 601             #pragma omp critical (MagickCore_EvaluateImages)
 602 #endif
 603             proceed=SetImageProgress(images,EvaluateImageTag,progress++,
 604               image->rows);
 605             if (proceed == MagickFalse)
 606               status=MagickFalse;
 607           }
 608       }
 609     }
 610   else
 611     {
 612 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 613       key=GetRandomSecretKey(random_info[0]);
 614       #pragma omp parallel for schedule(static) shared(progress,status) \
 615         magick_number_threads(image,images,image->rows,key == ~0UL)
 616 #endif
 617       for (y=0; y < (ssize_t) image->rows; y++)
 618       {
 619         CacheView
 620           *image_view;
 621
 622         const Image
 623           *next;
 624
 625         const int
 626           id = GetOpenMPThreadId();
 627
 628         register ssize_t
 629           i,
 630           x;
 631
 632         register PixelChannels
 633           *evaluate_pixel;
 634
 635         register Quantum
 636           *magick_restrict q;
 637
 638         ssize_t
 639           j;
 640
 641         if (status == MagickFalse)
 642           continue;
 643         q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,image->columns,1,
 644           exception);
 645         if (q == (Quantum *) NULL)
 646           {
 647             status=MagickFalse;
 648             continue;
 649           }
 650         evaluate_pixel=evaluate_pixels[id];
 651         for (j=0; j < (ssize_t) image->columns; j++)
 652           for (i=0; i < MaxPixelChannels; i++)
 653             evaluate_pixel[j].channel[i]=0.0;
 654         next=images;
 655         for (j=0; j < (ssize_t) number_images; j++)
 656         {
 657           register const Quantum
 658             *p;
 659
 660           image_view=AcquireVirtualCacheView(next,exception);
 661           p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,
 662             exception);
 663           if (p == (const Quantum *) NULL)
 664             {
 665               image_view=DestroyCacheView(image_view);
 666               break;
 667             }
 668           for (x=0; x < (ssize_t) image->columns; x++)
 669           {
 670             register ssize_t
 671               i;
 672
 673             for (i=0; i < (ssize_t) GetPixelChannels(next); i++)
 674             {
 675               PixelChannel channel = GetPixelChannelChannel(image,i);
 676               PixelTrait traits = GetPixelChannelTraits(next,channel);
 677               PixelTrait evaluate_traits=GetPixelChannelTraits(image,channel);
 678               if ((traits == UndefinedPixelTrait) ||
 679                   (evaluate_traits == UndefinedPixelTrait))
 680                 continue;
 681               if ((traits & UpdatePixelTrait) == 0)
 682                 continue;
 683               evaluate_pixel[x].channel[i]=ApplyEvaluateOperator(
 684                 random_info[id],GetPixelChannel(image,channel,p),j == 0 ?
 685                 AddEvaluateOperator : op,evaluate_pixel[x].channel[i]);
 686             }
 687             p+=GetPixelChannels(next);
 688           }
 689           image_view=DestroyCacheView(image_view);
 690           next=GetNextImageInList(next);
 691         }
 692         for (x=0; x < (ssize_t) image->columns; x++)
 693         {
 694           register ssize_t
 695              i;
 696
 697           switch (op)
 698           {
 699             case MeanEvaluateOperator:
 700             {
 701               for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 702                 evaluate_pixel[x].channel[i]/=(double) number_images;
 703               break;
 704             }
 705             case MultiplyEvaluateOperator:
 706             {
 707               for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 708               {
 709                 register ssize_t
 710                   j;
 711
 712                 for (j=0; j < (ssize_t) (number_images-1); j++)
 713                   evaluate_pixel[x].channel[i]*=QuantumScale;
 714               }
 715               break;
 716             }
 717             case RootMeanSquareEvaluateOperator:
 718             {
 719               for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 720                 evaluate_pixel[x].channel[i]=sqrt(evaluate_pixel[x].channel[i]/
 721                   number_images);
 722               break;
 723             }
 724             default:
 725               break;
 726           }
 727         }
 728         for (x=0; x < (ssize_t) image->columns; x++)
 729         {
 730           register ssize_t
 731             i;
 732
 733           for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 734           {
 735             PixelChannel channel = GetPixelChannelChannel(image,i);
 736             PixelTrait traits = GetPixelChannelTraits(image,channel);
 737             if (traits == UndefinedPixelTrait)
 738               continue;
 739             if ((traits & UpdatePixelTrait) == 0)
 740               continue;
 741             q[i]=ClampToQuantum(evaluate_pixel[x].channel[i]);
 742           }
 743           q+=GetPixelChannels(image);
 744         }
 745         if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
 746           status=MagickFalse;
 747         if (images->progress_monitor != (MagickProgressMonitor) NULL)
 748           {
 749             MagickBooleanType
 750               proceed;
 751
 752 #if   defined(MAGICKCORE_OPENMP_SUPPORT)
 753             #pragma omp critical (MagickCore_EvaluateImages)
 754 #endif
 755             proceed=SetImageProgress(images,EvaluateImageTag,progress++,
 756               image->rows);
 757             if (proceed == MagickFalse)
 758               status=MagickFalse;
 759           }
 760       }
 761     }
 762   evaluate_view=DestroyCacheView(evaluate_view);
 763   evaluate_pixels=DestroyPixelThreadSet(evaluate_pixels);
 764   random_info=DestroyRandomInfoThreadSet(random_info);
 765   if (status == MagickFalse)
 766     image=DestroyImage(image);
 767   return(image);
 768 }
 769
 770 MagickExport MagickBooleanType EvaluateImage(Image *image,
 771   const MagickEvaluateOperator op,const double value,ExceptionInfo *exception)
 772 {
 773   CacheView
 774     *image_view;
 775
 776   MagickBooleanType
 777     status;
 778
 779   MagickOffsetType
 780     progress;
 781
 782   RandomInfo
 783     **magick_restrict random_info;
 784
 785   ssize_t
 786     y;
 787
 788 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 789   unsigned long
 790     key;
 791 #endif
 792
 793   assert(image != (Image *) NULL);
 794   assert(image->signature == MagickCoreSignature);
 795   if (image->debug != MagickFalse)
 796     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 797   assert(exception != (ExceptionInfo *) NULL);
 798   assert(exception->signature == MagickCoreSignature);
 799   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
 800     return(MagickFalse);
 801   status=MagickTrue;
 802   progress=0;
 803   random_info=AcquireRandomInfoThreadSet();
 804   image_view=AcquireAuthenticCacheView(image,exception);
 805 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 806   key=GetRandomSecretKey(random_info[0]);
 807   #pragma omp parallel for schedule(static) shared(progress,status) \
 808     magick_number_threads(image,image,image->rows,key == ~0UL)
 809 #endif
 810   for (y=0; y < (ssize_t) image->rows; y++)
 811   {
 812     const int
 813       id = GetOpenMPThreadId();
 814
 815     register Quantum
 816       *magick_restrict q;
 817
 818     register ssize_t
 819       x;
 820
 821     if (status == MagickFalse)
 822       continue;
 823     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
 824     if (q == (Quantum *) NULL)
 825       {
 826         status=MagickFalse;
 827         continue;
 828       }
 829     for (x=0; x < (ssize_t) image->columns; x++)
 830     {
 831       double
 832         result;
 833
 834       register ssize_t
 835         i;
 836
 837       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 838       {
 839         PixelChannel channel = GetPixelChannelChannel(image,i);
 840         PixelTrait traits = GetPixelChannelTraits(image,channel);
 841         if (traits == UndefinedPixelTrait)
 842           continue;
 843         if ((traits & CopyPixelTrait) != 0)
 844           continue;
 845         if ((traits & UpdatePixelTrait) == 0)
 846           continue;
 847         result=ApplyEvaluateOperator(random_info[id],q[i],op,value);
 848         if (op == MeanEvaluateOperator)
 849           result/=2.0;
 850         q[i]=ClampToQuantum(result);
 851       }
 852       q+=GetPixelChannels(image);
 853     }
 854     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
 855       status=MagickFalse;
 856     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 857       {
 858         MagickBooleanType
 859           proceed;
 860
 861 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 862         #pragma omp critical (MagickCore_EvaluateImage)
 863 #endif
 864         proceed=SetImageProgress(image,EvaluateImageTag,progress++,image->rows);
 865         if (proceed == MagickFalse)
 866           status=MagickFalse;
 867       }
 868   }
 869   image_view=DestroyCacheView(image_view);
 870   random_info=DestroyRandomInfoThreadSet(random_info);
 871   return(status);
 872 }
 873 \f
 874 /*
 875 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 876 %                                                                             %
 877 %                                                                             %
 878 %                                                                             %
 879 %     F u n c t i o n I m a g e                                               %
 880 %                                                                             %
 881 %                                                                             %
 882 %                                                                             %
 883 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 884 %
 885 %  FunctionImage() applies a value to the image with an arithmetic, relational,
 886 %  or logical operator to an image. Use these operations to lighten or darken
 887 %  an image, to increase or decrease contrast in an image, or to produce the
 888 %  "negative" of an image.
 889 %
 890 %  The format of the FunctionImage method is:
 891 %
 892 %      MagickBooleanType FunctionImage(Image *image,
 893 %        const MagickFunction function,const ssize_t number_parameters,
 894 %        const double *parameters,ExceptionInfo *exception)
 895 %
 896 %  A description of each parameter follows:
 897 %
 898 %    o image: the image.
 899 %
 900 %    o function: A channel function.
 901 %
 902 %    o parameters: one or more parameters.
 903 %
 904 %    o exception: return any errors or warnings in this structure.
 905 %
 906 */
 907
 908 static Quantum ApplyFunction(Quantum pixel,const MagickFunction function,
 909   const size_t number_parameters,const double *parameters,
 910   ExceptionInfo *exception)
 911 {
 912   double
 913     result;
 914
 915   register ssize_t
 916     i;
 917
 918   (void) exception;
 919   result=0.0;
 920   switch (function)
 921   {
 922     case PolynomialFunction:
 923     {
 924       /*
 925         Polynomial: polynomial constants, highest to lowest order (e.g. c0*x^3+
 926         c1*x^2+c2*x+c3).
 927       */
 928       result=0.0;
 929       for (i=0; i < (ssize_t) number_parameters; i++)
 930         result=result*QuantumScale*pixel+parameters[i];
 931       result*=QuantumRange;
 932       break;
 933     }
 934     case SinusoidFunction:
 935     {
 936       double
 937         amplitude,
 938         bias,
 939         frequency,
 940         phase;
 941
 942       /*
 943         Sinusoid: frequency, phase, amplitude, bias.
 944       */
 945       frequency=(number_parameters >= 1) ? parameters[0] : 1.0;
 946       phase=(number_parameters >= 2) ? parameters[1] : 0.0;
 947       amplitude=(number_parameters >= 3) ? parameters[2] : 0.5;
 948       bias=(number_parameters >= 4) ? parameters[3] : 0.5;
 949       result=(double) (QuantumRange*(amplitude*sin((double) (2.0*
 950         MagickPI*(frequency*QuantumScale*pixel+phase/360.0)))+bias));
 951       break;
 952     }
 953     case ArcsinFunction:
 954     {
 955       double
 956         bias,
 957         center,
 958         range,
 959         width;
 960
 961       /*
 962         Arcsin (peged at range limits for invalid results): width, center,
 963         range, and bias.
 964       */
 965       width=(number_parameters >= 1) ? parameters[0] : 1.0;
 966       center=(number_parameters >= 2) ? parameters[1] : 0.5;
 967       range=(number_parameters >= 3) ? parameters[2] : 1.0;
 968       bias=(number_parameters >= 4) ? parameters[3] : 0.5;
 969       result=2.0/width*(QuantumScale*pixel-center);
 970       if ( result <= -1.0 )
 971         result=bias-range/2.0;
 972       else
 973         if (result >= 1.0)
 974           result=bias+range/2.0;
 975         else
 976           result=(double) (range/MagickPI*asin((double) result)+bias);
 977       result*=QuantumRange;
 978       break;
 979     }
 980     case ArctanFunction:
 981     {
 982       double
 983         center,
 984         bias,
 985         range,
 986         slope;
 987
 988       /*
 989         Arctan: slope, center, range, and bias.
 990       */
 991       slope=(number_parameters >= 1) ? parameters[0] : 1.0;
 992       center=(number_parameters >= 2) ? parameters[1] : 0.5;
 993       range=(number_parameters >= 3) ? parameters[2] : 1.0;
 994       bias=(number_parameters >= 4) ? parameters[3] : 0.5;
 995       result=(double) (MagickPI*slope*(QuantumScale*pixel-center));
 996       result=(double) (QuantumRange*(range/MagickPI*atan((double)
 997         result)+bias));
 998       break;
 999     }
1000     case UndefinedFunction:
1001       break;
1002   }
1003   return(ClampToQuantum(result));
1004 }
1005
1006 MagickExport MagickBooleanType FunctionImage(Image *image,
1007   const MagickFunction function,const size_t number_parameters,
1008   const double *parameters,ExceptionInfo *exception)
1009 {
1010 #define FunctionImageTag  "Function/Image "
1011
1012   CacheView
1013     *image_view;
1014
1015   MagickBooleanType
1016     status;
1017
1018   MagickOffsetType
1019     progress;
1020
1021   ssize_t
1022     y;
1023
1024   assert(image != (Image *) NULL);
1025   assert(image->signature == MagickCoreSignature);
1026   if (image->debug != MagickFalse)
1027     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1028   assert(exception != (ExceptionInfo *) NULL);
1029   assert(exception->signature == MagickCoreSignature);
1030 #if defined(MAGICKCORE_OPENCL_SUPPORT)
1031   if (AccelerateFunctionImage(image,function,number_parameters,parameters,
1032         exception) != MagickFalse)
1033     return(MagickTrue);
1034 #endif
1035   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
1036     return(MagickFalse);
1037   status=MagickTrue;
1038   progress=0;
1039   image_view=AcquireAuthenticCacheView(image,exception);
1040 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1041   #pragma omp parallel for schedule(static) shared(progress,status) \
1042     magick_number_threads(image,image,image->rows,1)
1043 #endif
1044   for (y=0; y < (ssize_t) image->rows; y++)
1045   {
1046     register Quantum
1047       *magick_restrict q;
1048
1049     register ssize_t
1050       x;
1051
1052     if (status == MagickFalse)
1053       continue;
1054     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
1055     if (q == (Quantum *) NULL)
1056       {
1057         status=MagickFalse;
1058         continue;
1059       }
1060     for (x=0; x < (ssize_t) image->columns; x++)
1061     {
1062       register ssize_t
1063         i;
1064
1065       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1066       {
1067         PixelChannel channel = GetPixelChannelChannel(image,i);
1068         PixelTrait traits = GetPixelChannelTraits(image,channel);
1069         if (traits == UndefinedPixelTrait)
1070           continue;
1071         if ((traits & UpdatePixelTrait) == 0)
1072           continue;
1073         q[i]=ApplyFunction(q[i],function,number_parameters,parameters,
1074           exception);
1075       }
1076       q+=GetPixelChannels(image);
1077     }
1078     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1079       status=MagickFalse;
1080     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1081       {
1082         MagickBooleanType
1083           proceed;
1084
1085 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1086         #pragma omp critical (MagickCore_FunctionImage)
1087 #endif
1088         proceed=SetImageProgress(image,FunctionImageTag,progress++,image->rows);
1089         if (proceed == MagickFalse)
1090           status=MagickFalse;
1091       }
1092   }
1093   image_view=DestroyCacheView(image_view);
1094   return(status);
1095 }
1096 \f
1097 /*
1098 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1099 %                                                                             %
1100 %                                                                             %
1101 %                                                                             %
1102 %   G e t I m a g e E n t r o p y                                             %
1103 %                                                                             %
1104 %                                                                             %
1105 %                                                                             %
1106 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1107 %
1108 %  GetImageEntropy() returns the entropy of one or more image channels.
1109 %
1110 %  The format of the GetImageEntropy method is:
1111 %
1112 %      MagickBooleanType GetImageEntropy(const Image *image,double *entropy,
1113 %        ExceptionInfo *exception)
1114 %
1115 %  A description of each parameter follows:
1116 %
1117 %    o image: the image.
1118 %
1119 %    o entropy: the average entropy of the selected channels.
1120 %
1121 %    o exception: return any errors or warnings in this structure.
1122 %
1123 */
1124 MagickExport MagickBooleanType GetImageEntropy(const Image *image,
1125   double *entropy,ExceptionInfo *exception)
1126 {
1127   ChannelStatistics
1128     *channel_statistics;
1129
1130   assert(image != (Image *) NULL);
1131   assert(image->signature == MagickCoreSignature);
1132   if (image->debug != MagickFalse)
1133     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1134   channel_statistics=GetImageStatistics(image,exception);
1135   if (channel_statistics == (ChannelStatistics *) NULL)
1136     return(MagickFalse);
1137   *entropy=channel_statistics[CompositePixelChannel].entropy;
1138   channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1139     channel_statistics);
1140   return(MagickTrue);
1141 }
1142 \f
1143 /*
1144 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1145 %                                                                             %
1146 %                                                                             %
1147 %                                                                             %
1148 %   G e t I m a g e E x t r e m a                                             %
1149 %                                                                             %
1150 %                                                                             %
1151 %                                                                             %
1152 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1153 %
1154 %  GetImageExtrema() returns the extrema of one or more image channels.
1155 %
1156 %  The format of the GetImageExtrema method is:
1157 %
1158 %      MagickBooleanType GetImageExtrema(const Image *image,size_t *minima,
1159 %        size_t *maxima,ExceptionInfo *exception)
1160 %
1161 %  A description of each parameter follows:
1162 %
1163 %    o image: the image.
1164 %
1165 %    o minima: the minimum value in the channel.
1166 %
1167 %    o maxima: the maximum value in the channel.
1168 %
1169 %    o exception: return any errors or warnings in this structure.
1170 %
1171 */
1172 MagickExport MagickBooleanType GetImageExtrema(const Image *image,
1173   size_t *minima,size_t *maxima,ExceptionInfo *exception)
1174 {
1175   double
1176     max,
1177     min;
1178
1179   MagickBooleanType
1180     status;
1181
1182   assert(image != (Image *) NULL);
1183   assert(image->signature == MagickCoreSignature);
1184   if (image->debug != MagickFalse)
1185     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1186   status=GetImageRange(image,&min,&max,exception);
1187   *minima=(size_t) ceil(min-0.5);
1188   *maxima=(size_t) floor(max+0.5);
1189   return(status);
1190 }
1191 \f
1192 /*
1193 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1194 %                                                                             %
1195 %                                                                             %
1196 %                                                                             %
1197 %   G e t I m a g e K u r t o s i s                                           %
1198 %                                                                             %
1199 %                                                                             %
1200 %                                                                             %
1201 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1202 %
1203 %  GetImageKurtosis() returns the kurtosis and skewness of one or more image
1204 %  channels.
1205 %
1206 %  The format of the GetImageKurtosis method is:
1207 %
1208 %      MagickBooleanType GetImageKurtosis(const Image *image,double *kurtosis,
1209 %        double *skewness,ExceptionInfo *exception)
1210 %
1211 %  A description of each parameter follows:
1212 %
1213 %    o image: the image.
1214 %
1215 %    o kurtosis: the kurtosis of the channel.
1216 %
1217 %    o skewness: the skewness of the channel.
1218 %
1219 %    o exception: return any errors or warnings in this structure.
1220 %
1221 */
1222 MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
1223   double *kurtosis,double *skewness,ExceptionInfo *exception)
1224 {
1225   ChannelStatistics
1226     *channel_statistics;
1227
1228   assert(image != (Image *) NULL);
1229   assert(image->signature == MagickCoreSignature);
1230   if (image->debug != MagickFalse)
1231     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1232   channel_statistics=GetImageStatistics(image,exception);
1233   if (channel_statistics == (ChannelStatistics *) NULL)
1234     return(MagickFalse);
1235   *kurtosis=channel_statistics[CompositePixelChannel].kurtosis;
1236   *skewness=channel_statistics[CompositePixelChannel].skewness;
1237   channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1238     channel_statistics);
1239   return(MagickTrue);
1240 }
1241 \f
1242 /*
1243 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1244 %                                                                             %
1245 %                                                                             %
1246 %                                                                             %
1247 %   G e t I m a g e M e a n                                                   %
1248 %                                                                             %
1249 %                                                                             %
1250 %                                                                             %
1251 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1252 %
1253 %  GetImageMean() returns the mean and standard deviation of one or more image
1254 %  channels.
1255 %
1256 %  The format of the GetImageMean method is:
1257 %
1258 %      MagickBooleanType GetImageMean(const Image *image,double *mean,
1259 %        double *standard_deviation,ExceptionInfo *exception)
1260 %
1261 %  A description of each parameter follows:
1262 %
1263 %    o image: the image.
1264 %
1265 %    o mean: the average value in the channel.
1266 %
1267 %    o standard_deviation: the standard deviation of the channel.
1268 %
1269 %    o exception: return any errors or warnings in this structure.
1270 %
1271 */
1272 MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean,
1273   double *standard_deviation,ExceptionInfo *exception)
1274 {
1275   ChannelStatistics
1276     *channel_statistics;
1277
1278   assert(image != (Image *) NULL);
1279   assert(image->signature == MagickCoreSignature);
1280   if (image->debug != MagickFalse)
1281     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1282   channel_statistics=GetImageStatistics(image,exception);
1283   if (channel_statistics == (ChannelStatistics *) NULL)
1284     return(MagickFalse);
1285   *mean=channel_statistics[CompositePixelChannel].mean;
1286   *standard_deviation=
1287     channel_statistics[CompositePixelChannel].standard_deviation;
1288   channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1289     channel_statistics);
1290   return(MagickTrue);
1291 }
1292 \f
1293 /*
1294 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1295 %                                                                             %
1296 %                                                                             %
1297 %                                                                             %
1298 %   G e t I m a g e M o m e n t s                                             %
1299 %                                                                             %
1300 %                                                                             %
1301 %                                                                             %
1302 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1303 %
1304 %  GetImageMoments() returns the normalized moments of one or more image
1305 %  channels.
1306 %
1307 %  The format of the GetImageMoments method is:
1308 %
1309 %      ChannelMoments *GetImageMoments(const Image *image,
1310 %        ExceptionInfo *exception)
1311 %
1312 %  A description of each parameter follows:
1313 %
1314 %    o image: the image.
1315 %
1316 %    o exception: return any errors or warnings in this structure.
1317 %
1318 */
1319
1320 static size_t GetImageChannels(const Image *image)
1321 {
1322   register ssize_t
1323     i;
1324
1325   size_t
1326     channels;
1327
1328   channels=0;
1329   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1330   {
1331     PixelChannel channel = GetPixelChannelChannel(image,i);
1332     PixelTrait traits = GetPixelChannelTraits(image,channel);
1333     if (traits == UndefinedPixelTrait)
1334       continue;
1335     if ((traits & UpdatePixelTrait) == 0)
1336       continue;
1337     channels++;
1338   }
1339   return((size_t) (channels == 0 ? 1 : channels));
1340 }
1341
1342 MagickExport ChannelMoments *GetImageMoments(const Image *image,
1343   ExceptionInfo *exception)
1344 {
1345 #define MaxNumberImageMoments  8
1346
1347   CacheView
1348     *image_view;
1349
1350   ChannelMoments
1351     *channel_moments;
1352
1353   double
1354     M00[MaxPixelChannels+1],
1355     M01[MaxPixelChannels+1],
1356     M02[MaxPixelChannels+1],
1357     M03[MaxPixelChannels+1],
1358     M10[MaxPixelChannels+1],
1359     M11[MaxPixelChannels+1],
1360     M12[MaxPixelChannels+1],
1361     M20[MaxPixelChannels+1],
1362     M21[MaxPixelChannels+1],
1363     M22[MaxPixelChannels+1],
1364     M30[MaxPixelChannels+1];
1365
1366   PointInfo
1367     centroid[MaxPixelChannels+1];
1368
1369   ssize_t
1370     channel,
1371     y;
1372
1373   assert(image != (Image *) NULL);
1374   assert(image->signature == MagickCoreSignature);
1375   if (image->debug != MagickFalse)
1376     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1377   channel_moments=(ChannelMoments *) AcquireQuantumMemory(MaxPixelChannels+1,
1378     sizeof(*channel_moments));
1379   if (channel_moments == (ChannelMoments *) NULL)
1380     return(channel_moments);
1381   (void) memset(channel_moments,0,(MaxPixelChannels+1)*
1382     sizeof(*channel_moments));
1383   (void) memset(centroid,0,sizeof(centroid));
1384   (void) memset(M00,0,sizeof(M00));
1385   (void) memset(M01,0,sizeof(M01));
1386   (void) memset(M02,0,sizeof(M02));
1387   (void) memset(M03,0,sizeof(M03));
1388   (void) memset(M10,0,sizeof(M10));
1389   (void) memset(M11,0,sizeof(M11));
1390   (void) memset(M12,0,sizeof(M12));
1391   (void) memset(M20,0,sizeof(M20));
1392   (void) memset(M21,0,sizeof(M21));
1393   (void) memset(M22,0,sizeof(M22));
1394   (void) memset(M30,0,sizeof(M30));
1395   image_view=AcquireVirtualCacheView(image,exception);
1396   for (y=0; y < (ssize_t) image->rows; y++)
1397   {
1398     register const Quantum
1399       *magick_restrict p;
1400
1401     register ssize_t
1402       x;
1403
1404     /*
1405       Compute center of mass (centroid).
1406     */
1407     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1408     if (p == (const Quantum *) NULL)
1409       break;
1410     for (x=0; x < (ssize_t) image->columns; x++)
1411     {
1412       register ssize_t
1413         i;
1414
1415       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1416       {
1417         PixelChannel channel = GetPixelChannelChannel(image,i);
1418         PixelTrait traits = GetPixelChannelTraits(image,channel);
1419         if (traits == UndefinedPixelTrait)
1420           continue;
1421         if ((traits & UpdatePixelTrait) == 0)
1422           continue;
1423         M00[channel]+=QuantumScale*p[i];
1424         M00[MaxPixelChannels]+=QuantumScale*p[i];
1425         M10[channel]+=x*QuantumScale*p[i];
1426         M10[MaxPixelChannels]+=x*QuantumScale*p[i];
1427         M01[channel]+=y*QuantumScale*p[i];
1428         M01[MaxPixelChannels]+=y*QuantumScale*p[i];
1429       }
1430       p+=GetPixelChannels(image);
1431     }
1432   }
1433   for (channel=0; channel <= MaxPixelChannels; channel++)
1434   {
1435     /*
1436        Compute center of mass (centroid).
1437     */
1438     if (M00[channel] < MagickEpsilon)
1439       {
1440         M00[channel]+=MagickEpsilon;
1441         centroid[channel].x=(double) image->columns/2.0;
1442         centroid[channel].y=(double) image->rows/2.0;
1443         continue;
1444       }
1445     M00[channel]+=MagickEpsilon;
1446     centroid[channel].x=M10[channel]/M00[channel];
1447     centroid[channel].y=M01[channel]/M00[channel];
1448   }
1449   for (y=0; y < (ssize_t) image->rows; y++)
1450   {
1451     register const Quantum
1452       *magick_restrict p;
1453
1454     register ssize_t
1455       x;
1456
1457     /*
1458       Compute the image moments.
1459     */
1460     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1461     if (p == (const Quantum *) NULL)
1462       break;
1463     for (x=0; x < (ssize_t) image->columns; x++)
1464     {
1465       register ssize_t
1466         i;
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) shared(status,initialize) \
1803     magick_number_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       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1835       {
1836         PixelChannel channel = GetPixelChannelChannel(image,i);
1837         PixelTrait traits = GetPixelChannelTraits(image,channel);
1838         if (traits == UndefinedPixelTrait)
1839           continue;
1840         if ((traits & UpdatePixelTrait) == 0)
1841           continue;
1842                                 if (row_initialize != MagickFalse)
1843           {
1844             row_minima=(double) p[i];
1845             row_maxima=(double) p[i];
1846             row_initialize=MagickFalse;
1847           }
1848         else
1849           {
1850             if ((double) p[i] < row_minima)
1851               row_minima=(double) p[i];
1852             if ((double) p[i] > row_maxima)
1853               row_maxima=(double) p[i];
1854          }
1855       }
1856       p+=GetPixelChannels(image);
1857     }
1858 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1859 #pragma omp critical (MagickCore_GetImageRange)
1860 #endif
1861     {
1862       if (initialize != MagickFalse)
1863         {
1864           *minima=row_minima;
1865           *maxima=row_maxima;
1866           initialize=MagickFalse;
1867         }
1868       else
1869         {
1870           if (row_minima < *minima)
1871             *minima=row_minima;
1872           if (row_maxima > *maxima)
1873             *maxima=row_maxima;
1874         }
1875     }
1876   }
1877   image_view=DestroyCacheView(image_view);
1878   return(status);
1879 }
1880 \f
1881 /*
1882 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1883 %                                                                             %
1884 %                                                                             %
1885 %                                                                             %
1886 %   G e t I m a g e S t a t i s t i c s                                       %
1887 %                                                                             %
1888 %                                                                             %
1889 %                                                                             %
1890 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1891 %
1892 %  GetImageStatistics() returns statistics for each channel in the image.  The
1893 %  statistics include the channel depth, its minima, maxima, mean, standard
1894 %  deviation, kurtosis and skewness.  You can access the red channel mean, for
1895 %  example, like this:
1896 %
1897 %      channel_statistics=GetImageStatistics(image,exception);
1898 %      red_mean=channel_statistics[RedPixelChannel].mean;
1899 %
1900 %  Use MagickRelinquishMemory() to free the statistics buffer.
1901 %
1902 %  The format of the GetImageStatistics method is:
1903 %
1904 %      ChannelStatistics *GetImageStatistics(const Image *image,
1905 %        ExceptionInfo *exception)
1906 %
1907 %  A description of each parameter follows:
1908 %
1909 %    o image: the image.
1910 %
1911 %    o exception: return any errors or warnings in this structure.
1912 %
1913 */
1914 MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
1915   ExceptionInfo *exception)
1916 {
1917   ChannelStatistics
1918     *channel_statistics;
1919
1920   double
1921     area,
1922     *histogram,
1923     standard_deviation;
1924
1925   MagickStatusType
1926     status;
1927
1928   QuantumAny
1929     range;
1930
1931   register ssize_t
1932     i;
1933
1934   size_t
1935     depth;
1936
1937   ssize_t
1938     y;
1939
1940   assert(image != (Image *) NULL);
1941   assert(image->signature == MagickCoreSignature);
1942   if (image->debug != MagickFalse)
1943     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1944   histogram=(double *) AcquireQuantumMemory(MaxMap+1UL,GetPixelChannels(image)*
1945     sizeof(*histogram));
1946   channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(
1947     MaxPixelChannels+1,sizeof(*channel_statistics));
1948   if ((channel_statistics == (ChannelStatistics *) NULL) ||
1949       (histogram == (double *) NULL))
1950     {
1951       if (histogram != (double *) NULL)
1952         histogram=(double *) RelinquishMagickMemory(histogram);
1953       if (channel_statistics != (ChannelStatistics *) NULL)
1954         channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1955           channel_statistics);
1956       return(channel_statistics);
1957     }
1958   (void) memset(channel_statistics,0,(MaxPixelChannels+1)*
1959     sizeof(*channel_statistics));
1960   for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
1961   {
1962     channel_statistics[i].depth=1;
1963     channel_statistics[i].maxima=(-MagickMaximumValue);
1964     channel_statistics[i].minima=MagickMaximumValue;
1965   }
1966   (void) memset(histogram,0,(MaxMap+1)*GetPixelChannels(image)*
1967     sizeof(*histogram));
1968   for (y=0; y < (ssize_t) image->rows; y++)
1969   {
1970     register const Quantum
1971       *magick_restrict p;
1972
1973     register ssize_t
1974       x;
1975
1976     /*
1977       Compute pixel statistics.
1978     */
1979     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1980     if (p == (const Quantum *) NULL)
1981       break;
1982     for (x=0; x < (ssize_t) image->columns; x++)
1983     {
1984       register ssize_t
1985         i;
1986
1987       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1988       {
1989         PixelChannel channel = GetPixelChannelChannel(image,i);
1990         PixelTrait traits = GetPixelChannelTraits(image,channel);
1991         if (traits == UndefinedPixelTrait)
1992           continue;
1993         if ((traits & UpdatePixelTrait) == 0)
1994           continue;
1995         if (channel_statistics[channel].depth != MAGICKCORE_QUANTUM_DEPTH)
1996           {
1997             depth=channel_statistics[channel].depth;
1998             range=GetQuantumRange(depth);
1999             status=p[i] != ScaleAnyToQuantum(ScaleQuantumToAny(p[i],range),
2000               range) ? MagickTrue : MagickFalse;
2001             if (status != MagickFalse)
2002               {
2003                 channel_statistics[channel].depth++;
2004                 i--;
2005                 continue;
2006               }
2007           }
2008         if ((double) p[i] < channel_statistics[channel].minima)
2009           channel_statistics[channel].minima=(double) p[i];
2010         if ((double) p[i] > channel_statistics[channel].maxima)
2011           channel_statistics[channel].maxima=(double) p[i];
2012         channel_statistics[channel].sum+=p[i];
2013         channel_statistics[channel].sum_squared+=(double) p[i]*p[i];
2014         channel_statistics[channel].sum_cubed+=(double) p[i]*p[i]*p[i];
2015         channel_statistics[channel].sum_fourth_power+=(double) p[i]*p[i]*p[i]*
2016           p[i];
2017         channel_statistics[channel].area++;
2018         if ((double) p[i] < channel_statistics[CompositePixelChannel].minima)
2019           channel_statistics[CompositePixelChannel].minima=(double) p[i];
2020         if ((double) p[i] > channel_statistics[CompositePixelChannel].maxima)
2021           channel_statistics[CompositePixelChannel].maxima=(double) p[i];
2022         histogram[GetPixelChannels(image)*ScaleQuantumToMap(
2023           ClampToQuantum((double) p[i]))+i]++;
2024         channel_statistics[CompositePixelChannel].sum+=(double) p[i];
2025         channel_statistics[CompositePixelChannel].sum_squared+=(double)
2026           p[i]*p[i];
2027         channel_statistics[CompositePixelChannel].sum_cubed+=(double)
2028           p[i]*p[i]*p[i];
2029         channel_statistics[CompositePixelChannel].sum_fourth_power+=(double)
2030           p[i]*p[i]*p[i]*p[i];
2031         channel_statistics[CompositePixelChannel].area++;
2032       }
2033       p+=GetPixelChannels(image);
2034     }
2035   }
2036   for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
2037   {
2038     /*
2039       Normalize pixel statistics.
2040     */
2041     area=PerceptibleReciprocal(channel_statistics[i].area);
2042     channel_statistics[i].sum*=area;
2043     channel_statistics[i].sum_squared*=area;
2044     channel_statistics[i].sum_cubed*=area;
2045     channel_statistics[i].sum_fourth_power*=area;
2046     channel_statistics[i].mean=channel_statistics[i].sum;
2047     channel_statistics[i].variance=channel_statistics[i].sum_squared;
2048     standard_deviation=sqrt(channel_statistics[i].variance-
2049       (channel_statistics[i].mean*channel_statistics[i].mean));
2050     standard_deviation=sqrt(PerceptibleReciprocal(channel_statistics[i].area-
2051       1.0)*channel_statistics[i].area*standard_deviation*standard_deviation);
2052     channel_statistics[i].standard_deviation=standard_deviation;
2053   }
2054   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2055   {
2056     double
2057       number_bins;
2058
2059     register ssize_t
2060       j;
2061
2062     /*
2063       Compute pixel entropy.
2064     */
2065     PixelChannel channel = GetPixelChannelChannel(image,i);
2066     number_bins=0.0;
2067     for (j=0; j <= (ssize_t) MaxMap; j++)
2068       if (histogram[GetPixelChannels(image)*j+i] > 0.0)
2069         number_bins++;
2070     area=PerceptibleReciprocal(channel_statistics[channel].area);
2071     for (j=0; j <= (ssize_t) MaxMap; j++)
2072     {
2073       double
2074         count;
2075
2076       count=area*histogram[GetPixelChannels(image)*j+i];
2077       channel_statistics[channel].entropy+=-count*MagickLog10(count)*
2078         PerceptibleReciprocal(MagickLog10(number_bins));
2079       channel_statistics[CompositePixelChannel].entropy+=-count*
2080         MagickLog10(count)*PerceptibleReciprocal(MagickLog10(number_bins))/
2081         GetPixelChannels(image);
2082     }
2083   }
2084   histogram=(double *) RelinquishMagickMemory(histogram);
2085   for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
2086   {
2087     /*
2088       Compute kurtosis & skewness statistics.
2089     */
2090     standard_deviation=PerceptibleReciprocal(
2091       channel_statistics[i].standard_deviation);
2092     channel_statistics[i].skewness=(channel_statistics[i].sum_cubed-3.0*
2093       channel_statistics[i].mean*channel_statistics[i].sum_squared+2.0*
2094       channel_statistics[i].mean*channel_statistics[i].mean*
2095       channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2096       standard_deviation);
2097     channel_statistics[i].kurtosis=(channel_statistics[i].sum_fourth_power-4.0*
2098       channel_statistics[i].mean*channel_statistics[i].sum_cubed+6.0*
2099       channel_statistics[i].mean*channel_statistics[i].mean*
2100       channel_statistics[i].sum_squared-3.0*channel_statistics[i].mean*
2101       channel_statistics[i].mean*1.0*channel_statistics[i].mean*
2102       channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2103       standard_deviation*standard_deviation)-3.0;
2104   }
2105   channel_statistics[CompositePixelChannel].mean=0.0;
2106   channel_statistics[CompositePixelChannel].standard_deviation=0.0;
2107   channel_statistics[CompositePixelChannel].entropy=0.0;
2108   for (i=0; i < (ssize_t) MaxPixelChannels; i++)
2109   {
2110     channel_statistics[CompositePixelChannel].mean+=
2111       channel_statistics[i].mean;
2112     channel_statistics[CompositePixelChannel].standard_deviation+=
2113       channel_statistics[i].standard_deviation;
2114     channel_statistics[CompositePixelChannel].entropy+=
2115       channel_statistics[i].entropy;
2116   }
2117   channel_statistics[CompositePixelChannel].mean/=(double)
2118     GetImageChannels(image);
2119   channel_statistics[CompositePixelChannel].standard_deviation/=(double)
2120     GetImageChannels(image);
2121   channel_statistics[CompositePixelChannel].entropy/=(double)
2122     GetImageChannels(image);
2123   if (y < (ssize_t) image->rows)
2124     channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
2125       channel_statistics);
2126   return(channel_statistics);
2127 }
2128 \f
2129 /*
2130 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2131 %                                                                             %
2132 %                                                                             %
2133 %                                                                             %
2134 %     P o l y n o m i a l I m a g e                                           %
2135 %                                                                             %
2136 %                                                                             %
2137 %                                                                             %
2138 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2139 %
2140 %  PolynomialImage() returns a new image where each pixel is the sum of the
2141 %  pixels in the image sequence after applying its corresponding terms
2142 %  (coefficient and degree pairs).
2143 %
2144 %  The format of the PolynomialImage method is:
2145 %
2146 %      Image *PolynomialImage(const Image *images,const size_t number_terms,
2147 %        const double *terms,ExceptionInfo *exception)
2148 %
2149 %  A description of each parameter follows:
2150 %
2151 %    o images: the image sequence.
2152 %
2153 %    o number_terms: the number of terms in the list.  The actual list length
2154 %      is 2 x number_terms + 1 (the constant).
2155 %
2156 %    o terms: the list of polynomial coefficients and degree pairs and a
2157 %      constant.
2158 %
2159 %    o exception: return any errors or warnings in this structure.
2160 %
2161 */
2162
2163 MagickExport Image *PolynomialImage(const Image *images,
2164   const size_t number_terms,const double *terms,ExceptionInfo *exception)
2165 {
2166 #define PolynomialImageTag  "Polynomial/Image"
2167
2168   CacheView
2169     *polynomial_view;
2170
2171   Image
2172     *image;
2173
2174   MagickBooleanType
2175     status;
2176
2177   MagickOffsetType
2178     progress;
2179
2180   PixelChannels
2181     **magick_restrict polynomial_pixels;
2182
2183   size_t
2184     number_images;
2185
2186   ssize_t
2187     y;
2188
2189   assert(images != (Image *) NULL);
2190   assert(images->signature == MagickCoreSignature);
2191   if (images->debug != MagickFalse)
2192     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
2193   assert(exception != (ExceptionInfo *) NULL);
2194   assert(exception->signature == MagickCoreSignature);
2195   image=AcquireImageCanvas(images,exception);
2196   if (image == (Image *) NULL)
2197     return((Image *) NULL);
2198   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
2199     {
2200       image=DestroyImage(image);
2201       return((Image *) NULL);
2202     }
2203   number_images=GetImageListLength(images);
2204   polynomial_pixels=AcquirePixelThreadSet(images);
2205   if (polynomial_pixels == (PixelChannels **) NULL)
2206     {
2207       image=DestroyImage(image);
2208       (void) ThrowMagickException(exception,GetMagickModule(),
2209         ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
2210       return((Image *) NULL);
2211     }
2212   /*
2213     Polynomial image pixels.
2214   */
2215   status=MagickTrue;
2216   progress=0;
2217   polynomial_view=AcquireAuthenticCacheView(image,exception);
2218 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2219   #pragma omp parallel for schedule(static) shared(progress,status) \
2220     magick_number_threads(image,image,image->rows,1)
2221 #endif
2222   for (y=0; y < (ssize_t) image->rows; y++)
2223   {
2224     CacheView
2225       *image_view;
2226
2227     const Image
2228       *next;
2229
2230     const int
2231       id = GetOpenMPThreadId();
2232
2233     register ssize_t
2234       i,
2235       x;
2236
2237     register PixelChannels
2238       *polynomial_pixel;
2239
2240     register Quantum
2241       *magick_restrict q;
2242
2243     ssize_t
2244       j;
2245
2246     if (status == MagickFalse)
2247       continue;
2248     q=QueueCacheViewAuthenticPixels(polynomial_view,0,y,image->columns,1,
2249       exception);
2250     if (q == (Quantum *) NULL)
2251       {
2252         status=MagickFalse;
2253         continue;
2254       }
2255     polynomial_pixel=polynomial_pixels[id];
2256     for (j=0; j < (ssize_t) image->columns; j++)
2257       for (i=0; i < MaxPixelChannels; i++)
2258         polynomial_pixel[j].channel[i]=0.0;
2259     next=images;
2260     for (j=0; j < (ssize_t) number_images; j++)
2261     {
2262       register const Quantum
2263         *p;
2264
2265       if (j >= (ssize_t) number_terms)
2266         continue;
2267       image_view=AcquireVirtualCacheView(next,exception);
2268       p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2269       if (p == (const Quantum *) NULL)
2270         {
2271           image_view=DestroyCacheView(image_view);
2272           break;
2273         }
2274       for (x=0; x < (ssize_t) image->columns; x++)
2275       {
2276         register ssize_t
2277           i;
2278
2279         for (i=0; i < (ssize_t) GetPixelChannels(next); i++)
2280         {
2281           MagickRealType
2282             coefficient,
2283             degree;
2284
2285           PixelChannel channel = GetPixelChannelChannel(image,i);
2286           PixelTrait traits = GetPixelChannelTraits(next,channel);
2287           PixelTrait polynomial_traits=GetPixelChannelTraits(image,channel);
2288           if ((traits == UndefinedPixelTrait) ||
2289               (polynomial_traits == UndefinedPixelTrait))
2290             continue;
2291           if ((traits & UpdatePixelTrait) == 0)
2292             continue;
2293           coefficient=(MagickRealType) terms[2*j];
2294           degree=(MagickRealType) terms[(j << 1)+1];
2295           polynomial_pixel[x].channel[i]+=coefficient*
2296             pow(QuantumScale*GetPixelChannel(image,channel,p),degree);
2297         }
2298         p+=GetPixelChannels(next);
2299       }
2300       image_view=DestroyCacheView(image_view);
2301       next=GetNextImageInList(next);
2302     }
2303     for (x=0; x < (ssize_t) image->columns; x++)
2304     {
2305       register ssize_t
2306         i;
2307
2308       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2309       {
2310         PixelChannel channel = GetPixelChannelChannel(image,i);
2311         PixelTrait traits = GetPixelChannelTraits(image,channel);
2312         if (traits == UndefinedPixelTrait)
2313           continue;
2314         if ((traits & UpdatePixelTrait) == 0)
2315           continue;
2316         q[i]=ClampToQuantum(QuantumRange*polynomial_pixel[x].channel[i]);
2317       }
2318       q+=GetPixelChannels(image);
2319     }
2320     if (SyncCacheViewAuthenticPixels(polynomial_view,exception) == MagickFalse)
2321       status=MagickFalse;
2322     if (images->progress_monitor != (MagickProgressMonitor) NULL)
2323       {
2324         MagickBooleanType
2325           proceed;
2326
2327 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2328         #pragma omp critical (MagickCore_PolynomialImages)
2329 #endif
2330         proceed=SetImageProgress(images,PolynomialImageTag,progress++,
2331           image->rows);
2332         if (proceed == MagickFalse)
2333           status=MagickFalse;
2334       }
2335   }
2336   polynomial_view=DestroyCacheView(polynomial_view);
2337   polynomial_pixels=DestroyPixelThreadSet(polynomial_pixels);
2338   if (status == MagickFalse)
2339     image=DestroyImage(image);
2340   return(image);
2341 }
2342 \f
2343 /*
2344 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2345 %                                                                             %
2346 %                                                                             %
2347 %                                                                             %
2348 %     S t a t i s t i c I m a g e                                             %
2349 %                                                                             %
2350 %                                                                             %
2351 %                                                                             %
2352 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2353 %
2354 %  StatisticImage() makes each pixel the min / max / median / mode / etc. of
2355 %  the neighborhood of the specified width and height.
2356 %
2357 %  The format of the StatisticImage method is:
2358 %
2359 %      Image *StatisticImage(const Image *image,const StatisticType type,
2360 %        const size_t width,const size_t height,ExceptionInfo *exception)
2361 %
2362 %  A description of each parameter follows:
2363 %
2364 %    o image: the image.
2365 %
2366 %    o type: the statistic type (median, mode, etc.).
2367 %
2368 %    o width: the width of the pixel neighborhood.
2369 %
2370 %    o height: the height of the pixel neighborhood.
2371 %
2372 %    o exception: return any errors or warnings in this structure.
2373 %
2374 */
2375
2376 typedef struct _SkipNode
2377 {
2378   size_t
2379     next[9],
2380     count,
2381     signature;
2382 } SkipNode;
2383
2384 typedef struct _SkipList
2385 {
2386   ssize_t
2387     level;
2388
2389   SkipNode
2390     *nodes;
2391 } SkipList;
2392
2393 typedef struct _PixelList
2394 {
2395   size_t
2396     length,
2397     seed;
2398
2399   SkipList
2400     skip_list;
2401
2402   size_t
2403     signature;
2404 } PixelList;
2405
2406 static PixelList *DestroyPixelList(PixelList *pixel_list)
2407 {
2408   if (pixel_list == (PixelList *) NULL)
2409     return((PixelList *) NULL);
2410   if (pixel_list->skip_list.nodes != (SkipNode *) NULL)
2411     pixel_list->skip_list.nodes=(SkipNode *) RelinquishAlignedMemory(
2412       pixel_list->skip_list.nodes);
2413   pixel_list=(PixelList *) RelinquishMagickMemory(pixel_list);
2414   return(pixel_list);
2415 }
2416
2417 static PixelList **DestroyPixelListThreadSet(PixelList **pixel_list)
2418 {
2419   register ssize_t
2420     i;
2421
2422   assert(pixel_list != (PixelList **) NULL);
2423   for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
2424     if (pixel_list[i] != (PixelList *) NULL)
2425       pixel_list[i]=DestroyPixelList(pixel_list[i]);
2426   pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list);
2427   return(pixel_list);
2428 }
2429
2430 static PixelList *AcquirePixelList(const size_t width,const size_t height)
2431 {
2432   PixelList
2433     *pixel_list;
2434
2435   pixel_list=(PixelList *) AcquireMagickMemory(sizeof(*pixel_list));
2436   if (pixel_list == (PixelList *) NULL)
2437     return(pixel_list);
2438   (void) memset((void *) pixel_list,0,sizeof(*pixel_list));
2439   pixel_list->length=width*height;
2440   pixel_list->skip_list.nodes=(SkipNode *) AcquireAlignedMemory(65537UL,
2441     sizeof(*pixel_list->skip_list.nodes));
2442   if (pixel_list->skip_list.nodes == (SkipNode *) NULL)
2443     return(DestroyPixelList(pixel_list));
2444   (void) memset(pixel_list->skip_list.nodes,0,65537UL*
2445     sizeof(*pixel_list->skip_list.nodes));
2446   pixel_list->signature=MagickCoreSignature;
2447   return(pixel_list);
2448 }
2449
2450 static PixelList **AcquirePixelListThreadSet(const size_t width,
2451   const size_t height)
2452 {
2453   PixelList
2454     **pixel_list;
2455
2456   register ssize_t
2457     i;
2458
2459   size_t
2460     number_threads;
2461
2462   number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
2463   pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
2464     sizeof(*pixel_list));
2465   if (pixel_list == (PixelList **) NULL)
2466     return((PixelList **) NULL);
2467   (void) memset(pixel_list,0,number_threads*sizeof(*pixel_list));
2468   for (i=0; i < (ssize_t) number_threads; i++)
2469   {
2470     pixel_list[i]=AcquirePixelList(width,height);
2471     if (pixel_list[i] == (PixelList *) NULL)
2472       return(DestroyPixelListThreadSet(pixel_list));
2473   }
2474   return(pixel_list);
2475 }
2476
2477 static void AddNodePixelList(PixelList *pixel_list,const size_t color)
2478 {
2479   register SkipList
2480     *p;
2481
2482   register ssize_t
2483     level;
2484
2485   size_t
2486     search,
2487     update[9];
2488
2489   /*
2490     Initialize the node.
2491   */
2492   p=(&pixel_list->skip_list);
2493   p->nodes[color].signature=pixel_list->signature;
2494   p->nodes[color].count=1;
2495   /*
2496     Determine where it belongs in the list.
2497   */
2498   search=65536UL;
2499   for (level=p->level; level >= 0; level--)
2500   {
2501     while (p->nodes[search].next[level] < color)
2502       search=p->nodes[search].next[level];
2503     update[level]=search;
2504   }
2505   /*
2506     Generate a pseudo-random level for this node.
2507   */
2508   for (level=0; ; level++)
2509   {
2510     pixel_list->seed=(pixel_list->seed*42893621L)+1L;
2511     if ((pixel_list->seed & 0x300) != 0x300)
2512       break;
2513   }
2514   if (level > 8)
2515     level=8;
2516   if (level > (p->level+2))
2517     level=p->level+2;
2518   /*
2519     If we're raising the list's level, link back to the root node.
2520   */
2521   while (level > p->level)
2522   {
2523     p->level++;
2524     update[p->level]=65536UL;
2525   }
2526   /*
2527     Link the node into the skip-list.
2528   */
2529   do
2530   {
2531     p->nodes[color].next[level]=p->nodes[update[level]].next[level];
2532     p->nodes[update[level]].next[level]=color;
2533   } while (level-- > 0);
2534 }
2535
2536 static inline void GetMaximumPixelList(PixelList *pixel_list,Quantum *pixel)
2537 {
2538   register SkipList
2539     *p;
2540
2541   size_t
2542     color,
2543     maximum;
2544
2545   ssize_t
2546     count;
2547
2548   /*
2549     Find the maximum value for each of the color.
2550   */
2551   p=(&pixel_list->skip_list);
2552   color=65536L;
2553   count=0;
2554   maximum=p->nodes[color].next[0];
2555   do
2556   {
2557     color=p->nodes[color].next[0];
2558     if (color > maximum)
2559       maximum=color;
2560     count+=p->nodes[color].count;
2561   } while (count < (ssize_t) pixel_list->length);
2562   *pixel=ScaleShortToQuantum((unsigned short) maximum);
2563 }
2564
2565 static inline void GetMeanPixelList(PixelList *pixel_list,Quantum *pixel)
2566 {
2567   double
2568     sum;
2569
2570   register SkipList
2571     *p;
2572
2573   size_t
2574     color;
2575
2576   ssize_t
2577     count;
2578
2579   /*
2580     Find the mean value for each of the color.
2581   */
2582   p=(&pixel_list->skip_list);
2583   color=65536L;
2584   count=0;
2585   sum=0.0;
2586   do
2587   {
2588     color=p->nodes[color].next[0];
2589     sum+=(double) p->nodes[color].count*color;
2590     count+=p->nodes[color].count;
2591   } while (count < (ssize_t) pixel_list->length);
2592   sum/=pixel_list->length;
2593   *pixel=ScaleShortToQuantum((unsigned short) sum);
2594 }
2595
2596 static inline void GetMedianPixelList(PixelList *pixel_list,Quantum *pixel)
2597 {
2598   register SkipList
2599     *p;
2600
2601   size_t
2602     color;
2603
2604   ssize_t
2605     count;
2606
2607   /*
2608     Find the median value for each of the color.
2609   */
2610   p=(&pixel_list->skip_list);
2611   color=65536L;
2612   count=0;
2613   do
2614   {
2615     color=p->nodes[color].next[0];
2616     count+=p->nodes[color].count;
2617   } while (count <= (ssize_t) (pixel_list->length >> 1));
2618   *pixel=ScaleShortToQuantum((unsigned short) color);
2619 }
2620
2621 static inline void GetMinimumPixelList(PixelList *pixel_list,Quantum *pixel)
2622 {
2623   register SkipList
2624     *p;
2625
2626   size_t
2627     color,
2628     minimum;
2629
2630   ssize_t
2631     count;
2632
2633   /*
2634     Find the minimum value for each of the color.
2635   */
2636   p=(&pixel_list->skip_list);
2637   count=0;
2638   color=65536UL;
2639   minimum=p->nodes[color].next[0];
2640   do
2641   {
2642     color=p->nodes[color].next[0];
2643     if (color < minimum)
2644       minimum=color;
2645     count+=p->nodes[color].count;
2646   } while (count < (ssize_t) pixel_list->length);
2647   *pixel=ScaleShortToQuantum((unsigned short) minimum);
2648 }
2649
2650 static inline void GetModePixelList(PixelList *pixel_list,Quantum *pixel)
2651 {
2652   register SkipList
2653     *p;
2654
2655   size_t
2656     color,
2657     max_count,
2658     mode;
2659
2660   ssize_t
2661     count;
2662
2663   /*
2664     Make each pixel the 'predominant color' of the specified neighborhood.
2665   */
2666   p=(&pixel_list->skip_list);
2667   color=65536L;
2668   mode=color;
2669   max_count=p->nodes[mode].count;
2670   count=0;
2671   do
2672   {
2673     color=p->nodes[color].next[0];
2674     if (p->nodes[color].count > max_count)
2675       {
2676         mode=color;
2677         max_count=p->nodes[mode].count;
2678       }
2679     count+=p->nodes[color].count;
2680   } while (count < (ssize_t) pixel_list->length);
2681   *pixel=ScaleShortToQuantum((unsigned short) mode);
2682 }
2683
2684 static inline void GetNonpeakPixelList(PixelList *pixel_list,Quantum *pixel)
2685 {
2686   register SkipList
2687     *p;
2688
2689   size_t
2690     color,
2691     next,
2692     previous;
2693
2694   ssize_t
2695     count;
2696
2697   /*
2698     Finds the non peak value for each of the colors.
2699   */
2700   p=(&pixel_list->skip_list);
2701   color=65536L;
2702   next=p->nodes[color].next[0];
2703   count=0;
2704   do
2705   {
2706     previous=color;
2707     color=next;
2708     next=p->nodes[color].next[0];
2709     count+=p->nodes[color].count;
2710   } while (count <= (ssize_t) (pixel_list->length >> 1));
2711   if ((previous == 65536UL) && (next != 65536UL))
2712     color=next;
2713   else
2714     if ((previous != 65536UL) && (next == 65536UL))
2715       color=previous;
2716   *pixel=ScaleShortToQuantum((unsigned short) color);
2717 }
2718
2719 static inline void GetRootMeanSquarePixelList(PixelList *pixel_list,
2720   Quantum *pixel)
2721 {
2722   double
2723     sum;
2724
2725   register SkipList
2726     *p;
2727
2728   size_t
2729     color;
2730
2731   ssize_t
2732     count;
2733
2734   /*
2735     Find the root mean square value for each of the color.
2736   */
2737   p=(&pixel_list->skip_list);
2738   color=65536L;
2739   count=0;
2740   sum=0.0;
2741   do
2742   {
2743     color=p->nodes[color].next[0];
2744     sum+=(double) (p->nodes[color].count*color*color);
2745     count+=p->nodes[color].count;
2746   } while (count < (ssize_t) pixel_list->length);
2747   sum/=pixel_list->length;
2748   *pixel=ScaleShortToQuantum((unsigned short) sqrt(sum));
2749 }
2750
2751 static inline void GetStandardDeviationPixelList(PixelList *pixel_list,
2752   Quantum *pixel)
2753 {
2754   double
2755     sum,
2756     sum_squared;
2757
2758   register SkipList
2759     *p;
2760
2761   size_t
2762     color;
2763
2764   ssize_t
2765     count;
2766
2767   /*
2768     Find the standard-deviation value for each of the color.
2769   */
2770   p=(&pixel_list->skip_list);
2771   color=65536L;
2772   count=0;
2773   sum=0.0;
2774   sum_squared=0.0;
2775   do
2776   {
2777     register ssize_t
2778       i;
2779
2780     color=p->nodes[color].next[0];
2781     sum+=(double) p->nodes[color].count*color;
2782     for (i=0; i < (ssize_t) p->nodes[color].count; i++)
2783       sum_squared+=((double) color)*((double) color);
2784     count+=p->nodes[color].count;
2785   } while (count < (ssize_t) pixel_list->length);
2786   sum/=pixel_list->length;
2787   sum_squared/=pixel_list->length;
2788   *pixel=ScaleShortToQuantum((unsigned short) sqrt(sum_squared-(sum*sum)));
2789 }
2790
2791 static inline void InsertPixelList(const Quantum pixel,PixelList *pixel_list)
2792 {
2793   size_t
2794     signature;
2795
2796   unsigned short
2797     index;
2798
2799   index=ScaleQuantumToShort(pixel);
2800   signature=pixel_list->skip_list.nodes[index].signature;
2801   if (signature == pixel_list->signature)
2802     {
2803       pixel_list->skip_list.nodes[index].count++;
2804       return;
2805     }
2806   AddNodePixelList(pixel_list,index);
2807 }
2808
2809 static void ResetPixelList(PixelList *pixel_list)
2810 {
2811   int
2812     level;
2813
2814   register SkipNode
2815     *root;
2816
2817   register SkipList
2818     *p;
2819
2820   /*
2821     Reset the skip-list.
2822   */
2823   p=(&pixel_list->skip_list);
2824   root=p->nodes+65536UL;
2825   p->level=0;
2826   for (level=0; level < 9; level++)
2827     root->next[level]=65536UL;
2828   pixel_list->seed=pixel_list->signature++;
2829 }
2830
2831 MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
2832   const size_t width,const size_t height,ExceptionInfo *exception)
2833 {
2834 #define StatisticImageTag  "Statistic/Image"
2835
2836   CacheView
2837     *image_view,
2838     *statistic_view;
2839
2840   Image
2841     *statistic_image;
2842
2843   MagickBooleanType
2844     status;
2845
2846   MagickOffsetType
2847     progress;
2848
2849   PixelList
2850     **magick_restrict pixel_list;
2851
2852   ssize_t
2853     center,
2854     y;
2855
2856   /*
2857     Initialize statistics image attributes.
2858   */
2859   assert(image != (Image *) NULL);
2860   assert(image->signature == MagickCoreSignature);
2861   if (image->debug != MagickFalse)
2862     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2863   assert(exception != (ExceptionInfo *) NULL);
2864   assert(exception->signature == MagickCoreSignature);
2865   statistic_image=CloneImage(image,0,0,MagickTrue,
2866     exception);
2867   if (statistic_image == (Image *) NULL)
2868     return((Image *) NULL);
2869   status=SetImageStorageClass(statistic_image,DirectClass,exception);
2870   if (status == MagickFalse)
2871     {
2872       statistic_image=DestroyImage(statistic_image);
2873       return((Image *) NULL);
2874     }
2875   pixel_list=AcquirePixelListThreadSet(MagickMax(width,1),MagickMax(height,1));
2876   if (pixel_list == (PixelList **) NULL)
2877     {
2878       statistic_image=DestroyImage(statistic_image);
2879       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2880     }
2881   /*
2882     Make each pixel the min / max / median / mode / etc. of the neighborhood.
2883   */
2884   center=(ssize_t) GetPixelChannels(image)*(image->columns+MagickMax(width,1))*
2885     (MagickMax(height,1)/2L)+GetPixelChannels(image)*(MagickMax(width,1)/2L);
2886   status=MagickTrue;
2887   progress=0;
2888   image_view=AcquireVirtualCacheView(image,exception);
2889   statistic_view=AcquireAuthenticCacheView(statistic_image,exception);
2890 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2891   #pragma omp parallel for schedule(static) shared(progress,status) \
2892     magick_number_threads(image,statistic_image,statistic_image->rows,1)
2893 #endif
2894   for (y=0; y < (ssize_t) statistic_image->rows; y++)
2895   {
2896     const int
2897       id = GetOpenMPThreadId();
2898
2899     register const Quantum
2900       *magick_restrict p;
2901
2902     register Quantum
2903       *magick_restrict q;
2904
2905     register ssize_t
2906       x;
2907
2908     if (status == MagickFalse)
2909       continue;
2910     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) MagickMax(width,1)/2L),y-
2911       (ssize_t) (MagickMax(height,1)/2L),image->columns+MagickMax(width,1),
2912       MagickMax(height,1),exception);
2913     q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns,      1,exception);
2914     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2915       {
2916         status=MagickFalse;
2917         continue;
2918       }
2919     for (x=0; x < (ssize_t) statistic_image->columns; x++)
2920     {
2921       register ssize_t
2922         i;
2923
2924       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2925       {
2926         Quantum
2927           pixel;
2928
2929         register const Quantum
2930           *magick_restrict pixels;
2931
2932         register ssize_t
2933           u;
2934
2935         ssize_t
2936           v;
2937
2938         PixelChannel channel = GetPixelChannelChannel(image,i);
2939         PixelTrait traits = GetPixelChannelTraits(image,channel);
2940         PixelTrait statistic_traits=GetPixelChannelTraits(statistic_image,
2941           channel);
2942         if ((traits == UndefinedPixelTrait) ||
2943             (statistic_traits == UndefinedPixelTrait))
2944           continue;
2945         if (((statistic_traits & CopyPixelTrait) != 0) ||
2946             (GetPixelWriteMask(image,p) <= (QuantumRange/2)))
2947           {
2948             SetPixelChannel(statistic_image,channel,p[center+i],q);
2949             continue;
2950           }
2951         if ((statistic_traits & UpdatePixelTrait) == 0)
2952           continue;
2953         pixels=p;
2954         ResetPixelList(pixel_list[id]);
2955         for (v=0; v < (ssize_t) MagickMax(height,1); v++)
2956         {
2957           for (u=0; u < (ssize_t) MagickMax(width,1); u++)
2958           {
2959             InsertPixelList(pixels[i],pixel_list[id]);
2960             pixels+=GetPixelChannels(image);
2961           }
2962           pixels+=GetPixelChannels(image)*image->columns;
2963         }
2964         switch (type)
2965         {
2966           case GradientStatistic:
2967           {
2968             double
2969               maximum,
2970               minimum;
2971
2972             GetMinimumPixelList(pixel_list[id],&pixel);
2973             minimum=(double) pixel;
2974             GetMaximumPixelList(pixel_list[id],&pixel);
2975             maximum=(double) pixel;
2976             pixel=ClampToQuantum(MagickAbsoluteValue(maximum-minimum));
2977             break;
2978           }
2979           case MaximumStatistic:
2980           {
2981             GetMaximumPixelList(pixel_list[id],&pixel);
2982             break;
2983           }
2984           case MeanStatistic:
2985           {
2986             GetMeanPixelList(pixel_list[id],&pixel);
2987             break;
2988           }
2989           case MedianStatistic:
2990           default:
2991           {
2992             GetMedianPixelList(pixel_list[id],&pixel);
2993             break;
2994           }
2995           case MinimumStatistic:
2996           {
2997             GetMinimumPixelList(pixel_list[id],&pixel);
2998             break;
2999           }
3000           case ModeStatistic:
3001           {
3002             GetModePixelList(pixel_list[id],&pixel);
3003             break;
3004           }
3005           case NonpeakStatistic:
3006           {
3007             GetNonpeakPixelList(pixel_list[id],&pixel);
3008             break;
3009           }
3010           case RootMeanSquareStatistic:
3011           {
3012             GetRootMeanSquarePixelList(pixel_list[id],&pixel);
3013             break;
3014           }
3015           case StandardDeviationStatistic:
3016           {
3017             GetStandardDeviationPixelList(pixel_list[id],&pixel);
3018             break;
3019           }
3020         }
3021         SetPixelChannel(statistic_image,channel,pixel,q);
3022       }
3023       p+=GetPixelChannels(image);
3024       q+=GetPixelChannels(statistic_image);
3025     }
3026     if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
3027       status=MagickFalse;
3028     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3029       {
3030         MagickBooleanType
3031           proceed;
3032
3033 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3034         #pragma omp critical (MagickCore_StatisticImage)
3035 #endif
3036         proceed=SetImageProgress(image,StatisticImageTag,progress++,
3037           image->rows);
3038         if (proceed == MagickFalse)
3039           status=MagickFalse;
3040       }
3041   }
3042   statistic_view=DestroyCacheView(statistic_view);
3043   image_view=DestroyCacheView(image_view);
3044   pixel_list=DestroyPixelListThreadSet(pixel_list);
3045   if (status == MagickFalse)
3046     statistic_image=DestroyImage(statistic_image);
3047   return(statistic_image);
3048 }