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