statistic.c revision 50908ac2139298d3f39d279d75c5868caea64316 
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 RootMeanSquareEvaluateOperator:
373    {
374      result=(double) (pixel*pixel+value);
375      break;
376    }
377    case SetEvaluateOperator:
378    {
379      result=value;
380      break;
381    }
382    case SineEvaluateOperator:
383    {
384      result=(double) (QuantumRange*(0.5*sin((double) (2.0*MagickPI*
385        QuantumScale*pixel*value))+0.5));
386      break;
387    }
388    case SubtractEvaluateOperator:
389    {
390      result=(double) (pixel-value);
391      break;
392    }
393    case SumEvaluateOperator:
394    {
395      result=(double) (pixel+value);
396      break;
397    }
398    case ThresholdEvaluateOperator:
399    {
400      result=(double) (((double) pixel <= value) ? 0 : QuantumRange);
401      break;
402    }
403    case ThresholdBlackEvaluateOperator:
404    {
405      result=(double) (((double) pixel <= value) ? 0 : pixel);
406      break;
407    }
408    case ThresholdWhiteEvaluateOperator:
409    {
410      result=(double) (((double) pixel > value) ? QuantumRange : pixel);
411      break;
412    }
413    case UniformNoiseEvaluateOperator:
414    {
415      result=(double) GenerateDifferentialNoise(random_info,pixel,UniformNoise,
416        value);
417      break;
418    }
419    case XorEvaluateOperator:
420    {
421      result=(double) ((size_t) pixel ^ (size_t) (value+0.5));
422      break;
423    }
424  }
425  return(result);
426}
427
428MagickExport Image *EvaluateImages(const Image *images,
429  const MagickEvaluateOperator op,ExceptionInfo *exception)
430{
431#define EvaluateImageTag  "Evaluate/Image"
432
433  CacheView
434    *evaluate_view;
435
436  Image
437    *image;
438
439  MagickBooleanType
440    status;
441
442  MagickOffsetType
443    progress;
444
445  PixelChannels
446    **restrict evaluate_pixels;
447
448  RandomInfo
449    **restrict random_info;
450
451  size_t
452    number_images;
453
454  ssize_t
455    y;
456
457#if defined(MAGICKCORE_OPENMP_SUPPORT)
458  unsigned long
459    key;
460#endif
461
462  assert(images != (Image *) NULL);
463  assert(images->signature == MagickSignature);
464  if (images->debug != MagickFalse)
465    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
466  assert(exception != (ExceptionInfo *) NULL);
467  assert(exception->signature == MagickSignature);
468  image=CloneImage(images,images->columns,images->rows,MagickTrue,
469    exception);
470  if (image == (Image *) NULL)
471    return((Image *) NULL);
472  if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
473    {
474      image=DestroyImage(image);
475      return((Image *) NULL);
476    }
477  number_images=GetImageListLength(images);
478  evaluate_pixels=AcquirePixelThreadSet(images,number_images);
479  if (evaluate_pixels == (PixelChannels **) NULL)
480    {
481      image=DestroyImage(image);
482      (void) ThrowMagickException(exception,GetMagickModule(),
483        ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
484      return((Image *) NULL);
485    }
486  /*
487    Evaluate image pixels.
488  */
489  status=MagickTrue;
490  progress=0;
491  random_info=AcquireRandomInfoThreadSet();
492#if defined(MAGICKCORE_OPENMP_SUPPORT)
493  key=GetRandomSecretKey(random_info[0]);
494#endif
495  evaluate_view=AcquireAuthenticCacheView(image,exception);
496  if (op == MedianEvaluateOperator)
497    {
498#if defined(MAGICKCORE_OPENMP_SUPPORT)
499      #pragma omp parallel for schedule(static,4) shared(progress,status) \
500        magick_threads(image,images,image->rows,key == ~0UL)
501#endif
502      for (y=0; y < (ssize_t) image->rows; y++)
503      {
504        CacheView
505          *image_view;
506
507        const Image
508          *next;
509
510        const int
511          id = GetOpenMPThreadId();
512
513        register PixelChannels
514          *evaluate_pixel;
515
516        register Quantum
517          *restrict q;
518
519        register ssize_t
520          x;
521
522        if (status == MagickFalse)
523          continue;
524        q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,image->columns,1,
525          exception);
526        if (q == (Quantum *) NULL)
527          {
528            status=MagickFalse;
529            continue;
530          }
531        evaluate_pixel=evaluate_pixels[id];
532        for (x=0; x < (ssize_t) image->columns; x++)
533        {
534          register ssize_t
535            j,
536            k;
537
538          for (j=0; j < (ssize_t) number_images; j++)
539            for (k=0; k < MaxPixelChannels; k++)
540              evaluate_pixel[j].channel[k]=0.0;
541          next=images;
542          for (j=0; j < (ssize_t) number_images; j++)
543          {
544            register const Quantum
545              *p;
546
547            register ssize_t
548              i;
549
550            image_view=AcquireVirtualCacheView(next,exception);
551            p=GetCacheViewVirtualPixels(image_view,x,y,1,1,exception);
552            if (p == (const Quantum *) NULL)
553              {
554                image_view=DestroyCacheView(image_view);
555                break;
556              }
557            for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
558            {
559              PixelChannel channel=GetPixelChannelChannel(image,i);
560              PixelTrait evaluate_traits=GetPixelChannelTraits(image,channel);
561              PixelTrait traits=GetPixelChannelTraits(next,channel);
562              if ((traits == UndefinedPixelTrait) ||
563                  (evaluate_traits == UndefinedPixelTrait))
564                continue;
565              if ((evaluate_traits & UpdatePixelTrait) == 0)
566                continue;
567              evaluate_pixel[j].channel[i]=ApplyEvaluateOperator(
568                random_info[id],GetPixelChannel(image,channel,p),op,
569                evaluate_pixel[j].channel[i]);
570            }
571            image_view=DestroyCacheView(image_view);
572            next=GetNextImageInList(next);
573          }
574          qsort((void *) evaluate_pixel,number_images,sizeof(*evaluate_pixel),
575            IntensityCompare);
576          for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
577            q[k]=ClampToQuantum(evaluate_pixel[j/2].channel[k]);
578          q+=GetPixelChannels(image);
579        }
580        if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
581          status=MagickFalse;
582        if (images->progress_monitor != (MagickProgressMonitor) NULL)
583          {
584            MagickBooleanType
585              proceed;
586
587#if   defined(MAGICKCORE_OPENMP_SUPPORT)
588            #pragma omp critical (MagickCore_EvaluateImages)
589#endif
590            proceed=SetImageProgress(images,EvaluateImageTag,progress++,
591              image->rows);
592            if (proceed == MagickFalse)
593              status=MagickFalse;
594          }
595      }
596    }
597  else
598    {
599#if defined(MAGICKCORE_OPENMP_SUPPORT)
600      #pragma omp parallel for schedule(static,4) shared(progress,status) \
601        magick_threads(image,images,image->rows,key == ~0UL)
602#endif
603      for (y=0; y < (ssize_t) image->rows; y++)
604      {
605        CacheView
606          *image_view;
607
608        const Image
609          *next;
610
611        const int
612          id = GetOpenMPThreadId();
613
614        register ssize_t
615          i,
616          x;
617
618        register PixelChannels
619          *evaluate_pixel;
620
621        register Quantum
622          *restrict q;
623
624        ssize_t
625          j;
626
627        if (status == MagickFalse)
628          continue;
629        q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,image->columns,1,
630          exception);
631        if (q == (Quantum *) NULL)
632          {
633            status=MagickFalse;
634            continue;
635          }
636        evaluate_pixel=evaluate_pixels[id];
637        for (j=0; j < (ssize_t) image->columns; j++)
638          for (i=0; i < MaxPixelChannels; i++)
639            evaluate_pixel[j].channel[i]=0.0;
640        next=images;
641        for (j=0; j < (ssize_t) number_images; j++)
642        {
643          register const Quantum
644            *p;
645
646          image_view=AcquireVirtualCacheView(next,exception);
647          p=GetCacheViewVirtualPixels(image_view,0,y,next->columns,1,exception);
648          if (p == (const Quantum *) NULL)
649            {
650              image_view=DestroyCacheView(image_view);
651              break;
652            }
653          for (x=0; x < (ssize_t) next->columns; x++)
654          {
655            register ssize_t
656              i;
657
658            if (GetPixelReadMask(next,p) == 0)
659              {
660                p+=GetPixelChannels(next);
661                continue;
662              }
663            for (i=0; i < (ssize_t) GetPixelChannels(next); i++)
664            {
665              PixelChannel channel=GetPixelChannelChannel(image,i);
666              PixelTrait  traits=GetPixelChannelTraits(next,channel);
667              PixelTrait  evaluate_traits=GetPixelChannelTraits(image,channel);
668              if ((traits == UndefinedPixelTrait) ||
669                  (evaluate_traits == UndefinedPixelTrait))
670                continue;
671              if ((traits & UpdatePixelTrait) == 0)
672                continue;
673              evaluate_pixel[x].channel[i]=ApplyEvaluateOperator(
674                random_info[id],GetPixelChannel(image,channel,p),j == 0 ?
675                AddEvaluateOperator : op,evaluate_pixel[x].channel[i]);
676            }
677            p+=GetPixelChannels(next);
678          }
679          image_view=DestroyCacheView(image_view);
680          next=GetNextImageInList(next);
681        }
682        for (x=0; x < (ssize_t) image->columns; x++)
683        {
684          register ssize_t
685             i;
686
687          switch (op)
688          {
689            case MeanEvaluateOperator:
690            {
691              for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
692                evaluate_pixel[x].channel[i]/=(double) number_images;
693              break;
694            }
695            case MultiplyEvaluateOperator:
696            {
697              for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
698              {
699                register ssize_t
700                  j;
701
702                for (j=0; j < (ssize_t) (number_images-1); j++)
703                  evaluate_pixel[x].channel[i]*=QuantumScale;
704              }
705              break;
706            }
707            case RootMeanSquareEvaluateOperator:
708            {
709              for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
710                evaluate_pixel[x].channel[i]=sqrt(evaluate_pixel[x].channel[i]/
711                  number_images);
712              break;
713            }
714            default:
715              break;
716          }
717        }
718        for (x=0; x < (ssize_t) image->columns; x++)
719        {
720          register ssize_t
721            i;
722
723          if (GetPixelReadMask(image,q) == 0)
724            {
725              q+=GetPixelChannels(image);
726              continue;
727            }
728          for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
729          {
730            PixelChannel channel=GetPixelChannelChannel(image,i);
731            PixelTrait traits=GetPixelChannelTraits(image,channel);
732            if (traits == UndefinedPixelTrait)
733              continue;
734            if ((traits & UpdatePixelTrait) == 0)
735              continue;
736            q[i]=ClampToQuantum(evaluate_pixel[x].channel[i]);
737          }
738          q+=GetPixelChannels(image);
739        }
740        if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
741          status=MagickFalse;
742        if (images->progress_monitor != (MagickProgressMonitor) NULL)
743          {
744            MagickBooleanType
745              proceed;
746
747#if   defined(MAGICKCORE_OPENMP_SUPPORT)
748            #pragma omp critical (MagickCore_EvaluateImages)
749#endif
750            proceed=SetImageProgress(images,EvaluateImageTag,progress++,
751              image->rows);
752            if (proceed == MagickFalse)
753              status=MagickFalse;
754          }
755      }
756    }
757  evaluate_view=DestroyCacheView(evaluate_view);
758  evaluate_pixels=DestroyPixelThreadSet(evaluate_pixels);
759  random_info=DestroyRandomInfoThreadSet(random_info);
760  if (status == MagickFalse)
761    image=DestroyImage(image);
762  return(image);
763}
764
765MagickExport MagickBooleanType EvaluateImage(Image *image,
766  const MagickEvaluateOperator op,const double value,ExceptionInfo *exception)
767{
768  CacheView
769    *image_view;
770
771  MagickBooleanType
772    status;
773
774  MagickOffsetType
775    progress;
776
777  RandomInfo
778    **restrict random_info;
779
780  ssize_t
781    y;
782
783#if defined(MAGICKCORE_OPENMP_SUPPORT)
784  unsigned long
785    key;
786#endif
787
788  assert(image != (Image *) NULL);
789  assert(image->signature == MagickSignature);
790  if (image->debug != MagickFalse)
791    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
792  assert(exception != (ExceptionInfo *) NULL);
793  assert(exception->signature == MagickSignature);
794  if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
795    return(MagickFalse);
796  status=MagickTrue;
797  progress=0;
798  random_info=AcquireRandomInfoThreadSet();
799#if defined(MAGICKCORE_OPENMP_SUPPORT)
800  key=GetRandomSecretKey(random_info[0]);
801#endif
802  image_view=AcquireAuthenticCacheView(image,exception);
803#if defined(MAGICKCORE_OPENMP_SUPPORT)
804  #pragma omp parallel for schedule(static,4) shared(progress,status) \
805    magick_threads(image,image,image->rows,key == ~0UL)
806#endif
807  for (y=0; y < (ssize_t) image->rows; y++)
808  {
809    const int
810      id = GetOpenMPThreadId();
811
812    register Quantum
813      *restrict q;
814
815    register ssize_t
816      x;
817
818    if (status == MagickFalse)
819      continue;
820    q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
821    if (q == (Quantum *) NULL)
822      {
823        status=MagickFalse;
824        continue;
825      }
826    for (x=0; x < (ssize_t) image->columns; x++)
827    {
828      double
829        result;
830
831      register ssize_t
832        i;
833
834      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
835      {
836        PixelChannel channel=GetPixelChannelChannel(image,i);
837        PixelTrait traits=GetPixelChannelTraits(image,channel);
838        if (traits == UndefinedPixelTrait)
839          continue;
840        if (((traits & CopyPixelTrait) != 0) ||
841            (GetPixelReadMask(image,q) == 0))
842          continue;
843        result=ApplyEvaluateOperator(random_info[id],q[i],op,value);
844        if (op == MeanEvaluateOperator)
845          result/=2.0;
846        q[i]=ClampToQuantum(result);
847      }
848      q+=GetPixelChannels(image);
849    }
850    if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
851      status=MagickFalse;
852    if (image->progress_monitor != (MagickProgressMonitor) NULL)
853      {
854        MagickBooleanType
855          proceed;
856
857#if defined(MAGICKCORE_OPENMP_SUPPORT)
858        #pragma omp critical (MagickCore_EvaluateImage)
859#endif
860        proceed=SetImageProgress(image,EvaluateImageTag,progress++,image->rows);
861        if (proceed == MagickFalse)
862          status=MagickFalse;
863      }
864  }
865  image_view=DestroyCacheView(image_view);
866  random_info=DestroyRandomInfoThreadSet(random_info);
867  return(status);
868}
869
870/*
871%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
872%                                                                             %
873%                                                                             %
874%                                                                             %
875%     F u n c t i o n I m a g e                                               %
876%                                                                             %
877%                                                                             %
878%                                                                             %
879%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
880%
881%  FunctionImage() applies a value to the image with an arithmetic, relational,
882%  or logical operator to an image. Use these operations to lighten or darken
883%  an image, to increase or decrease contrast in an image, or to produce the
884%  "negative" of an image.
885%
886%  The format of the FunctionImage method is:
887%
888%      MagickBooleanType FunctionImage(Image *image,
889%        const MagickFunction function,const ssize_t number_parameters,
890%        const double *parameters,ExceptionInfo *exception)
891%
892%  A description of each parameter follows:
893%
894%    o image: the image.
895%
896%    o function: A channel function.
897%
898%    o parameters: one or more parameters.
899%
900%    o exception: return any errors or warnings in this structure.
901%
902*/
903
904static Quantum ApplyFunction(Quantum pixel,const MagickFunction function,
905  const size_t number_parameters,const double *parameters,
906  ExceptionInfo *exception)
907{
908  double
909    result;
910
911  register ssize_t
912    i;
913
914  (void) exception;
915  result=0.0;
916  switch (function)
917  {
918    case PolynomialFunction:
919    {
920      /*
921        Polynomial: polynomial constants, highest to lowest order (e.g. c0*x^3+
922        c1*x^2+c2*x+c3).
923      */
924      result=0.0;
925      for (i=0; i < (ssize_t) number_parameters; i++)
926        result=result*QuantumScale*pixel+parameters[i];
927      result*=QuantumRange;
928      break;
929    }
930    case SinusoidFunction:
931    {
932      double
933        amplitude,
934        bias,
935        frequency,
936        phase;
937
938      /*
939        Sinusoid: frequency, phase, amplitude, bias.
940      */
941      frequency=(number_parameters >= 1) ? parameters[0] : 1.0;
942      phase=(number_parameters >= 2) ? parameters[1] : 0.0;
943      amplitude=(number_parameters >= 3) ? parameters[2] : 0.5;
944      bias=(number_parameters >= 4) ? parameters[3] : 0.5;
945      result=(double) (QuantumRange*(amplitude*sin((double) (2.0*
946        MagickPI*(frequency*QuantumScale*pixel+phase/360.0)))+bias));
947      break;
948    }
949    case ArcsinFunction:
950    {
951      double
952        bias,
953        center,
954        range,
955        width;
956
957      /*
958        Arcsin (peged at range limits for invalid results): width, center,
959        range, and bias.
960      */
961      width=(number_parameters >= 1) ? parameters[0] : 1.0;
962      center=(number_parameters >= 2) ? parameters[1] : 0.5;
963      range=(number_parameters >= 3) ? parameters[2] : 1.0;
964      bias=(number_parameters >= 4) ? parameters[3] : 0.5;
965      result=2.0/width*(QuantumScale*pixel-center);
966      if ( result <= -1.0 )
967        result=bias-range/2.0;
968      else
969        if (result >= 1.0)
970          result=bias+range/2.0;
971        else
972          result=(double) (range/MagickPI*asin((double) result)+bias);
973      result*=QuantumRange;
974      break;
975    }
976    case ArctanFunction:
977    {
978      double
979        center,
980        bias,
981        range,
982        slope;
983
984      /*
985        Arctan: slope, center, range, and bias.
986      */
987      slope=(number_parameters >= 1) ? parameters[0] : 1.0;
988      center=(number_parameters >= 2) ? parameters[1] : 0.5;
989      range=(number_parameters >= 3) ? parameters[2] : 1.0;
990      bias=(number_parameters >= 4) ? parameters[3] : 0.5;
991      result=(double) (MagickPI*slope*(QuantumScale*pixel-center));
992      result=(double) (QuantumRange*(range/MagickPI*atan((double)
993        result)+bias));
994      break;
995    }
996    case UndefinedFunction:
997      break;
998  }
999  return(ClampToQuantum(result));
1000}
1001
1002MagickExport MagickBooleanType FunctionImage(Image *image,
1003  const MagickFunction function,const size_t number_parameters,
1004  const double *parameters,ExceptionInfo *exception)
1005{
1006#define FunctionImageTag  "Function/Image "
1007
1008  CacheView
1009    *image_view;
1010
1011  MagickBooleanType
1012    status;
1013
1014  MagickOffsetType
1015    progress;
1016
1017  ssize_t
1018    y;
1019
1020  assert(image != (Image *) NULL);
1021  assert(image->signature == MagickSignature);
1022  if (image->debug != MagickFalse)
1023    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1024  assert(exception != (ExceptionInfo *) NULL);
1025  assert(exception->signature == MagickSignature);
1026  if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
1027    return(MagickFalse);
1028  status=MagickTrue;
1029  progress=0;
1030  image_view=AcquireAuthenticCacheView(image,exception);
1031#if defined(MAGICKCORE_OPENMP_SUPPORT)
1032  #pragma omp parallel for schedule(static,4) shared(progress,status) \
1033    magick_threads(image,image,image->rows,1)
1034#endif
1035  for (y=0; y < (ssize_t) image->rows; y++)
1036  {
1037    register Quantum
1038      *restrict q;
1039
1040    register ssize_t
1041      x;
1042
1043    if (status == MagickFalse)
1044      continue;
1045    q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
1046    if (q == (Quantum *) NULL)
1047      {
1048        status=MagickFalse;
1049        continue;
1050      }
1051    for (x=0; x < (ssize_t) image->columns; x++)
1052    {
1053      register ssize_t
1054        i;
1055
1056      if (GetPixelReadMask(image,q) == 0)
1057        {
1058          q+=GetPixelChannels(image);
1059          continue;
1060        }
1061      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1062      {
1063        PixelChannel channel=GetPixelChannelChannel(image,i);
1064        PixelTrait traits=GetPixelChannelTraits(image,channel);
1065        if (traits == UndefinedPixelTrait)
1066          continue;
1067        if ((traits & UpdatePixelTrait) == 0)
1068          continue;
1069        q[i]=ApplyFunction(q[i],function,number_parameters,parameters,
1070          exception);
1071      }
1072      q+=GetPixelChannels(image);
1073    }
1074    if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1075      status=MagickFalse;
1076    if (image->progress_monitor != (MagickProgressMonitor) NULL)
1077      {
1078        MagickBooleanType
1079          proceed;
1080
1081#if defined(MAGICKCORE_OPENMP_SUPPORT)
1082        #pragma omp critical (MagickCore_FunctionImage)
1083#endif
1084        proceed=SetImageProgress(image,FunctionImageTag,progress++,image->rows);
1085        if (proceed == MagickFalse)
1086          status=MagickFalse;
1087      }
1088  }
1089  image_view=DestroyCacheView(image_view);
1090  return(status);
1091}
1092
1093/*
1094%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1095%                                                                             %
1096%                                                                             %
1097%                                                                             %
1098%   G e t I m a g e E n t r o p y                                             %
1099%                                                                             %
1100%                                                                             %
1101%                                                                             %
1102%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1103%
1104%  GetImageEntropy() returns the entropy of one or more image channels.
1105%
1106%  The format of the GetImageEntropy method is:
1107%
1108%      MagickBooleanType GetImageEntropy(const Image *image,double *entropy,
1109%        ExceptionInfo *exception)
1110%
1111%  A description of each parameter follows:
1112%
1113%    o image: the image.
1114%
1115%    o entropy: the average entropy of the selected channels.
1116%
1117%    o exception: return any errors or warnings in this structure.
1118%
1119*/
1120MagickExport MagickBooleanType GetImageEntropy(const Image *image,
1121  double *entropy,ExceptionInfo *exception)
1122{
1123  double
1124    area;
1125
1126  ChannelStatistics
1127    *channel_statistics;
1128
1129  register ssize_t
1130    i;
1131
1132  assert(image != (Image *) NULL);
1133  assert(image->signature == MagickSignature);
1134  if (image->debug != MagickFalse)
1135    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1136  channel_statistics=GetImageStatistics(image,exception);
1137  if (channel_statistics == (ChannelStatistics *) NULL)
1138    return(MagickFalse);
1139  area=0.0;
1140  channel_statistics[CompositePixelChannel].entropy=0.0;
1141  channel_statistics[CompositePixelChannel].standard_deviation=0.0;
1142  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1143  {
1144    PixelChannel channel=GetPixelChannelChannel(image,i);
1145    PixelTrait traits=GetPixelChannelTraits(image,channel);
1146    if (traits == UndefinedPixelTrait)
1147      continue;
1148    if ((traits & UpdatePixelTrait) == 0)
1149      continue;
1150    channel_statistics[CompositePixelChannel].entropy+=
1151      channel_statistics[i].entropy;
1152    area++;
1153  }
1154  channel_statistics[CompositePixelChannel].entropy/=area;
1155  channel_statistics[CompositePixelChannel].standard_deviation=
1156    sqrt(channel_statistics[CompositePixelChannel].standard_deviation/area);
1157  *entropy=channel_statistics[CompositePixelChannel].entropy;
1158  channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1159    channel_statistics);
1160  return(MagickTrue);
1161}
1162
1163/*
1164%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1165%                                                                             %
1166%                                                                             %
1167%                                                                             %
1168%   G e t I m a g e E x t r e m a                                             %
1169%                                                                             %
1170%                                                                             %
1171%                                                                             %
1172%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1173%
1174%  GetImageExtrema() returns the extrema of one or more image channels.
1175%
1176%  The format of the GetImageExtrema method is:
1177%
1178%      MagickBooleanType GetImageExtrema(const Image *image,size_t *minima,
1179%        size_t *maxima,ExceptionInfo *exception)
1180%
1181%  A description of each parameter follows:
1182%
1183%    o image: the image.
1184%
1185%    o minima: the minimum value in the channel.
1186%
1187%    o maxima: the maximum value in the channel.
1188%
1189%    o exception: return any errors or warnings in this structure.
1190%
1191*/
1192MagickExport MagickBooleanType GetImageExtrema(const Image *image,
1193  size_t *minima,size_t *maxima,ExceptionInfo *exception)
1194{
1195  double
1196    max,
1197    min;
1198
1199  MagickBooleanType
1200    status;
1201
1202  assert(image != (Image *) NULL);
1203  assert(image->signature == MagickSignature);
1204  if (image->debug != MagickFalse)
1205    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1206  status=GetImageRange(image,&min,&max,exception);
1207  *minima=(size_t) ceil(min-0.5);
1208  *maxima=(size_t) floor(max+0.5);
1209  return(status);
1210}
1211
1212/*
1213%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1214%                                                                             %
1215%                                                                             %
1216%                                                                             %
1217%   G e t I m a g e K u r t o s i s                                           %
1218%                                                                             %
1219%                                                                             %
1220%                                                                             %
1221%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1222%
1223%  GetImageKurtosis() returns the kurtosis and skewness of one or more image
1224%  channels.
1225%
1226%  The format of the GetImageKurtosis method is:
1227%
1228%      MagickBooleanType GetImageKurtosis(const Image *image,double *kurtosis,
1229%        double *skewness,ExceptionInfo *exception)
1230%
1231%  A description of each parameter follows:
1232%
1233%    o image: the image.
1234%
1235%    o kurtosis: the kurtosis of the channel.
1236%
1237%    o skewness: the skewness of the channel.
1238%
1239%    o exception: return any errors or warnings in this structure.
1240%
1241*/
1242MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
1243  double *kurtosis,double *skewness,ExceptionInfo *exception)
1244{
1245  CacheView
1246    *image_view;
1247
1248  double
1249    area,
1250    mean,
1251    standard_deviation,
1252    sum_squares,
1253    sum_cubes,
1254    sum_fourth_power;
1255
1256  MagickBooleanType
1257    status;
1258
1259  ssize_t
1260    y;
1261
1262  assert(image != (Image *) NULL);
1263  assert(image->signature == MagickSignature);
1264  if (image->debug != MagickFalse)
1265    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1266  status=MagickTrue;
1267  *kurtosis=0.0;
1268  *skewness=0.0;
1269  area=0.0;
1270  mean=0.0;
1271  standard_deviation=0.0;
1272  sum_squares=0.0;
1273  sum_cubes=0.0;
1274  sum_fourth_power=0.0;
1275  image_view=AcquireVirtualCacheView(image,exception);
1276#if defined(MAGICKCORE_OPENMP_SUPPORT)
1277  #pragma omp parallel for schedule(static,4) shared(status) \
1278    magick_threads(image,image,image->rows,1)
1279#endif
1280  for (y=0; y < (ssize_t) image->rows; y++)
1281  {
1282    register const Quantum
1283      *restrict p;
1284
1285    register ssize_t
1286      x;
1287
1288    if (status == MagickFalse)
1289      continue;
1290    p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1291    if (p == (const Quantum *) NULL)
1292      {
1293        status=MagickFalse;
1294        continue;
1295      }
1296    for (x=0; x < (ssize_t) image->columns; x++)
1297    {
1298      register ssize_t
1299        i;
1300
1301      if (GetPixelReadMask(image,p) == 0)
1302        {
1303          p+=GetPixelChannels(image);
1304          continue;
1305        }
1306      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1307      {
1308        PixelChannel channel=GetPixelChannelChannel(image,i);
1309        PixelTrait traits=GetPixelChannelTraits(image,channel);
1310        if (traits == UndefinedPixelTrait)
1311          continue;
1312        if ((traits & UpdatePixelTrait) == 0)
1313          continue;
1314#if defined(MAGICKCORE_OPENMP_SUPPORT)
1315        #pragma omp critical (MagickCore_GetImageKurtosis)
1316#endif
1317        {
1318          mean+=p[i];
1319          sum_squares+=(double) p[i]*p[i];
1320          sum_cubes+=(double) p[i]*p[i]*p[i];
1321          sum_fourth_power+=(double) p[i]*p[i]*p[i]*p[i];
1322          area++;
1323        }
1324      }
1325      p+=GetPixelChannels(image);
1326    }
1327  }
1328  image_view=DestroyCacheView(image_view);
1329  if (area != 0.0)
1330    {
1331      mean/=area;
1332      sum_squares/=area;
1333      sum_cubes/=area;
1334      sum_fourth_power/=area;
1335    }
1336  standard_deviation=sqrt(sum_squares-(mean*mean));
1337  if (standard_deviation != 0.0)
1338    {
1339      *kurtosis=sum_fourth_power-4.0*mean*sum_cubes+6.0*mean*mean*sum_squares-
1340        3.0*mean*mean*mean*mean;
1341      *kurtosis/=standard_deviation*standard_deviation*standard_deviation*
1342        standard_deviation;
1343      *kurtosis-=3.0;
1344      *skewness=sum_cubes-3.0*mean*sum_squares+2.0*mean*mean*mean;
1345      *skewness/=standard_deviation*standard_deviation*standard_deviation;
1346    }
1347  return(status);
1348}
1349
1350/*
1351%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1352%                                                                             %
1353%                                                                             %
1354%                                                                             %
1355%   G e t I m a g e M e a n                                                   %
1356%                                                                             %
1357%                                                                             %
1358%                                                                             %
1359%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1360%
1361%  GetImageMean() returns the mean and standard deviation of one or more image
1362%  channels.
1363%
1364%  The format of the GetImageMean method is:
1365%
1366%      MagickBooleanType GetImageMean(const Image *image,double *mean,
1367%        double *standard_deviation,ExceptionInfo *exception)
1368%
1369%  A description of each parameter follows:
1370%
1371%    o image: the image.
1372%
1373%    o mean: the average value in the channel.
1374%
1375%    o standard_deviation: the standard deviation of the channel.
1376%
1377%    o exception: return any errors or warnings in this structure.
1378%
1379*/
1380MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean,
1381  double *standard_deviation,ExceptionInfo *exception)
1382{
1383  double
1384    area;
1385
1386  ChannelStatistics
1387    *channel_statistics;
1388
1389  register ssize_t
1390    i;
1391
1392  assert(image != (Image *) NULL);
1393  assert(image->signature == MagickSignature);
1394  if (image->debug != MagickFalse)
1395    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1396  channel_statistics=GetImageStatistics(image,exception);
1397  if (channel_statistics == (ChannelStatistics *) NULL)
1398    return(MagickFalse);
1399  area=0.0;
1400  channel_statistics[CompositePixelChannel].mean=0.0;
1401  channel_statistics[CompositePixelChannel].standard_deviation=0.0;
1402  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1403  {
1404    PixelChannel channel=GetPixelChannelChannel(image,i);
1405    PixelTrait traits=GetPixelChannelTraits(image,channel);
1406    if (traits == UndefinedPixelTrait)
1407      continue;
1408    if ((traits & UpdatePixelTrait) == 0)
1409      continue;
1410    channel_statistics[CompositePixelChannel].mean+=channel_statistics[i].mean;
1411    channel_statistics[CompositePixelChannel].standard_deviation+=
1412      channel_statistics[i].variance-channel_statistics[i].mean*
1413      channel_statistics[i].mean;
1414    area++;
1415  }
1416  channel_statistics[CompositePixelChannel].mean/=area;
1417  channel_statistics[CompositePixelChannel].standard_deviation=
1418    sqrt(channel_statistics[CompositePixelChannel].standard_deviation/area);
1419  *mean=channel_statistics[CompositePixelChannel].mean;
1420  *standard_deviation=
1421    channel_statistics[CompositePixelChannel].standard_deviation;
1422  channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1423    channel_statistics);
1424  return(MagickTrue);
1425}
1426
1427/*
1428%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1429%                                                                             %
1430%                                                                             %
1431%                                                                             %
1432%   G e t I m a g e M o m e n t s                                             %
1433%                                                                             %
1434%                                                                             %
1435%                                                                             %
1436%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1437%
1438%  GetImageMoments() returns the normalized moments of one or more image
1439%  channels.
1440%
1441%  The format of the GetImageMoments method is:
1442%
1443%      ChannelMoments *GetImageMoments(const Image *image,
1444%        ExceptionInfo *exception)
1445%
1446%  A description of each parameter follows:
1447%
1448%    o image: the image.
1449%
1450%    o exception: return any errors or warnings in this structure.
1451%
1452*/
1453
1454static size_t GetImageChannels(const Image *image)
1455{
1456  register ssize_t
1457    i;
1458
1459  size_t
1460    channels;
1461
1462  channels=0;
1463  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1464  {
1465    PixelChannel channel=GetPixelChannelChannel(image,i);
1466    PixelTrait traits=GetPixelChannelTraits(image,channel);
1467    if ((traits & UpdatePixelTrait) != 0)
1468      channels++;
1469  }
1470  return((size_t) (channels == 0 ? 1 : channels));
1471}
1472
1473MagickExport ChannelMoments *GetImageMoments(const Image *image,
1474  ExceptionInfo *exception)
1475{
1476#define MaxNumberImageMoments  8
1477
1478  CacheView
1479    *image_view;
1480
1481  ChannelMoments
1482    *channel_moments;
1483
1484  double
1485    M00[MaxPixelChannels+1],
1486    M01[MaxPixelChannels+1],
1487    M02[MaxPixelChannels+1],
1488    M03[MaxPixelChannels+1],
1489    M10[MaxPixelChannels+1],
1490    M11[MaxPixelChannels+1],
1491    M12[MaxPixelChannels+1],
1492    M20[MaxPixelChannels+1],
1493    M21[MaxPixelChannels+1],
1494    M22[MaxPixelChannels+1],
1495    M30[MaxPixelChannels+1];
1496
1497  PointInfo
1498    centroid[MaxPixelChannels+1];
1499
1500  ssize_t
1501    channel,
1502    y;
1503
1504  assert(image != (Image *) NULL);
1505  assert(image->signature == MagickSignature);
1506  if (image->debug != MagickFalse)
1507    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1508  channel_moments=(ChannelMoments *) AcquireQuantumMemory(MaxPixelChannels+1,
1509    sizeof(*channel_moments));
1510  if (channel_moments == (ChannelMoments *) NULL)
1511    return(channel_moments);
1512  (void) ResetMagickMemory(channel_moments,0,(MaxPixelChannels+1)*
1513    sizeof(*channel_moments));
1514  (void) ResetMagickMemory(centroid,0,sizeof(centroid));
1515  (void) ResetMagickMemory(M00,0,sizeof(M00));
1516  (void) ResetMagickMemory(M01,0,sizeof(M01));
1517  (void) ResetMagickMemory(M02,0,sizeof(M02));
1518  (void) ResetMagickMemory(M03,0,sizeof(M03));
1519  (void) ResetMagickMemory(M10,0,sizeof(M10));
1520  (void) ResetMagickMemory(M11,0,sizeof(M11));
1521  (void) ResetMagickMemory(M12,0,sizeof(M12));
1522  (void) ResetMagickMemory(M20,0,sizeof(M20));
1523  (void) ResetMagickMemory(M21,0,sizeof(M21));
1524  (void) ResetMagickMemory(M22,0,sizeof(M22));
1525  (void) ResetMagickMemory(M30,0,sizeof(M30));
1526  image_view=AcquireVirtualCacheView(image,exception);
1527  for (y=0; y < (ssize_t) image->rows; y++)
1528  {
1529    register const Quantum
1530      *restrict p;
1531
1532    register ssize_t
1533      x;
1534
1535    /*
1536      Compute center of mass (centroid).
1537    */
1538    p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1539    if (p == (const Quantum *) NULL)
1540      break;
1541    for (x=0; x < (ssize_t) image->columns; x++)
1542    {
1543      register ssize_t
1544        i;
1545
1546      if (GetPixelReadMask(image,p) == 0)
1547        {
1548          p+=GetPixelChannels(image);
1549          continue;
1550        }
1551      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1552      {
1553        PixelChannel channel=GetPixelChannelChannel(image,i);
1554        PixelTrait traits=GetPixelChannelTraits(image,channel);
1555        if (traits == UndefinedPixelTrait)
1556          continue;
1557        if ((traits & UpdatePixelTrait) == 0)
1558          continue;
1559        M00[channel]+=QuantumScale*p[i];
1560        M00[MaxPixelChannels]+=QuantumScale*p[i];
1561        M10[channel]+=x*QuantumScale*p[i];
1562        M10[MaxPixelChannels]+=x*QuantumScale*p[i];
1563        M01[channel]+=y*QuantumScale*p[i];
1564        M01[MaxPixelChannels]+=y*QuantumScale*p[i];
1565      }
1566      p+=GetPixelChannels(image);
1567    }
1568  }
1569  for (channel=0; channel <= MaxPixelChannels; channel++)
1570  {
1571    /*
1572       Compute center of mass (centroid).
1573    */
1574    if (M00[channel] < MagickEpsilon)
1575      {
1576        M00[channel]+=MagickEpsilon;
1577        centroid[channel].x=(double) image->columns/2.0;
1578        centroid[channel].y=(double) image->rows/2.0;
1579        continue;
1580      }
1581    M00[channel]+=MagickEpsilon;
1582    centroid[channel].x=M10[channel]/M00[channel];
1583    centroid[channel].y=M01[channel]/M00[channel];
1584  }
1585  for (y=0; y < (ssize_t) image->rows; y++)
1586  {
1587    register const Quantum
1588      *restrict p;
1589
1590    register ssize_t
1591      x;
1592
1593    /*
1594      Compute the image moments.
1595    */
1596    p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1597    if (p == (const Quantum *) NULL)
1598      break;
1599    for (x=0; x < (ssize_t) image->columns; x++)
1600    {
1601      register ssize_t
1602        i;
1603
1604      if (GetPixelReadMask(image,p) == 0)
1605        {
1606          p+=GetPixelChannels(image);
1607          continue;
1608        }
1609      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1610      {
1611        PixelChannel channel=GetPixelChannelChannel(image,i);
1612        PixelTrait traits=GetPixelChannelTraits(image,channel);
1613        if (traits == UndefinedPixelTrait)
1614          continue;
1615        if ((traits & UpdatePixelTrait) == 0)
1616          continue;
1617        M11[channel]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
1618          QuantumScale*p[i];
1619        M11[MaxPixelChannels]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
1620          QuantumScale*p[i];
1621        M20[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1622          QuantumScale*p[i];
1623        M20[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1624          QuantumScale*p[i];
1625        M02[channel]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
1626          QuantumScale*p[i];
1627        M02[MaxPixelChannels]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
1628          QuantumScale*p[i];
1629        M21[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1630          (y-centroid[channel].y)*QuantumScale*p[i];
1631        M21[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1632          (y-centroid[channel].y)*QuantumScale*p[i];
1633        M12[channel]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
1634          (y-centroid[channel].y)*QuantumScale*p[i];
1635        M12[MaxPixelChannels]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
1636          (y-centroid[channel].y)*QuantumScale*p[i];
1637        M22[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1638          (y-centroid[channel].y)*(y-centroid[channel].y)*QuantumScale*p[i];
1639        M22[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1640          (y-centroid[channel].y)*(y-centroid[channel].y)*QuantumScale*p[i];
1641        M30[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1642          (x-centroid[channel].x)*QuantumScale*p[i];
1643        M30[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1644          (x-centroid[channel].x)*QuantumScale*p[i];
1645        M03[channel]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
1646          (y-centroid[channel].y)*QuantumScale*p[i];
1647        M03[MaxPixelChannels]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
1648          (y-centroid[channel].y)*QuantumScale*p[i];
1649      }
1650      p+=GetPixelChannels(image);
1651    }
1652  }
1653  M00[MaxPixelChannels]/=GetImageChannels(image);
1654  M01[MaxPixelChannels]/=GetImageChannels(image);
1655  M02[MaxPixelChannels]/=GetImageChannels(image);
1656  M03[MaxPixelChannels]/=GetImageChannels(image);
1657  M10[MaxPixelChannels]/=GetImageChannels(image);
1658  M11[MaxPixelChannels]/=GetImageChannels(image);
1659  M12[MaxPixelChannels]/=GetImageChannels(image);
1660  M20[MaxPixelChannels]/=GetImageChannels(image);
1661  M21[MaxPixelChannels]/=GetImageChannels(image);
1662  M22[MaxPixelChannels]/=GetImageChannels(image);
1663  M30[MaxPixelChannels]/=GetImageChannels(image);
1664  for (channel=0; channel <= MaxPixelChannels; channel++)
1665  {
1666    /*
1667      Compute elliptical angle, major and minor axes, eccentricity, & intensity.
1668    */
1669    channel_moments[channel].centroid=centroid[channel];
1670    channel_moments[channel].ellipse_axis.x=sqrt((2.0/M00[channel])*
1671      ((M20[channel]+M02[channel])+sqrt(4.0*M11[channel]*M11[channel]+
1672      (M20[channel]-M02[channel])*(M20[channel]-M02[channel]))));
1673    channel_moments[channel].ellipse_axis.y=sqrt((2.0/M00[channel])*
1674      ((M20[channel]+M02[channel])-sqrt(4.0*M11[channel]*M11[channel]+
1675      (M20[channel]-M02[channel])*(M20[channel]-M02[channel]))));
1676    channel_moments[channel].ellipse_angle=RadiansToDegrees(0.5*atan(2.0*
1677      M11[channel]/(M20[channel]-M02[channel]+MagickEpsilon)));
1678    channel_moments[channel].ellipse_eccentricity=sqrt(1.0-(
1679      channel_moments[channel].ellipse_axis.y/
1680      (channel_moments[channel].ellipse_axis.x+MagickEpsilon)));
1681    channel_moments[channel].ellipse_intensity=M00[channel]/
1682      (MagickPI*channel_moments[channel].ellipse_axis.x*
1683      channel_moments[channel].ellipse_axis.y+MagickEpsilon);
1684  }
1685  for (channel=0; channel <= MaxPixelChannels; channel++)
1686  {
1687    /*
1688      Normalize image moments.
1689    */
1690    M10[channel]=0.0;
1691    M01[channel]=0.0;
1692    M11[channel]/=pow(M00[channel],1.0+(1.0+1.0)/2.0);
1693    M20[channel]/=pow(M00[channel],1.0+(2.0+0.0)/2.0);
1694    M02[channel]/=pow(M00[channel],1.0+(0.0+2.0)/2.0);
1695    M21[channel]/=pow(M00[channel],1.0+(2.0+1.0)/2.0);
1696    M12[channel]/=pow(M00[channel],1.0+(1.0+2.0)/2.0);
1697    M22[channel]/=pow(M00[channel],1.0+(2.0+2.0)/2.0);
1698    M30[channel]/=pow(M00[channel],1.0+(3.0+0.0)/2.0);
1699    M03[channel]/=pow(M00[channel],1.0+(0.0+3.0)/2.0);
1700    M00[channel]=1.0;
1701  }
1702  image_view=DestroyCacheView(image_view);
1703  for (channel=0; channel <= MaxPixelChannels; channel++)
1704  {
1705    /*
1706      Compute Hu invariant moments.
1707    */
1708    channel_moments[channel].invariant[0]=M20[channel]+M02[channel];
1709    channel_moments[channel].invariant[1]=(M20[channel]-M02[channel])*
1710      (M20[channel]-M02[channel])+4.0*M11[channel]*M11[channel];
1711    channel_moments[channel].invariant[2]=(M30[channel]-3.0*M12[channel])*
1712      (M30[channel]-3.0*M12[channel])+(3.0*M21[channel]-M03[channel])*
1713      (3.0*M21[channel]-M03[channel]);
1714    channel_moments[channel].invariant[3]=(M30[channel]+M12[channel])*
1715      (M30[channel]+M12[channel])+(M21[channel]+M03[channel])*
1716      (M21[channel]+M03[channel]);
1717    channel_moments[channel].invariant[4]=(M30[channel]-3.0*M12[channel])*
1718      (M30[channel]+M12[channel])*((M30[channel]+M12[channel])*
1719      (M30[channel]+M12[channel])-3.0*(M21[channel]+M03[channel])*
1720      (M21[channel]+M03[channel]))+(3.0*M21[channel]-M03[channel])*
1721      (M21[channel]+M03[channel])*(3.0*(M30[channel]+M12[channel])*
1722      (M30[channel]+M12[channel])-(M21[channel]+M03[channel])*
1723      (M21[channel]+M03[channel]));
1724    channel_moments[channel].invariant[5]=(M20[channel]-M02[channel])*
1725      ((M30[channel]+M12[channel])*(M30[channel]+M12[channel])-
1726      (M21[channel]+M03[channel])*(M21[channel]+M03[channel]))+
1727      4.0*M11[channel]*(M30[channel]+M12[channel])*(M21[channel]+M03[channel]);
1728    channel_moments[channel].invariant[6]=(3.0*M21[channel]-M03[channel])*
1729      (M30[channel]+M12[channel])*((M30[channel]+M12[channel])*
1730      (M30[channel]+M12[channel])-3.0*(M21[channel]+M03[channel])*
1731      (M21[channel]+M03[channel]))-(M30[channel]-3*M12[channel])*
1732      (M21[channel]+M03[channel])*(3.0*(M30[channel]+M12[channel])*
1733      (M30[channel]+M12[channel])-(M21[channel]+M03[channel])*
1734      (M21[channel]+M03[channel]));
1735    channel_moments[channel].invariant[7]=M11[channel]*((M30[channel]+
1736      M12[channel])*(M30[channel]+M12[channel])-(M03[channel]+M21[channel])*
1737      (M03[channel]+M21[channel]))-(M20[channel]-M02[channel])*
1738      (M30[channel]+M12[channel])*(M03[channel]+M21[channel]);
1739  }
1740  if (y < (ssize_t) image->rows)
1741    channel_moments=(ChannelMoments *) RelinquishMagickMemory(channel_moments);
1742  return(channel_moments);
1743}
1744
1745/*
1746%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1747%                                                                             %
1748%                                                                             %
1749%                                                                             %
1750%   G e t I m a g e C h a n n e l P e r c e p t u a l H a s h                 %
1751%                                                                             %
1752%                                                                             %
1753%                                                                             %
1754%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1755%
1756%  GetImagePerceptualHash() returns the perceptual hash of one or more
1757%  image channels.
1758%
1759%  The format of the GetImagePerceptualHash method is:
1760%
1761%      ChannelPerceptualHash *GetImagePerceptualHash(const Image *image,
1762%        ExceptionInfo *exception)
1763%
1764%  A description of each parameter follows:
1765%
1766%    o image: the image.
1767%
1768%    o exception: return any errors or warnings in this structure.
1769%
1770*/
1771
1772static inline double MagickLog10(const double x)
1773{
1774#define Log10Epsilon  (1.0e-11)
1775
1776 if (fabs(x) < Log10Epsilon)
1777   return(log10(Log10Epsilon));
1778 return(log10(fabs(x)));
1779}
1780
1781MagickExport ChannelPerceptualHash *GetImagePerceptualHash(
1782  const Image *image,ExceptionInfo *exception)
1783{
1784  ChannelMoments
1785    *moments;
1786
1787  ChannelPerceptualHash
1788    *perceptual_hash;
1789
1790  Image
1791    *hash_image;
1792
1793  MagickBooleanType
1794    status;
1795
1796  register ssize_t
1797    i;
1798
1799  ssize_t
1800    channel;
1801
1802  /*
1803    Blur then transform to sRGB colorspace.
1804  */
1805  hash_image=BlurImage(image,0.0,1.0,exception);
1806  if (hash_image == (Image *) NULL)
1807    return((ChannelPerceptualHash *) NULL);
1808  hash_image->depth=8;
1809  status=TransformImageColorspace(hash_image,sRGBColorspace,exception);
1810  if (status == MagickFalse)
1811    return((ChannelPerceptualHash *) NULL);
1812  moments=GetImageMoments(hash_image,exception);
1813  hash_image=DestroyImage(hash_image);
1814  if (moments == (ChannelMoments *) NULL)
1815    return((ChannelPerceptualHash *) NULL);
1816  perceptual_hash=(ChannelPerceptualHash *) AcquireQuantumMemory(
1817    CompositeChannels+1UL,sizeof(*perceptual_hash));
1818  if (perceptual_hash == (ChannelPerceptualHash *) NULL)
1819    return((ChannelPerceptualHash *) NULL);
1820  for (channel=0; channel <= MaxPixelChannels; channel++)
1821    for (i=0; i < MaximumNumberOfImageMoments; i++)
1822      perceptual_hash[channel].srgb_hu_phash[i]=
1823        (-MagickLog10(moments[channel].invariant[i]));
1824  moments=(ChannelMoments *) RelinquishMagickMemory(moments);
1825  /*
1826    Blur then transform to HCLp colorspace.
1827  */
1828  hash_image=BlurImage(image,0.0,1.0,exception);
1829  if (hash_image == (Image *) NULL)
1830    {
1831      perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1832        perceptual_hash);
1833      return((ChannelPerceptualHash *) NULL);
1834    }
1835  hash_image->depth=8;
1836  status=TransformImageColorspace(hash_image,HCLpColorspace,exception);
1837  if (status == MagickFalse)
1838    {
1839      perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1840        perceptual_hash);
1841      return((ChannelPerceptualHash *) NULL);
1842    }
1843  moments=GetImageMoments(hash_image,exception);
1844  hash_image=DestroyImage(hash_image);
1845  if (moments == (ChannelMoments *) NULL)
1846    {
1847      perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1848        perceptual_hash);
1849      return((ChannelPerceptualHash *) NULL);
1850    }
1851  for (channel=0; channel <= MaxPixelChannels; channel++)
1852    for (i=0; i < MaximumNumberOfImageMoments; i++)
1853      perceptual_hash[channel].hclp_hu_phash[i]=
1854        (-MagickLog10(moments[channel].invariant[i]));
1855  moments=(ChannelMoments *) RelinquishMagickMemory(moments);
1856  return(perceptual_hash);
1857}
1858
1859/*
1860%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1861%                                                                             %
1862%                                                                             %
1863%                                                                             %
1864%   G e t I m a g e R a n g e                                                 %
1865%                                                                             %
1866%                                                                             %
1867%                                                                             %
1868%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1869%
1870%  GetImageRange() returns the range of one or more image channels.
1871%
1872%  The format of the GetImageRange method is:
1873%
1874%      MagickBooleanType GetImageRange(const Image *image,double *minima,
1875%        double *maxima,ExceptionInfo *exception)
1876%
1877%  A description of each parameter follows:
1878%
1879%    o image: the image.
1880%
1881%    o minima: the minimum value in the channel.
1882%
1883%    o maxima: the maximum value in the channel.
1884%
1885%    o exception: return any errors or warnings in this structure.
1886%
1887*/
1888MagickExport MagickBooleanType GetImageRange(const Image *image,double *minima,
1889  double *maxima,ExceptionInfo *exception)
1890{
1891  CacheView
1892    *image_view;
1893
1894  MagickBooleanType
1895    initialize,
1896    status;
1897
1898  ssize_t
1899    y;
1900
1901  assert(image != (Image *) NULL);
1902  assert(image->signature == MagickSignature);
1903  if (image->debug != MagickFalse)
1904    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1905  status=MagickTrue;
1906  initialize=MagickTrue;
1907  *maxima=0.0;
1908  *minima=0.0;
1909  image_view=AcquireVirtualCacheView(image,exception);
1910#if defined(MAGICKCORE_OPENMP_SUPPORT)
1911  #pragma omp parallel for schedule(static,4) shared(status,initialize) \
1912    magick_threads(image,image,image->rows,1)
1913#endif
1914  for (y=0; y < (ssize_t) image->rows; y++)
1915  {
1916    register const Quantum
1917      *restrict p;
1918
1919    register ssize_t
1920      x;
1921
1922    if (status == MagickFalse)
1923      continue;
1924    p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1925    if (p == (const Quantum *) NULL)
1926      {
1927        status=MagickFalse;
1928        continue;
1929      }
1930    for (x=0; x < (ssize_t) image->columns; x++)
1931    {
1932      register ssize_t
1933        i;
1934
1935      if (GetPixelReadMask(image,p) == 0)
1936        {
1937          p+=GetPixelChannels(image);
1938          continue;
1939        }
1940      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1941      {
1942        PixelChannel channel=GetPixelChannelChannel(image,i);
1943        PixelTrait traits=GetPixelChannelTraits(image,channel);
1944        if (traits == UndefinedPixelTrait)
1945          continue;
1946        if ((traits & UpdatePixelTrait) == 0)
1947          continue;
1948#if defined(MAGICKCORE_OPENMP_SUPPORT)
1949        #pragma omp critical (MagickCore_GetImageRange)
1950#endif
1951        {
1952          if (initialize != MagickFalse)
1953            {
1954              *minima=(double) p[i];
1955              *maxima=(double) p[i];
1956              initialize=MagickFalse;
1957            }
1958          else
1959            {
1960              if ((double) p[i] < *minima)
1961                *minima=(double) p[i];
1962              if ((double) p[i] > *maxima)
1963                *maxima=(double) p[i];
1964           }
1965        }
1966      }
1967      p+=GetPixelChannels(image);
1968    }
1969  }
1970  image_view=DestroyCacheView(image_view);
1971  return(status);
1972}
1973
1974/*
1975%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1976%                                                                             %
1977%                                                                             %
1978%                                                                             %
1979%   G e t I m a g e S t a t i s t i c s                                       %
1980%                                                                             %
1981%                                                                             %
1982%                                                                             %
1983%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1984%
1985%  GetImageStatistics() returns statistics for each channel in the image.  The
1986%  statistics include the channel depth, its minima, maxima, mean, standard
1987%  deviation, kurtosis and skewness.  You can access the red channel mean, for
1988%  example, like this:
1989%
1990%      channel_statistics=GetImageStatistics(image,exception);
1991%      red_mean=channel_statistics[RedPixelChannel].mean;
1992%
1993%  Use MagickRelinquishMemory() to free the statistics buffer.
1994%
1995%  The format of the GetImageStatistics method is:
1996%
1997%      ChannelStatistics *GetImageStatistics(const Image *image,
1998%        ExceptionInfo *exception)
1999%
2000%  A description of each parameter follows:
2001%
2002%    o image: the image.
2003%
2004%    o exception: return any errors or warnings in this structure.
2005%
2006*/
2007MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
2008  ExceptionInfo *exception)
2009{
2010  ChannelStatistics
2011    *channel_statistics;
2012
2013  double
2014    *histogram;
2015
2016  MagickStatusType
2017    status;
2018
2019  QuantumAny
2020    range;
2021
2022  register ssize_t
2023    i;
2024
2025  size_t
2026    channels,
2027    depth;
2028
2029  ssize_t
2030    y;
2031
2032  assert(image != (Image *) NULL);
2033  assert(image->signature == MagickSignature);
2034  if (image->debug != MagickFalse)
2035    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2036  histogram=(double *) AcquireQuantumMemory(MaxMap+1UL,GetPixelChannels(image)*
2037    sizeof(*histogram));
2038  channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(
2039    MaxPixelChannels+1,sizeof(*channel_statistics));
2040  if ((channel_statistics == (ChannelStatistics *) NULL) ||
2041      (histogram == (double *) NULL))
2042    {
2043      if (histogram != (double *) NULL)
2044        histogram=(double *) RelinquishMagickMemory(histogram);
2045      if (channel_statistics != (ChannelStatistics *) NULL)
2046        channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
2047          channel_statistics);
2048      return(channel_statistics);
2049    }
2050  (void) ResetMagickMemory(channel_statistics,0,(MaxPixelChannels+1)*
2051    sizeof(*channel_statistics));
2052  for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
2053  {
2054    channel_statistics[i].depth=1;
2055    channel_statistics[i].maxima=(-MagickMaximumValue);
2056    channel_statistics[i].minima=MagickMaximumValue;
2057  }
2058  (void) ResetMagickMemory(histogram,0,(MaxMap+1)*GetPixelChannels(image)*
2059    sizeof(*histogram));
2060  for (y=0; y < (ssize_t) image->rows; y++)
2061  {
2062    register const Quantum
2063      *restrict p;
2064
2065    register ssize_t
2066      x;
2067
2068    p=GetVirtualPixels(image,0,y,image->columns,1,exception);
2069    if (p == (const Quantum *) NULL)
2070      break;
2071    for (x=0; x < (ssize_t) image->columns; x++)
2072    {
2073      register ssize_t
2074        i;
2075
2076      if (GetPixelReadMask(image,p) == 0)
2077        {
2078          p+=GetPixelChannels(image);
2079          continue;
2080        }
2081      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2082      {
2083        PixelChannel channel=GetPixelChannelChannel(image,i);
2084        PixelTrait traits=GetPixelChannelTraits(image,channel);
2085        if (traits == UndefinedPixelTrait)
2086          continue;
2087        if (channel_statistics[channel].depth != MAGICKCORE_QUANTUM_DEPTH)
2088          {
2089            depth=channel_statistics[channel].depth;
2090            range=GetQuantumRange(depth);
2091            status=p[i] != ScaleAnyToQuantum(ScaleQuantumToAny(p[i],range),
2092              range) ? MagickTrue : MagickFalse;
2093            if (status != MagickFalse)
2094              {
2095                channel_statistics[channel].depth++;
2096                i--;
2097                continue;
2098              }
2099          }
2100        if ((double) p[i] < channel_statistics[channel].minima)
2101          channel_statistics[channel].minima=(double) p[i];
2102        if ((double) p[i] > channel_statistics[channel].maxima)
2103          channel_statistics[channel].maxima=(double) p[i];
2104        channel_statistics[channel].sum+=p[i];
2105        channel_statistics[channel].sum_squared+=(double) p[i]*p[i];
2106        channel_statistics[channel].sum_cubed+=(double) p[i]*p[i]*p[i];
2107        channel_statistics[channel].sum_fourth_power+=(double) p[i]*p[i]*p[i]*
2108          p[i];
2109        channel_statistics[channel].area++;
2110        histogram[GetPixelChannels(image)*ScaleQuantumToMap(
2111          ClampToQuantum((double) p[i]))+i]++;
2112      }
2113      p+=GetPixelChannels(image);
2114    }
2115  }
2116  for (i=0; i < (ssize_t) MaxPixelChannels; i++)
2117  {
2118    double
2119      area;
2120
2121    register ssize_t
2122      j;
2123
2124    area=PerceptibleReciprocal(channel_statistics[i].area);
2125    channel_statistics[i].sum*=area;
2126    channel_statistics[i].sum_squared*=area;
2127    channel_statistics[i].sum_cubed*=area;
2128    channel_statistics[i].sum_fourth_power*=area;
2129    channel_statistics[i].mean=channel_statistics[i].sum;
2130    channel_statistics[i].variance=channel_statistics[i].sum_squared;
2131    channel_statistics[i].standard_deviation=sqrt(
2132      channel_statistics[i].variance-(channel_statistics[i].mean*
2133      channel_statistics[i].mean));
2134    for (j=0; j < (ssize_t) (MaxMap+1U); j++)
2135    {
2136      double
2137        count;
2138
2139      count=histogram[GetPixelChannels(image)*j+i]*area;
2140      channel_statistics[i].entropy+=-count*MagickLog10(count)/
2141        MagickLog10(MaxMap+1.0);
2142    }
2143  }
2144  for (i=0; i < (ssize_t) MaxPixelChannels; i++)
2145  {
2146    channel_statistics[CompositePixelChannel].area+=channel_statistics[i].area;
2147    channel_statistics[CompositePixelChannel].minima=MagickMin(
2148      channel_statistics[CompositePixelChannel].minima,
2149      channel_statistics[i].minima);
2150    channel_statistics[CompositePixelChannel].maxima=EvaluateMax(
2151      channel_statistics[CompositePixelChannel].maxima,
2152      channel_statistics[i].maxima);
2153    channel_statistics[CompositePixelChannel].sum+=channel_statistics[i].sum;
2154    channel_statistics[CompositePixelChannel].sum_squared+=
2155      channel_statistics[i].sum_squared;
2156    channel_statistics[CompositePixelChannel].sum_cubed+=
2157      channel_statistics[i].sum_cubed;
2158    channel_statistics[CompositePixelChannel].sum_fourth_power+=
2159      channel_statistics[i].sum_fourth_power;
2160    channel_statistics[CompositePixelChannel].mean+=channel_statistics[i].mean;
2161    channel_statistics[CompositePixelChannel].variance+=
2162      channel_statistics[i].variance-channel_statistics[i].mean*
2163      channel_statistics[i].mean;
2164    channel_statistics[CompositePixelChannel].standard_deviation+=
2165      channel_statistics[i].variance-channel_statistics[i].mean*
2166      channel_statistics[i].mean;
2167    if (channel_statistics[i].entropy > MagickEpsilon)
2168      channel_statistics[CompositePixelChannel].entropy+=
2169        channel_statistics[i].entropy;
2170  }
2171  channels=GetImageChannels(image);
2172  channel_statistics[CompositePixelChannel].area/=channels;
2173  channel_statistics[CompositePixelChannel].sum/=channels;
2174  channel_statistics[CompositePixelChannel].sum_squared/=channels;
2175  channel_statistics[CompositePixelChannel].sum_cubed/=channels;
2176  channel_statistics[CompositePixelChannel].sum_fourth_power/=channels;
2177  channel_statistics[CompositePixelChannel].mean/=channels;
2178  channel_statistics[CompositePixelChannel].variance/=channels;
2179  channel_statistics[CompositePixelChannel].standard_deviation=
2180    sqrt(channel_statistics[CompositePixelChannel].standard_deviation/channels);
2181  channel_statistics[CompositePixelChannel].kurtosis/=channels;
2182  channel_statistics[CompositePixelChannel].skewness/=channels;
2183  channel_statistics[CompositePixelChannel].entropy/=channels;
2184  for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
2185  {
2186    double
2187      standard_deviation;
2188
2189    if (channel_statistics[i].standard_deviation == 0.0)
2190      continue;
2191    standard_deviation=PerceptibleReciprocal(
2192      channel_statistics[i].standard_deviation);
2193    channel_statistics[i].skewness=(channel_statistics[i].sum_cubed-3.0*
2194      channel_statistics[i].mean*channel_statistics[i].sum_squared+2.0*
2195      channel_statistics[i].mean*channel_statistics[i].mean*
2196      channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2197      standard_deviation);
2198    channel_statistics[i].kurtosis=(channel_statistics[i].sum_fourth_power-4.0*
2199      channel_statistics[i].mean*channel_statistics[i].sum_cubed+6.0*
2200      channel_statistics[i].mean*channel_statistics[i].mean*
2201      channel_statistics[i].sum_squared-3.0*channel_statistics[i].mean*
2202      channel_statistics[i].mean*1.0*channel_statistics[i].mean*
2203      channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2204      standard_deviation*standard_deviation)-3.0;
2205  }
2206  histogram=(double *) RelinquishMagickMemory(histogram);
2207  if (y < (ssize_t) image->rows)
2208    channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
2209      channel_statistics);
2210  return(channel_statistics);
2211}
2212
2213/*
2214%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2215%                                                                             %
2216%                                                                             %
2217%                                                                             %
2218%     P o l y n o m i a l I m a g e                                           %
2219%                                                                             %
2220%                                                                             %
2221%                                                                             %
2222%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2223%
2224%  PolynomialImage() returns a new image where each pixel is the sum of the
2225%  pixels in the image sequence after applying its corresponding terms
2226%  (coefficient and degree pairs).
2227%
2228%  The format of the PolynomialImage method is:
2229%
2230%      Image *PolynomialImage(const Image *images,const size_t number_terms,
2231%        const double *terms,ExceptionInfo *exception)
2232%
2233%  A description of each parameter follows:
2234%
2235%    o images: the image sequence.
2236%
2237%    o number_terms: the number of terms in the list.  The actual list length
2238%      is 2 x number_terms + 1 (the constant).
2239%
2240%    o terms: the list of polynomial coefficients and degree pairs and a
2241%      constant.
2242%
2243%    o exception: return any errors or warnings in this structure.
2244%
2245*/
2246
2247MagickExport Image *PolynomialImage(const Image *images,
2248  const size_t number_terms,const double *terms,ExceptionInfo *exception)
2249{
2250#define PolynomialImageTag  "Polynomial/Image"
2251
2252  CacheView
2253    *polynomial_view;
2254
2255  Image
2256    *image;
2257
2258  MagickBooleanType
2259    status;
2260
2261  MagickOffsetType
2262    progress;
2263
2264  PixelChannels
2265    **restrict polynomial_pixels;
2266
2267  size_t
2268    number_images;
2269
2270  ssize_t
2271    y;
2272
2273  assert(images != (Image *) NULL);
2274  assert(images->signature == MagickSignature);
2275  if (images->debug != MagickFalse)
2276    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
2277  assert(exception != (ExceptionInfo *) NULL);
2278  assert(exception->signature == MagickSignature);
2279  image=CloneImage(images,images->columns,images->rows,MagickTrue,
2280    exception);
2281  if (image == (Image *) NULL)
2282    return((Image *) NULL);
2283  if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
2284    {
2285      image=DestroyImage(image);
2286      return((Image *) NULL);
2287    }
2288  number_images=GetImageListLength(images);
2289  polynomial_pixels=AcquirePixelThreadSet(images,number_images);
2290  if (polynomial_pixels == (PixelChannels **) NULL)
2291    {
2292      image=DestroyImage(image);
2293      (void) ThrowMagickException(exception,GetMagickModule(),
2294        ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
2295      return((Image *) NULL);
2296    }
2297  /*
2298    Polynomial image pixels.
2299  */
2300  status=MagickTrue;
2301  progress=0;
2302  polynomial_view=AcquireAuthenticCacheView(image,exception);
2303#if defined(MAGICKCORE_OPENMP_SUPPORT)
2304  #pragma omp parallel for schedule(static,4) shared(progress,status) \
2305    magick_threads(image,image,image->rows,1)
2306#endif
2307  for (y=0; y < (ssize_t) image->rows; y++)
2308  {
2309    CacheView
2310      *image_view;
2311
2312    const Image
2313      *next;
2314
2315    const int
2316      id = GetOpenMPThreadId();
2317
2318    register ssize_t
2319      i,
2320      x;
2321
2322    register PixelChannels
2323      *polynomial_pixel;
2324
2325    register Quantum
2326      *restrict q;
2327
2328    ssize_t
2329      j;
2330
2331    if (status == MagickFalse)
2332      continue;
2333    q=QueueCacheViewAuthenticPixels(polynomial_view,0,y,image->columns,1,
2334      exception);
2335    if (q == (Quantum *) NULL)
2336      {
2337        status=MagickFalse;
2338        continue;
2339      }
2340    polynomial_pixel=polynomial_pixels[id];
2341    for (j=0; j < (ssize_t) image->columns; j++)
2342      for (i=0; i < MaxPixelChannels; i++)
2343        polynomial_pixel[j].channel[i]=0.0;
2344    next=images;
2345    for (j=0; j < (ssize_t) number_images; j++)
2346    {
2347      register const Quantum
2348        *p;
2349
2350      if (j >= (ssize_t) number_terms)
2351        continue;
2352      image_view=AcquireVirtualCacheView(next,exception);
2353      p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2354      if (p == (const Quantum *) NULL)
2355        {
2356          image_view=DestroyCacheView(image_view);
2357          break;
2358        }
2359      for (x=0; x < (ssize_t) image->columns; x++)
2360      {
2361        register ssize_t
2362          i;
2363
2364        if (GetPixelReadMask(next,p) == 0)
2365          {
2366            p+=GetPixelChannels(next);
2367            continue;
2368          }
2369        for (i=0; i < (ssize_t) GetPixelChannels(next); i++)
2370        {
2371          MagickRealType
2372            coefficient,
2373            degree;
2374
2375          PixelChannel channel=GetPixelChannelChannel(image,i);
2376          PixelTrait traits=GetPixelChannelTraits(next,channel);
2377          PixelTrait polynomial_traits=GetPixelChannelTraits(image,channel);
2378          if ((traits == UndefinedPixelTrait) ||
2379              (polynomial_traits == UndefinedPixelTrait))
2380            continue;
2381          if ((traits & UpdatePixelTrait) == 0)
2382            continue;
2383          coefficient=(MagickRealType) terms[2*i];
2384          degree=(MagickRealType) terms[(i << 1)+1];
2385          polynomial_pixel[x].channel[i]+=coefficient*
2386            pow(QuantumScale*GetPixelChannel(image,channel,p),degree);
2387        }
2388        p+=GetPixelChannels(next);
2389      }
2390      image_view=DestroyCacheView(image_view);
2391      next=GetNextImageInList(next);
2392    }
2393    for (x=0; x < (ssize_t) image->columns; x++)
2394    {
2395      register ssize_t
2396        i;
2397
2398      if (GetPixelReadMask(image,q) == 0)
2399        {
2400          q+=GetPixelChannels(image);
2401          continue;
2402        }
2403      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2404      {
2405        PixelChannel channel=GetPixelChannelChannel(image,i);
2406        PixelTrait traits=GetPixelChannelTraits(image,channel);
2407        if (traits == UndefinedPixelTrait)
2408          continue;
2409        if ((traits & UpdatePixelTrait) == 0)
2410          continue;
2411        q[i]=ClampToQuantum(QuantumRange*polynomial_pixel[x].channel[i]);
2412      }
2413      q+=GetPixelChannels(image);
2414    }
2415    if (SyncCacheViewAuthenticPixels(polynomial_view,exception) == MagickFalse)
2416      status=MagickFalse;
2417    if (images->progress_monitor != (MagickProgressMonitor) NULL)
2418      {
2419        MagickBooleanType
2420          proceed;
2421
2422#if defined(MAGICKCORE_OPENMP_SUPPORT)
2423        #pragma omp critical (MagickCore_PolynomialImages)
2424#endif
2425        proceed=SetImageProgress(images,PolynomialImageTag,progress++,
2426          image->rows);
2427        if (proceed == MagickFalse)
2428          status=MagickFalse;
2429      }
2430  }
2431  polynomial_view=DestroyCacheView(polynomial_view);
2432  polynomial_pixels=DestroyPixelThreadSet(polynomial_pixels);
2433  if (status == MagickFalse)
2434    image=DestroyImage(image);
2435  return(image);
2436}
2437
2438/*
2439%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2440%                                                                             %
2441%                                                                             %
2442%                                                                             %
2443%     S t a t i s t i c I m a g e                                             %
2444%                                                                             %
2445%                                                                             %
2446%                                                                             %
2447%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2448%
2449%  StatisticImage() makes each pixel the min / max / median / mode / etc. of
2450%  the neighborhood of the specified width and height.
2451%
2452%  The format of the StatisticImage method is:
2453%
2454%      Image *StatisticImage(const Image *image,const StatisticType type,
2455%        const size_t width,const size_t height,ExceptionInfo *exception)
2456%
2457%  A description of each parameter follows:
2458%
2459%    o image: the image.
2460%
2461%    o type: the statistic type (median, mode, etc.).
2462%
2463%    o width: the width of the pixel neighborhood.
2464%
2465%    o height: the height of the pixel neighborhood.
2466%
2467%    o exception: return any errors or warnings in this structure.
2468%
2469*/
2470
2471typedef struct _SkipNode
2472{
2473  size_t
2474    next[9],
2475    count,
2476    signature;
2477} SkipNode;
2478
2479typedef struct _SkipList
2480{
2481  ssize_t
2482    level;
2483
2484  SkipNode
2485    *nodes;
2486} SkipList;
2487
2488typedef struct _PixelList
2489{
2490  size_t
2491    length,
2492    seed;
2493
2494  SkipList
2495    skip_list;
2496
2497  size_t
2498    signature;
2499} PixelList;
2500
2501static PixelList *DestroyPixelList(PixelList *pixel_list)
2502{
2503  if (pixel_list == (PixelList *) NULL)
2504    return((PixelList *) NULL);
2505  if (pixel_list->skip_list.nodes != (SkipNode *) NULL)
2506    pixel_list->skip_list.nodes=(SkipNode *) RelinquishMagickMemory(
2507      pixel_list->skip_list.nodes);
2508  pixel_list=(PixelList *) RelinquishMagickMemory(pixel_list);
2509  return(pixel_list);
2510}
2511
2512static PixelList **DestroyPixelListThreadSet(PixelList **pixel_list)
2513{
2514  register ssize_t
2515    i;
2516
2517  assert(pixel_list != (PixelList **) NULL);
2518  for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
2519    if (pixel_list[i] != (PixelList *) NULL)
2520      pixel_list[i]=DestroyPixelList(pixel_list[i]);
2521  pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list);
2522  return(pixel_list);
2523}
2524
2525static PixelList *AcquirePixelList(const size_t width,const size_t height)
2526{
2527  PixelList
2528    *pixel_list;
2529
2530  pixel_list=(PixelList *) AcquireMagickMemory(sizeof(*pixel_list));
2531  if (pixel_list == (PixelList *) NULL)
2532    return(pixel_list);
2533  (void) ResetMagickMemory((void *) pixel_list,0,sizeof(*pixel_list));
2534  pixel_list->length=width*height;
2535  pixel_list->skip_list.nodes=(SkipNode *) AcquireQuantumMemory(65537UL,
2536    sizeof(*pixel_list->skip_list.nodes));
2537  if (pixel_list->skip_list.nodes == (SkipNode *) NULL)
2538    return(DestroyPixelList(pixel_list));
2539  (void) ResetMagickMemory(pixel_list->skip_list.nodes,0,65537UL*
2540    sizeof(*pixel_list->skip_list.nodes));
2541  pixel_list->signature=MagickSignature;
2542  return(pixel_list);
2543}
2544
2545static PixelList **AcquirePixelListThreadSet(const size_t width,
2546  const size_t height)
2547{
2548  PixelList
2549    **pixel_list;
2550
2551  register ssize_t
2552    i;
2553
2554  size_t
2555    number_threads;
2556
2557  number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
2558  pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
2559    sizeof(*pixel_list));
2560  if (pixel_list == (PixelList **) NULL)
2561    return((PixelList **) NULL);
2562  (void) ResetMagickMemory(pixel_list,0,number_threads*sizeof(*pixel_list));
2563  for (i=0; i < (ssize_t) number_threads; i++)
2564  {
2565    pixel_list[i]=AcquirePixelList(width,height);
2566    if (pixel_list[i] == (PixelList *) NULL)
2567      return(DestroyPixelListThreadSet(pixel_list));
2568  }
2569  return(pixel_list);
2570}
2571
2572static void AddNodePixelList(PixelList *pixel_list,const size_t color)
2573{
2574  register SkipList
2575    *p;
2576
2577  register ssize_t
2578    level;
2579
2580  size_t
2581    search,
2582    update[9];
2583
2584  /*
2585    Initialize the node.
2586  */
2587  p=(&pixel_list->skip_list);
2588  p->nodes[color].signature=pixel_list->signature;
2589  p->nodes[color].count=1;
2590  /*
2591    Determine where it belongs in the list.
2592  */
2593  search=65536UL;
2594  for (level=p->level; level >= 0; level--)
2595  {
2596    while (p->nodes[search].next[level] < color)
2597      search=p->nodes[search].next[level];
2598    update[level]=search;
2599  }
2600  /*
2601    Generate a pseudo-random level for this node.
2602  */
2603  for (level=0; ; level++)
2604  {
2605    pixel_list->seed=(pixel_list->seed*42893621L)+1L;
2606    if ((pixel_list->seed & 0x300) != 0x300)
2607      break;
2608  }
2609  if (level > 8)
2610    level=8;
2611  if (level > (p->level+2))
2612    level=p->level+2;
2613  /*
2614    If we're raising the list's level, link back to the root node.
2615  */
2616  while (level > p->level)
2617  {
2618    p->level++;
2619    update[p->level]=65536UL;
2620  }
2621  /*
2622    Link the node into the skip-list.
2623  */
2624  do
2625  {
2626    p->nodes[color].next[level]=p->nodes[update[level]].next[level];
2627    p->nodes[update[level]].next[level]=color;
2628  } while (level-- > 0);
2629}
2630
2631static inline void GetMaximumPixelList(PixelList *pixel_list,Quantum *pixel)
2632{
2633  register SkipList
2634    *p;
2635
2636  size_t
2637    color,
2638    maximum;
2639
2640  ssize_t
2641    count;
2642
2643  /*
2644    Find the maximum value for each of the color.
2645  */
2646  p=(&pixel_list->skip_list);
2647  color=65536L;
2648  count=0;
2649  maximum=p->nodes[color].next[0];
2650  do
2651  {
2652    color=p->nodes[color].next[0];
2653    if (color > maximum)
2654      maximum=color;
2655    count+=p->nodes[color].count;
2656  } while (count < (ssize_t) pixel_list->length);
2657  *pixel=ScaleShortToQuantum((unsigned short) maximum);
2658}
2659
2660static inline void GetMeanPixelList(PixelList *pixel_list,Quantum *pixel)
2661{
2662  double
2663    sum;
2664
2665  register SkipList
2666    *p;
2667
2668  size_t
2669    color;
2670
2671  ssize_t
2672    count;
2673
2674  /*
2675    Find the mean value for each of the color.
2676  */
2677  p=(&pixel_list->skip_list);
2678  color=65536L;
2679  count=0;
2680  sum=0.0;
2681  do
2682  {
2683    color=p->nodes[color].next[0];
2684    sum+=(double) p->nodes[color].count*color;
2685    count+=p->nodes[color].count;
2686  } while (count < (ssize_t) pixel_list->length);
2687  sum/=pixel_list->length;
2688  *pixel=ScaleShortToQuantum((unsigned short) sum);
2689}
2690
2691static inline void GetMedianPixelList(PixelList *pixel_list,Quantum *pixel)
2692{
2693  register SkipList
2694    *p;
2695
2696  size_t
2697    color;
2698
2699  ssize_t
2700    count;
2701
2702  /*
2703    Find the median value for each of the color.
2704  */
2705  p=(&pixel_list->skip_list);
2706  color=65536L;
2707  count=0;
2708  do
2709  {
2710    color=p->nodes[color].next[0];
2711    count+=p->nodes[color].count;
2712  } while (count <= (ssize_t) (pixel_list->length >> 1));
2713  *pixel=ScaleShortToQuantum((unsigned short) color);
2714}
2715
2716static inline void GetMinimumPixelList(PixelList *pixel_list,Quantum *pixel)
2717{
2718  register SkipList
2719    *p;
2720
2721  size_t
2722    color,
2723    minimum;
2724
2725  ssize_t
2726    count;
2727
2728  /*
2729    Find the minimum value for each of the color.
2730  */
2731  p=(&pixel_list->skip_list);
2732  count=0;
2733  color=65536UL;
2734  minimum=p->nodes[color].next[0];
2735  do
2736  {
2737    color=p->nodes[color].next[0];
2738    if (color < minimum)
2739      minimum=color;
2740    count+=p->nodes[color].count;
2741  } while (count < (ssize_t) pixel_list->length);
2742  *pixel=ScaleShortToQuantum((unsigned short) minimum);
2743}
2744
2745static inline void GetModePixelList(PixelList *pixel_list,Quantum *pixel)
2746{
2747  register SkipList
2748    *p;
2749
2750  size_t
2751    color,
2752    max_count,
2753    mode;
2754
2755  ssize_t
2756    count;
2757
2758  /*
2759    Make each pixel the 'predominant color' of the specified neighborhood.
2760  */
2761  p=(&pixel_list->skip_list);
2762  color=65536L;
2763  mode=color;
2764  max_count=p->nodes[mode].count;
2765  count=0;
2766  do
2767  {
2768    color=p->nodes[color].next[0];
2769    if (p->nodes[color].count > max_count)
2770      {
2771        mode=color;
2772        max_count=p->nodes[mode].count;
2773      }
2774    count+=p->nodes[color].count;
2775  } while (count < (ssize_t) pixel_list->length);
2776  *pixel=ScaleShortToQuantum((unsigned short) mode);
2777}
2778
2779static inline void GetNonpeakPixelList(PixelList *pixel_list,Quantum *pixel)
2780{
2781  register SkipList
2782    *p;
2783
2784  size_t
2785    color,
2786    next,
2787    previous;
2788
2789  ssize_t
2790    count;
2791
2792  /*
2793    Finds the non peak value for each of the colors.
2794  */
2795  p=(&pixel_list->skip_list);
2796  color=65536L;
2797  next=p->nodes[color].next[0];
2798  count=0;
2799  do
2800  {
2801    previous=color;
2802    color=next;
2803    next=p->nodes[color].next[0];
2804    count+=p->nodes[color].count;
2805  } while (count <= (ssize_t) (pixel_list->length >> 1));
2806  if ((previous == 65536UL) && (next != 65536UL))
2807    color=next;
2808  else
2809    if ((previous != 65536UL) && (next == 65536UL))
2810      color=previous;
2811  *pixel=ScaleShortToQuantum((unsigned short) color);
2812}
2813
2814static inline void GetRootMeanSquarePixelList(PixelList *pixel_list,
2815  Quantum *pixel)
2816{
2817  double
2818    sum;
2819
2820  register SkipList
2821    *p;
2822
2823  size_t
2824    color;
2825
2826  ssize_t
2827    count;
2828
2829  /*
2830    Find the root mean square value for each of the color.
2831  */
2832  p=(&pixel_list->skip_list);
2833  color=65536L;
2834  count=0;
2835  sum=0.0;
2836  do
2837  {
2838    color=p->nodes[color].next[0];
2839    sum+=(double) (p->nodes[color].count*color*color);
2840    count+=p->nodes[color].count;
2841  } while (count < (ssize_t) pixel_list->length);
2842  sum/=pixel_list->length;
2843  *pixel=ScaleShortToQuantum((unsigned short) sqrt(sum));
2844}
2845
2846static inline void GetStandardDeviationPixelList(PixelList *pixel_list,
2847  Quantum *pixel)
2848{
2849  double
2850    sum,
2851    sum_squared;
2852
2853  register SkipList
2854    *p;
2855
2856  size_t
2857    color;
2858
2859  ssize_t
2860    count;
2861
2862  /*
2863    Find the standard-deviation value for each of the color.
2864  */
2865  p=(&pixel_list->skip_list);
2866  color=65536L;
2867  count=0;
2868  sum=0.0;
2869  sum_squared=0.0;
2870  do
2871  {
2872    register ssize_t
2873      i;
2874
2875    color=p->nodes[color].next[0];
2876    sum+=(double) p->nodes[color].count*color;
2877    for (i=0; i < (ssize_t) p->nodes[color].count; i++)
2878      sum_squared+=((double) color)*((double) color);
2879    count+=p->nodes[color].count;
2880  } while (count < (ssize_t) pixel_list->length);
2881  sum/=pixel_list->length;
2882  sum_squared/=pixel_list->length;
2883  *pixel=ScaleShortToQuantum((unsigned short) sqrt(sum_squared-(sum*sum)));
2884}
2885
2886static inline void InsertPixelList(const Quantum pixel,PixelList *pixel_list)
2887{
2888  size_t
2889    signature;
2890
2891  unsigned short
2892    index;
2893
2894  index=ScaleQuantumToShort(pixel);
2895  signature=pixel_list->skip_list.nodes[index].signature;
2896  if (signature == pixel_list->signature)
2897    {
2898      pixel_list->skip_list.nodes[index].count++;
2899      return;
2900    }
2901  AddNodePixelList(pixel_list,index);
2902}
2903
2904static inline double MagickAbsoluteValue(const double x)
2905{
2906  if (x < 0)
2907    return(-x);
2908  return(x);
2909}
2910
2911static inline size_t MagickMax(const size_t x,const size_t y)
2912{
2913  if (x > y)
2914    return(x);
2915  return(y);
2916}
2917
2918static void ResetPixelList(PixelList *pixel_list)
2919{
2920  int
2921    level;
2922
2923  register SkipNode
2924    *root;
2925
2926  register SkipList
2927    *p;
2928
2929  /*
2930    Reset the skip-list.
2931  */
2932  p=(&pixel_list->skip_list);
2933  root=p->nodes+65536UL;
2934  p->level=0;
2935  for (level=0; level < 9; level++)
2936    root->next[level]=65536UL;
2937  pixel_list->seed=pixel_list->signature++;
2938}
2939
2940MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
2941  const size_t width,const size_t height,ExceptionInfo *exception)
2942{
2943#define StatisticImageTag  "Statistic/Image"
2944
2945  CacheView
2946    *image_view,
2947    *statistic_view;
2948
2949  Image
2950    *statistic_image;
2951
2952  MagickBooleanType
2953    status;
2954
2955  MagickOffsetType
2956    progress;
2957
2958  PixelList
2959    **restrict pixel_list;
2960
2961  ssize_t
2962    center,
2963    y;
2964
2965  /*
2966    Initialize statistics image attributes.
2967  */
2968  assert(image != (Image *) NULL);
2969  assert(image->signature == MagickSignature);
2970  if (image->debug != MagickFalse)
2971    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2972  assert(exception != (ExceptionInfo *) NULL);
2973  assert(exception->signature == MagickSignature);
2974  statistic_image=CloneImage(image,image->columns,image->rows,MagickTrue,
2975    exception);
2976  if (statistic_image == (Image *) NULL)
2977    return((Image *) NULL);
2978  status=SetImageStorageClass(statistic_image,DirectClass,exception);
2979  if (status == MagickFalse)
2980    {
2981      statistic_image=DestroyImage(statistic_image);
2982      return((Image *) NULL);
2983    }
2984  pixel_list=AcquirePixelListThreadSet(MagickMax(width,1),MagickMax(height,1));
2985  if (pixel_list == (PixelList **) NULL)
2986    {
2987      statistic_image=DestroyImage(statistic_image);
2988      ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2989    }
2990  /*
2991    Make each pixel the min / max / median / mode / etc. of the neighborhood.
2992  */
2993  center=(ssize_t) GetPixelChannels(image)*(image->columns+MagickMax(width,1))*
2994    (MagickMax(height,1)/2L)+GetPixelChannels(image)*(MagickMax(width,1)/2L);
2995  status=MagickTrue;
2996  progress=0;
2997  image_view=AcquireVirtualCacheView(image,exception);
2998  statistic_view=AcquireAuthenticCacheView(statistic_image,exception);
2999#if defined(MAGICKCORE_OPENMP_SUPPORT)
3000  #pragma omp parallel for schedule(static,4) shared(progress,status) \
3001    magick_threads(image,statistic_image,statistic_image->rows,1)
3002#endif
3003  for (y=0; y < (ssize_t) statistic_image->rows; y++)
3004  {
3005    const int
3006      id = GetOpenMPThreadId();
3007
3008    register const Quantum
3009      *restrict p;
3010
3011    register Quantum
3012      *restrict q;
3013
3014    register ssize_t
3015      x;
3016
3017    if (status == MagickFalse)
3018      continue;
3019    p=GetCacheViewVirtualPixels(image_view,-((ssize_t) MagickMax(width,1)/2L),y-
3020      (ssize_t) (MagickMax(height,1)/2L),image->columns+MagickMax(width,1),
3021      MagickMax(height,1),exception);
3022    q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns,      1,exception);
3023    if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3024      {
3025        status=MagickFalse;
3026        continue;
3027      }
3028    for (x=0; x < (ssize_t) statistic_image->columns; x++)
3029    {
3030      register ssize_t
3031        i;
3032
3033      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3034      {
3035        Quantum
3036          pixel;
3037
3038        register const Quantum
3039          *restrict pixels;
3040
3041        register ssize_t
3042          u;
3043
3044        ssize_t
3045          v;
3046
3047        PixelChannel channel=GetPixelChannelChannel(image,i);
3048        PixelTrait traits=GetPixelChannelTraits(image,channel);
3049        PixelTrait statistic_traits=GetPixelChannelTraits(statistic_image,
3050          channel);
3051        if ((traits == UndefinedPixelTrait) ||
3052            (statistic_traits == UndefinedPixelTrait))
3053          continue;
3054        if (((statistic_traits & CopyPixelTrait) != 0) ||
3055            (GetPixelReadMask(image,p) == 0))
3056          {
3057            SetPixelChannel(statistic_image,channel,p[center+i],q);
3058            continue;
3059          }
3060        pixels=p;
3061        ResetPixelList(pixel_list[id]);
3062        for (v=0; v < (ssize_t) MagickMax(height,1); v++)
3063        {
3064          for (u=0; u < (ssize_t) MagickMax(width,1); u++)
3065          {
3066            InsertPixelList(pixels[i],pixel_list[id]);
3067            pixels+=GetPixelChannels(image);
3068          }
3069          pixels+=(image->columns-1)*GetPixelChannels(image);
3070        }
3071        switch (type)
3072        {
3073          case GradientStatistic:
3074          {
3075            double
3076              maximum,
3077              minimum;
3078
3079            GetMinimumPixelList(pixel_list[id],&pixel);
3080            minimum=(double) pixel;
3081            GetMaximumPixelList(pixel_list[id],&pixel);
3082            maximum=(double) pixel;
3083            pixel=ClampToQuantum(MagickAbsoluteValue(maximum-minimum));
3084            break;
3085          }
3086          case MaximumStatistic:
3087          {
3088            GetMaximumPixelList(pixel_list[id],&pixel);
3089            break;
3090          }
3091          case MeanStatistic:
3092          {
3093            GetMeanPixelList(pixel_list[id],&pixel);
3094            break;
3095          }
3096          case MedianStatistic:
3097          default:
3098          {
3099            GetMedianPixelList(pixel_list[id],&pixel);
3100            break;
3101          }
3102          case MinimumStatistic:
3103          {
3104            GetMinimumPixelList(pixel_list[id],&pixel);
3105            break;
3106          }
3107          case ModeStatistic:
3108          {
3109            GetModePixelList(pixel_list[id],&pixel);
3110            break;
3111          }
3112          case NonpeakStatistic:
3113          {
3114            GetNonpeakPixelList(pixel_list[id],&pixel);
3115            break;
3116          }
3117          case RootMeanSquareStatistic:
3118          {
3119            GetRootMeanSquarePixelList(pixel_list[id],&pixel);
3120            break;
3121          }
3122          case StandardDeviationStatistic:
3123          {
3124            GetStandardDeviationPixelList(pixel_list[id],&pixel);
3125            break;
3126          }
3127        }
3128        SetPixelChannel(statistic_image,channel,pixel,q);
3129      }
3130      p+=GetPixelChannels(image);
3131      q+=GetPixelChannels(statistic_image);
3132    }
3133    if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
3134      status=MagickFalse;
3135    if (image->progress_monitor != (MagickProgressMonitor) NULL)
3136      {
3137        MagickBooleanType
3138          proceed;
3139
3140#if defined(MAGICKCORE_OPENMP_SUPPORT)
3141        #pragma omp critical (MagickCore_StatisticImage)
3142#endif
3143        proceed=SetImageProgress(image,StatisticImageTag,progress++,
3144          image->rows);
3145        if (proceed == MagickFalse)
3146          status=MagickFalse;
3147      }
3148  }
3149  statistic_view=DestroyCacheView(statistic_view);
3150  image_view=DestroyCacheView(image_view);
3151  pixel_list=DestroyPixelListThreadSet(pixel_list);
3152  if (status == MagickFalse)
3153    statistic_image=DestroyImage(statistic_image);
3154  return(statistic_image);
3155}
3156