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