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