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