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  if (GetPixelReadMask(image,p) <= (QuantumRange/2))
1992  {
1993  p+=GetPixelChannels(image);
1994  continue;
1995  }
1996  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1997  {
1998  PixelChannel channel = GetPixelChannelChannel(image,i);
1999  PixelTrait traits = GetPixelChannelTraits(image,channel);
2000  if (traits == UndefinedPixelTrait)
2001  continue;
2002  if ((traits & UpdatePixelTrait) == 0)
2003  continue;
2004  if (channel_statistics[channel].depth != MAGICKCORE_QUANTUM_DEPTH)
2005  {
2006  depth=channel_statistics[channel].depth;
2007  range=GetQuantumRange(depth);
2008  status=p[i] != ScaleAnyToQuantum(ScaleQuantumToAny(p[i],range),
2009  range) ? MagickTrue : MagickFalse;
2010  if (status != MagickFalse)
2011  {
2012  channel_statistics[channel].depth++;
2013  i--;
2014  continue;
2015  }
2016  }
2017  if ((double) p[i] < channel_statistics[channel].minima)
2018  channel_statistics[channel].minima=(double) p[i];
2019  if ((double) p[i] > channel_statistics[channel].maxima)
2020  channel_statistics[channel].maxima=(double) p[i];
2021  channel_statistics[channel].sum+=p[i];
2022  channel_statistics[channel].sum_squared+=(double) p[i]*p[i];
2023  channel_statistics[channel].sum_cubed+=(double) p[i]*p[i]*p[i];
2024  channel_statistics[channel].sum_fourth_power+=(double) p[i]*p[i]*p[i]*
2025  p[i];
2026  channel_statistics[channel].area++;
2027  if ((double) p[i] < channel_statistics[CompositePixelChannel].minima)
2028  channel_statistics[CompositePixelChannel].minima=(double) p[i];
2029  if ((double) p[i] > channel_statistics[CompositePixelChannel].maxima)
2030  channel_statistics[CompositePixelChannel].maxima=(double) p[i];
2031  histogram[GetPixelChannels(image)*ScaleQuantumToMap(
2032  ClampToQuantum((double) p[i]))+i]++;
2033  channel_statistics[CompositePixelChannel].sum+=(double) p[i];
2034  channel_statistics[CompositePixelChannel].sum_squared+=(double)
2035  p[i]*p[i];
2036  channel_statistics[CompositePixelChannel].sum_cubed+=(double)
2037  p[i]*p[i]*p[i];
2038  channel_statistics[CompositePixelChannel].sum_fourth_power+=(double)
2039  p[i]*p[i]*p[i]*p[i];
2040  channel_statistics[CompositePixelChannel].area++;
2041  }
2042  p+=GetPixelChannels(image);
2043  }
2044  }
2045  for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
2046  {
2047  /*
2048  Normalize pixel statistics.
2049  */
2050  area=PerceptibleReciprocal(channel_statistics[i].area);
2051  channel_statistics[i].sum*=area;
2052  channel_statistics[i].sum_squared*=area;
2053  channel_statistics[i].sum_cubed*=area;
2054  channel_statistics[i].sum_fourth_power*=area;
2055  channel_statistics[i].mean=channel_statistics[i].sum;
2056  channel_statistics[i].variance=channel_statistics[i].sum_squared;
2057  standard_deviation=sqrt(channel_statistics[i].variance-
2058  (channel_statistics[i].mean*channel_statistics[i].mean));
2059  standard_deviation=sqrt(PerceptibleReciprocal(channel_statistics[i].area-
2060  1.0)*channel_statistics[i].area*standard_deviation*standard_deviation);
2061  channel_statistics[i].standard_deviation=standard_deviation;
2062  }
2063  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2064  {
2065  double
2066  number_bins;
2067 
2068  register ssize_t
2069  j;
2070 
2071  /*
2072  Compute pixel entropy.
2073  */
2074  PixelChannel channel = GetPixelChannelChannel(image,i);
2075  number_bins=0.0;
2076  for (j=0; j <= (ssize_t) MaxMap; j++)
2077  if (histogram[GetPixelChannels(image)*j+i] > 0.0)
2078  number_bins++;
2079  area=PerceptibleReciprocal(channel_statistics[channel].area);
2080  for (j=0; j <= (ssize_t) MaxMap; j++)
2081  {
2082  double
2083  count;
2084 
2085  count=area*histogram[GetPixelChannels(image)*j+i];
2086  channel_statistics[channel].entropy+=-count*MagickLog10(count)*
2087  PerceptibleReciprocal(MagickLog10(number_bins));
2088  channel_statistics[CompositePixelChannel].entropy+=-count*
2089  MagickLog10(count)*PerceptibleReciprocal(MagickLog10(number_bins))/
2090  GetPixelChannels(image);
2091  }
2092  }
2093  histogram=(double *) RelinquishMagickMemory(histogram);
2094  for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
2095  {
2096  /*
2097  Compute kurtosis & skewness statistics.
2098  */
2099  standard_deviation=PerceptibleReciprocal(
2100  channel_statistics[i].standard_deviation);
2101  channel_statistics[i].skewness=(channel_statistics[i].sum_cubed-3.0*
2102  channel_statistics[i].mean*channel_statistics[i].sum_squared+2.0*
2103  channel_statistics[i].mean*channel_statistics[i].mean*
2104  channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2105  standard_deviation);
2106  channel_statistics[i].kurtosis=(channel_statistics[i].sum_fourth_power-4.0*
2107  channel_statistics[i].mean*channel_statistics[i].sum_cubed+6.0*
2108  channel_statistics[i].mean*channel_statistics[i].mean*
2109  channel_statistics[i].sum_squared-3.0*channel_statistics[i].mean*
2110  channel_statistics[i].mean*1.0*channel_statistics[i].mean*
2111  channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2112  standard_deviation*standard_deviation)-3.0;
2113  }
2114  channel_statistics[CompositePixelChannel].mean=0.0;
2115  channel_statistics[CompositePixelChannel].standard_deviation=0.0;
2116  channel_statistics[CompositePixelChannel].entropy=0.0;
2117  for (i=0; i < (ssize_t) MaxPixelChannels; i++)
2118  {
2119  channel_statistics[CompositePixelChannel].mean+=
2120  channel_statistics[i].mean;
2121  channel_statistics[CompositePixelChannel].standard_deviation+=
2122  channel_statistics[i].standard_deviation;
2123  channel_statistics[CompositePixelChannel].entropy+=
2124  channel_statistics[i].entropy;
2125  }
2126  channel_statistics[CompositePixelChannel].mean/=(double)
2127  GetImageChannels(image);
2128  channel_statistics[CompositePixelChannel].standard_deviation/=(double)
2129  GetImageChannels(image);
2130  channel_statistics[CompositePixelChannel].entropy/=(double)
2131  GetImageChannels(image);
2132  if (y < (ssize_t) image->rows)
2133  channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
2134  channel_statistics);
2135  return(channel_statistics);
2136 }
2137 
2138 /*
2139 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2140 % %
2141 % %
2142 % %
2143 % P o l y n o m i a l I m a g e %
2144 % %
2145 % %
2146 % %
2147 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2148 %
2149 % PolynomialImage() returns a new image where each pixel is the sum of the
2150 % pixels in the image sequence after applying its corresponding terms
2151 % (coefficient and degree pairs).
2152 %
2153 % The format of the PolynomialImage method is:
2154 %
2155 % Image *PolynomialImage(const Image *images,const size_t number_terms,
2156 % const double *terms,ExceptionInfo *exception)
2157 %
2158 % A description of each parameter follows:
2159 %
2160 % o images: the image sequence.
2161 %
2162 % o number_terms: the number of terms in the list. The actual list length
2163 % is 2 x number_terms + 1 (the constant).
2164 %
2165 % o terms: the list of polynomial coefficients and degree pairs and a
2166 % constant.
2167 %
2168 % o exception: return any errors or warnings in this structure.
2169 %
2170 */
2171 
2173  const size_t number_terms,const double *terms,ExceptionInfo *exception)
2174 {
2175 #define PolynomialImageTag "Polynomial/Image"
2176 
2177  CacheView
2178  *polynomial_view;
2179 
2180  Image
2181  *image;
2182 
2184  status;
2185 
2187  progress;
2188 
2190  **magick_restrict polynomial_pixels;
2191 
2192  size_t
2193  number_images;
2194 
2195  ssize_t
2196  y;
2197 
2198  assert(images != (Image *) NULL);
2199  assert(images->signature == MagickCoreSignature);
2200  if (images->debug != MagickFalse)
2201  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
2202  assert(exception != (ExceptionInfo *) NULL);
2203  assert(exception->signature == MagickCoreSignature);
2204  image=AcquireImageCanvas(images,exception);
2205  if (image == (Image *) NULL)
2206  return((Image *) NULL);
2207  if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
2208  {
2209  image=DestroyImage(image);
2210  return((Image *) NULL);
2211  }
2212  number_images=GetImageListLength(images);
2213  polynomial_pixels=AcquirePixelThreadSet(images);
2214  if (polynomial_pixels == (PixelChannels **) NULL)
2215  {
2216  image=DestroyImage(image);
2217  (void) ThrowMagickException(exception,GetMagickModule(),
2218  ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
2219  return((Image *) NULL);
2220  }
2221  /*
2222  Polynomial image pixels.
2223  */
2224  status=MagickTrue;
2225  progress=0;
2226  polynomial_view=AcquireAuthenticCacheView(image,exception);
2227 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2228  #pragma omp parallel for schedule(static) shared(progress,status) \
2229  magick_number_threads(image,image,image->rows,1)
2230 #endif
2231  for (y=0; y < (ssize_t) image->rows; y++)
2232  {
2233  CacheView
2234  *image_view;
2235 
2236  const Image
2237  *next;
2238 
2239  const int
2240  id = GetOpenMPThreadId();
2241 
2242  register ssize_t
2243  i,
2244  x;
2245 
2246  register PixelChannels
2247  *polynomial_pixel;
2248 
2249  register Quantum
2250  *magick_restrict q;
2251 
2252  ssize_t
2253  j;
2254 
2255  if (status == MagickFalse)
2256  continue;
2257  q=QueueCacheViewAuthenticPixels(polynomial_view,0,y,image->columns,1,
2258  exception);
2259  if (q == (Quantum *) NULL)
2260  {
2261  status=MagickFalse;
2262  continue;
2263  }
2264  polynomial_pixel=polynomial_pixels[id];
2265  for (j=0; j < (ssize_t) image->columns; j++)
2266  for (i=0; i < MaxPixelChannels; i++)
2267  polynomial_pixel[j].channel[i]=0.0;
2268  next=images;
2269  for (j=0; j < (ssize_t) number_images; j++)
2270  {
2271  register const Quantum
2272  *p;
2273 
2274  if (j >= (ssize_t) number_terms)
2275  continue;
2276  image_view=AcquireVirtualCacheView(next,exception);
2277  p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2278  if (p == (const Quantum *) NULL)
2279  {
2280  image_view=DestroyCacheView(image_view);
2281  break;
2282  }
2283  for (x=0; x < (ssize_t) image->columns; x++)
2284  {
2285  register ssize_t
2286  i;
2287 
2288  for (i=0; i < (ssize_t) GetPixelChannels(next); i++)
2289  {
2291  coefficient,
2292  degree;
2293 
2294  PixelChannel channel = GetPixelChannelChannel(image,i);
2295  PixelTrait traits = GetPixelChannelTraits(next,channel);
2296  PixelTrait polynomial_traits=GetPixelChannelTraits(image,channel);
2297  if ((traits == UndefinedPixelTrait) ||
2298  (polynomial_traits == UndefinedPixelTrait))
2299  continue;
2300  if ((traits & UpdatePixelTrait) == 0)
2301  continue;
2302  coefficient=(MagickRealType) terms[2*j];
2303  degree=(MagickRealType) terms[(j << 1)+1];
2304  polynomial_pixel[x].channel[i]+=coefficient*
2305  pow(QuantumScale*GetPixelChannel(image,channel,p),degree);
2306  }
2307  p+=GetPixelChannels(next);
2308  }
2309  image_view=DestroyCacheView(image_view);
2310  next=GetNextImageInList(next);
2311  }
2312  for (x=0; x < (ssize_t) image->columns; x++)
2313  {
2314  register ssize_t
2315  i;
2316 
2317  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2318  {
2319  PixelChannel channel = GetPixelChannelChannel(image,i);
2320  PixelTrait traits = GetPixelChannelTraits(image,channel);
2321  if (traits == UndefinedPixelTrait)
2322  continue;
2323  if ((traits & UpdatePixelTrait) == 0)
2324  continue;
2325  q[i]=ClampToQuantum(QuantumRange*polynomial_pixel[x].channel[i]);
2326  }
2327  q+=GetPixelChannels(image);
2328  }
2329  if (SyncCacheViewAuthenticPixels(polynomial_view,exception) == MagickFalse)
2330  status=MagickFalse;
2331  if (images->progress_monitor != (MagickProgressMonitor) NULL)
2332  {
2334  proceed;
2335 
2336 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2337  #pragma omp atomic
2338 #endif
2339  progress++;
2340  proceed=SetImageProgress(images,PolynomialImageTag,progress,
2341  image->rows);
2342  if (proceed == MagickFalse)
2343  status=MagickFalse;
2344  }
2345  }
2346  polynomial_view=DestroyCacheView(polynomial_view);
2347  polynomial_pixels=DestroyPixelThreadSet(polynomial_pixels);
2348  if (status == MagickFalse)
2349  image=DestroyImage(image);
2350  return(image);
2351 }
2352 
2353 /*
2354 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2355 % %
2356 % %
2357 % %
2358 % S t a t i s t i c I m a g e %
2359 % %
2360 % %
2361 % %
2362 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2363 %
2364 % StatisticImage() makes each pixel the min / max / median / mode / etc. of
2365 % the neighborhood of the specified width and height.
2366 %
2367 % The format of the StatisticImage method is:
2368 %
2369 % Image *StatisticImage(const Image *image,const StatisticType type,
2370 % const size_t width,const size_t height,ExceptionInfo *exception)
2371 %
2372 % A description of each parameter follows:
2373 %
2374 % o image: the image.
2375 %
2376 % o type: the statistic type (median, mode, etc.).
2377 %
2378 % o width: the width of the pixel neighborhood.
2379 %
2380 % o height: the height of the pixel neighborhood.
2381 %
2382 % o exception: return any errors or warnings in this structure.
2383 %
2384 */
2385 
2386 typedef struct _SkipNode
2387 {
2388  size_t
2389  next[9],
2390  count,
2391  signature;
2392 } SkipNode;
2393 
2394 typedef struct _SkipList
2395 {
2396  ssize_t
2398 
2399  SkipNode
2401 } SkipList;
2402 
2403 typedef struct _PixelList
2404 {
2405  size_t
2407  seed;
2408 
2409  SkipList
2411 
2412  size_t
2414 } PixelList;
2415 
2417 {
2418  if (pixel_list == (PixelList *) NULL)
2419  return((PixelList *) NULL);
2420  if (pixel_list->skip_list.nodes != (SkipNode *) NULL)
2422  pixel_list->skip_list.nodes);
2423  pixel_list=(PixelList *) RelinquishMagickMemory(pixel_list);
2424  return(pixel_list);
2425 }
2426 
2428 {
2429  register ssize_t
2430  i;
2431 
2432  assert(pixel_list != (PixelList **) NULL);
2433  for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
2434  if (pixel_list[i] != (PixelList *) NULL)
2435  pixel_list[i]=DestroyPixelList(pixel_list[i]);
2436  pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list);
2437  return(pixel_list);
2438 }
2439 
2440 static PixelList *AcquirePixelList(const size_t width,const size_t height)
2441 {
2442  PixelList
2443  *pixel_list;
2444 
2445  pixel_list=(PixelList *) AcquireMagickMemory(sizeof(*pixel_list));
2446  if (pixel_list == (PixelList *) NULL)
2447  return(pixel_list);
2448  (void) memset((void *) pixel_list,0,sizeof(*pixel_list));
2449  pixel_list->length=width*height;
2450  pixel_list->skip_list.nodes=(SkipNode *) AcquireAlignedMemory(65537UL,
2451  sizeof(*pixel_list->skip_list.nodes));
2452  if (pixel_list->skip_list.nodes == (SkipNode *) NULL)
2453  return(DestroyPixelList(pixel_list));
2454  (void) memset(pixel_list->skip_list.nodes,0,65537UL*
2455  sizeof(*pixel_list->skip_list.nodes));
2456  pixel_list->signature=MagickCoreSignature;
2457  return(pixel_list);
2458 }
2459 
2460 static PixelList **AcquirePixelListThreadSet(const size_t width,
2461  const size_t height)
2462 {
2463  PixelList
2464  **pixel_list;
2465 
2466  register ssize_t
2467  i;
2468 
2469  size_t
2470  number_threads;
2471 
2472  number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
2473  pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
2474  sizeof(*pixel_list));
2475  if (pixel_list == (PixelList **) NULL)
2476  return((PixelList **) NULL);
2477  (void) memset(pixel_list,0,number_threads*sizeof(*pixel_list));
2478  for (i=0; i < (ssize_t) number_threads; i++)
2479  {
2480  pixel_list[i]=AcquirePixelList(width,height);
2481  if (pixel_list[i] == (PixelList *) NULL)
2482  return(DestroyPixelListThreadSet(pixel_list));
2483  }
2484  return(pixel_list);
2485 }
2486 
2487 static void AddNodePixelList(PixelList *pixel_list,const size_t color)
2488 {
2489  register SkipList
2490  *p;
2491 
2492  register ssize_t
2493  level;
2494 
2495  size_t
2496  search,
2497  update[9];
2498 
2499  /*
2500  Initialize the node.
2501  */
2502  p=(&pixel_list->skip_list);
2503  p->nodes[color].signature=pixel_list->signature;
2504  p->nodes[color].count=1;
2505  /*
2506  Determine where it belongs in the list.
2507  */
2508  search=65536UL;
2509  for (level=p->level; level >= 0; level--)
2510  {
2511  while (p->nodes[search].next[level] < color)
2512  search=p->nodes[search].next[level];
2513  update[level]=search;
2514  }
2515  /*
2516  Generate a pseudo-random level for this node.
2517  */
2518  for (level=0; ; level++)
2519  {
2520  pixel_list->seed=(pixel_list->seed*42893621L)+1L;
2521  if ((pixel_list->seed & 0x300) != 0x300)
2522  break;
2523  }
2524  if (level > 8)
2525  level=8;
2526  if (level > (p->level+2))
2527  level=p->level+2;
2528  /*
2529  If we're raising the list's level, link back to the root node.
2530  */
2531  while (level > p->level)
2532  {
2533  p->level++;
2534  update[p->level]=65536UL;
2535  }
2536  /*
2537  Link the node into the skip-list.
2538  */
2539  do
2540  {
2541  p->nodes[color].next[level]=p->nodes[update[level]].next[level];
2542  p->nodes[update[level]].next[level]=color;
2543  } while (level-- > 0);
2544 }
2545 
2546 static inline void GetMaximumPixelList(PixelList *pixel_list,Quantum *pixel)
2547 {
2548  register SkipList
2549  *p;
2550 
2551  size_t
2552  color,
2553  maximum;
2554 
2555  ssize_t
2556  count;
2557 
2558  /*
2559  Find the maximum value for each of the color.
2560  */
2561  p=(&pixel_list->skip_list);
2562  color=65536L;
2563  count=0;
2564  maximum=p->nodes[color].next[0];
2565  do
2566  {
2567  color=p->nodes[color].next[0];
2568  if (color > maximum)
2569  maximum=color;
2570  count+=p->nodes[color].count;
2571  } while (count < (ssize_t) pixel_list->length);
2572  *pixel=ScaleShortToQuantum((unsigned short) maximum);
2573 }
2574 
2575 static inline void GetMeanPixelList(PixelList *pixel_list,Quantum *pixel)
2576 {
2577  double
2578  sum;
2579 
2580  register SkipList
2581  *p;
2582 
2583  size_t
2584  color;
2585 
2586  ssize_t
2587  count;
2588 
2589  /*
2590  Find the mean value for each of the color.
2591  */
2592  p=(&pixel_list->skip_list);
2593  color=65536L;
2594  count=0;
2595  sum=0.0;
2596  do
2597  {
2598  color=p->nodes[color].next[0];
2599  sum+=(double) p->nodes[color].count*color;
2600  count+=p->nodes[color].count;
2601  } while (count < (ssize_t) pixel_list->length);
2602  sum/=pixel_list->length;
2603  *pixel=ScaleShortToQuantum((unsigned short) sum);
2604 }
2605 
2606 static inline void GetMedianPixelList(PixelList *pixel_list,Quantum *pixel)
2607 {
2608  register SkipList
2609  *p;
2610 
2611  size_t
2612  color;
2613 
2614  ssize_t
2615  count;
2616 
2617  /*
2618  Find the median value for each of the color.
2619  */
2620  p=(&pixel_list->skip_list);
2621  color=65536L;
2622  count=0;
2623  do
2624  {
2625  color=p->nodes[color].next[0];
2626  count+=p->nodes[color].count;
2627  } while (count <= (ssize_t) (pixel_list->length >> 1));
2628  *pixel=ScaleShortToQuantum((unsigned short) color);
2629 }
2630 
2631 static inline void GetMinimumPixelList(PixelList *pixel_list,Quantum *pixel)
2632 {
2633  register SkipList
2634  *p;
2635 
2636  size_t
2637  color,
2638  minimum;
2639 
2640  ssize_t
2641  count;
2642 
2643  /*
2644  Find the minimum value for each of the color.
2645  */
2646  p=(&pixel_list->skip_list);
2647  count=0;
2648  color=65536UL;
2649  minimum=p->nodes[color].next[0];
2650  do
2651  {
2652  color=p->nodes[color].next[0];
2653  if (color < minimum)
2654  minimum=color;
2655  count+=p->nodes[color].count;
2656  } while (count < (ssize_t) pixel_list->length);
2657  *pixel=ScaleShortToQuantum((unsigned short) minimum);
2658 }
2659 
2660 static inline void GetModePixelList(PixelList *pixel_list,Quantum *pixel)
2661 {
2662  register SkipList
2663  *p;
2664 
2665  size_t
2666  color,
2667  max_count,
2668  mode;
2669 
2670  ssize_t
2671  count;
2672 
2673  /*
2674  Make each pixel the 'predominant color' of the specified neighborhood.
2675  */
2676  p=(&pixel_list->skip_list);
2677  color=65536L;
2678  mode=color;
2679  max_count=p->nodes[mode].count;
2680  count=0;
2681  do
2682  {
2683  color=p->nodes[color].next[0];
2684  if (p->nodes[color].count > max_count)
2685  {
2686  mode=color;
2687  max_count=p->nodes[mode].count;
2688  }
2689  count+=p->nodes[color].count;
2690  } while (count < (ssize_t) pixel_list->length);
2691  *pixel=ScaleShortToQuantum((unsigned short) mode);
2692 }
2693 
2694 static inline void GetNonpeakPixelList(PixelList *pixel_list,Quantum *pixel)
2695 {
2696  register SkipList
2697  *p;
2698 
2699  size_t
2700  color,
2701  next,
2702  previous;
2703 
2704  ssize_t
2705  count;
2706 
2707  /*
2708  Finds the non peak value for each of the colors.
2709  */
2710  p=(&pixel_list->skip_list);
2711  color=65536L;
2712  next=p->nodes[color].next[0];
2713  count=0;
2714  do
2715  {
2716  previous=color;
2717  color=next;
2718  next=p->nodes[color].next[0];
2719  count+=p->nodes[color].count;
2720  } while (count <= (ssize_t) (pixel_list->length >> 1));
2721  if ((previous == 65536UL) && (next != 65536UL))
2722  color=next;
2723  else
2724  if ((previous != 65536UL) && (next == 65536UL))
2725  color=previous;
2726  *pixel=ScaleShortToQuantum((unsigned short) color);
2727 }
2728 
2729 static inline void GetRootMeanSquarePixelList(PixelList *pixel_list,
2730  Quantum *pixel)
2731 {
2732  double
2733  sum;
2734 
2735  register SkipList
2736  *p;
2737 
2738  size_t
2739  color;
2740 
2741  ssize_t
2742  count;
2743 
2744  /*
2745  Find the root mean square value for each of the color.
2746  */
2747  p=(&pixel_list->skip_list);
2748  color=65536L;
2749  count=0;
2750  sum=0.0;
2751  do
2752  {
2753  color=p->nodes[color].next[0];
2754  sum+=(double) (p->nodes[color].count*color*color);
2755  count+=p->nodes[color].count;
2756  } while (count < (ssize_t) pixel_list->length);
2757  sum/=pixel_list->length;
2758  *pixel=ScaleShortToQuantum((unsigned short) sqrt(sum));
2759 }
2760 
2761 static inline void GetStandardDeviationPixelList(PixelList *pixel_list,
2762  Quantum *pixel)
2763 {
2764  double
2765  sum,
2766  sum_squared;
2767 
2768  register SkipList
2769  *p;
2770 
2771  size_t
2772  color;
2773 
2774  ssize_t
2775  count;
2776 
2777  /*
2778  Find the standard-deviation value for each of the color.
2779  */
2780  p=(&pixel_list->skip_list);
2781  color=65536L;
2782  count=0;
2783  sum=0.0;
2784  sum_squared=0.0;
2785  do
2786  {
2787  register ssize_t
2788  i;
2789 
2790  color=p->nodes[color].next[0];
2791  sum+=(double) p->nodes[color].count*color;
2792  for (i=0; i < (ssize_t) p->nodes[color].count; i++)
2793  sum_squared+=((double) color)*((double) color);
2794  count+=p->nodes[color].count;
2795  } while (count < (ssize_t) pixel_list->length);
2796  sum/=pixel_list->length;
2797  sum_squared/=pixel_list->length;
2798  *pixel=ScaleShortToQuantum((unsigned short) sqrt(sum_squared-(sum*sum)));
2799 }
2800 
2801 static inline void InsertPixelList(const Quantum pixel,PixelList *pixel_list)
2802 {
2803  size_t
2804  signature;
2805 
2806  unsigned short
2807  index;
2808 
2809  index=ScaleQuantumToShort(pixel);
2810  signature=pixel_list->skip_list.nodes[index].signature;
2811  if (signature == pixel_list->signature)
2812  {
2813  pixel_list->skip_list.nodes[index].count++;
2814  return;
2815  }
2816  AddNodePixelList(pixel_list,index);
2817 }
2818 
2819 static void ResetPixelList(PixelList *pixel_list)
2820 {
2821  int
2822  level;
2823 
2824  register SkipNode
2825  *root;
2826 
2827  register SkipList
2828  *p;
2829 
2830  /*
2831  Reset the skip-list.
2832  */
2833  p=(&pixel_list->skip_list);
2834  root=p->nodes+65536UL;
2835  p->level=0;
2836  for (level=0; level < 9; level++)
2837  root->next[level]=65536UL;
2838  pixel_list->seed=pixel_list->signature++;
2839 }
2840 
2842  const size_t width,const size_t height,ExceptionInfo *exception)
2843 {
2844 #define StatisticImageTag "Statistic/Image"
2845 
2846  CacheView
2847  *image_view,
2848  *statistic_view;
2849 
2850  Image
2851  *statistic_image;
2852 
2854  status;
2855 
2857  progress;
2858 
2859  PixelList
2860  **magick_restrict pixel_list;
2861 
2862  ssize_t
2863  center,
2864  y;
2865 
2866  /*
2867  Initialize statistics image attributes.
2868  */
2869  assert(image != (Image *) NULL);
2870  assert(image->signature == MagickCoreSignature);
2871  if (image->debug != MagickFalse)
2872  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2873  assert(exception != (ExceptionInfo *) NULL);
2874  assert(exception->signature == MagickCoreSignature);
2875  statistic_image=CloneImage(image,0,0,MagickTrue,
2876  exception);
2877  if (statistic_image == (Image *) NULL)
2878  return((Image *) NULL);
2879  status=SetImageStorageClass(statistic_image,DirectClass,exception);
2880  if (status == MagickFalse)
2881  {
2882  statistic_image=DestroyImage(statistic_image);
2883  return((Image *) NULL);
2884  }
2885  pixel_list=AcquirePixelListThreadSet(MagickMax(width,1),MagickMax(height,1));
2886  if (pixel_list == (PixelList **) NULL)
2887  {
2888  statistic_image=DestroyImage(statistic_image);
2889  ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2890  }
2891  /*
2892  Make each pixel the min / max / median / mode / etc. of the neighborhood.
2893  */
2894  center=(ssize_t) GetPixelChannels(image)*(image->columns+MagickMax(width,1))*
2895  (MagickMax(height,1)/2L)+GetPixelChannels(image)*(MagickMax(width,1)/2L);
2896  status=MagickTrue;
2897  progress=0;
2898  image_view=AcquireVirtualCacheView(image,exception);
2899  statistic_view=AcquireAuthenticCacheView(statistic_image,exception);
2900 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2901  #pragma omp parallel for schedule(static) shared(progress,status) \
2902  magick_number_threads(image,statistic_image,statistic_image->rows,1)
2903 #endif
2904  for (y=0; y < (ssize_t) statistic_image->rows; y++)
2905  {
2906  const int
2907  id = GetOpenMPThreadId();
2908 
2909  register const Quantum
2910  *magick_restrict p;
2911 
2912  register Quantum
2913  *magick_restrict q;
2914 
2915  register ssize_t
2916  x;
2917 
2918  if (status == MagickFalse)
2919  continue;
2920  p=GetCacheViewVirtualPixels(image_view,-((ssize_t) MagickMax(width,1)/2L),y-
2921  (ssize_t) (MagickMax(height,1)/2L),image->columns+MagickMax(width,1),
2922  MagickMax(height,1),exception);
2923  q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns, 1,exception);
2924  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2925  {
2926  status=MagickFalse;
2927  continue;
2928  }
2929  for (x=0; x < (ssize_t) statistic_image->columns; x++)
2930  {
2931  register ssize_t
2932  i;
2933 
2934  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2935  {
2936  Quantum
2937  pixel;
2938 
2939  register const Quantum
2940  *magick_restrict pixels;
2941 
2942  register ssize_t
2943  u;
2944 
2945  ssize_t
2946  v;
2947 
2948  PixelChannel channel = GetPixelChannelChannel(image,i);
2949  PixelTrait traits = GetPixelChannelTraits(image,channel);
2950  PixelTrait statistic_traits=GetPixelChannelTraits(statistic_image,
2951  channel);
2952  if ((traits == UndefinedPixelTrait) ||
2953  (statistic_traits == UndefinedPixelTrait))
2954  continue;
2955  if (((statistic_traits & CopyPixelTrait) != 0) ||
2956  (GetPixelWriteMask(image,p) <= (QuantumRange/2)))
2957  {
2958  SetPixelChannel(statistic_image,channel,p[center+i],q);
2959  continue;
2960  }
2961  if ((statistic_traits & UpdatePixelTrait) == 0)
2962  continue;
2963  pixels=p;
2964  ResetPixelList(pixel_list[id]);
2965  for (v=0; v < (ssize_t) MagickMax(height,1); v++)
2966  {
2967  for (u=0; u < (ssize_t) MagickMax(width,1); u++)
2968  {
2969  InsertPixelList(pixels[i],pixel_list[id]);
2970  pixels+=GetPixelChannels(image);
2971  }
2972  pixels+=GetPixelChannels(image)*image->columns;
2973  }
2974  switch (type)
2975  {
2976  case GradientStatistic:
2977  {
2978  double
2979  maximum,
2980  minimum;
2981 
2982  GetMinimumPixelList(pixel_list[id],&pixel);
2983  minimum=(double) pixel;
2984  GetMaximumPixelList(pixel_list[id],&pixel);
2985  maximum=(double) pixel;
2986  pixel=ClampToQuantum(MagickAbsoluteValue(maximum-minimum));
2987  break;
2988  }
2989  case MaximumStatistic:
2990  {
2991  GetMaximumPixelList(pixel_list[id],&pixel);
2992  break;
2993  }
2994  case MeanStatistic:
2995  {
2996  GetMeanPixelList(pixel_list[id],&pixel);
2997  break;
2998  }
2999  case MedianStatistic:
3000  default:
3001  {
3002  GetMedianPixelList(pixel_list[id],&pixel);
3003  break;
3004  }
3005  case MinimumStatistic:
3006  {
3007  GetMinimumPixelList(pixel_list[id],&pixel);
3008  break;
3009  }
3010  case ModeStatistic:
3011  {
3012  GetModePixelList(pixel_list[id],&pixel);
3013  break;
3014  }
3015  case NonpeakStatistic:
3016  {
3017  GetNonpeakPixelList(pixel_list[id],&pixel);
3018  break;
3019  }
3021  {
3022  GetRootMeanSquarePixelList(pixel_list[id],&pixel);
3023  break;
3024  }
3026  {
3027  GetStandardDeviationPixelList(pixel_list[id],&pixel);
3028  break;
3029  }
3030  }
3031  SetPixelChannel(statistic_image,channel,pixel,q);
3032  }
3033  p+=GetPixelChannels(image);
3034  q+=GetPixelChannels(statistic_image);
3035  }
3036  if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
3037  status=MagickFalse;
3038  if (image->progress_monitor != (MagickProgressMonitor) NULL)
3039  {
3041  proceed;
3042 
3043 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3044  #pragma omp atomic
3045 #endif
3046  progress++;
3047  proceed=SetImageProgress(image,StatisticImageTag,progress,image->rows);
3048  if (proceed == MagickFalse)
3049  status=MagickFalse;
3050  }
3051  }
3052  statistic_view=DestroyCacheView(statistic_view);
3053  image_view=DestroyCacheView(image_view);
3054  pixel_list=DestroyPixelListThreadSet(pixel_list);
3055  if (status == MagickFalse)
3056  statistic_image=DestroyImage(statistic_image);
3057  return(statistic_image);
3058 }
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:2416
static void GetMaximumPixelList(PixelList *pixel_list, Quantum *pixel)
Definition: statistic.c:2546
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
static Quantum GetPixelReadMask(const Image *magick_restrict image, const Quantum *magick_restrict pixel)
#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:2389
#define MaximumNumberOfImageMoments
Definition: statistic.h:25
static void GetMeanPixelList(PixelList *pixel_list, Quantum *pixel)
Definition: statistic.c:2575
#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:2406
#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:2400
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:2694
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:2389
size_t signature
Definition: statistic.c:2413
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:3229
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:2427
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:1145
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:790
static void GetMinimumPixelList(PixelList *pixel_list, Quantum *pixel)
Definition: statistic.c:2631
static void InsertPixelList(const Quantum pixel, PixelList *pixel_list)
Definition: statistic.c:2801
struct _Image * next
Definition: image.h:348
static PixelList * AcquirePixelList(const size_t width, const size_t height)
Definition: statistic.c:2440
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:2616
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:2406
ssize_t level
Definition: statistic.c:2397
#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:2606
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:2660
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:2410
static void GetStandardDeviationPixelList(PixelList *pixel_list, Quantum *pixel)
Definition: statistic.c:2761
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:113
size_t next[9]
Definition: statistic.c:2389
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:2487
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:2819
MagickExport Image * PolynomialImage(const Image *images, const size_t number_terms, const double *terms, ExceptionInfo *exception)
Definition: statistic.c:2172
static void GetRootMeanSquarePixelList(PixelList *pixel_list, Quantum *pixel)
Definition: statistic.c:2729
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:2460
MagickExport Image * DestroyImage(Image *image)
Definition: image.c:1180
MagickExport Image * CloneImage(const Image *image, const size_t columns, const size_t rows, const MagickBooleanType detach, ExceptionInfo *exception)
Definition: image.c:795
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:2841
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