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