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-2012 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/segment.h"
  85 #include "MagickCore/semaphore.h"
  86 #include "MagickCore/signature-private.h"
  87 #include "MagickCore/statistic.h"
  88 #include "MagickCore/string_.h"
  89 #include "MagickCore/thread-private.h"
  90 #include "MagickCore/timer.h"
  91 #include "MagickCore/utility.h"
  92 #include "MagickCore/version.h"
  93 \f
  94 /*
  95 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  96 %                                                                             %
  97 %                                                                             %
  98 %                                                                             %
  99 %     E v a l u a t e I m a g e                                               %
 100 %                                                                             %
 101 %                                                                             %
 102 %                                                                             %
 103 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 104 %
 105 %  EvaluateImage() applies a value to the image with an arithmetic, relational,
 106 %  or logical operator to an image. Use these operations to lighten or darken
 107 %  an image, to increase or decrease contrast in an image, or to produce the
 108 %  "negative" of an image.
 109 %
 110 %  The format of the EvaluateImage method is:
 111 %
 112 %      MagickBooleanType EvaluateImage(Image *image,
 113 %        const MagickEvaluateOperator op,const double value,
 114 %        ExceptionInfo *exception)
 115 %      MagickBooleanType EvaluateImages(Image *images,
 116 %        const MagickEvaluateOperator op,const double value,
 117 %        ExceptionInfo *exception)
 118 %
 119 %  A description of each parameter follows:
 120 %
 121 %    o image: the image.
 122 %
 123 %    o op: A channel op.
 124 %
 125 %    o value: A value value.
 126 %
 127 %    o exception: return any errors or warnings in this structure.
 128 %
 129 */
 130
 131 typedef struct _PixelChannels
 132 {
 133   MagickRealType
 134     channel[CompositePixelChannel];
 135 } PixelChannels;
 136
 137 static PixelChannels **DestroyPixelThreadSet(PixelChannels **pixels)
 138 {
 139   register ssize_t
 140     i;
 141
 142   assert(pixels != (PixelChannels **) NULL);
 143   for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
 144     if (pixels[i] != (PixelChannels *) NULL)
 145       pixels[i]=(PixelChannels *) RelinquishMagickMemory(pixels[i]);
 146   pixels=(PixelChannels **) RelinquishMagickMemory(pixels);
 147   return(pixels);
 148 }
 149
 150 static PixelChannels **AcquirePixelThreadSet(const Image *image,
 151   const size_t number_images)
 152 {
 153   register ssize_t
 154     i;
 155
 156   PixelChannels
 157     **pixels;
 158
 159   size_t
 160     length,
 161     number_threads;
 162
 163   number_threads=GetOpenMPMaximumThreads();
 164   pixels=(PixelChannels **) AcquireQuantumMemory(number_threads,
 165     sizeof(*pixels));
 166   if (pixels == (PixelChannels **) NULL)
 167     return((PixelChannels **) NULL);
 168   (void) ResetMagickMemory(pixels,0,number_threads*sizeof(*pixels));
 169   for (i=0; i < (ssize_t) number_threads; i++)
 170   {
 171     register ssize_t
 172       j;
 173
 174     length=image->columns;
 175     if (length < number_images)
 176       length=number_images;
 177     pixels[i]=(PixelChannels *) AcquireQuantumMemory(length,sizeof(**pixels));
 178     if (pixels[i] == (PixelChannels *) NULL)
 179       return(DestroyPixelThreadSet(pixels));
 180     for (j=0; j < (ssize_t) length; j++)
 181     {
 182       register ssize_t
 183         k;
 184
 185       for (k=0; k < MaxPixelChannels; k++)
 186         pixels[i][j].channel[k]=0.0;
 187     }
 188   }
 189   return(pixels);
 190 }
 191
 192 static inline double EvaluateMax(const double x,const double y)
 193 {
 194   if (x > y)
 195     return(x);
 196   return(y);
 197 }
 198
 199 #if defined(__cplusplus) || defined(c_plusplus)
 200 extern "C" {
 201 #endif
 202
 203 static int IntensityCompare(const void *x,const void *y)
 204 {
 205   const PixelChannels
 206     *color_1,
 207     *color_2;
 208
 209   MagickRealType
 210     distance;
 211
 212   register ssize_t
 213     i;
 214
 215   color_1=(const PixelChannels *) x;
 216   color_2=(const PixelChannels *) y;
 217   distance=0.0;
 218   for (i=0; i < MaxPixelChannels; i++)
 219     distance+=color_1->channel[i]-(MagickRealType) color_2->channel[i];
 220   return(distance < 0 ? -1 : distance > 0 ? 1 : 0);
 221 }
 222
 223 #if defined(__cplusplus) || defined(c_plusplus)
 224 }
 225 #endif
 226
 227 static inline double MagickMin(const double x,const double y)
 228 {
 229   if (x < y)
 230     return(x);
 231   return(y);
 232 }
 233
 234 static MagickRealType ApplyEvaluateOperator(RandomInfo *random_info,
 235   Quantum pixel,const MagickEvaluateOperator op,const MagickRealType value)
 236 {
 237   MagickRealType
 238     result;
 239
 240   result=0.0;
 241   switch (op)
 242   {
 243     case UndefinedEvaluateOperator:
 244       break;
 245     case AbsEvaluateOperator:
 246     {
 247       result=(MagickRealType) fabs((double) (pixel+value));
 248       break;
 249     }
 250     case AddEvaluateOperator:
 251     {
 252       result=(MagickRealType) (pixel+value);
 253       break;
 254     }
 255     case AddModulusEvaluateOperator:
 256     {
 257       /*
 258         This returns a 'floored modulus' of the addition which is a positive
 259         result.  It differs from % or fmod() that returns a 'truncated modulus'
 260         result, where floor() is replaced by trunc() and could return a
 261         negative result (which is clipped).
 262       */
 263       result=pixel+value;
 264       result-=(QuantumRange+1.0)*floor((double) result/(QuantumRange+1.0));
 265       break;
 266     }
 267     case AndEvaluateOperator:
 268     {
 269       result=(MagickRealType) ((size_t) pixel & (size_t) (value+0.5));
 270       break;
 271     }
 272     case CosineEvaluateOperator:
 273     {
 274       result=(MagickRealType) (QuantumRange*(0.5*cos((double) (2.0*MagickPI*
 275         QuantumScale*pixel*value))+0.5));
 276       break;
 277     }
 278     case DivideEvaluateOperator:
 279     {
 280       result=pixel/(value == 0.0 ? 1.0 : value);
 281       break;
 282     }
 283     case ExponentialEvaluateOperator:
 284     {
 285       result=(MagickRealType) (QuantumRange*exp((double) (value*QuantumScale*
 286         pixel)));
 287       break;
 288     }
 289     case GaussianNoiseEvaluateOperator:
 290     {
 291       result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
 292         GaussianNoise,value);
 293       break;
 294     }
 295     case ImpulseNoiseEvaluateOperator:
 296     {
 297       result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
 298         ImpulseNoise,value);
 299       break;
 300     }
 301     case LaplacianNoiseEvaluateOperator:
 302     {
 303       result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
 304         LaplacianNoise,value);
 305       break;
 306     }
 307     case LeftShiftEvaluateOperator:
 308     {
 309       result=(MagickRealType) ((size_t) pixel << (size_t) (value+0.5));
 310       break;
 311     }
 312     case LogEvaluateOperator:
 313     {
 314       result=(MagickRealType) (QuantumRange*log((double) (QuantumScale*value*
 315         pixel+1.0))/log((double) (value+1.0)));
 316       break;
 317     }
 318     case MaxEvaluateOperator:
 319     {
 320       result=(MagickRealType) EvaluateMax((double) pixel,value);
 321       break;
 322     }
 323     case MeanEvaluateOperator:
 324     {
 325       result=(MagickRealType) (pixel+value);
 326       break;
 327     }
 328     case MedianEvaluateOperator:
 329     {
 330       result=(MagickRealType) (pixel+value);
 331       break;
 332     }
 333     case MinEvaluateOperator:
 334     {
 335       result=(MagickRealType) MagickMin((double) pixel,value);
 336       break;
 337     }
 338     case MultiplicativeNoiseEvaluateOperator:
 339     {
 340       result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
 341         MultiplicativeGaussianNoise,value);
 342       break;
 343     }
 344     case MultiplyEvaluateOperator:
 345     {
 346       result=(MagickRealType) (value*pixel);
 347       break;
 348     }
 349     case OrEvaluateOperator:
 350     {
 351       result=(MagickRealType) ((size_t) pixel | (size_t) (value+0.5));
 352       break;
 353     }
 354     case PoissonNoiseEvaluateOperator:
 355     {
 356       result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
 357         PoissonNoise,value);
 358       break;
 359     }
 360     case PowEvaluateOperator:
 361     {
 362       result=(MagickRealType) (QuantumRange*pow((double) (QuantumScale*pixel),
 363         (double) value));
 364       break;
 365     }
 366     case RightShiftEvaluateOperator:
 367     {
 368       result=(MagickRealType) ((size_t) pixel >> (size_t) (value+0.5));
 369       break;
 370     }
 371     case SetEvaluateOperator:
 372     {
 373       result=value;
 374       break;
 375     }
 376     case SineEvaluateOperator:
 377     {
 378       result=(MagickRealType) (QuantumRange*(0.5*sin((double) (2.0*MagickPI*
 379         QuantumScale*pixel*value))+0.5));
 380       break;
 381     }
 382     case SubtractEvaluateOperator:
 383     {
 384       result=(MagickRealType) (pixel-value);
 385       break;
 386     }
 387     case ThresholdEvaluateOperator:
 388     {
 389       result=(MagickRealType) (((MagickRealType) pixel <= value) ? 0 :
 390         QuantumRange);
 391       break;
 392     }
 393     case ThresholdBlackEvaluateOperator:
 394     {
 395       result=(MagickRealType) (((MagickRealType) pixel <= value) ? 0 : pixel);
 396       break;
 397     }
 398     case ThresholdWhiteEvaluateOperator:
 399     {
 400       result=(MagickRealType) (((MagickRealType) pixel > value) ? QuantumRange :
 401         pixel);
 402       break;
 403     }
 404     case UniformNoiseEvaluateOperator:
 405     {
 406       result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
 407         UniformNoise,value);
 408       break;
 409     }
 410     case XorEvaluateOperator:
 411     {
 412       result=(MagickRealType) ((size_t) pixel ^ (size_t) (value+0.5));
 413       break;
 414     }
 415   }
 416   return(result);
 417 }
 418
 419 MagickExport Image *EvaluateImages(const Image *images,
 420   const MagickEvaluateOperator op,ExceptionInfo *exception)
 421 {
 422 #define EvaluateImageTag  "Evaluate/Image"
 423
 424   CacheView
 425     *evaluate_view;
 426
 427   const Image
 428     *next;
 429
 430   Image
 431     *evaluate_image;
 432
 433   MagickBooleanType
 434     status;
 435
 436   MagickOffsetType
 437     progress;
 438
 439   PixelChannels
 440     **restrict evaluate_pixels;
 441
 442   RandomInfo
 443     **restrict random_info;
 444
 445   size_t
 446     number_images;
 447
 448   ssize_t
 449     y;
 450
 451   /*
 452     Ensure the image are the same size.
 453   */
 454   assert(images != (Image *) NULL);
 455   assert(images->signature == MagickSignature);
 456   if (images->debug != MagickFalse)
 457     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
 458   assert(exception != (ExceptionInfo *) NULL);
 459   assert(exception->signature == MagickSignature);
 460   for (next=images; next != (Image *) NULL; next=GetNextImageInList(next))
 461     if ((next->columns != images->columns) || (next->rows != images->rows))
 462       {
 463         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
 464           "ImageWidthsOrHeightsDiffer","`%s'",images->filename);
 465         return((Image *) NULL);
 466       }
 467   /*
 468     Initialize evaluate next attributes.
 469   */
 470   evaluate_image=CloneImage(images,images->columns,images->rows,MagickTrue,
 471     exception);
 472   if (evaluate_image == (Image *) NULL)
 473     return((Image *) NULL);
 474   if (SetImageStorageClass(evaluate_image,DirectClass,exception) == MagickFalse)
 475     {
 476       evaluate_image=DestroyImage(evaluate_image);
 477       return((Image *) NULL);
 478     }
 479   number_images=GetImageListLength(images);
 480   evaluate_pixels=AcquirePixelThreadSet(images,number_images);
 481   if (evaluate_pixels == (PixelChannels **) NULL)
 482     {
 483       evaluate_image=DestroyImage(evaluate_image);
 484       (void) ThrowMagickException(exception,GetMagickModule(),
 485         ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
 486       return((Image *) NULL);
 487     }
 488   /*
 489     Evaluate image pixels.
 490   */
 491   status=MagickTrue;
 492   progress=0;
 493   random_info=AcquireRandomInfoThreadSet();
 494   evaluate_view=AcquireCacheView(evaluate_image);
 495   if (op == MedianEvaluateOperator)
 496     {
 497 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 498       #pragma omp parallel for schedule(dynamic) shared(progress,status)
 499 #endif
 500       for (y=0; y < (ssize_t) evaluate_image->rows; y++)
 501       {
 502         CacheView
 503           *image_view;
 504
 505         const Image
 506           *next;
 507
 508         const int
 509           id = GetOpenMPThreadId();
 510
 511         register PixelChannels
 512           *evaluate_pixel;
 513
 514         register Quantum
 515           *restrict q;
 516
 517         register ssize_t
 518           x;
 519
 520         if (status == MagickFalse)
 521           continue;
 522         q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,
 523           evaluate_image->columns,1,exception);
 524         if (q == (Quantum *) NULL)
 525           {
 526             status=MagickFalse;
 527             continue;
 528           }
 529         evaluate_pixel=evaluate_pixels[id];
 530         for (x=0; x < (ssize_t) evaluate_image->columns; x++)
 531         {
 532           register ssize_t
 533             j,
 534             k;
 535
 536           for (j=0; j < (ssize_t) number_images; j++)
 537             for (k=0; k < MaxPixelChannels; k++)
 538               evaluate_pixel[j].channel[k]=0.0;
 539           next=images;
 540           for (j=0; j < (ssize_t) number_images; j++)
 541           {
 542             register const Quantum
 543               *p;
 544
 545             register ssize_t
 546               i;
 547
 548             image_view=AcquireCacheView(next);
 549             p=GetCacheViewVirtualPixels(image_view,x,y,1,1,exception);
 550             if (p == (const Quantum *) NULL)
 551               {
 552                 image_view=DestroyCacheView(image_view);
 553                 break;
 554               }
 555             for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++)
 556             {
 557               PixelChannel
 558                 channel;
 559
 560               PixelTrait
 561                 evaluate_traits,
 562                 traits;
 563
 564               channel=GetPixelChannelMapChannel(evaluate_image,i);
 565               evaluate_traits=GetPixelChannelMapTraits(evaluate_image,channel);
 566               traits=GetPixelChannelMapTraits(next,channel);
 567               if ((traits == UndefinedPixelTrait) ||
 568                   (evaluate_traits == UndefinedPixelTrait))
 569                 continue;
 570               if ((evaluate_traits & UpdatePixelTrait) == 0)
 571                 continue;
 572               evaluate_pixel[j].channel[i]=ApplyEvaluateOperator(
 573                 random_info[id],GetPixelChannel(evaluate_image,channel,p),op,
 574                 evaluate_pixel[j].channel[i]);
 575             }
 576             image_view=DestroyCacheView(image_view);
 577             next=GetNextImageInList(next);
 578           }
 579           qsort((void *) evaluate_pixel,number_images,sizeof(*evaluate_pixel),
 580             IntensityCompare);
 581           for (k=0; k < (ssize_t) GetPixelChannels(evaluate_image); k++)
 582             q[k]=ClampToQuantum(evaluate_pixel[j/2].channel[k]);
 583           q+=GetPixelChannels(evaluate_image);
 584         }
 585         if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
 586           status=MagickFalse;
 587         if (images->progress_monitor != (MagickProgressMonitor) NULL)
 588           {
 589             MagickBooleanType
 590               proceed;
 591
 592 #if   defined(MAGICKCORE_OPENMP_SUPPORT)
 593             #pragma omp critical (MagickCore_EvaluateImages)
 594 #endif
 595             proceed=SetImageProgress(images,EvaluateImageTag,progress++,
 596               evaluate_image->rows);
 597             if (proceed == MagickFalse)
 598               status=MagickFalse;
 599           }
 600       }
 601     }
 602   else
 603     {
 604 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 605       #pragma omp parallel for schedule(dynamic) shared(progress,status)
 606 #endif
 607       for (y=0; y < (ssize_t) evaluate_image->rows; y++)
 608       {
 609         CacheView
 610           *image_view;
 611
 612         const Image
 613           *next;
 614
 615         const int
 616           id = GetOpenMPThreadId();
 617
 618         register ssize_t
 619           i,
 620           x;
 621
 622         register PixelChannels
 623           *evaluate_pixel;
 624
 625         register Quantum
 626           *restrict q;
 627
 628         ssize_t
 629           j;
 630
 631         if (status == MagickFalse)
 632           continue;
 633         q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,
 634           evaluate_image->columns,1,exception);
 635         if (q == (Quantum *) NULL)
 636           {
 637             status=MagickFalse;
 638             continue;
 639           }
 640         evaluate_pixel=evaluate_pixels[id];
 641         for (j=0; j < (ssize_t) evaluate_image->columns; j++)
 642           for (i=0; i < MaxPixelChannels; i++)
 643             evaluate_pixel[j].channel[i]=0.0;
 644         next=images;
 645         for (j=0; j < (ssize_t) number_images; j++)
 646         {
 647           register const Quantum
 648             *p;
 649
 650           image_view=AcquireCacheView(next);
 651           p=GetCacheViewVirtualPixels(image_view,0,y,next->columns,1,exception);
 652           if (p == (const Quantum *) NULL)
 653             {
 654               image_view=DestroyCacheView(image_view);
 655               break;
 656             }
 657           for (x=0; x < (ssize_t) next->columns; x++)
 658           {
 659             register ssize_t
 660               i;
 661
 662             for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++)
 663             {
 664               PixelChannel
 665                 channel;
 666
 667               PixelTrait
 668                 evaluate_traits,
 669                 traits;
 670
 671               channel=GetPixelChannelMapChannel(evaluate_image,i);
 672               traits=GetPixelChannelMapTraits(next,channel);
 673               evaluate_traits=GetPixelChannelMapTraits(evaluate_image,channel);
 674               if ((traits == UndefinedPixelTrait) ||
 675                   (evaluate_traits == UndefinedPixelTrait))
 676                 continue;
 677               if ((traits & UpdatePixelTrait) == 0)
 678                 continue;
 679               evaluate_pixel[x].channel[i]=ApplyEvaluateOperator(
 680                 random_info[id],GetPixelChannel(evaluate_image,channel,p),j ==
 681                 0 ? AddEvaluateOperator : op,evaluate_pixel[x].channel[i]);
 682             }
 683             p+=GetPixelChannels(next);
 684           }
 685           image_view=DestroyCacheView(image_view);
 686           next=GetNextImageInList(next);
 687         }
 688         for (x=0; x < (ssize_t) evaluate_image->columns; x++)
 689         {
 690           register ssize_t
 691              i;
 692
 693           switch (op)
 694           {
 695             case MeanEvaluateOperator:
 696             {
 697               for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++)
 698                 evaluate_pixel[x].channel[i]/=(MagickRealType) number_images;
 699               break;
 700             }
 701             case MultiplyEvaluateOperator:
 702             {
 703               for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++)
 704               {
 705                 register ssize_t
 706                   j;
 707
 708                 for (j=0; j < (ssize_t) (number_images-1); j++)
 709                   evaluate_pixel[x].channel[i]*=QuantumScale;
 710               }
 711               break;
 712             }
 713             default:
 714               break;
 715           }
 716         }
 717         for (x=0; x < (ssize_t) evaluate_image->columns; x++)
 718         {
 719           register ssize_t
 720             i;
 721
 722           for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++)
 723           {
 724             PixelChannel
 725               channel;
 726
 727             PixelTrait
 728               traits;
 729
 730             channel=GetPixelChannelMapChannel(evaluate_image,i);
 731             traits=GetPixelChannelMapTraits(evaluate_image,channel);
 732             if (traits == UndefinedPixelTrait)
 733               continue;
 734             if ((traits & UpdatePixelTrait) == 0)
 735               continue;
 736             q[i]=ClampToQuantum(evaluate_pixel[x].channel[i]);
 737           }
 738           q+=GetPixelChannels(evaluate_image);
 739         }
 740         if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
 741           status=MagickFalse;
 742         if (images->progress_monitor != (MagickProgressMonitor) NULL)
 743           {
 744             MagickBooleanType
 745               proceed;
 746
 747 #if   defined(MAGICKCORE_OPENMP_SUPPORT)
 748             #pragma omp critical (MagickCore_EvaluateImages)
 749 #endif
 750             proceed=SetImageProgress(images,EvaluateImageTag,progress++,
 751               evaluate_image->rows);
 752             if (proceed == MagickFalse)
 753               status=MagickFalse;
 754           }
 755       }
 756     }
 757   evaluate_view=DestroyCacheView(evaluate_view);
 758   evaluate_pixels=DestroyPixelThreadSet(evaluate_pixels);
 759   random_info=DestroyRandomInfoThreadSet(random_info);
 760   if (status == MagickFalse)
 761     evaluate_image=DestroyImage(evaluate_image);
 762   return(evaluate_image);
 763 }
 764
 765 MagickExport MagickBooleanType EvaluateImage(Image *image,
 766   const MagickEvaluateOperator op,const double value,ExceptionInfo *exception)
 767 {
 768   CacheView
 769     *image_view;
 770
 771   MagickBooleanType
 772     status;
 773
 774   MagickOffsetType
 775     progress;
 776
 777   RandomInfo
 778     **restrict random_info;
 779
 780   ssize_t
 781     y;
 782
 783   assert(image != (Image *) NULL);
 784   assert(image->signature == MagickSignature);
 785   if (image->debug != MagickFalse)
 786     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 787   assert(exception != (ExceptionInfo *) NULL);
 788   assert(exception->signature == MagickSignature);
 789   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
 790     return(MagickFalse);
 791   status=MagickTrue;
 792   progress=0;
 793   random_info=AcquireRandomInfoThreadSet();
 794   image_view=AcquireCacheView(image);
 795 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 796   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
 797 #endif
 798   for (y=0; y < (ssize_t) image->rows; y++)
 799   {
 800     const int
 801       id = GetOpenMPThreadId();
 802
 803     register Quantum
 804       *restrict q;
 805
 806     register ssize_t
 807       x;
 808
 809     if (status == MagickFalse)
 810       continue;
 811     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
 812     if (q == (Quantum *) NULL)
 813       {
 814         status=MagickFalse;
 815         continue;
 816       }
 817     for (x=0; x < (ssize_t) image->columns; x++)
 818     {
 819       register ssize_t
 820         i;
 821
 822       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 823       {
 824         PixelChannel
 825           channel;
 826
 827         PixelTrait
 828           traits;
 829
 830         channel=GetPixelChannelMapChannel(image,i);
 831         traits=GetPixelChannelMapTraits(image,channel);
 832         if (traits == UndefinedPixelTrait)
 833           continue;
 834         q[i]=ClampToQuantum(ApplyEvaluateOperator(random_info[id],q[i],op,
 835           value));
 836       }
 837       q+=GetPixelChannels(image);
 838     }
 839     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
 840       status=MagickFalse;
 841     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 842       {
 843         MagickBooleanType
 844           proceed;
 845
 846 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 847   #pragma omp critical (MagickCore_EvaluateImage)
 848 #endif
 849         proceed=SetImageProgress(image,EvaluateImageTag,progress++,image->rows);
 850         if (proceed == MagickFalse)
 851           status=MagickFalse;
 852       }
 853   }
 854   image_view=DestroyCacheView(image_view);
 855   random_info=DestroyRandomInfoThreadSet(random_info);
 856   return(status);
 857 }
 858 \f
 859 /*
 860 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 861 %                                                                             %
 862 %                                                                             %
 863 %                                                                             %
 864 %     F u n c t i o n I m a g e                                               %
 865 %                                                                             %
 866 %                                                                             %
 867 %                                                                             %
 868 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 869 %
 870 %  FunctionImage() applies a value to the image with an arithmetic, relational,
 871 %  or logical operator to an image. Use these operations to lighten or darken
 872 %  an image, to increase or decrease contrast in an image, or to produce the
 873 %  "negative" of an image.
 874 %
 875 %  The format of the FunctionImage method is:
 876 %
 877 %      MagickBooleanType FunctionImage(Image *image,
 878 %        const MagickFunction function,const ssize_t number_parameters,
 879 %        const double *parameters,ExceptionInfo *exception)
 880 %
 881 %  A description of each parameter follows:
 882 %
 883 %    o image: the image.
 884 %
 885 %    o function: A channel function.
 886 %
 887 %    o parameters: one or more parameters.
 888 %
 889 %    o exception: return any errors or warnings in this structure.
 890 %
 891 */
 892
 893 static Quantum ApplyFunction(Quantum pixel,const MagickFunction function,
 894   const size_t number_parameters,const double *parameters,
 895   ExceptionInfo *exception)
 896 {
 897   MagickRealType
 898     result;
 899
 900   register ssize_t
 901     i;
 902
 903   (void) exception;
 904   result=0.0;
 905   switch (function)
 906   {
 907     case PolynomialFunction:
 908     {
 909       /*
 910         Polynomial: polynomial constants, highest to lowest order (e.g. c0*x^3+
 911         c1*x^2+c2*x+c3).
 912       */
 913       result=0.0;
 914       for (i=0; i < (ssize_t) number_parameters; i++)
 915         result=result*QuantumScale*pixel+parameters[i];
 916       result*=QuantumRange;
 917       break;
 918     }
 919     case SinusoidFunction:
 920     {
 921       MagickRealType
 922         amplitude,
 923         bias,
 924         frequency,
 925         phase;
 926
 927       /*
 928         Sinusoid: frequency, phase, amplitude, bias.
 929       */
 930       frequency=(number_parameters >= 1) ? parameters[0] : 1.0;
 931       phase=(number_parameters >= 2) ? parameters[1] : 0.0;
 932       amplitude=(number_parameters >= 3) ? parameters[2] : 0.5;
 933       bias=(number_parameters >= 4) ? parameters[3] : 0.5;
 934       result=(MagickRealType) (QuantumRange*(amplitude*sin((double) (2.0*
 935         MagickPI*(frequency*QuantumScale*pixel+phase/360.0)))+bias));
 936       break;
 937     }
 938     case ArcsinFunction:
 939     {
 940       MagickRealType
 941         bias,
 942         center,
 943         range,
 944         width;
 945
 946       /*
 947         Arcsin (peged at range limits for invalid results): width, center,
 948         range, and bias.
 949       */
 950       width=(number_parameters >= 1) ? parameters[0] : 1.0;
 951       center=(number_parameters >= 2) ? parameters[1] : 0.5;
 952       range=(number_parameters >= 3) ? parameters[2] : 1.0;
 953       bias=(number_parameters >= 4) ? parameters[3] : 0.5;
 954       result=2.0/width*(QuantumScale*pixel-center);
 955       if ( result <= -1.0 )
 956         result=bias-range/2.0;
 957       else
 958         if (result >= 1.0)
 959           result=bias+range/2.0;
 960         else
 961           result=(MagickRealType) (range/MagickPI*asin((double) result)+bias);
 962       result*=QuantumRange;
 963       break;
 964     }
 965     case ArctanFunction:
 966     {
 967       MagickRealType
 968         center,
 969         bias,
 970         range,
 971         slope;
 972
 973       /*
 974         Arctan: slope, center, range, and bias.
 975       */
 976       slope=(number_parameters >= 1) ? parameters[0] : 1.0;
 977       center=(number_parameters >= 2) ? parameters[1] : 0.5;
 978       range=(number_parameters >= 3) ? parameters[2] : 1.0;
 979       bias=(number_parameters >= 4) ? parameters[3] : 0.5;
 980       result=(MagickRealType) (MagickPI*slope*(QuantumScale*pixel-center));
 981       result=(MagickRealType) (QuantumRange*(range/MagickPI*atan((double)
 982         result)+bias));
 983       break;
 984     }
 985     case UndefinedFunction:
 986       break;
 987   }
 988   return(ClampToQuantum(result));
 989 }
 990
 991 MagickExport MagickBooleanType FunctionImage(Image *image,
 992   const MagickFunction function,const size_t number_parameters,
 993   const double *parameters,ExceptionInfo *exception)
 994 {
 995 #define FunctionImageTag  "Function/Image "
 996
 997   CacheView
 998     *image_view;
 999
1000   MagickBooleanType
1001     status;
1002
1003   MagickOffsetType
1004     progress;
1005
1006   ssize_t
1007     y;
1008
1009   assert(image != (Image *) NULL);
1010   assert(image->signature == MagickSignature);
1011   if (image->debug != MagickFalse)
1012     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1013   assert(exception != (ExceptionInfo *) NULL);
1014   assert(exception->signature == MagickSignature);
1015   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
1016     return(MagickFalse);
1017   status=MagickTrue;
1018   progress=0;
1019   image_view=AcquireCacheView(image);
1020 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1021   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1022 #endif
1023   for (y=0; y < (ssize_t) image->rows; y++)
1024   {
1025     register Quantum
1026       *restrict q;
1027
1028     register ssize_t
1029       x;
1030
1031     if (status == MagickFalse)
1032       continue;
1033     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
1034     if (q == (Quantum *) NULL)
1035       {
1036         status=MagickFalse;
1037         continue;
1038       }
1039     for (x=0; x < (ssize_t) image->columns; x++)
1040     {
1041       register ssize_t
1042         i;
1043
1044       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1045       {
1046         PixelChannel
1047           channel;
1048
1049         PixelTrait
1050           traits;
1051
1052         channel=GetPixelChannelMapChannel(image,i);
1053         traits=GetPixelChannelMapTraits(image,channel);
1054         if (traits == UndefinedPixelTrait)
1055           continue;
1056         if ((traits & UpdatePixelTrait) == 0)
1057           continue;
1058         q[i]=ApplyFunction(q[i],function,number_parameters,parameters,
1059           exception);
1060       }
1061       q+=GetPixelChannels(image);
1062     }
1063     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1064       status=MagickFalse;
1065     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1066       {
1067         MagickBooleanType
1068           proceed;
1069
1070 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1071   #pragma omp critical (MagickCore_FunctionImage)
1072 #endif
1073         proceed=SetImageProgress(image,FunctionImageTag,progress++,image->rows);
1074         if (proceed == MagickFalse)
1075           status=MagickFalse;
1076       }
1077   }
1078   image_view=DestroyCacheView(image_view);
1079   return(status);
1080 }
1081 \f
1082 /*
1083 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1084 %                                                                             %
1085 %                                                                             %
1086 %                                                                             %
1087 %   G e t I m a g e E x t r e m a                                             %
1088 %                                                                             %
1089 %                                                                             %
1090 %                                                                             %
1091 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1092 %
1093 %  GetImageExtrema() returns the extrema of one or more image channels.
1094 %
1095 %  The format of the GetImageExtrema method is:
1096 %
1097 %      MagickBooleanType GetImageExtrema(const Image *image,size_t *minima,
1098 %        size_t *maxima,ExceptionInfo *exception)
1099 %
1100 %  A description of each parameter follows:
1101 %
1102 %    o image: the image.
1103 %
1104 %    o minima: the minimum value in the channel.
1105 %
1106 %    o maxima: the maximum value in the channel.
1107 %
1108 %    o exception: return any errors or warnings in this structure.
1109 %
1110 */
1111 MagickExport MagickBooleanType GetImageExtrema(const Image *image,
1112   size_t *minima,size_t *maxima,ExceptionInfo *exception)
1113 {
1114   double
1115     max,
1116     min;
1117
1118   MagickBooleanType
1119     status;
1120
1121   assert(image != (Image *) NULL);
1122   assert(image->signature == MagickSignature);
1123   if (image->debug != MagickFalse)
1124     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1125   status=GetImageRange(image,&min,&max,exception);
1126   *minima=(size_t) ceil(min-0.5);
1127   *maxima=(size_t) floor(max+0.5);
1128   return(status);
1129 }
1130 \f
1131 /*
1132 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1133 %                                                                             %
1134 %                                                                             %
1135 %                                                                             %
1136 %   G e t I m a g e M e a n                                                   %
1137 %                                                                             %
1138 %                                                                             %
1139 %                                                                             %
1140 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1141 %
1142 %  GetImageMean() returns the mean and standard deviation of one or more
1143 %  image channels.
1144 %
1145 %  The format of the GetImageMean method is:
1146 %
1147 %      MagickBooleanType GetImageMean(const Image *image,double *mean,
1148 %        double *standard_deviation,ExceptionInfo *exception)
1149 %
1150 %  A description of each parameter follows:
1151 %
1152 %    o image: the image.
1153 %
1154 %    o mean: the average value in the channel.
1155 %
1156 %    o standard_deviation: the standard deviation of the channel.
1157 %
1158 %    o exception: return any errors or warnings in this structure.
1159 %
1160 */
1161 MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean,
1162   double *standard_deviation,ExceptionInfo *exception)
1163 {
1164   ChannelStatistics
1165     *channel_statistics;
1166
1167   register ssize_t
1168     i;
1169
1170   size_t
1171     area;
1172
1173   assert(image != (Image *) NULL);
1174   assert(image->signature == MagickSignature);
1175   if (image->debug != MagickFalse)
1176     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1177   channel_statistics=GetImageStatistics(image,exception);
1178   if (channel_statistics == (ChannelStatistics *) NULL)
1179     return(MagickFalse);
1180   area=0;
1181   channel_statistics[CompositePixelChannel].mean=0.0;
1182   channel_statistics[CompositePixelChannel].standard_deviation=0.0;
1183   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1184   {
1185     PixelChannel
1186       channel;
1187
1188     PixelTrait
1189       traits;
1190
1191     channel=GetPixelChannelMapChannel(image,i);
1192     traits=GetPixelChannelMapTraits(image,channel);
1193     if (traits == UndefinedPixelTrait)
1194       continue;
1195     if ((traits & UpdatePixelTrait) == 0)
1196       continue;
1197     channel_statistics[CompositePixelChannel].mean+=channel_statistics[i].mean;
1198     channel_statistics[CompositePixelChannel].standard_deviation+=
1199       channel_statistics[i].variance-channel_statistics[i].mean*
1200       channel_statistics[i].mean;
1201     area++;
1202   }
1203   channel_statistics[CompositePixelChannel].mean/=area;
1204   channel_statistics[CompositePixelChannel].standard_deviation=
1205     sqrt(channel_statistics[CompositePixelChannel].standard_deviation/area);
1206   *mean=channel_statistics[CompositePixelChannel].mean;
1207   *standard_deviation=
1208     channel_statistics[CompositePixelChannel].standard_deviation;
1209   channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1210     channel_statistics);
1211   return(MagickTrue);
1212 }
1213 \f
1214 /*
1215 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1216 %                                                                             %
1217 %                                                                             %
1218 %                                                                             %
1219 %   G e t I m a g e K u r t o s i s                                           %
1220 %                                                                             %
1221 %                                                                             %
1222 %                                                                             %
1223 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1224 %
1225 %  GetImageKurtosis() returns the kurtosis and skewness of one or more
1226 %  image channels.
1227 %
1228 %  The format of the GetImageKurtosis method is:
1229 %
1230 %      MagickBooleanType GetImageKurtosis(const Image *image,double *kurtosis,
1231 %        double *skewness,ExceptionInfo *exception)
1232 %
1233 %  A description of each parameter follows:
1234 %
1235 %    o image: the image.
1236 %
1237 %    o kurtosis: the kurtosis of the channel.
1238 %
1239 %    o skewness: the skewness of the channel.
1240 %
1241 %    o exception: return any errors or warnings in this structure.
1242 %
1243 */
1244 MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
1245   double *kurtosis,double *skewness,ExceptionInfo *exception)
1246 {
1247   CacheView
1248     *image_view;
1249
1250   double
1251     area,
1252     mean,
1253     standard_deviation,
1254     sum_squares,
1255     sum_cubes,
1256     sum_fourth_power;
1257
1258   MagickBooleanType
1259     status;
1260
1261   ssize_t
1262     y;
1263
1264   assert(image != (Image *) NULL);
1265   assert(image->signature == MagickSignature);
1266   if (image->debug != MagickFalse)
1267     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1268   status=MagickTrue;
1269   *kurtosis=0.0;
1270   *skewness=0.0;
1271   area=0.0;
1272   mean=0.0;
1273   standard_deviation=0.0;
1274   sum_squares=0.0;
1275   sum_cubes=0.0;
1276   sum_fourth_power=0.0;
1277   image_view=AcquireCacheView(image);
1278 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1279   #pragma omp parallel for schedule(dynamic) shared(status) omp_throttle(1)
1280 #endif
1281   for (y=0; y < (ssize_t) image->rows; y++)
1282   {
1283     register const Quantum
1284       *restrict p;
1285
1286     register ssize_t
1287       x;
1288
1289     if (status == MagickFalse)
1290       continue;
1291     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1292     if (p == (const Quantum *) NULL)
1293       {
1294         status=MagickFalse;
1295         continue;
1296       }
1297     for (x=0; x < (ssize_t) image->columns; x++)
1298     {
1299       register ssize_t
1300         i;
1301
1302       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1303       {
1304         PixelChannel
1305           channel;
1306
1307         PixelTrait
1308           traits;
1309
1310         channel=GetPixelChannelMapChannel(image,i);
1311         traits=GetPixelChannelMapTraits(image,channel);
1312         if (traits == UndefinedPixelTrait)
1313           continue;
1314         if ((traits & UpdatePixelTrait) == 0)
1315           continue;
1316 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1317         #pragma omp critical (MagickCore_GetImageKurtosis)
1318 #endif
1319         {
1320           mean+=p[i];
1321           sum_squares+=(double) p[i]*p[i];
1322           sum_cubes+=(double) p[i]*p[i]*p[i];
1323           sum_fourth_power+=(double) p[i]*p[i]*p[i]*p[i];
1324           area++;
1325         }
1326       }
1327       p+=GetPixelChannels(image);
1328     }
1329   }
1330   image_view=DestroyCacheView(image_view);
1331   if (area != 0.0)
1332     {
1333       mean/=area;
1334       sum_squares/=area;
1335       sum_cubes/=area;
1336       sum_fourth_power/=area;
1337     }
1338   standard_deviation=sqrt(sum_squares-(mean*mean));
1339   if (standard_deviation != 0.0)
1340     {
1341       *kurtosis=sum_fourth_power-4.0*mean*sum_cubes+6.0*mean*mean*sum_squares-
1342         3.0*mean*mean*mean*mean;
1343       *kurtosis/=standard_deviation*standard_deviation*standard_deviation*
1344         standard_deviation;
1345       *kurtosis-=3.0;
1346       *skewness=sum_cubes-3.0*mean*sum_squares+2.0*mean*mean*mean;
1347       *skewness/=standard_deviation*standard_deviation*standard_deviation;
1348     }
1349   return(status);
1350 }
1351 \f
1352 /*
1353 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1354 %                                                                             %
1355 %                                                                             %
1356 %                                                                             %
1357 %   G e t I m a g e R a n g e                                                 %
1358 %                                                                             %
1359 %                                                                             %
1360 %                                                                             %
1361 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1362 %
1363 %  GetImageRange() returns the range of one or more image channels.
1364 %
1365 %  The format of the GetImageRange method is:
1366 %
1367 %      MagickBooleanType GetImageRange(const Image *image,double *minima,
1368 %        double *maxima,ExceptionInfo *exception)
1369 %
1370 %  A description of each parameter follows:
1371 %
1372 %    o image: the image.
1373 %
1374 %    o minima: the minimum value in the channel.
1375 %
1376 %    o maxima: the maximum value in the channel.
1377 %
1378 %    o exception: return any errors or warnings in this structure.
1379 %
1380 */
1381 MagickExport MagickBooleanType GetImageRange(const Image *image,double *minima,
1382   double *maxima,ExceptionInfo *exception)
1383 {
1384   CacheView
1385     *image_view;
1386
1387   MagickBooleanType
1388     status;
1389
1390   ssize_t
1391     y;
1392
1393   assert(image != (Image *) NULL);
1394   assert(image->signature == MagickSignature);
1395   if (image->debug != MagickFalse)
1396     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1397   status=MagickTrue;
1398   *maxima=(-MagickHuge);
1399   *minima=MagickHuge;
1400   image_view=AcquireCacheView(image);
1401 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1402   #pragma omp parallel for schedule(dynamic) shared(status) omp_throttle(1)
1403 #endif
1404   for (y=0; y < (ssize_t) image->rows; y++)
1405   {
1406     register const Quantum
1407       *restrict p;
1408
1409     register ssize_t
1410       x;
1411
1412     if (status == MagickFalse)
1413       continue;
1414     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1415     if (p == (const Quantum *) NULL)
1416       {
1417         status=MagickFalse;
1418         continue;
1419       }
1420     for (x=0; x < (ssize_t) image->columns; x++)
1421     {
1422       register ssize_t
1423         i;
1424
1425       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1426       {
1427         PixelChannel
1428           channel;
1429
1430         PixelTrait
1431           traits;
1432
1433         channel=GetPixelChannelMapChannel(image,i);
1434         traits=GetPixelChannelMapTraits(image,channel);
1435         if (traits == UndefinedPixelTrait)
1436           continue;
1437         if ((traits & UpdatePixelTrait) == 0)
1438           continue;
1439 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1440         #pragma omp critical (MagickCore_GetImageRange)
1441 #endif
1442         {
1443           if (p[i] < *minima)
1444             *minima=(double) p[i];
1445           if (p[i] > *maxima)
1446             *maxima=(double) p[i];
1447         }
1448       }
1449       p+=GetPixelChannels(image);
1450     }
1451   }
1452   image_view=DestroyCacheView(image_view);
1453   return(status);
1454 }
1455 \f
1456 /*
1457 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1458 %                                                                             %
1459 %                                                                             %
1460 %                                                                             %
1461 %   G e t I m a g e S t a t i s t i c s                                       %
1462 %                                                                             %
1463 %                                                                             %
1464 %                                                                             %
1465 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1466 %
1467 %  GetImageStatistics() returns statistics for each channel in the
1468 %  image.  The statistics include the channel depth, its minima, maxima, mean,
1469 %  standard deviation, kurtosis and skewness.  You can access the red channel
1470 %  mean, for example, like this:
1471 %
1472 %      channel_statistics=GetImageStatistics(image,exception);
1473 %      red_mean=channel_statistics[RedPixelChannel].mean;
1474 %
1475 %  Use MagickRelinquishMemory() to free the statistics buffer.
1476 %
1477 %  The format of the GetImageStatistics method is:
1478 %
1479 %      ChannelStatistics *GetImageStatistics(const Image *image,
1480 %        ExceptionInfo *exception)
1481 %
1482 %  A description of each parameter follows:
1483 %
1484 %    o image: the image.
1485 %
1486 %    o exception: return any errors or warnings in this structure.
1487 %
1488 */
1489
1490 static size_t GetImageChannels(const Image *image)
1491 {
1492   register ssize_t
1493     i;
1494
1495   size_t
1496     channels;
1497
1498   channels=0;
1499   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1500   {
1501     PixelChannel
1502       channel;
1503
1504     PixelTrait
1505       traits;
1506
1507     channel=GetPixelChannelMapChannel(image,i);
1508     traits=GetPixelChannelMapTraits(image,channel);
1509     if ((traits & UpdatePixelTrait) != 0)
1510       channels++;
1511   }
1512   return(channels);
1513 }
1514
1515 MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
1516   ExceptionInfo *exception)
1517 {
1518   ChannelStatistics
1519     *channel_statistics;
1520
1521   double
1522     area;
1523
1524   MagickStatusType
1525     status;
1526
1527   QuantumAny
1528     range;
1529
1530   register ssize_t
1531     i;
1532
1533   size_t
1534     channels,
1535     depth;
1536
1537   ssize_t
1538     y;
1539
1540   assert(image != (Image *) NULL);
1541   assert(image->signature == MagickSignature);
1542   if (image->debug != MagickFalse)
1543     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1544   channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(
1545     MaxPixelChannels+1,sizeof(*channel_statistics));
1546   if (channel_statistics == (ChannelStatistics *) NULL)
1547     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1548   (void) ResetMagickMemory(channel_statistics,0,(MaxPixelChannels+1)*
1549     sizeof(*channel_statistics));
1550   for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
1551   {
1552     channel_statistics[i].depth=1;
1553     channel_statistics[i].maxima=(-MagickHuge);
1554     channel_statistics[i].minima=MagickHuge;
1555   }
1556   for (y=0; y < (ssize_t) image->rows; y++)
1557   {
1558     register const Quantum
1559       *restrict p;
1560
1561     register ssize_t
1562       x;
1563
1564     p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1565     if (p == (const Quantum *) NULL)
1566       break;
1567     for (x=0; x < (ssize_t) image->columns; x++)
1568     {
1569       register ssize_t
1570         i;
1571
1572       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1573       {
1574         PixelChannel
1575           channel;
1576
1577         PixelTrait
1578           traits;
1579
1580         channel=GetPixelChannelMapChannel(image,i);
1581         traits=GetPixelChannelMapTraits(image,channel);
1582         if (traits == UndefinedPixelTrait)
1583           continue;
1584         if (channel_statistics[i].depth != MAGICKCORE_QUANTUM_DEPTH)
1585           {
1586             depth=channel_statistics[i].depth;
1587             range=GetQuantumRange(depth);
1588             status=p[i] != ScaleAnyToQuantum(ScaleQuantumToAny(p[i],range),
1589               range) ? MagickTrue : MagickFalse;
1590             if (status != MagickFalse)
1591               {
1592                 channel_statistics[channel].depth++;
1593                 continue;
1594               }
1595           }
1596         if ((double) p[i] < channel_statistics[channel].minima)
1597           channel_statistics[channel].minima=(double) p[i];
1598         if ((double) p[i] > channel_statistics[channel].maxima)
1599           channel_statistics[channel].maxima=(double) p[i];
1600         channel_statistics[channel].sum+=p[i];
1601         channel_statistics[channel].sum_squared+=(double) p[i]*p[i];
1602         channel_statistics[channel].sum_cubed+=(double) p[i]*p[i]*p[i];
1603         channel_statistics[channel].sum_fourth_power+=(double) p[i]*p[i]*p[i]*
1604           p[i];
1605       }
1606       p+=GetPixelChannels(image);
1607     }
1608   }
1609   area=(double) image->columns*image->rows;
1610   for (i=0; i < (ssize_t) MaxPixelChannels; i++)
1611   {
1612     channel_statistics[i].sum/=area;
1613     channel_statistics[i].sum_squared/=area;
1614     channel_statistics[i].sum_cubed/=area;
1615     channel_statistics[i].sum_fourth_power/=area;
1616     channel_statistics[i].mean=channel_statistics[i].sum;
1617     channel_statistics[i].variance=channel_statistics[i].sum_squared;
1618     channel_statistics[i].standard_deviation=sqrt(
1619       channel_statistics[i].variance-(channel_statistics[i].mean*
1620       channel_statistics[i].mean));
1621   }
1622   for (i=0; i < (ssize_t) MaxPixelChannels; i++)
1623   {
1624     channel_statistics[CompositePixelChannel].depth=(size_t) EvaluateMax(
1625       (double) channel_statistics[CompositePixelChannel].depth,(double)
1626       channel_statistics[i].depth);
1627     channel_statistics[CompositePixelChannel].minima=MagickMin(
1628       channel_statistics[CompositePixelChannel].minima,
1629       channel_statistics[i].minima);
1630     channel_statistics[CompositePixelChannel].maxima=EvaluateMax(
1631       channel_statistics[CompositePixelChannel].maxima,
1632       channel_statistics[i].maxima);
1633     channel_statistics[CompositePixelChannel].sum+=channel_statistics[i].sum;
1634     channel_statistics[CompositePixelChannel].sum_squared+=
1635       channel_statistics[i].sum_squared;
1636     channel_statistics[CompositePixelChannel].sum_cubed+=
1637       channel_statistics[i].sum_cubed;
1638     channel_statistics[CompositePixelChannel].sum_fourth_power+=
1639       channel_statistics[i].sum_fourth_power;
1640     channel_statistics[CompositePixelChannel].mean+=channel_statistics[i].mean;
1641     channel_statistics[CompositePixelChannel].variance+=
1642       channel_statistics[i].variance-channel_statistics[i].mean*
1643       channel_statistics[i].mean;
1644     channel_statistics[CompositePixelChannel].standard_deviation+=
1645       channel_statistics[i].variance-channel_statistics[i].mean*
1646       channel_statistics[i].mean;
1647   }
1648   channels=GetImageChannels(image);
1649   channel_statistics[CompositePixelChannel].sum/=channels;
1650   channel_statistics[CompositePixelChannel].sum_squared/=channels;
1651   channel_statistics[CompositePixelChannel].sum_cubed/=channels;
1652   channel_statistics[CompositePixelChannel].sum_fourth_power/=channels;
1653   channel_statistics[CompositePixelChannel].mean/=channels;
1654   channel_statistics[CompositePixelChannel].variance/=channels;
1655   channel_statistics[CompositePixelChannel].standard_deviation=
1656     sqrt(channel_statistics[CompositePixelChannel].standard_deviation/channels);
1657   channel_statistics[CompositePixelChannel].kurtosis/=channels;
1658   channel_statistics[CompositePixelChannel].skewness/=channels;
1659   for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
1660   {
1661     if (channel_statistics[i].standard_deviation == 0.0)
1662       continue;
1663     channel_statistics[i].skewness=(channel_statistics[i].sum_cubed-3.0*
1664       channel_statistics[i].mean*channel_statistics[i].sum_squared+2.0*
1665       channel_statistics[i].mean*channel_statistics[i].mean*
1666       channel_statistics[i].mean)/(channel_statistics[i].standard_deviation*
1667       channel_statistics[i].standard_deviation*
1668       channel_statistics[i].standard_deviation);
1669     channel_statistics[i].kurtosis=(channel_statistics[i].sum_fourth_power-4.0*
1670       channel_statistics[i].mean*channel_statistics[i].sum_cubed+6.0*
1671       channel_statistics[i].mean*channel_statistics[i].mean*
1672       channel_statistics[i].sum_squared-3.0*channel_statistics[i].mean*
1673       channel_statistics[i].mean*1.0*channel_statistics[i].mean*
1674       channel_statistics[i].mean)/(channel_statistics[i].standard_deviation*
1675       channel_statistics[i].standard_deviation*
1676       channel_statistics[i].standard_deviation*
1677       channel_statistics[i].standard_deviation)-3.0;
1678   }
1679   return(channel_statistics);
1680 }
1681 \f
1682 /*
1683 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1684 %                                                                             %
1685 %                                                                             %
1686 %                                                                             %
1687 %     S t a t i s t i c I m a g e                                             %
1688 %                                                                             %
1689 %                                                                             %
1690 %                                                                             %
1691 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1692 %
1693 %  StatisticImage() makes each pixel the min / max / median / mode / etc. of
1694 %  the neighborhood of the specified width and height.
1695 %
1696 %  The format of the StatisticImage method is:
1697 %
1698 %      Image *StatisticImage(const Image *image,const StatisticType type,
1699 %        const size_t width,const size_t height,ExceptionInfo *exception)
1700 %
1701 %  A description of each parameter follows:
1702 %
1703 %    o image: the image.
1704 %
1705 %    o type: the statistic type (median, mode, etc.).
1706 %
1707 %    o width: the width of the pixel neighborhood.
1708 %
1709 %    o height: the height of the pixel neighborhood.
1710 %
1711 %    o exception: return any errors or warnings in this structure.
1712 %
1713 */
1714
1715 typedef struct _SkipNode
1716 {
1717   size_t
1718     next[9],
1719     count,
1720     signature;
1721 } SkipNode;
1722
1723 typedef struct _SkipList
1724 {
1725   ssize_t
1726     level;
1727
1728   SkipNode
1729     *nodes;
1730 } SkipList;
1731
1732 typedef struct _PixelList
1733 {
1734   size_t
1735     length,
1736     seed;
1737
1738   SkipList
1739     skip_list;
1740
1741   size_t
1742     signature;
1743 } PixelList;
1744
1745 static PixelList *DestroyPixelList(PixelList *pixel_list)
1746 {
1747   if (pixel_list == (PixelList *) NULL)
1748     return((PixelList *) NULL);
1749   if (pixel_list->skip_list.nodes != (SkipNode *) NULL)
1750     pixel_list->skip_list.nodes=(SkipNode *) RelinquishMagickMemory(
1751       pixel_list->skip_list.nodes);
1752   pixel_list=(PixelList *) RelinquishMagickMemory(pixel_list);
1753   return(pixel_list);
1754 }
1755
1756 static PixelList **DestroyPixelListThreadSet(PixelList **pixel_list)
1757 {
1758   register ssize_t
1759     i;
1760
1761   assert(pixel_list != (PixelList **) NULL);
1762   for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
1763     if (pixel_list[i] != (PixelList *) NULL)
1764       pixel_list[i]=DestroyPixelList(pixel_list[i]);
1765   pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list);
1766   return(pixel_list);
1767 }
1768
1769 static PixelList *AcquirePixelList(const size_t width,const size_t height)
1770 {
1771   PixelList
1772     *pixel_list;
1773
1774   pixel_list=(PixelList *) AcquireMagickMemory(sizeof(*pixel_list));
1775   if (pixel_list == (PixelList *) NULL)
1776     return(pixel_list);
1777   (void) ResetMagickMemory((void *) pixel_list,0,sizeof(*pixel_list));
1778   pixel_list->length=width*height;
1779   pixel_list->skip_list.nodes=(SkipNode *) AcquireQuantumMemory(65537UL,
1780     sizeof(*pixel_list->skip_list.nodes));
1781   if (pixel_list->skip_list.nodes == (SkipNode *) NULL)
1782     return(DestroyPixelList(pixel_list));
1783   (void) ResetMagickMemory(pixel_list->skip_list.nodes,0,65537UL*
1784     sizeof(*pixel_list->skip_list.nodes));
1785   pixel_list->signature=MagickSignature;
1786   return(pixel_list);
1787 }
1788
1789 static PixelList **AcquirePixelListThreadSet(const size_t width,
1790   const size_t height)
1791 {
1792   PixelList
1793     **pixel_list;
1794
1795   register ssize_t
1796     i;
1797
1798   size_t
1799     number_threads;
1800
1801   number_threads=GetOpenMPMaximumThreads();
1802   pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
1803     sizeof(*pixel_list));
1804   if (pixel_list == (PixelList **) NULL)
1805     return((PixelList **) NULL);
1806   (void) ResetMagickMemory(pixel_list,0,number_threads*sizeof(*pixel_list));
1807   for (i=0; i < (ssize_t) number_threads; i++)
1808   {
1809     pixel_list[i]=AcquirePixelList(width,height);
1810     if (pixel_list[i] == (PixelList *) NULL)
1811       return(DestroyPixelListThreadSet(pixel_list));
1812   }
1813   return(pixel_list);
1814 }
1815
1816 static void AddNodePixelList(PixelList *pixel_list,const size_t color)
1817 {
1818   register SkipList
1819     *p;
1820
1821   register ssize_t
1822     level;
1823
1824   size_t
1825     search,
1826     update[9];
1827
1828   /*
1829     Initialize the node.
1830   */
1831   p=(&pixel_list->skip_list);
1832   p->nodes[color].signature=pixel_list->signature;
1833   p->nodes[color].count=1;
1834   /*
1835     Determine where it belongs in the list.
1836   */
1837   search=65536UL;
1838   for (level=p->level; level >= 0; level--)
1839   {
1840     while (p->nodes[search].next[level] < color)
1841       search=p->nodes[search].next[level];
1842     update[level]=search;
1843   }
1844   /*
1845     Generate a pseudo-random level for this node.
1846   */
1847   for (level=0; ; level++)
1848   {
1849     pixel_list->seed=(pixel_list->seed*42893621L)+1L;
1850     if ((pixel_list->seed & 0x300) != 0x300)
1851       break;
1852   }
1853   if (level > 8)
1854     level=8;
1855   if (level > (p->level+2))
1856     level=p->level+2;
1857   /*
1858     If we're raising the list's level, link back to the root node.
1859   */
1860   while (level > p->level)
1861   {
1862     p->level++;
1863     update[p->level]=65536UL;
1864   }
1865   /*
1866     Link the node into the skip-list.
1867   */
1868   do
1869   {
1870     p->nodes[color].next[level]=p->nodes[update[level]].next[level];
1871     p->nodes[update[level]].next[level]=color;
1872   } while (level-- > 0);
1873 }
1874
1875 static inline void GetMaximumPixelList(PixelList *pixel_list,Quantum *pixel)
1876 {
1877   register SkipList
1878     *p;
1879
1880   size_t
1881     color,
1882     maximum;
1883
1884   ssize_t
1885     count;
1886
1887   /*
1888     Find the maximum value for each of the color.
1889   */
1890   p=(&pixel_list->skip_list);
1891   color=65536L;
1892   count=0;
1893   maximum=p->nodes[color].next[0];
1894   do
1895   {
1896     color=p->nodes[color].next[0];
1897     if (color > maximum)
1898       maximum=color;
1899     count+=p->nodes[color].count;
1900   } while (count < (ssize_t) pixel_list->length);
1901   *pixel=ScaleShortToQuantum((unsigned short) maximum);
1902 }
1903
1904 static inline void GetMeanPixelList(PixelList *pixel_list,Quantum *pixel)
1905 {
1906   MagickRealType
1907     sum;
1908
1909   register SkipList
1910     *p;
1911
1912   size_t
1913     color;
1914
1915   ssize_t
1916     count;
1917
1918   /*
1919     Find the mean value for each of the color.
1920   */
1921   p=(&pixel_list->skip_list);
1922   color=65536L;
1923   count=0;
1924   sum=0.0;
1925   do
1926   {
1927     color=p->nodes[color].next[0];
1928     sum+=(MagickRealType) p->nodes[color].count*color;
1929     count+=p->nodes[color].count;
1930   } while (count < (ssize_t) pixel_list->length);
1931   sum/=pixel_list->length;
1932   *pixel=ScaleShortToQuantum((unsigned short) sum);
1933 }
1934
1935 static inline void GetMedianPixelList(PixelList *pixel_list,Quantum *pixel)
1936 {
1937   register SkipList
1938     *p;
1939
1940   size_t
1941     color;
1942
1943   ssize_t
1944     count;
1945
1946   /*
1947     Find the median value for each of the color.
1948   */
1949   p=(&pixel_list->skip_list);
1950   color=65536L;
1951   count=0;
1952   do
1953   {
1954     color=p->nodes[color].next[0];
1955     count+=p->nodes[color].count;
1956   } while (count <= (ssize_t) (pixel_list->length >> 1));
1957   *pixel=ScaleShortToQuantum((unsigned short) color);
1958 }
1959
1960 static inline void GetMinimumPixelList(PixelList *pixel_list,Quantum *pixel)
1961 {
1962   register SkipList
1963     *p;
1964
1965   size_t
1966     color,
1967     minimum;
1968
1969   ssize_t
1970     count;
1971
1972   /*
1973     Find the minimum value for each of the color.
1974   */
1975   p=(&pixel_list->skip_list);
1976   count=0;
1977   color=65536UL;
1978   minimum=p->nodes[color].next[0];
1979   do
1980   {
1981     color=p->nodes[color].next[0];
1982     if (color < minimum)
1983       minimum=color;
1984     count+=p->nodes[color].count;
1985   } while (count < (ssize_t) pixel_list->length);
1986   *pixel=ScaleShortToQuantum((unsigned short) minimum);
1987 }
1988
1989 static inline void GetModePixelList(PixelList *pixel_list,Quantum *pixel)
1990 {
1991   register SkipList
1992     *p;
1993
1994   size_t
1995     color,
1996     max_count,
1997     mode;
1998
1999   ssize_t
2000     count;
2001
2002   /*
2003     Make each pixel the 'predominant color' of the specified neighborhood.
2004   */
2005   p=(&pixel_list->skip_list);
2006   color=65536L;
2007   mode=color;
2008   max_count=p->nodes[mode].count;
2009   count=0;
2010   do
2011   {
2012     color=p->nodes[color].next[0];
2013     if (p->nodes[color].count > max_count)
2014       {
2015         mode=color;
2016         max_count=p->nodes[mode].count;
2017       }
2018     count+=p->nodes[color].count;
2019   } while (count < (ssize_t) pixel_list->length);
2020   *pixel=ScaleShortToQuantum((unsigned short) mode);
2021 }
2022
2023 static inline void GetNonpeakPixelList(PixelList *pixel_list,Quantum *pixel)
2024 {
2025   register SkipList
2026     *p;
2027
2028   size_t
2029     color,
2030     next,
2031     previous;
2032
2033   ssize_t
2034     count;
2035
2036   /*
2037     Finds the non peak value for each of the colors.
2038   */
2039   p=(&pixel_list->skip_list);
2040   color=65536L;
2041   next=p->nodes[color].next[0];
2042   count=0;
2043   do
2044   {
2045     previous=color;
2046     color=next;
2047     next=p->nodes[color].next[0];
2048     count+=p->nodes[color].count;
2049   } while (count <= (ssize_t) (pixel_list->length >> 1));
2050   if ((previous == 65536UL) && (next != 65536UL))
2051     color=next;
2052   else
2053     if ((previous != 65536UL) && (next == 65536UL))
2054       color=previous;
2055   *pixel=ScaleShortToQuantum((unsigned short) color);
2056 }
2057
2058 static inline void GetStandardDeviationPixelList(PixelList *pixel_list,
2059   Quantum *pixel)
2060 {
2061   MagickRealType
2062     sum,
2063     sum_squared;
2064
2065   register SkipList
2066     *p;
2067
2068   size_t
2069     color;
2070
2071   ssize_t
2072     count;
2073
2074   /*
2075     Find the standard-deviation value for each of the color.
2076   */
2077   p=(&pixel_list->skip_list);
2078   color=65536L;
2079   count=0;
2080   sum=0.0;
2081   sum_squared=0.0;
2082   do
2083   {
2084     register ssize_t
2085       i;
2086
2087     color=p->nodes[color].next[0];
2088     sum+=(MagickRealType) p->nodes[color].count*color;
2089     for (i=0; i < (ssize_t) p->nodes[color].count; i++)
2090       sum_squared+=((MagickRealType) color)*((MagickRealType) color);
2091     count+=p->nodes[color].count;
2092   } while (count < (ssize_t) pixel_list->length);
2093   sum/=pixel_list->length;
2094   sum_squared/=pixel_list->length;
2095   *pixel=ScaleShortToQuantum((unsigned short) sqrt(sum_squared-(sum*sum)));
2096 }
2097
2098 static inline void InsertPixelList(const Image *image,const Quantum pixel,
2099   PixelList *pixel_list)
2100 {
2101   size_t
2102     signature;
2103
2104   unsigned short
2105     index;
2106
2107   index=ScaleQuantumToShort(pixel);
2108   signature=pixel_list->skip_list.nodes[index].signature;
2109   if (signature == pixel_list->signature)
2110     {
2111       pixel_list->skip_list.nodes[index].count++;
2112       return;
2113     }
2114   AddNodePixelList(pixel_list,index);
2115 }
2116
2117 static inline MagickRealType MagickAbsoluteValue(const MagickRealType x)
2118 {
2119   if (x < 0)
2120     return(-x);
2121   return(x);
2122 }
2123
2124 static inline size_t MagickMax(const size_t x,const size_t y)
2125 {
2126   if (x > y)
2127     return(x);
2128   return(y);
2129 }
2130
2131 static void ResetPixelList(PixelList *pixel_list)
2132 {
2133   int
2134     level;
2135
2136   register SkipNode
2137     *root;
2138
2139   register SkipList
2140     *p;
2141
2142   /*
2143     Reset the skip-list.
2144   */
2145   p=(&pixel_list->skip_list);
2146   root=p->nodes+65536UL;
2147   p->level=0;
2148   for (level=0; level < 9; level++)
2149     root->next[level]=65536UL;
2150   pixel_list->seed=pixel_list->signature++;
2151 }
2152
2153 MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
2154   const size_t width,const size_t height,ExceptionInfo *exception)
2155 {
2156 #define StatisticImageTag  "Statistic/Image"
2157
2158   CacheView
2159     *image_view,
2160     *statistic_view;
2161
2162   Image
2163     *statistic_image;
2164
2165   MagickBooleanType
2166     status;
2167
2168   MagickOffsetType
2169     progress;
2170
2171   PixelList
2172     **restrict pixel_list;
2173
2174   ssize_t
2175     center,
2176     y;
2177
2178   /*
2179     Initialize statistics image attributes.
2180   */
2181   assert(image != (Image *) NULL);
2182   assert(image->signature == MagickSignature);
2183   if (image->debug != MagickFalse)
2184     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2185   assert(exception != (ExceptionInfo *) NULL);
2186   assert(exception->signature == MagickSignature);
2187   statistic_image=CloneImage(image,image->columns,image->rows,MagickTrue,
2188     exception);
2189   if (statistic_image == (Image *) NULL)
2190     return((Image *) NULL);
2191   status=SetImageStorageClass(statistic_image,DirectClass,exception);
2192   if (status == MagickFalse)
2193     {
2194       statistic_image=DestroyImage(statistic_image);
2195       return((Image *) NULL);
2196     }
2197   pixel_list=AcquirePixelListThreadSet(MagickMax(width,1),MagickMax(height,1));
2198   if (pixel_list == (PixelList **) NULL)
2199     {
2200       statistic_image=DestroyImage(statistic_image);
2201       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2202     }
2203   /*
2204     Make each pixel the min / max / median / mode / etc. of the neighborhood.
2205   */
2206   center=(ssize_t) GetPixelChannels(image)*(image->columns+MagickMax(width,1))*
2207     (MagickMax(height,1)/2L)+GetPixelChannels(image)*(MagickMax(width,1)/2L);
2208   status=MagickTrue;
2209   progress=0;
2210   image_view=AcquireCacheView(image);
2211   statistic_view=AcquireCacheView(statistic_image);
2212 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2213   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2214 #endif
2215   for (y=0; y < (ssize_t) statistic_image->rows; y++)
2216   {
2217     const int
2218       id = GetOpenMPThreadId();
2219
2220     register const Quantum
2221       *restrict p;
2222
2223     register Quantum
2224       *restrict q;
2225
2226     register ssize_t
2227       x;
2228
2229     if (status == MagickFalse)
2230       continue;
2231     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) MagickMax(width,1)/2L),y-
2232       (ssize_t) (MagickMax(height,1)/2L),image->columns+MagickMax(width,1),
2233       MagickMax(height,1),exception);
2234     q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns,      1,exception);
2235     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2236       {
2237         status=MagickFalse;
2238         continue;
2239       }
2240     for (x=0; x < (ssize_t) statistic_image->columns; x++)
2241     {
2242       register ssize_t
2243         i;
2244
2245       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2246       {
2247         PixelChannel
2248           channel;
2249
2250         PixelTrait
2251           statistic_traits,
2252           traits;
2253
2254         Quantum
2255           pixel;
2256
2257         register const Quantum
2258           *restrict pixels;
2259
2260         register ssize_t
2261           u;
2262
2263         ssize_t
2264           v;
2265
2266         channel=GetPixelChannelMapChannel(image,i);
2267         traits=GetPixelChannelMapTraits(image,channel);
2268         statistic_traits=GetPixelChannelMapTraits(statistic_image,channel);
2269         if ((traits == UndefinedPixelTrait) ||
2270             (statistic_traits == UndefinedPixelTrait))
2271           continue;
2272         if ((statistic_traits & CopyPixelTrait) != 0)
2273           {
2274             SetPixelChannel(statistic_image,channel,p[center+i],q);
2275             continue;
2276           }
2277         pixels=p;
2278         ResetPixelList(pixel_list[id]);
2279         for (v=0; v < (ssize_t) MagickMax(height,1); v++)
2280         {
2281           for (u=0; u < (ssize_t) MagickMax(width,1); u++)
2282           {
2283             InsertPixelList(image,pixels[i],pixel_list[id]);
2284             pixels+=GetPixelChannels(image);
2285           }
2286           pixels+=image->columns*GetPixelChannels(image);
2287         }
2288         switch (type)
2289         {
2290           case GradientStatistic:
2291           {
2292             MagickRealType
2293               maximum,
2294               minimum;
2295
2296             GetMinimumPixelList(pixel_list[id],&pixel);
2297             minimum=(MagickRealType) pixel;
2298             GetMaximumPixelList(pixel_list[id],&pixel);
2299             maximum=(MagickRealType) pixel;
2300             pixel=ClampToQuantum(MagickAbsoluteValue(maximum-minimum));
2301             break;
2302           }
2303           case MaximumStatistic:
2304           {
2305             GetMaximumPixelList(pixel_list[id],&pixel);
2306             break;
2307           }
2308           case MeanStatistic:
2309           {
2310             GetMeanPixelList(pixel_list[id],&pixel);
2311             break;
2312           }
2313           case MedianStatistic:
2314           default:
2315           {
2316             GetMedianPixelList(pixel_list[id],&pixel);
2317             break;
2318           }
2319           case MinimumStatistic:
2320           {
2321             GetMinimumPixelList(pixel_list[id],&pixel);
2322             break;
2323           }
2324           case ModeStatistic:
2325           {
2326             GetModePixelList(pixel_list[id],&pixel);
2327             break;
2328           }
2329           case NonpeakStatistic:
2330           {
2331             GetNonpeakPixelList(pixel_list[id],&pixel);
2332             break;
2333           }
2334           case StandardDeviationStatistic:
2335           {
2336             GetStandardDeviationPixelList(pixel_list[id],&pixel);
2337             break;
2338           }
2339         }
2340         SetPixelChannel(statistic_image,channel,pixel,q);
2341       }
2342       p+=GetPixelChannels(image);
2343       q+=GetPixelChannels(statistic_image);
2344     }
2345     if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
2346       status=MagickFalse;
2347     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2348       {
2349         MagickBooleanType
2350           proceed;
2351
2352 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2353   #pragma omp critical (MagickCore_StatisticImage)
2354 #endif
2355         proceed=SetImageProgress(image,StatisticImageTag,progress++,
2356           image->rows);
2357         if (proceed == MagickFalse)
2358           status=MagickFalse;
2359       }
2360   }
2361   statistic_view=DestroyCacheView(statistic_view);
2362   image_view=DestroyCacheView(image_view);
2363   pixel_list=DestroyPixelListThreadSet(pixel_list);
2364   return(statistic_image);
2365 }