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