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