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 %                                John Cristy                                  %
  17 %                                 July 1992                                   %
  18 %                                                                             %
  19 %                                                                             %
  20 %  Copyright 1999-2013 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 (GetPixelMask(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 (GetPixelMask(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             (GetPixelMask(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 (GetPixelMask(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 M e a n                                                   %
1131 %                                                                             %
1132 %                                                                             %
1133 %                                                                             %
1134 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1135 %
1136 %  GetImageMean() returns the mean and standard deviation of one or more image
1137 %  channels.
1138 %
1139 %  The format of the GetImageMean method is:
1140 %
1141 %      MagickBooleanType GetImageMean(const Image *image,double *mean,
1142 %        double *standard_deviation,ExceptionInfo *exception)
1143 %
1144 %  A description of each parameter follows:
1145 %
1146 %    o image: the image.
1147 %
1148 %    o mean: the average value in the channel.
1149 %
1150 %    o standard_deviation: the standard deviation of the channel.
1151 %
1152 %    o exception: return any errors or warnings in this structure.
1153 %
1154 */
1155 MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean,
1156   double *standard_deviation,ExceptionInfo *exception)
1157 {
1158   double
1159     area;
1160
1161   ChannelStatistics
1162     *channel_statistics;
1163
1164   register ssize_t
1165     i;
1166
1167   assert(image != (Image *) NULL);
1168   assert(image->signature == MagickSignature);
1169   if (image->debug != MagickFalse)
1170     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1171   channel_statistics=GetImageStatistics(image,exception);
1172   if (channel_statistics == (ChannelStatistics *) NULL)
1173     return(MagickFalse);
1174   area=0.0;
1175   channel_statistics[CompositePixelChannel].mean=0.0;
1176   channel_statistics[CompositePixelChannel].standard_deviation=0.0;
1177   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1178   {
1179     PixelChannel channel=GetPixelChannelChannel(image,i);
1180     PixelTrait traits=GetPixelChannelTraits(image,channel);
1181     if (traits == UndefinedPixelTrait)
1182       continue;
1183     if ((traits & UpdatePixelTrait) == 0)
1184       continue;
1185     channel_statistics[CompositePixelChannel].mean+=channel_statistics[i].mean;
1186     channel_statistics[CompositePixelChannel].standard_deviation+=
1187       channel_statistics[i].variance-channel_statistics[i].mean*
1188       channel_statistics[i].mean;
1189     area++;
1190   }
1191   channel_statistics[CompositePixelChannel].mean/=area;
1192   channel_statistics[CompositePixelChannel].standard_deviation=
1193     sqrt(channel_statistics[CompositePixelChannel].standard_deviation/area);
1194   *mean=channel_statistics[CompositePixelChannel].mean;
1195   *standard_deviation=
1196     channel_statistics[CompositePixelChannel].standard_deviation;
1197   channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1198     channel_statistics);
1199   return(MagickTrue);
1200 }
1201 \f
1202 /*
1203 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1204 %                                                                             %
1205 %                                                                             %
1206 %                                                                             %
1207 %   G e t I m a g e K u r t o s i s                                           %
1208 %                                                                             %
1209 %                                                                             %
1210 %                                                                             %
1211 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1212 %
1213 %  GetImageKurtosis() returns the kurtosis and skewness of one or more
1214 %  image channels.
1215 %
1216 %  The format of the GetImageKurtosis method is:
1217 %
1218 %      MagickBooleanType GetImageKurtosis(const Image *image,double *kurtosis,
1219 %        double *skewness,ExceptionInfo *exception)
1220 %
1221 %  A description of each parameter follows:
1222 %
1223 %    o image: the image.
1224 %
1225 %    o kurtosis: the kurtosis of the channel.
1226 %
1227 %    o skewness: the skewness of the channel.
1228 %
1229 %    o exception: return any errors or warnings in this structure.
1230 %
1231 */
1232 MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
1233   double *kurtosis,double *skewness,ExceptionInfo *exception)
1234 {
1235   CacheView
1236     *image_view;
1237
1238   double
1239     area,
1240     mean,
1241     standard_deviation,
1242     sum_squares,
1243     sum_cubes,
1244     sum_fourth_power;
1245
1246   MagickBooleanType
1247     status;
1248
1249   ssize_t
1250     y;
1251
1252   assert(image != (Image *) NULL);
1253   assert(image->signature == MagickSignature);
1254   if (image->debug != MagickFalse)
1255     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1256   status=MagickTrue;
1257   *kurtosis=0.0;
1258   *skewness=0.0;
1259   area=0.0;
1260   mean=0.0;
1261   standard_deviation=0.0;
1262   sum_squares=0.0;
1263   sum_cubes=0.0;
1264   sum_fourth_power=0.0;
1265   image_view=AcquireVirtualCacheView(image,exception);
1266 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1267   #pragma omp parallel for schedule(static,4) shared(status) \
1268     magick_threads(image,image,image->rows,1)
1269 #endif
1270   for (y=0; y < (ssize_t) image->rows; y++)
1271   {
1272     register const Quantum
1273       *restrict p;
1274
1275     register ssize_t
1276       x;
1277
1278     if (status == MagickFalse)
1279       continue;
1280     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1281     if (p == (const Quantum *) NULL)
1282       {
1283         status=MagickFalse;
1284         continue;
1285       }
1286     for (x=0; x < (ssize_t) image->columns; x++)
1287     {
1288       register ssize_t
1289         i;
1290
1291       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1292       {
1293         PixelChannel channel=GetPixelChannelChannel(image,i);
1294         PixelTrait traits=GetPixelChannelTraits(image,channel);
1295         if (traits == UndefinedPixelTrait)
1296           continue;
1297         if (((traits & UpdatePixelTrait) == 0) ||
1298             (GetPixelMask(image,p) != 0))
1299           continue;
1300 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1301         #pragma omp critical (MagickCore_GetImageKurtosis)
1302 #endif
1303         {
1304           mean+=p[i];
1305           sum_squares+=(double) p[i]*p[i];
1306           sum_cubes+=(double) p[i]*p[i]*p[i];
1307           sum_fourth_power+=(double) p[i]*p[i]*p[i]*p[i];
1308           area++;
1309         }
1310       }
1311       p+=GetPixelChannels(image);
1312     }
1313   }
1314   image_view=DestroyCacheView(image_view);
1315   if (area != 0.0)
1316     {
1317       mean/=area;
1318       sum_squares/=area;
1319       sum_cubes/=area;
1320       sum_fourth_power/=area;
1321     }
1322   standard_deviation=sqrt(sum_squares-(mean*mean));
1323   if (standard_deviation != 0.0)
1324     {
1325       *kurtosis=sum_fourth_power-4.0*mean*sum_cubes+6.0*mean*mean*sum_squares-
1326         3.0*mean*mean*mean*mean;
1327       *kurtosis/=standard_deviation*standard_deviation*standard_deviation*
1328         standard_deviation;
1329       *kurtosis-=3.0;
1330       *skewness=sum_cubes-3.0*mean*sum_squares+2.0*mean*mean*mean;
1331       *skewness/=standard_deviation*standard_deviation*standard_deviation;
1332     }
1333   return(status);
1334 }
1335 \f
1336 /*
1337 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1338 %                                                                             %
1339 %                                                                             %
1340 %                                                                             %
1341 %   G e t I m a g e R a n g e                                                 %
1342 %                                                                             %
1343 %                                                                             %
1344 %                                                                             %
1345 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1346 %
1347 %  GetImageRange() returns the range of one or more image channels.
1348 %
1349 %  The format of the GetImageRange method is:
1350 %
1351 %      MagickBooleanType GetImageRange(const Image *image,double *minima,
1352 %        double *maxima,ExceptionInfo *exception)
1353 %
1354 %  A description of each parameter follows:
1355 %
1356 %    o image: the image.
1357 %
1358 %    o minima: the minimum value in the channel.
1359 %
1360 %    o maxima: the maximum value in the channel.
1361 %
1362 %    o exception: return any errors or warnings in this structure.
1363 %
1364 */
1365 MagickExport MagickBooleanType GetImageRange(const Image *image,double *minima,
1366   double *maxima,ExceptionInfo *exception)
1367 {
1368   CacheView
1369     *image_view;
1370
1371   MagickBooleanType
1372     initialize,
1373     status;
1374
1375   ssize_t
1376     y;
1377
1378   assert(image != (Image *) NULL);
1379   assert(image->signature == MagickSignature);
1380   if (image->debug != MagickFalse)
1381     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1382   status=MagickTrue;
1383   initialize=MagickTrue;
1384   *maxima=0.0;
1385   *minima=0.0;
1386   image_view=AcquireVirtualCacheView(image,exception);
1387 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1388   #pragma omp parallel for schedule(static,4) shared(status,initialize) \
1389     magick_threads(image,image,image->rows,1)
1390 #endif
1391   for (y=0; y < (ssize_t) image->rows; y++)
1392   {
1393     register const Quantum
1394       *restrict p;
1395
1396     register ssize_t
1397       x;
1398
1399     if (status == MagickFalse)
1400       continue;
1401     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1402     if (p == (const Quantum *) NULL)
1403       {
1404         status=MagickFalse;
1405         continue;
1406       }
1407     for (x=0; x < (ssize_t) image->columns; x++)
1408     {
1409       register ssize_t
1410         i;
1411
1412       if (GetPixelMask(image,p) != 0)
1413         {
1414           p+=GetPixelChannels(image);
1415           continue;
1416         }
1417       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1418       {
1419         PixelChannel channel=GetPixelChannelChannel(image,i);
1420         PixelTrait traits=GetPixelChannelTraits(image,channel);
1421         if (traits == UndefinedPixelTrait)
1422           continue;
1423         if ((traits & UpdatePixelTrait) == 0)
1424           continue;
1425 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1426         #pragma omp critical (MagickCore_GetImageRange)
1427 #endif
1428         {
1429           if (initialize != MagickFalse)
1430             {
1431               *minima=(double) p[i];
1432               *maxima=(double) p[i];
1433               initialize=MagickFalse;
1434             }
1435           else
1436             {
1437               if ((double) p[i] < *minima)
1438                 *minima=(double) p[i];
1439               if ((double) p[i] > *maxima)
1440                 *maxima=(double) p[i];
1441            }
1442         }
1443       }
1444       p+=GetPixelChannels(image);
1445     }
1446   }
1447   image_view=DestroyCacheView(image_view);
1448   return(status);
1449 }
1450 \f
1451 /*
1452 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1453 %                                                                             %
1454 %                                                                             %
1455 %                                                                             %
1456 %   G e t I m a g e S t a t i s t i c s                                       %
1457 %                                                                             %
1458 %                                                                             %
1459 %                                                                             %
1460 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1461 %
1462 %  GetImageStatistics() returns statistics for each channel in the image.  The
1463 %  statistics include the channel depth, its minima, maxima, mean, standard
1464 %  deviation, kurtosis and skewness.  You can access the red channel mean, for
1465 %  example, like this:
1466 %
1467 %      channel_statistics=GetImageStatistics(image,exception);
1468 %      red_mean=channel_statistics[RedPixelChannel].mean;
1469 %
1470 %  Use MagickRelinquishMemory() to free the statistics buffer.
1471 %
1472 %  The format of the GetImageStatistics method is:
1473 %
1474 %      ChannelStatistics *GetImageStatistics(const Image *image,
1475 %        ExceptionInfo *exception)
1476 %
1477 %  A description of each parameter follows:
1478 %
1479 %    o image: the image.
1480 %
1481 %    o exception: return any errors or warnings in this structure.
1482 %
1483 */
1484
1485 static size_t GetImageChannels(const Image *image)
1486 {
1487   register ssize_t
1488     i;
1489
1490   size_t
1491     channels;
1492
1493   channels=0;
1494   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1495   {
1496     PixelChannel channel=GetPixelChannelChannel(image,i);
1497     PixelTrait traits=GetPixelChannelTraits(image,channel);
1498     if ((traits & UpdatePixelTrait) != 0)
1499       channels++;
1500   }
1501   return(channels);
1502 }
1503
1504 MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
1505   ExceptionInfo *exception)
1506 {
1507   ChannelStatistics
1508     *channel_statistics;
1509
1510   MagickStatusType
1511     status;
1512
1513   QuantumAny
1514     range;
1515
1516   register ssize_t
1517     i;
1518
1519   size_t
1520     channels,
1521     depth;
1522
1523   ssize_t
1524     y;
1525
1526   assert(image != (Image *) NULL);
1527   assert(image->signature == MagickSignature);
1528   if (image->debug != MagickFalse)
1529     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1530   channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(
1531     MaxPixelChannels+1,sizeof(*channel_statistics));
1532   if (channel_statistics == (ChannelStatistics *) NULL)
1533     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1534   (void) ResetMagickMemory(channel_statistics,0,(MaxPixelChannels+1)*
1535     sizeof(*channel_statistics));
1536   for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
1537   {
1538     channel_statistics[i].depth=1;
1539     channel_statistics[i].maxima=(-MagickHuge);
1540     channel_statistics[i].minima=MagickHuge;
1541   }
1542   for (y=0; y < (ssize_t) image->rows; y++)
1543   {
1544     register const Quantum
1545       *restrict p;
1546
1547     register ssize_t
1548       x;
1549
1550     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1551     if (p == (const Quantum *) NULL)
1552       break;
1553     for (x=0; x < (ssize_t) image->columns; x++)
1554     {
1555       register ssize_t
1556         i;
1557
1558       if (GetPixelMask(image,p) != 0)
1559         {
1560           p+=GetPixelChannels(image);
1561           continue;
1562         }
1563       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1564       {
1565         PixelChannel channel=GetPixelChannelChannel(image,i);
1566         PixelTrait traits=GetPixelChannelTraits(image,channel);
1567         if (traits == UndefinedPixelTrait)
1568           continue;
1569         if (channel_statistics[channel].depth != MAGICKCORE_QUANTUM_DEPTH)
1570           {
1571             depth=channel_statistics[channel].depth;
1572             range=GetQuantumRange(depth);
1573             status=p[i] != ScaleAnyToQuantum(ScaleQuantumToAny(p[i],range),
1574               range) ? MagickTrue : MagickFalse;
1575             if (status != MagickFalse)
1576               {
1577                 channel_statistics[channel].depth++;
1578                 i--;
1579                 continue;
1580               }
1581           }
1582         if ((double) p[i] < channel_statistics[channel].minima)
1583           channel_statistics[channel].minima=(double) p[i];
1584         if ((double) p[i] > channel_statistics[channel].maxima)
1585           channel_statistics[channel].maxima=(double) p[i];
1586         channel_statistics[channel].sum+=p[i];
1587         channel_statistics[channel].sum_squared+=(double) p[i]*p[i];
1588         channel_statistics[channel].sum_cubed+=(double) p[i]*p[i]*p[i];
1589         channel_statistics[channel].sum_fourth_power+=(double) p[i]*p[i]*p[i]*
1590           p[i];
1591         channel_statistics[channel].area++;
1592       }
1593       p+=GetPixelChannels(image);
1594     }
1595   }
1596   for (i=0; i < (ssize_t) MaxPixelChannels; i++)
1597   {
1598     double
1599       area;
1600
1601     area=PerceptibleReciprocal(channel_statistics[i].area);
1602     channel_statistics[i].sum*=area;
1603     channel_statistics[i].sum_squared*=area;
1604     channel_statistics[i].sum_cubed*=area;
1605     channel_statistics[i].sum_fourth_power*=area;
1606     channel_statistics[i].mean=channel_statistics[i].sum;
1607     channel_statistics[i].variance=channel_statistics[i].sum_squared;
1608     channel_statistics[i].standard_deviation=sqrt(
1609       channel_statistics[i].variance-(channel_statistics[i].mean*
1610       channel_statistics[i].mean));
1611   }
1612   for (i=0; i < (ssize_t) MaxPixelChannels; i++)
1613   {
1614     channel_statistics[CompositePixelChannel].depth=(size_t) EvaluateMax(
1615       (double) channel_statistics[CompositePixelChannel].depth,(double)
1616       channel_statistics[i].depth);
1617     channel_statistics[CompositePixelChannel].minima=MagickMin(
1618       channel_statistics[CompositePixelChannel].minima,
1619       channel_statistics[i].minima);
1620     channel_statistics[CompositePixelChannel].maxima=EvaluateMax(
1621       channel_statistics[CompositePixelChannel].maxima,
1622       channel_statistics[i].maxima);
1623     channel_statistics[CompositePixelChannel].sum+=channel_statistics[i].sum;
1624     channel_statistics[CompositePixelChannel].sum_squared+=
1625       channel_statistics[i].sum_squared;
1626     channel_statistics[CompositePixelChannel].sum_cubed+=
1627       channel_statistics[i].sum_cubed;
1628     channel_statistics[CompositePixelChannel].sum_fourth_power+=
1629       channel_statistics[i].sum_fourth_power;
1630     channel_statistics[CompositePixelChannel].mean+=channel_statistics[i].mean;
1631     channel_statistics[CompositePixelChannel].variance+=
1632       channel_statistics[i].variance-channel_statistics[i].mean*
1633       channel_statistics[i].mean;
1634     channel_statistics[CompositePixelChannel].standard_deviation+=
1635       channel_statistics[i].variance-channel_statistics[i].mean*
1636       channel_statistics[i].mean;
1637   }
1638   channels=GetImageChannels(image);
1639   channel_statistics[CompositePixelChannel].sum/=channels;
1640   channel_statistics[CompositePixelChannel].sum_squared/=channels;
1641   channel_statistics[CompositePixelChannel].sum_cubed/=channels;
1642   channel_statistics[CompositePixelChannel].sum_fourth_power/=channels;
1643   channel_statistics[CompositePixelChannel].mean/=channels;
1644   channel_statistics[CompositePixelChannel].variance/=channels;
1645   channel_statistics[CompositePixelChannel].standard_deviation=
1646     sqrt(channel_statistics[CompositePixelChannel].standard_deviation/channels);
1647   channel_statistics[CompositePixelChannel].kurtosis/=channels;
1648   channel_statistics[CompositePixelChannel].skewness/=channels;
1649   for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
1650   {
1651     double
1652       standard_deviation;
1653
1654     standard_deviation=PerceptibleReciprocal(
1655       channel_statistics[i].standard_deviation);
1656     channel_statistics[i].skewness=(channel_statistics[i].sum_cubed-3.0*
1657       channel_statistics[i].mean*channel_statistics[i].sum_squared+2.0*
1658       channel_statistics[i].mean*channel_statistics[i].mean*
1659       channel_statistics[i].mean)*(standard_deviation*standard_deviation*
1660       standard_deviation);
1661     channel_statistics[i].kurtosis=(channel_statistics[i].sum_fourth_power-4.0*
1662       channel_statistics[i].mean*channel_statistics[i].sum_cubed+6.0*
1663       channel_statistics[i].mean*channel_statistics[i].mean*
1664       channel_statistics[i].sum_squared-3.0*channel_statistics[i].mean*
1665       channel_statistics[i].mean*1.0*channel_statistics[i].mean*
1666       channel_statistics[i].mean)*(standard_deviation*standard_deviation*
1667       standard_deviation*standard_deviation)-3.0;
1668   }
1669   return(channel_statistics);
1670 }
1671 \f
1672 /*
1673 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1674 %                                                                             %
1675 %                                                                             %
1676 %                                                                             %
1677 %     P o l y n o m i a l I m a g e                                           %
1678 %                                                                             %
1679 %                                                                             %
1680 %                                                                             %
1681 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1682 %
1683 %  PolynomialImage() returns a new image where each pixel is the sum of the
1684 %  pixels in the image sequence after applying its corresponding terms
1685 %  (coefficient and degree pairs).
1686 %
1687 %  The format of the PolynomialImage method is:
1688 %
1689 %      Image *PolynomialImage(const Image *images,const size_t number_terms,
1690 %        const double *terms,ExceptionInfo *exception)
1691 %
1692 %  A description of each parameter follows:
1693 %
1694 %    o images: the image sequence.
1695 %
1696 %    o number_terms: the number of terms in the list.  The actual list length
1697 %      is 2 x number_terms + 1 (the constant).
1698 %
1699 %    o terms: the list of polynomial coefficients and degree pairs and a
1700 %      constant.
1701 %
1702 %    o exception: return any errors or warnings in this structure.
1703 %
1704 */
1705
1706 MagickExport Image *PolynomialImage(const Image *images,
1707   const size_t number_terms,const double *terms,ExceptionInfo *exception)
1708 {
1709 #define PolynomialImageTag  "Polynomial/Image"
1710
1711   CacheView
1712     *polynomial_view;
1713
1714   Image
1715     *image;
1716
1717   MagickBooleanType
1718     status;
1719
1720   MagickOffsetType
1721     progress;
1722
1723   PixelChannels
1724     **restrict polynomial_pixels;
1725
1726   size_t
1727     number_images;
1728
1729   ssize_t
1730     y;
1731
1732   assert(images != (Image *) NULL);
1733   assert(images->signature == MagickSignature);
1734   if (images->debug != MagickFalse)
1735     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
1736   assert(exception != (ExceptionInfo *) NULL);
1737   assert(exception->signature == MagickSignature);
1738   image=CloneImage(images,images->columns,images->rows,MagickTrue,
1739     exception);
1740   if (image == (Image *) NULL)
1741     return((Image *) NULL);
1742   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
1743     {
1744       image=DestroyImage(image);
1745       return((Image *) NULL);
1746     }
1747   number_images=GetImageListLength(images);
1748   polynomial_pixels=AcquirePixelThreadSet(images,number_images);
1749   if (polynomial_pixels == (PixelChannels **) NULL)
1750     {
1751       image=DestroyImage(image);
1752       (void) ThrowMagickException(exception,GetMagickModule(),
1753         ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
1754       return((Image *) NULL);
1755     }
1756   /*
1757     Polynomial image pixels.
1758   */
1759   status=MagickTrue;
1760   progress=0;
1761   polynomial_view=AcquireAuthenticCacheView(image,exception);
1762 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1763   #pragma omp parallel for schedule(static,4) shared(progress,status) \
1764     magick_threads(image,image,image->rows,1)
1765 #endif
1766   for (y=0; y < (ssize_t) image->rows; y++)
1767   {
1768     CacheView
1769       *image_view;
1770
1771     const Image
1772       *next;
1773
1774     const int
1775       id = GetOpenMPThreadId();
1776
1777     register ssize_t
1778       i,
1779       x;
1780
1781     register PixelChannels
1782       *polynomial_pixel;
1783
1784     register Quantum
1785       *restrict q;
1786
1787     ssize_t
1788       j;
1789
1790     if (status == MagickFalse)
1791       continue;
1792     q=QueueCacheViewAuthenticPixels(polynomial_view,0,y,image->columns,1,
1793       exception);
1794     if (q == (Quantum *) NULL)
1795       {
1796         status=MagickFalse;
1797         continue;
1798       }
1799     polynomial_pixel=polynomial_pixels[id];
1800     for (j=0; j < (ssize_t) image->columns; j++)
1801       for (i=0; i < MaxPixelChannels; i++)
1802         polynomial_pixel[j].channel[i]=0.0;
1803     next=images;
1804     for (j=0; j < (ssize_t) number_images; j++)
1805     {
1806       register const Quantum
1807         *p;
1808
1809       if (j >= (ssize_t) number_terms)
1810         continue;
1811       image_view=AcquireVirtualCacheView(next,exception);
1812       p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1813       if (p == (const Quantum *) NULL)
1814         {
1815           image_view=DestroyCacheView(image_view);
1816           break;
1817         }
1818       for (x=0; x < (ssize_t) image->columns; x++)
1819       {
1820         register ssize_t
1821           i;
1822
1823         if (GetPixelMask(next,p) != 0)
1824           {
1825             p+=GetPixelChannels(next);
1826             continue;
1827           }
1828         for (i=0; i < (ssize_t) GetPixelChannels(next); i++)
1829         {
1830           MagickRealType
1831             coefficient,
1832             degree;
1833
1834           PixelChannel channel=GetPixelChannelChannel(image,i);
1835           PixelTrait traits=GetPixelChannelTraits(next,channel);
1836           PixelTrait polynomial_traits=GetPixelChannelTraits(image,channel);
1837           if ((traits == UndefinedPixelTrait) ||
1838               (polynomial_traits == UndefinedPixelTrait))
1839             continue;
1840           if ((traits & UpdatePixelTrait) == 0)
1841             continue;
1842           coefficient=(MagickRealType) terms[2*i];
1843           degree=(MagickRealType) terms[(i << 1)+1];
1844           polynomial_pixel[x].channel[i]+=coefficient*
1845             pow(QuantumScale*GetPixelChannel(image,channel,p),degree);
1846         }
1847         p+=GetPixelChannels(next);
1848       }
1849       image_view=DestroyCacheView(image_view);
1850       next=GetNextImageInList(next);
1851     }
1852     for (x=0; x < (ssize_t) image->columns; x++)
1853     {
1854       register ssize_t
1855         i;
1856
1857       if (GetPixelMask(image,q) != 0)
1858         {
1859           q+=GetPixelChannels(image);
1860           continue;
1861         }
1862       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1863       {
1864         PixelChannel channel=GetPixelChannelChannel(image,i);
1865         PixelTrait traits=GetPixelChannelTraits(image,channel);
1866         if (traits == UndefinedPixelTrait)
1867           continue;
1868         if ((traits & UpdatePixelTrait) == 0)
1869           continue;
1870         q[i]=ClampToQuantum(QuantumRange*polynomial_pixel[x].channel[i]);
1871       }
1872       q+=GetPixelChannels(image);
1873     }
1874     if (SyncCacheViewAuthenticPixels(polynomial_view,exception) == MagickFalse)
1875       status=MagickFalse;
1876     if (images->progress_monitor != (MagickProgressMonitor) NULL)
1877       {
1878         MagickBooleanType
1879           proceed;
1880
1881 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1882         #pragma omp critical (MagickCore_PolynomialImages)
1883 #endif
1884         proceed=SetImageProgress(images,PolynomialImageTag,progress++,
1885           image->rows);
1886         if (proceed == MagickFalse)
1887           status=MagickFalse;
1888       }
1889   }
1890   polynomial_view=DestroyCacheView(polynomial_view);
1891   polynomial_pixels=DestroyPixelThreadSet(polynomial_pixels);
1892   if (status == MagickFalse)
1893     image=DestroyImage(image);
1894   return(image);
1895 }
1896 \f
1897 /*
1898 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1899 %                                                                             %
1900 %                                                                             %
1901 %                                                                             %
1902 %     S t a t i s t i c I m a g e                                             %
1903 %                                                                             %
1904 %                                                                             %
1905 %                                                                             %
1906 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1907 %
1908 %  StatisticImage() makes each pixel the min / max / median / mode / etc. of
1909 %  the neighborhood of the specified width and height.
1910 %
1911 %  The format of the StatisticImage method is:
1912 %
1913 %      Image *StatisticImage(const Image *image,const StatisticType type,
1914 %        const size_t width,const size_t height,ExceptionInfo *exception)
1915 %
1916 %  A description of each parameter follows:
1917 %
1918 %    o image: the image.
1919 %
1920 %    o type: the statistic type (median, mode, etc.).
1921 %
1922 %    o width: the width of the pixel neighborhood.
1923 %
1924 %    o height: the height of the pixel neighborhood.
1925 %
1926 %    o exception: return any errors or warnings in this structure.
1927 %
1928 */
1929
1930 typedef struct _SkipNode
1931 {
1932   size_t
1933     next[9],
1934     count,
1935     signature;
1936 } SkipNode;
1937
1938 typedef struct _SkipList
1939 {
1940   ssize_t
1941     level;
1942
1943   SkipNode
1944     *nodes;
1945 } SkipList;
1946
1947 typedef struct _PixelList
1948 {
1949   size_t
1950     length,
1951     seed;
1952
1953   SkipList
1954     skip_list;
1955
1956   size_t
1957     signature;
1958 } PixelList;
1959
1960 static PixelList *DestroyPixelList(PixelList *pixel_list)
1961 {
1962   if (pixel_list == (PixelList *) NULL)
1963     return((PixelList *) NULL);
1964   if (pixel_list->skip_list.nodes != (SkipNode *) NULL)
1965     pixel_list->skip_list.nodes=(SkipNode *) RelinquishMagickMemory(
1966       pixel_list->skip_list.nodes);
1967   pixel_list=(PixelList *) RelinquishMagickMemory(pixel_list);
1968   return(pixel_list);
1969 }
1970
1971 static PixelList **DestroyPixelListThreadSet(PixelList **pixel_list)
1972 {
1973   register ssize_t
1974     i;
1975
1976   assert(pixel_list != (PixelList **) NULL);
1977   for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
1978     if (pixel_list[i] != (PixelList *) NULL)
1979       pixel_list[i]=DestroyPixelList(pixel_list[i]);
1980   pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list);
1981   return(pixel_list);
1982 }
1983
1984 static PixelList *AcquirePixelList(const size_t width,const size_t height)
1985 {
1986   PixelList
1987     *pixel_list;
1988
1989   pixel_list=(PixelList *) AcquireMagickMemory(sizeof(*pixel_list));
1990   if (pixel_list == (PixelList *) NULL)
1991     return(pixel_list);
1992   (void) ResetMagickMemory((void *) pixel_list,0,sizeof(*pixel_list));
1993   pixel_list->length=width*height;
1994   pixel_list->skip_list.nodes=(SkipNode *) AcquireQuantumMemory(65537UL,
1995     sizeof(*pixel_list->skip_list.nodes));
1996   if (pixel_list->skip_list.nodes == (SkipNode *) NULL)
1997     return(DestroyPixelList(pixel_list));
1998   (void) ResetMagickMemory(pixel_list->skip_list.nodes,0,65537UL*
1999     sizeof(*pixel_list->skip_list.nodes));
2000   pixel_list->signature=MagickSignature;
2001   return(pixel_list);
2002 }
2003
2004 static PixelList **AcquirePixelListThreadSet(const size_t width,
2005   const size_t height)
2006 {
2007   PixelList
2008     **pixel_list;
2009
2010   register ssize_t
2011     i;
2012
2013   size_t
2014     number_threads;
2015
2016   number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
2017   pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
2018     sizeof(*pixel_list));
2019   if (pixel_list == (PixelList **) NULL)
2020     return((PixelList **) NULL);
2021   (void) ResetMagickMemory(pixel_list,0,number_threads*sizeof(*pixel_list));
2022   for (i=0; i < (ssize_t) number_threads; i++)
2023   {
2024     pixel_list[i]=AcquirePixelList(width,height);
2025     if (pixel_list[i] == (PixelList *) NULL)
2026       return(DestroyPixelListThreadSet(pixel_list));
2027   }
2028   return(pixel_list);
2029 }
2030
2031 static void AddNodePixelList(PixelList *pixel_list,const size_t color)
2032 {
2033   register SkipList
2034     *p;
2035
2036   register ssize_t
2037     level;
2038
2039   size_t
2040     search,
2041     update[9];
2042
2043   /*
2044     Initialize the node.
2045   */
2046   p=(&pixel_list->skip_list);
2047   p->nodes[color].signature=pixel_list->signature;
2048   p->nodes[color].count=1;
2049   /*
2050     Determine where it belongs in the list.
2051   */
2052   search=65536UL;
2053   for (level=p->level; level >= 0; level--)
2054   {
2055     while (p->nodes[search].next[level] < color)
2056       search=p->nodes[search].next[level];
2057     update[level]=search;
2058   }
2059   /*
2060     Generate a pseudo-random level for this node.
2061   */
2062   for (level=0; ; level++)
2063   {
2064     pixel_list->seed=(pixel_list->seed*42893621L)+1L;
2065     if ((pixel_list->seed & 0x300) != 0x300)
2066       break;
2067   }
2068   if (level > 8)
2069     level=8;
2070   if (level > (p->level+2))
2071     level=p->level+2;
2072   /*
2073     If we're raising the list's level, link back to the root node.
2074   */
2075   while (level > p->level)
2076   {
2077     p->level++;
2078     update[p->level]=65536UL;
2079   }
2080   /*
2081     Link the node into the skip-list.
2082   */
2083   do
2084   {
2085     p->nodes[color].next[level]=p->nodes[update[level]].next[level];
2086     p->nodes[update[level]].next[level]=color;
2087   } while (level-- > 0);
2088 }
2089
2090 static inline void GetMaximumPixelList(PixelList *pixel_list,Quantum *pixel)
2091 {
2092   register SkipList
2093     *p;
2094
2095   size_t
2096     color,
2097     maximum;
2098
2099   ssize_t
2100     count;
2101
2102   /*
2103     Find the maximum value for each of the color.
2104   */
2105   p=(&pixel_list->skip_list);
2106   color=65536L;
2107   count=0;
2108   maximum=p->nodes[color].next[0];
2109   do
2110   {
2111     color=p->nodes[color].next[0];
2112     if (color > maximum)
2113       maximum=color;
2114     count+=p->nodes[color].count;
2115   } while (count < (ssize_t) pixel_list->length);
2116   *pixel=ScaleShortToQuantum((unsigned short) maximum);
2117 }
2118
2119 static inline void GetMeanPixelList(PixelList *pixel_list,Quantum *pixel)
2120 {
2121   double
2122     sum;
2123
2124   register SkipList
2125     *p;
2126
2127   size_t
2128     color;
2129
2130   ssize_t
2131     count;
2132
2133   /*
2134     Find the mean value for each of the color.
2135   */
2136   p=(&pixel_list->skip_list);
2137   color=65536L;
2138   count=0;
2139   sum=0.0;
2140   do
2141   {
2142     color=p->nodes[color].next[0];
2143     sum+=(double) p->nodes[color].count*color;
2144     count+=p->nodes[color].count;
2145   } while (count < (ssize_t) pixel_list->length);
2146   sum/=pixel_list->length;
2147   *pixel=ScaleShortToQuantum((unsigned short) sum);
2148 }
2149
2150 static inline void GetMedianPixelList(PixelList *pixel_list,Quantum *pixel)
2151 {
2152   register SkipList
2153     *p;
2154
2155   size_t
2156     color;
2157
2158   ssize_t
2159     count;
2160
2161   /*
2162     Find the median value for each of the color.
2163   */
2164   p=(&pixel_list->skip_list);
2165   color=65536L;
2166   count=0;
2167   do
2168   {
2169     color=p->nodes[color].next[0];
2170     count+=p->nodes[color].count;
2171   } while (count <= (ssize_t) (pixel_list->length >> 1));
2172   *pixel=ScaleShortToQuantum((unsigned short) color);
2173 }
2174
2175 static inline void GetMinimumPixelList(PixelList *pixel_list,Quantum *pixel)
2176 {
2177   register SkipList
2178     *p;
2179
2180   size_t
2181     color,
2182     minimum;
2183
2184   ssize_t
2185     count;
2186
2187   /*
2188     Find the minimum value for each of the color.
2189   */
2190   p=(&pixel_list->skip_list);
2191   count=0;
2192   color=65536UL;
2193   minimum=p->nodes[color].next[0];
2194   do
2195   {
2196     color=p->nodes[color].next[0];
2197     if (color < minimum)
2198       minimum=color;
2199     count+=p->nodes[color].count;
2200   } while (count < (ssize_t) pixel_list->length);
2201   *pixel=ScaleShortToQuantum((unsigned short) minimum);
2202 }
2203
2204 static inline void GetModePixelList(PixelList *pixel_list,Quantum *pixel)
2205 {
2206   register SkipList
2207     *p;
2208
2209   size_t
2210     color,
2211     max_count,
2212     mode;
2213
2214   ssize_t
2215     count;
2216
2217   /*
2218     Make each pixel the 'predominant color' of the specified neighborhood.
2219   */
2220   p=(&pixel_list->skip_list);
2221   color=65536L;
2222   mode=color;
2223   max_count=p->nodes[mode].count;
2224   count=0;
2225   do
2226   {
2227     color=p->nodes[color].next[0];
2228     if (p->nodes[color].count > max_count)
2229       {
2230         mode=color;
2231         max_count=p->nodes[mode].count;
2232       }
2233     count+=p->nodes[color].count;
2234   } while (count < (ssize_t) pixel_list->length);
2235   *pixel=ScaleShortToQuantum((unsigned short) mode);
2236 }
2237
2238 static inline void GetNonpeakPixelList(PixelList *pixel_list,Quantum *pixel)
2239 {
2240   register SkipList
2241     *p;
2242
2243   size_t
2244     color,
2245     next,
2246     previous;
2247
2248   ssize_t
2249     count;
2250
2251   /*
2252     Finds the non peak value for each of the colors.
2253   */
2254   p=(&pixel_list->skip_list);
2255   color=65536L;
2256   next=p->nodes[color].next[0];
2257   count=0;
2258   do
2259   {
2260     previous=color;
2261     color=next;
2262     next=p->nodes[color].next[0];
2263     count+=p->nodes[color].count;
2264   } while (count <= (ssize_t) (pixel_list->length >> 1));
2265   if ((previous == 65536UL) && (next != 65536UL))
2266     color=next;
2267   else
2268     if ((previous != 65536UL) && (next == 65536UL))
2269       color=previous;
2270   *pixel=ScaleShortToQuantum((unsigned short) color);
2271 }
2272
2273 static inline void GetStandardDeviationPixelList(PixelList *pixel_list,
2274   Quantum *pixel)
2275 {
2276   double
2277     sum,
2278     sum_squared;
2279
2280   register SkipList
2281     *p;
2282
2283   size_t
2284     color;
2285
2286   ssize_t
2287     count;
2288
2289   /*
2290     Find the standard-deviation value for each of the color.
2291   */
2292   p=(&pixel_list->skip_list);
2293   color=65536L;
2294   count=0;
2295   sum=0.0;
2296   sum_squared=0.0;
2297   do
2298   {
2299     register ssize_t
2300       i;
2301
2302     color=p->nodes[color].next[0];
2303     sum+=(double) p->nodes[color].count*color;
2304     for (i=0; i < (ssize_t) p->nodes[color].count; i++)
2305       sum_squared+=((double) color)*((double) color);
2306     count+=p->nodes[color].count;
2307   } while (count < (ssize_t) pixel_list->length);
2308   sum/=pixel_list->length;
2309   sum_squared/=pixel_list->length;
2310   *pixel=ScaleShortToQuantum((unsigned short) sqrt(sum_squared-(sum*sum)));
2311 }
2312
2313 static inline void InsertPixelList(const Image *image,const Quantum pixel,
2314   PixelList *pixel_list)
2315 {
2316   size_t
2317     signature;
2318
2319   unsigned short
2320     index;
2321
2322   index=ScaleQuantumToShort(pixel);
2323   signature=pixel_list->skip_list.nodes[index].signature;
2324   if (signature == pixel_list->signature)
2325     {
2326       pixel_list->skip_list.nodes[index].count++;
2327       return;
2328     }
2329   AddNodePixelList(pixel_list,index);
2330 }
2331
2332 static inline double MagickAbsoluteValue(const double x)
2333 {
2334   if (x < 0)
2335     return(-x);
2336   return(x);
2337 }
2338
2339 static inline size_t MagickMax(const size_t x,const size_t y)
2340 {
2341   if (x > y)
2342     return(x);
2343   return(y);
2344 }
2345
2346 static void ResetPixelList(PixelList *pixel_list)
2347 {
2348   int
2349     level;
2350
2351   register SkipNode
2352     *root;
2353
2354   register SkipList
2355     *p;
2356
2357   /*
2358     Reset the skip-list.
2359   */
2360   p=(&pixel_list->skip_list);
2361   root=p->nodes+65536UL;
2362   p->level=0;
2363   for (level=0; level < 9; level++)
2364     root->next[level]=65536UL;
2365   pixel_list->seed=pixel_list->signature++;
2366 }
2367
2368 MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
2369   const size_t width,const size_t height,ExceptionInfo *exception)
2370 {
2371 #define StatisticImageTag  "Statistic/Image"
2372
2373   CacheView
2374     *image_view,
2375     *statistic_view;
2376
2377   Image
2378     *statistic_image;
2379
2380   MagickBooleanType
2381     status;
2382
2383   MagickOffsetType
2384     progress;
2385
2386   PixelList
2387     **restrict pixel_list;
2388
2389   ssize_t
2390     center,
2391     y;
2392
2393   /*
2394     Initialize statistics image attributes.
2395   */
2396   assert(image != (Image *) NULL);
2397   assert(image->signature == MagickSignature);
2398   if (image->debug != MagickFalse)
2399     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2400   assert(exception != (ExceptionInfo *) NULL);
2401   assert(exception->signature == MagickSignature);
2402   statistic_image=CloneImage(image,image->columns,image->rows,MagickTrue,
2403     exception);
2404   if (statistic_image == (Image *) NULL)
2405     return((Image *) NULL);
2406   status=SetImageStorageClass(statistic_image,DirectClass,exception);
2407   if (status == MagickFalse)
2408     {
2409       statistic_image=DestroyImage(statistic_image);
2410       return((Image *) NULL);
2411     }
2412   pixel_list=AcquirePixelListThreadSet(MagickMax(width,1),MagickMax(height,1));
2413   if (pixel_list == (PixelList **) NULL)
2414     {
2415       statistic_image=DestroyImage(statistic_image);
2416       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2417     }
2418   /*
2419     Make each pixel the min / max / median / mode / etc. of the neighborhood.
2420   */
2421   center=(ssize_t) GetPixelChannels(image)*(image->columns+MagickMax(width,1))*
2422     (MagickMax(height,1)/2L)+GetPixelChannels(image)*(MagickMax(width,1)/2L);
2423   status=MagickTrue;
2424   progress=0;
2425   image_view=AcquireVirtualCacheView(image,exception);
2426   statistic_view=AcquireAuthenticCacheView(statistic_image,exception);
2427 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2428   #pragma omp parallel for schedule(static,4) shared(progress,status) \
2429     magick_threads(image,statistic_image,statistic_image->rows,1)
2430 #endif
2431   for (y=0; y < (ssize_t) statistic_image->rows; y++)
2432   {
2433     const int
2434       id = GetOpenMPThreadId();
2435
2436     register const Quantum
2437       *restrict p;
2438
2439     register Quantum
2440       *restrict q;
2441
2442     register ssize_t
2443       x;
2444
2445     if (status == MagickFalse)
2446       continue;
2447     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) MagickMax(width,1)/2L),y-
2448       (ssize_t) (MagickMax(height,1)/2L),image->columns+MagickMax(width,1),
2449       MagickMax(height,1),exception);
2450     q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns,      1,exception);
2451     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2452       {
2453         status=MagickFalse;
2454         continue;
2455       }
2456     for (x=0; x < (ssize_t) statistic_image->columns; x++)
2457     {
2458       register ssize_t
2459         i;
2460
2461       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2462       {
2463         Quantum
2464           pixel;
2465
2466         register const Quantum
2467           *restrict pixels;
2468
2469         register ssize_t
2470           u;
2471
2472         ssize_t
2473           v;
2474
2475         PixelChannel channel=GetPixelChannelChannel(image,i);
2476         PixelTrait traits=GetPixelChannelTraits(image,channel);
2477         PixelTrait statistic_traits=GetPixelChannelTraits(statistic_image,
2478           channel);
2479         if ((traits == UndefinedPixelTrait) ||
2480             (statistic_traits == UndefinedPixelTrait))
2481           continue;
2482         if (((statistic_traits & CopyPixelTrait) != 0) ||
2483             (GetPixelMask(image,p) != 0))
2484           {
2485             SetPixelChannel(statistic_image,channel,p[center+i],q);
2486             continue;
2487           }
2488         pixels=p;
2489         ResetPixelList(pixel_list[id]);
2490         for (v=0; v < (ssize_t) MagickMax(height,1); v++)
2491         {
2492           for (u=0; u < (ssize_t) MagickMax(width,1); u++)
2493           {
2494             InsertPixelList(image,pixels[i],pixel_list[id]);
2495             pixels+=GetPixelChannels(image);
2496           }
2497           pixels+=image->columns*GetPixelChannels(image);
2498         }
2499         switch (type)
2500         {
2501           case GradientStatistic:
2502           {
2503             double
2504               maximum,
2505               minimum;
2506
2507             GetMinimumPixelList(pixel_list[id],&pixel);
2508             minimum=(double) pixel;
2509             GetMaximumPixelList(pixel_list[id],&pixel);
2510             maximum=(double) pixel;
2511             pixel=ClampToQuantum(MagickAbsoluteValue(maximum-minimum));
2512             break;
2513           }
2514           case MaximumStatistic:
2515           {
2516             GetMaximumPixelList(pixel_list[id],&pixel);
2517             break;
2518           }
2519           case MeanStatistic:
2520           {
2521             GetMeanPixelList(pixel_list[id],&pixel);
2522             break;
2523           }
2524           case MedianStatistic:
2525           default:
2526           {
2527             GetMedianPixelList(pixel_list[id],&pixel);
2528             break;
2529           }
2530           case MinimumStatistic:
2531           {
2532             GetMinimumPixelList(pixel_list[id],&pixel);
2533             break;
2534           }
2535           case ModeStatistic:
2536           {
2537             GetModePixelList(pixel_list[id],&pixel);
2538             break;
2539           }
2540           case NonpeakStatistic:
2541           {
2542             GetNonpeakPixelList(pixel_list[id],&pixel);
2543             break;
2544           }
2545           case StandardDeviationStatistic:
2546           {
2547             GetStandardDeviationPixelList(pixel_list[id],&pixel);
2548             break;
2549           }
2550         }
2551         SetPixelChannel(statistic_image,channel,pixel,q);
2552       }
2553       p+=GetPixelChannels(image);
2554       q+=GetPixelChannels(statistic_image);
2555     }
2556     if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
2557       status=MagickFalse;
2558     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2559       {
2560         MagickBooleanType
2561           proceed;
2562
2563 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2564         #pragma omp critical (MagickCore_StatisticImage)
2565 #endif
2566         proceed=SetImageProgress(image,StatisticImageTag,progress++,
2567           image->rows);
2568         if (proceed == MagickFalse)
2569           status=MagickFalse;
2570       }
2571   }
2572   statistic_view=DestroyCacheView(statistic_view);
2573   image_view=DestroyCacheView(image_view);
2574   pixel_list=DestroyPixelListThreadSet(pixel_list);
2575   if (status == MagickFalse)
2576     statistic_image=DestroyImage(statistic_image);
2577   return(statistic_image);
2578 }