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