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