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