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