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