granicus.if.org Git - imagemagick/blob - MagickCore/statistic.c

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