MagickCore  6.9.12-62
Convert, Edit, Or Compose Bitmap Images
 All Data Structures
statistic.c
1 /*
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % %
4 % %
5 % %
6 % SSSSS TTTTT AAA TTTTT IIIII SSSSS TTTTT IIIII CCCC %
7 % SS T A A T I SS T I C %
8 % SSS T AAAAA T I SSS T I C %
9 % SS T A A T I SS T I C %
10 % SSSSS T A A T IIIII SSSSS T IIIII CCCC %
11 % %
12 % %
13 % MagickCore Image Statistical Methods %
14 % %
15 % Software Design %
16 % Cristy %
17 % July 1992 %
18 % %
19 % %
20 % Copyright 1999-2021 ImageMagick Studio LLC, a non-profit organization %
21 % dedicated to making software imaging solutions freely available. %
22 % %
23 % You may not use this file except in compliance with the License. You may %
24 % obtain a copy of the License at %
25 % %
26 % https://imagemagick.org/script/license.php %
27 % %
28 % Unless required by applicable law or agreed to in writing, software %
29 % distributed under the License is distributed on an "AS IS" BASIS, %
30 % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. %
31 % See the License for the specific language governing permissions and %
32 % limitations under the License. %
33 % %
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35 %
36 %
37 %
38 */
39 
40 
41 /*
42  Include declarations.
43 */
44 #include "magick/studio.h"
45 #include "magick/accelerate-private.h"
46 #include "magick/animate.h"
47 #include "magick/animate.h"
48 #include "magick/blob.h"
49 #include "magick/blob-private.h"
50 #include "magick/cache.h"
51 #include "magick/cache-private.h"
52 #include "magick/cache-view.h"
53 #include "magick/client.h"
54 #include "magick/color.h"
55 #include "magick/color-private.h"
56 #include "magick/colorspace.h"
57 #include "magick/colorspace-private.h"
58 #include "magick/composite.h"
59 #include "magick/composite-private.h"
60 #include "magick/compress.h"
61 #include "magick/constitute.h"
62 #include "magick/deprecate.h"
63 #include "magick/display.h"
64 #include "magick/draw.h"
65 #include "magick/enhance.h"
66 #include "magick/exception.h"
67 #include "magick/exception-private.h"
68 #include "magick/gem.h"
69 #include "magick/geometry.h"
70 #include "magick/list.h"
71 #include "magick/image-private.h"
72 #include "magick/magic.h"
73 #include "magick/magick.h"
74 #include "magick/memory_.h"
75 #include "magick/module.h"
76 #include "magick/monitor.h"
77 #include "magick/monitor-private.h"
78 #include "magick/option.h"
79 #include "magick/paint.h"
80 #include "magick/pixel-private.h"
81 #include "magick/profile.h"
82 #include "magick/property.h"
83 #include "magick/quantize.h"
84 #include "magick/random_.h"
85 #include "magick/random-private.h"
86 #include "magick/resource_.h"
87 #include "magick/segment.h"
88 #include "magick/semaphore.h"
89 #include "magick/signature-private.h"
90 #include "magick/statistic.h"
91 #include "magick/string_.h"
92 #include "magick/thread-private.h"
93 #include "magick/timer.h"
94 #include "magick/utility.h"
95 #include "magick/version.h"
96 
97 /*
98 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
99 % %
100 % %
101 % %
102 % E v a l u a t e I m a g e %
103 % %
104 % %
105 % %
106 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
107 %
108 % EvaluateImage() applies a value to the image with an arithmetic, relational,
109 % or logical operator to an image. Use these operations to lighten or darken
110 % an image, to increase or decrease contrast in an image, or to produce the
111 % "negative" of an image.
112 %
113 % The format of the EvaluateImageChannel method is:
114 %
115 % MagickBooleanType EvaluateImage(Image *image,
116 % const MagickEvaluateOperator op,const double value,
117 % ExceptionInfo *exception)
118 % MagickBooleanType EvaluateImages(Image *images,
119 % const MagickEvaluateOperator op,const double value,
120 % ExceptionInfo *exception)
121 % MagickBooleanType EvaluateImageChannel(Image *image,
122 % const ChannelType channel,const MagickEvaluateOperator op,
123 % const double value,ExceptionInfo *exception)
124 %
125 % A description of each parameter follows:
126 %
127 % o image: the image.
128 %
129 % o channel: the channel.
130 %
131 % o op: A channel op.
132 %
133 % o value: A value value.
134 %
135 % o exception: return any errors or warnings in this structure.
136 %
137 */
138 
139 static MagickPixelPacket **DestroyPixelTLS(const Image *images,
140  MagickPixelPacket **pixels)
141 {
142  ssize_t
143  i;
144 
145  size_t
146  rows;
147 
148  assert(pixels != (MagickPixelPacket **) NULL);
149  rows=MagickMax(GetImageListLength(images),
150  (size_t) GetMagickResourceLimit(ThreadResource));
151  for (i=0; i < (ssize_t) rows; i++)
152  if (pixels[i] != (MagickPixelPacket *) NULL)
153  pixels[i]=(MagickPixelPacket *) RelinquishMagickMemory(pixels[i]);
154  pixels=(MagickPixelPacket **) RelinquishMagickMemory(pixels);
155  return(pixels);
156 }
157 
158 static MagickPixelPacket **AcquirePixelTLS(const Image *images)
159 {
160  const Image
161  *next;
162 
164  **pixels;
165 
166  ssize_t
167  i,
168  j;
169 
170  size_t
171  columns,
172  rows;
173 
174  rows=MagickMax(GetImageListLength(images),
175  (size_t) GetMagickResourceLimit(ThreadResource));
176  pixels=(MagickPixelPacket **) AcquireQuantumMemory(rows,sizeof(*pixels));
177  if (pixels == (MagickPixelPacket **) NULL)
178  return((MagickPixelPacket **) NULL);
179  (void) memset(pixels,0,rows*sizeof(*pixels));
180  columns=GetImageListLength(images);
181  for (next=images; next != (Image *) NULL; next=next->next)
182  columns=MagickMax(next->columns,columns);
183  for (i=0; i < (ssize_t) rows; i++)
184  {
185  pixels[i]=(MagickPixelPacket *) AcquireQuantumMemory(columns,
186  sizeof(**pixels));
187  if (pixels[i] == (MagickPixelPacket *) NULL)
188  return(DestroyPixelTLS(images,pixels));
189  for (j=0; j < (ssize_t) columns; j++)
190  GetMagickPixelPacket(images,&pixels[i][j]);
191  }
192  return(pixels);
193 }
194 
195 static inline double EvaluateMax(const double x,const double y)
196 {
197  if (x > y)
198  return(x);
199  return(y);
200 }
201 
202 #if defined(__cplusplus) || defined(c_plusplus)
203 extern "C" {
204 #endif
205 
206 static int IntensityCompare(const void *x,const void *y)
207 {
208  const MagickPixelPacket
209  *color_1,
210  *color_2;
211 
212  int
213  intensity;
214 
215  color_1=(const MagickPixelPacket *) x;
216  color_2=(const MagickPixelPacket *) y;
217  intensity=(int) MagickPixelIntensity(color_2)-(int)
218  MagickPixelIntensity(color_1);
219  return(intensity);
220 }
221 
222 #if defined(__cplusplus) || defined(c_plusplus)
223 }
224 #endif
225 
226 static MagickRealType ApplyEvaluateOperator(RandomInfo *random_info,
227  const Quantum pixel,const MagickEvaluateOperator op,
228  const MagickRealType value)
229 {
230  MagickRealType
231  result;
232 
233  ssize_t
234  i;
235 
236  result=0.0;
237  switch (op)
238  {
239  case UndefinedEvaluateOperator:
240  break;
241  case AbsEvaluateOperator:
242  {
243  result=(MagickRealType) fabs((double) (pixel+value));
244  break;
245  }
246  case AddEvaluateOperator:
247  {
248  result=(MagickRealType) (pixel+value);
249  break;
250  }
251  case AddModulusEvaluateOperator:
252  {
253  /*
254  This returns a 'floored modulus' of the addition which is a
255  positive result. It differs from % or fmod() which returns a
256  'truncated modulus' result, where floor() is replaced by trunc()
257  and could return a negative result (which is clipped).
258  */
259  result=pixel+value;
260  result-=(QuantumRange+1.0)*floor((double) result/(QuantumRange+1.0));
261  break;
262  }
263  case AndEvaluateOperator:
264  {
265  result=(MagickRealType) ((ssize_t) pixel & (ssize_t) (value+0.5));
266  break;
267  }
268  case CosineEvaluateOperator:
269  {
270  result=(MagickRealType) (QuantumRange*(0.5*cos((double) (2.0*MagickPI*
271  QuantumScale*pixel*value))+0.5));
272  break;
273  }
274  case DivideEvaluateOperator:
275  {
276  result=pixel/(value == 0.0 ? 1.0 : value);
277  break;
278  }
279  case ExponentialEvaluateOperator:
280  {
281  result=(MagickRealType) (QuantumRange*exp((double) (value*QuantumScale*
282  pixel)));
283  break;
284  }
285  case GaussianNoiseEvaluateOperator:
286  {
287  result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
288  GaussianNoise,value);
289  break;
290  }
291  case ImpulseNoiseEvaluateOperator:
292  {
293  result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
294  ImpulseNoise,value);
295  break;
296  }
297  case InverseLogEvaluateOperator:
298  {
299  result=(QuantumRange*pow((value+1.0),QuantumScale*pixel)-1.0)*
300  PerceptibleReciprocal(value);
301  break;
302  }
303  case LaplacianNoiseEvaluateOperator:
304  {
305  result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
306  LaplacianNoise,value);
307  break;
308  }
309  case LeftShiftEvaluateOperator:
310  {
311  result=(double) pixel;
312  for (i=0; i < (ssize_t) value; i++)
313  result*=2.0;
314  break;
315  }
316  case LogEvaluateOperator:
317  {
318  if ((QuantumScale*pixel) >= MagickEpsilon)
319  result=(MagickRealType) (QuantumRange*log((double) (QuantumScale*value*
320  pixel+1.0))/log((double) (value+1.0)));
321  break;
322  }
323  case MaxEvaluateOperator:
324  {
325  result=(MagickRealType) EvaluateMax((double) pixel,value);
326  break;
327  }
328  case MeanEvaluateOperator:
329  {
330  result=(MagickRealType) (pixel+value);
331  break;
332  }
333  case MedianEvaluateOperator:
334  {
335  result=(MagickRealType) (pixel+value);
336  break;
337  }
338  case MinEvaluateOperator:
339  {
340  result=(MagickRealType) MagickMin((double) pixel,value);
341  break;
342  }
343  case MultiplicativeNoiseEvaluateOperator:
344  {
345  result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
346  MultiplicativeGaussianNoise,value);
347  break;
348  }
349  case MultiplyEvaluateOperator:
350  {
351  result=(MagickRealType) (value*pixel);
352  break;
353  }
354  case OrEvaluateOperator:
355  {
356  result=(MagickRealType) ((ssize_t) pixel | (ssize_t) (value+0.5));
357  break;
358  }
359  case PoissonNoiseEvaluateOperator:
360  {
361  result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
362  PoissonNoise,value);
363  break;
364  }
365  case PowEvaluateOperator:
366  {
367  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 
2110 static inline double MagickLog10(const double x)
2111 {
2112 #define Log10Epsilon (1.0e-11)
2113 
2114  if (fabs(x) < Log10Epsilon)
2115  return(log10(Log10Epsilon));
2116  return(log10(fabs(x)));
2117 }
2118 
2119 MagickExport ChannelPerceptualHash *GetImageChannelPerceptualHash(
2120  const Image *image,ExceptionInfo *exception)
2121 {
2123  *moments;
2124 
2126  *perceptual_hash;
2127 
2128  Image
2129  *hash_image;
2130 
2131  MagickBooleanType
2132  status;
2133 
2134  ssize_t
2135  i;
2136 
2137  ssize_t
2138  channel;
2139 
2140  /*
2141  Blur then transform to sRGB colorspace.
2142  */
2143  hash_image=BlurImage(image,0.0,1.0,exception);
2144  if (hash_image == (Image *) NULL)
2145  return((ChannelPerceptualHash *) NULL);
2146  hash_image->depth=8;
2147  status=TransformImageColorspace(hash_image,sRGBColorspace);
2148  if (status == MagickFalse)
2149  return((ChannelPerceptualHash *) NULL);
2150  moments=GetImageChannelMoments(hash_image,exception);
2151  hash_image=DestroyImage(hash_image);
2152  if (moments == (ChannelMoments *) NULL)
2153  return((ChannelPerceptualHash *) NULL);
2154  perceptual_hash=(ChannelPerceptualHash *) AcquireQuantumMemory(
2155  CompositeChannels+1UL,sizeof(*perceptual_hash));
2156  if (perceptual_hash == (ChannelPerceptualHash *) NULL)
2157  return((ChannelPerceptualHash *) NULL);
2158  for (channel=0; channel <= CompositeChannels; channel++)
2159  for (i=0; i < MaximumNumberOfImageMoments; i++)
2160  perceptual_hash[channel].P[i]=(-MagickLog10(moments[channel].I[i]));
2161  moments=(ChannelMoments *) RelinquishMagickMemory(moments);
2162  /*
2163  Blur then transform to HCLp colorspace.
2164  */
2165  hash_image=BlurImage(image,0.0,1.0,exception);
2166  if (hash_image == (Image *) NULL)
2167  {
2168  perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
2169  perceptual_hash);
2170  return((ChannelPerceptualHash *) NULL);
2171  }
2172  hash_image->depth=8;
2173  status=TransformImageColorspace(hash_image,HCLpColorspace);
2174  if (status == MagickFalse)
2175  {
2176  perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
2177  perceptual_hash);
2178  return((ChannelPerceptualHash *) NULL);
2179  }
2180  moments=GetImageChannelMoments(hash_image,exception);
2181  hash_image=DestroyImage(hash_image);
2182  if (moments == (ChannelMoments *) NULL)
2183  {
2184  perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
2185  perceptual_hash);
2186  return((ChannelPerceptualHash *) NULL);
2187  }
2188  for (channel=0; channel <= CompositeChannels; channel++)
2189  for (i=0; i < MaximumNumberOfImageMoments; i++)
2190  perceptual_hash[channel].Q[i]=(-MagickLog10(moments[channel].I[i]));
2191  moments=(ChannelMoments *) RelinquishMagickMemory(moments);
2192  return(perceptual_hash);
2193 }
2194 
2195 /*
2196 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2197 % %
2198 % %
2199 % %
2200 % G e t I m a g e C h a n n e l R a n g e %
2201 % %
2202 % %
2203 % %
2204 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2205 %
2206 % GetImageChannelRange() returns the range of one or more image channels.
2207 %
2208 % The format of the GetImageChannelRange method is:
2209 %
2210 % MagickBooleanType GetImageChannelRange(const Image *image,
2211 % const ChannelType channel,double *minima,double *maxima,
2212 % ExceptionInfo *exception)
2213 %
2214 % A description of each parameter follows:
2215 %
2216 % o image: the image.
2217 %
2218 % o channel: the channel.
2219 %
2220 % o minima: the minimum value in the channel.
2221 %
2222 % o maxima: the maximum value in the channel.
2223 %
2224 % o exception: return any errors or warnings in this structure.
2225 %
2226 */
2227 
2228 MagickExport MagickBooleanType GetImageRange(const Image *image,
2229  double *minima,double *maxima,ExceptionInfo *exception)
2230 {
2231  return(GetImageChannelRange(image,CompositeChannels,minima,maxima,exception));
2232 }
2233 
2234 MagickExport MagickBooleanType GetImageChannelRange(const Image *image,
2235  const ChannelType channel,double *minima,double *maxima,
2236  ExceptionInfo *exception)
2237 {
2239  pixel;
2240 
2241  ssize_t
2242  y;
2243 
2244  assert(image != (Image *) NULL);
2245  assert(image->signature == MagickCoreSignature);
2246  if (IsEventLogging() != MagickFalse)
2247  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2248  *maxima=(-MagickMaximumValue);
2249  *minima=MagickMaximumValue;
2250  GetMagickPixelPacket(image,&pixel);
2251  for (y=0; y < (ssize_t) image->rows; y++)
2252  {
2253  const IndexPacket
2254  *magick_restrict indexes;
2255 
2256  const PixelPacket
2257  *magick_restrict p;
2258 
2259  ssize_t
2260  x;
2261 
2262  p=GetVirtualPixels(image,0,y,image->columns,1,exception);
2263  if (p == (const PixelPacket *) NULL)
2264  break;
2265  indexes=GetVirtualIndexQueue(image);
2266  for (x=0; x < (ssize_t) image->columns; x++)
2267  {
2268  SetMagickPixelPacket(image,p,indexes+x,&pixel);
2269  if ((channel & RedChannel) != 0)
2270  {
2271  if (pixel.red < *minima)
2272  *minima=(double) pixel.red;
2273  if (pixel.red > *maxima)
2274  *maxima=(double) pixel.red;
2275  }
2276  if ((channel & GreenChannel) != 0)
2277  {
2278  if (pixel.green < *minima)
2279  *minima=(double) pixel.green;
2280  if (pixel.green > *maxima)
2281  *maxima=(double) pixel.green;
2282  }
2283  if ((channel & BlueChannel) != 0)
2284  {
2285  if (pixel.blue < *minima)
2286  *minima=(double) pixel.blue;
2287  if (pixel.blue > *maxima)
2288  *maxima=(double) pixel.blue;
2289  }
2290  if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
2291  {
2292  if ((QuantumRange-pixel.opacity) < *minima)
2293  *minima=(double) (QuantumRange-pixel.opacity);
2294  if ((QuantumRange-pixel.opacity) > *maxima)
2295  *maxima=(double) (QuantumRange-pixel.opacity);
2296  }
2297  if (((channel & IndexChannel) != 0) &&
2298  (image->colorspace == CMYKColorspace))
2299  {
2300  if ((double) pixel.index < *minima)
2301  *minima=(double) pixel.index;
2302  if ((double) pixel.index > *maxima)
2303  *maxima=(double) pixel.index;
2304  }
2305  p++;
2306  }
2307  }
2308  return(y == (ssize_t) image->rows ? MagickTrue : MagickFalse);
2309 }
2310 
2311 /*
2312 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2313 % %
2314 % %
2315 % %
2316 % 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 %
2317 % %
2318 % %
2319 % %
2320 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2321 %
2322 % GetImageChannelStatistics() returns statistics for each channel in the
2323 % image. The statistics include the channel depth, its minima, maxima, mean,
2324 % standard deviation, kurtosis and skewness. You can access the red channel
2325 % mean, for example, like this:
2326 %
2327 % channel_statistics=GetImageChannelStatistics(image,exception);
2328 % red_mean=channel_statistics[RedChannel].mean;
2329 %
2330 % Use MagickRelinquishMemory() to free the statistics buffer.
2331 %
2332 % The format of the GetImageChannelStatistics method is:
2333 %
2334 % ChannelStatistics *GetImageChannelStatistics(const Image *image,
2335 % ExceptionInfo *exception)
2336 %
2337 % A description of each parameter follows:
2338 %
2339 % o image: the image.
2340 %
2341 % o exception: return any errors or warnings in this structure.
2342 %
2343 */
2344 MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
2345  ExceptionInfo *exception)
2346 {
2348  *channel_statistics;
2349 
2350  double
2351  area,
2352  standard_deviation;
2353 
2355  number_bins,
2356  *histogram;
2357 
2358  QuantumAny
2359  range;
2360 
2361  ssize_t
2362  i;
2363 
2364  size_t
2365  channels,
2366  depth,
2367  length;
2368 
2369  ssize_t
2370  y;
2371 
2372  assert(image != (Image *) NULL);
2373  assert(image->signature == MagickCoreSignature);
2374  if (IsEventLogging() != MagickFalse)
2375  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2376  length=CompositeChannels+1UL;
2377  channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(length,
2378  sizeof(*channel_statistics));
2379  histogram=(MagickPixelPacket *) AcquireQuantumMemory(MaxMap+1U,
2380  sizeof(*histogram));
2381  if ((channel_statistics == (ChannelStatistics *) NULL) ||
2382  (histogram == (MagickPixelPacket *) NULL))
2383  {
2384  if (histogram != (MagickPixelPacket *) NULL)
2385  histogram=(MagickPixelPacket *) RelinquishMagickMemory(histogram);
2386  if (channel_statistics != (ChannelStatistics *) NULL)
2387  channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
2388  channel_statistics);
2389  return(channel_statistics);
2390  }
2391  (void) memset(channel_statistics,0,length*
2392  sizeof(*channel_statistics));
2393  for (i=0; i <= (ssize_t) CompositeChannels; i++)
2394  {
2395  channel_statistics[i].depth=1;
2396  channel_statistics[i].maxima=(-MagickMaximumValue);
2397  channel_statistics[i].minima=MagickMaximumValue;
2398  }
2399  (void) memset(histogram,0,(MaxMap+1U)*sizeof(*histogram));
2400  (void) memset(&number_bins,0,sizeof(number_bins));
2401  for (y=0; y < (ssize_t) image->rows; y++)
2402  {
2403  const IndexPacket
2404  *magick_restrict indexes;
2405 
2406  const PixelPacket
2407  *magick_restrict p;
2408 
2409  ssize_t
2410  x;
2411 
2412  /*
2413  Compute pixel statistics.
2414  */
2415  p=GetVirtualPixels(image,0,y,image->columns,1,exception);
2416  if (p == (const PixelPacket *) NULL)
2417  break;
2418  indexes=GetVirtualIndexQueue(image);
2419  for (x=0; x < (ssize_t) image->columns; )
2420  {
2421  if (channel_statistics[RedChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2422  {
2423  depth=channel_statistics[RedChannel].depth;
2424  range=GetQuantumRange(depth);
2425  if (IsPixelAtDepth(GetPixelRed(p),range) == MagickFalse)
2426  {
2427  channel_statistics[RedChannel].depth++;
2428  continue;
2429  }
2430  }
2431  if (channel_statistics[GreenChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2432  {
2433  depth=channel_statistics[GreenChannel].depth;
2434  range=GetQuantumRange(depth);
2435  if (IsPixelAtDepth(GetPixelGreen(p),range) == MagickFalse)
2436  {
2437  channel_statistics[GreenChannel].depth++;
2438  continue;
2439  }
2440  }
2441  if (channel_statistics[BlueChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2442  {
2443  depth=channel_statistics[BlueChannel].depth;
2444  range=GetQuantumRange(depth);
2445  if (IsPixelAtDepth(GetPixelBlue(p),range) == MagickFalse)
2446  {
2447  channel_statistics[BlueChannel].depth++;
2448  continue;
2449  }
2450  }
2451  if (image->matte != MagickFalse)
2452  {
2453  if (channel_statistics[OpacityChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2454  {
2455  depth=channel_statistics[OpacityChannel].depth;
2456  range=GetQuantumRange(depth);
2457  if (IsPixelAtDepth(GetPixelAlpha(p),range) == MagickFalse)
2458  {
2459  channel_statistics[OpacityChannel].depth++;
2460  continue;
2461  }
2462  }
2463  }
2464  if (image->colorspace == CMYKColorspace)
2465  {
2466  if (channel_statistics[BlackChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2467  {
2468  depth=channel_statistics[BlackChannel].depth;
2469  range=GetQuantumRange(depth);
2470  if (IsPixelAtDepth(GetPixelIndex(indexes+x),range) == MagickFalse)
2471  {
2472  channel_statistics[BlackChannel].depth++;
2473  continue;
2474  }
2475  }
2476  }
2477  if ((double) GetPixelRed(p) < channel_statistics[RedChannel].minima)
2478  channel_statistics[RedChannel].minima=(double) GetPixelRed(p);
2479  if ((double) GetPixelRed(p) > channel_statistics[RedChannel].maxima)
2480  channel_statistics[RedChannel].maxima=(double) GetPixelRed(p);
2481  channel_statistics[RedChannel].sum+=GetPixelRed(p);
2482  channel_statistics[RedChannel].sum_squared+=(double) GetPixelRed(p)*
2483  GetPixelRed(p);
2484  channel_statistics[RedChannel].sum_cubed+=(double)
2485  GetPixelRed(p)*GetPixelRed(p)*GetPixelRed(p);
2486  channel_statistics[RedChannel].sum_fourth_power+=(double)
2487  GetPixelRed(p)*GetPixelRed(p)*GetPixelRed(p)*GetPixelRed(p);
2488  if ((double) GetPixelGreen(p) < channel_statistics[GreenChannel].minima)
2489  channel_statistics[GreenChannel].minima=(double) GetPixelGreen(p);
2490  if ((double) GetPixelGreen(p) > channel_statistics[GreenChannel].maxima)
2491  channel_statistics[GreenChannel].maxima=(double) GetPixelGreen(p);
2492  channel_statistics[GreenChannel].sum+=GetPixelGreen(p);
2493  channel_statistics[GreenChannel].sum_squared+=(double) GetPixelGreen(p)*
2494  GetPixelGreen(p);
2495  channel_statistics[GreenChannel].sum_cubed+=(double) GetPixelGreen(p)*
2496  GetPixelGreen(p)*GetPixelGreen(p);
2497  channel_statistics[GreenChannel].sum_fourth_power+=(double)
2498  GetPixelGreen(p)*GetPixelGreen(p)*GetPixelGreen(p)*GetPixelGreen(p);
2499  if ((double) GetPixelBlue(p) < channel_statistics[BlueChannel].minima)
2500  channel_statistics[BlueChannel].minima=(double) GetPixelBlue(p);
2501  if ((double) GetPixelBlue(p) > channel_statistics[BlueChannel].maxima)
2502  channel_statistics[BlueChannel].maxima=(double) GetPixelBlue(p);
2503  channel_statistics[BlueChannel].sum+=GetPixelBlue(p);
2504  channel_statistics[BlueChannel].sum_squared+=(double) GetPixelBlue(p)*
2505  GetPixelBlue(p);
2506  channel_statistics[BlueChannel].sum_cubed+=(double) GetPixelBlue(p)*
2507  GetPixelBlue(p)*GetPixelBlue(p);
2508  channel_statistics[BlueChannel].sum_fourth_power+=(double)
2509  GetPixelBlue(p)*GetPixelBlue(p)*GetPixelBlue(p)*GetPixelBlue(p);
2510  histogram[ScaleQuantumToMap(GetPixelRed(p))].red++;
2511  histogram[ScaleQuantumToMap(GetPixelGreen(p))].green++;
2512  histogram[ScaleQuantumToMap(GetPixelBlue(p))].blue++;
2513  if (image->matte != MagickFalse)
2514  {
2515  if ((double) GetPixelAlpha(p) < channel_statistics[OpacityChannel].minima)
2516  channel_statistics[OpacityChannel].minima=(double) GetPixelAlpha(p);
2517  if ((double) GetPixelAlpha(p) > channel_statistics[OpacityChannel].maxima)
2518  channel_statistics[OpacityChannel].maxima=(double) GetPixelAlpha(p);
2519  channel_statistics[OpacityChannel].sum+=GetPixelAlpha(p);
2520  channel_statistics[OpacityChannel].sum_squared+=(double)
2521  GetPixelAlpha(p)*GetPixelAlpha(p);
2522  channel_statistics[OpacityChannel].sum_cubed+=(double)
2523  GetPixelAlpha(p)*GetPixelAlpha(p)*GetPixelAlpha(p);
2524  channel_statistics[OpacityChannel].sum_fourth_power+=(double)
2525  GetPixelAlpha(p)*GetPixelAlpha(p)*GetPixelAlpha(p)*GetPixelAlpha(p);
2526  histogram[ScaleQuantumToMap(GetPixelAlpha(p))].opacity++;
2527  }
2528  if (image->colorspace == CMYKColorspace)
2529  {
2530  if ((double) GetPixelIndex(indexes+x) < channel_statistics[BlackChannel].minima)
2531  channel_statistics[BlackChannel].minima=(double)
2532  GetPixelIndex(indexes+x);
2533  if ((double) GetPixelIndex(indexes+x) > channel_statistics[BlackChannel].maxima)
2534  channel_statistics[BlackChannel].maxima=(double)
2535  GetPixelIndex(indexes+x);
2536  channel_statistics[BlackChannel].sum+=GetPixelIndex(indexes+x);
2537  channel_statistics[BlackChannel].sum_squared+=(double)
2538  GetPixelIndex(indexes+x)*GetPixelIndex(indexes+x);
2539  channel_statistics[BlackChannel].sum_cubed+=(double)
2540  GetPixelIndex(indexes+x)*GetPixelIndex(indexes+x)*
2541  GetPixelIndex(indexes+x);
2542  channel_statistics[BlackChannel].sum_fourth_power+=(double)
2543  GetPixelIndex(indexes+x)*GetPixelIndex(indexes+x)*
2544  GetPixelIndex(indexes+x)*GetPixelIndex(indexes+x);
2545  histogram[ScaleQuantumToMap(GetPixelIndex(indexes+x))].index++;
2546  }
2547  x++;
2548  p++;
2549  }
2550  }
2551  for (i=0; i < (ssize_t) CompositeChannels; i++)
2552  {
2553  double
2554  area,
2555  mean,
2556  standard_deviation;
2557 
2558  /*
2559  Normalize pixel statistics.
2560  */
2561  area=PerceptibleReciprocal((double) image->columns*image->rows);
2562  mean=channel_statistics[i].sum*area;
2563  channel_statistics[i].sum=mean;
2564  channel_statistics[i].sum_squared*=area;
2565  channel_statistics[i].sum_cubed*=area;
2566  channel_statistics[i].sum_fourth_power*=area;
2567  channel_statistics[i].mean=mean;
2568  channel_statistics[i].variance=channel_statistics[i].sum_squared;
2569  standard_deviation=sqrt(channel_statistics[i].variance-(mean*mean));
2570  area=PerceptibleReciprocal((double) image->columns*image->rows-1.0)*
2571  ((double) image->columns*image->rows);
2572  standard_deviation=sqrt(area*standard_deviation*standard_deviation);
2573  channel_statistics[i].standard_deviation=standard_deviation;
2574  }
2575  for (i=0; i < (ssize_t) (MaxMap+1U); i++)
2576  {
2577  if (histogram[i].red > 0.0)
2578  number_bins.red++;
2579  if (histogram[i].green > 0.0)
2580  number_bins.green++;
2581  if (histogram[i].blue > 0.0)
2582  number_bins.blue++;
2583  if ((image->matte != MagickFalse) && (histogram[i].opacity > 0.0))
2584  number_bins.opacity++;
2585  if ((image->colorspace == CMYKColorspace) && (histogram[i].index > 0.0))
2586  number_bins.index++;
2587  }
2588  area=PerceptibleReciprocal((double) image->columns*image->rows);
2589  for (i=0; i < (ssize_t) (MaxMap+1U); i++)
2590  {
2591  /*
2592  Compute pixel entropy.
2593  */
2594  histogram[i].red*=area;
2595  channel_statistics[RedChannel].entropy+=-histogram[i].red*
2596  MagickLog10(histogram[i].red)*
2597  PerceptibleReciprocal(MagickLog10((double) number_bins.red));
2598  histogram[i].green*=area;
2599  channel_statistics[GreenChannel].entropy+=-histogram[i].green*
2600  MagickLog10(histogram[i].green)*
2601  PerceptibleReciprocal(MagickLog10((double) number_bins.green));
2602  histogram[i].blue*=area;
2603  channel_statistics[BlueChannel].entropy+=-histogram[i].blue*
2604  MagickLog10(histogram[i].blue)*
2605  PerceptibleReciprocal(MagickLog10((double) number_bins.blue));
2606  if (image->matte != MagickFalse)
2607  {
2608  histogram[i].opacity*=area;
2609  channel_statistics[OpacityChannel].entropy+=-histogram[i].opacity*
2610  MagickLog10(histogram[i].opacity)*
2611  PerceptibleReciprocal(MagickLog10((double) number_bins.opacity));
2612  }
2613  if (image->colorspace == CMYKColorspace)
2614  {
2615  histogram[i].index*=area;
2616  channel_statistics[IndexChannel].entropy+=-histogram[i].index*
2617  MagickLog10(histogram[i].index)*
2618  PerceptibleReciprocal(MagickLog10((double) number_bins.index));
2619  }
2620  }
2621  /*
2622  Compute overall statistics.
2623  */
2624  for (i=0; i < (ssize_t) CompositeChannels; i++)
2625  {
2626  channel_statistics[CompositeChannels].depth=(size_t) EvaluateMax((double)
2627  channel_statistics[CompositeChannels].depth,(double)
2628  channel_statistics[i].depth);
2629  channel_statistics[CompositeChannels].minima=MagickMin(
2630  channel_statistics[CompositeChannels].minima,
2631  channel_statistics[i].minima);
2632  channel_statistics[CompositeChannels].maxima=EvaluateMax(
2633  channel_statistics[CompositeChannels].maxima,
2634  channel_statistics[i].maxima);
2635  channel_statistics[CompositeChannels].sum+=channel_statistics[i].sum;
2636  channel_statistics[CompositeChannels].sum_squared+=
2637  channel_statistics[i].sum_squared;
2638  channel_statistics[CompositeChannels].sum_cubed+=
2639  channel_statistics[i].sum_cubed;
2640  channel_statistics[CompositeChannels].sum_fourth_power+=
2641  channel_statistics[i].sum_fourth_power;
2642  channel_statistics[CompositeChannels].mean+=channel_statistics[i].mean;
2643  channel_statistics[CompositeChannels].variance+=
2644  channel_statistics[i].variance-channel_statistics[i].mean*
2645  channel_statistics[i].mean;
2646  standard_deviation=sqrt(channel_statistics[i].variance-
2647  (channel_statistics[i].mean*channel_statistics[i].mean));
2648  area=PerceptibleReciprocal((double) image->columns*image->rows-1.0)*
2649  ((double) image->columns*image->rows);
2650  standard_deviation=sqrt(area*standard_deviation*standard_deviation);
2651  channel_statistics[CompositeChannels].standard_deviation=standard_deviation;
2652  channel_statistics[CompositeChannels].entropy+=
2653  channel_statistics[i].entropy;
2654  }
2655  channels=3;
2656  if (image->matte != MagickFalse)
2657  channels++;
2658  if (image->colorspace == CMYKColorspace)
2659  channels++;
2660  channel_statistics[CompositeChannels].sum/=channels;
2661  channel_statistics[CompositeChannels].sum_squared/=channels;
2662  channel_statistics[CompositeChannels].sum_cubed/=channels;
2663  channel_statistics[CompositeChannels].sum_fourth_power/=channels;
2664  channel_statistics[CompositeChannels].mean/=channels;
2665  channel_statistics[CompositeChannels].kurtosis/=channels;
2666  channel_statistics[CompositeChannels].skewness/=channels;
2667  channel_statistics[CompositeChannels].entropy/=channels;
2668  i=CompositeChannels;
2669  area=PerceptibleReciprocal((double) channels*image->columns*image->rows);
2670  channel_statistics[i].variance=channel_statistics[i].sum_squared;
2671  channel_statistics[i].mean=channel_statistics[i].sum;
2672  standard_deviation=sqrt(channel_statistics[i].variance-
2673  (channel_statistics[i].mean*channel_statistics[i].mean));
2674  standard_deviation=sqrt(PerceptibleReciprocal((double) channels*
2675  image->columns*image->rows-1.0)*channels*image->columns*image->rows*
2676  standard_deviation*standard_deviation);
2677  channel_statistics[i].standard_deviation=standard_deviation;
2678  for (i=0; i <= (ssize_t) CompositeChannels; i++)
2679  {
2680  /*
2681  Compute kurtosis & skewness statistics.
2682  */
2683  standard_deviation=PerceptibleReciprocal(
2684  channel_statistics[i].standard_deviation);
2685  channel_statistics[i].skewness=(channel_statistics[i].sum_cubed-3.0*
2686  channel_statistics[i].mean*channel_statistics[i].sum_squared+2.0*
2687  channel_statistics[i].mean*channel_statistics[i].mean*
2688  channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2689  standard_deviation);
2690  channel_statistics[i].kurtosis=(channel_statistics[i].sum_fourth_power-4.0*
2691  channel_statistics[i].mean*channel_statistics[i].sum_cubed+6.0*
2692  channel_statistics[i].mean*channel_statistics[i].mean*
2693  channel_statistics[i].sum_squared-3.0*channel_statistics[i].mean*
2694  channel_statistics[i].mean*1.0*channel_statistics[i].mean*
2695  channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2696  standard_deviation*standard_deviation)-3.0;
2697  }
2698  channel_statistics[CompositeChannels].mean=0.0;
2699  channel_statistics[CompositeChannels].standard_deviation=0.0;
2700  for (i=0; i < (ssize_t) CompositeChannels; i++)
2701  {
2702  channel_statistics[CompositeChannels].mean+=
2703  channel_statistics[i].mean;
2704  channel_statistics[CompositeChannels].standard_deviation+=
2705  channel_statistics[i].standard_deviation;
2706  }
2707  channel_statistics[CompositeChannels].mean/=(double) channels;
2708  channel_statistics[CompositeChannels].standard_deviation/=(double) channels;
2709  histogram=(MagickPixelPacket *) RelinquishMagickMemory(histogram);
2710  if (y < (ssize_t) image->rows)
2711  channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
2712  channel_statistics);
2713  return(channel_statistics);
2714 }
2715 
2716 /*
2717 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2718 % %
2719 % %
2720 % %
2721 % P o l y n o m i a l I m a g e %
2722 % %
2723 % %
2724 % %
2725 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2726 %
2727 % PolynomialImage() returns a new image where each pixel is the sum of the
2728 % pixels in the image sequence after applying its corresponding terms
2729 % (coefficient and degree pairs).
2730 %
2731 % The format of the PolynomialImage method is:
2732 %
2733 % Image *PolynomialImage(const Image *images,const size_t number_terms,
2734 % const double *terms,ExceptionInfo *exception)
2735 % Image *PolynomialImageChannel(const Image *images,
2736 % const size_t number_terms,const ChannelType channel,
2737 % const double *terms,ExceptionInfo *exception)
2738 %
2739 % A description of each parameter follows:
2740 %
2741 % o images: the image sequence.
2742 %
2743 % o channel: the channel.
2744 %
2745 % o number_terms: the number of terms in the list. The actual list length
2746 % is 2 x number_terms + 1 (the constant).
2747 %
2748 % o terms: the list of polynomial coefficients and degree pairs and a
2749 % constant.
2750 %
2751 % o exception: return any errors or warnings in this structure.
2752 %
2753 */
2754 MagickExport Image *PolynomialImage(const Image *images,
2755  const size_t number_terms,const double *terms,ExceptionInfo *exception)
2756 {
2757  Image
2758  *polynomial_image;
2759 
2760  polynomial_image=PolynomialImageChannel(images,DefaultChannels,number_terms,
2761  terms,exception);
2762  return(polynomial_image);
2763 }
2764 
2765 MagickExport Image *PolynomialImageChannel(const Image *images,
2766  const ChannelType channel,const size_t number_terms,const double *terms,
2767  ExceptionInfo *exception)
2768 {
2769 #define PolynomialImageTag "Polynomial/Image"
2770 
2771  CacheView
2772  *polynomial_view;
2773 
2774  Image
2775  *image;
2776 
2777  MagickBooleanType
2778  status;
2779 
2780  MagickOffsetType
2781  progress;
2782 
2784  **magick_restrict polynomial_pixels,
2785  zero;
2786 
2787  ssize_t
2788  y;
2789 
2790  assert(images != (Image *) NULL);
2791  assert(images->signature == MagickCoreSignature);
2792  if (IsEventLogging() != MagickFalse)
2793  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
2794  assert(exception != (ExceptionInfo *) NULL);
2795  assert(exception->signature == MagickCoreSignature);
2796  image=AcquireImageCanvas(images,exception);
2797  if (image == (Image *) NULL)
2798  return((Image *) NULL);
2799  if (SetImageStorageClass(image,DirectClass) == MagickFalse)
2800  {
2801  InheritException(exception,&image->exception);
2802  image=DestroyImage(image);
2803  return((Image *) NULL);
2804  }
2805  polynomial_pixels=AcquirePixelTLS(images);
2806  if (polynomial_pixels == (MagickPixelPacket **) NULL)
2807  {
2808  image=DestroyImage(image);
2809  (void) ThrowMagickException(exception,GetMagickModule(),
2810  ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
2811  return((Image *) NULL);
2812  }
2813  /*
2814  Polynomial image pixels.
2815  */
2816  status=MagickTrue;
2817  progress=0;
2818  GetMagickPixelPacket(images,&zero);
2819  polynomial_view=AcquireAuthenticCacheView(image,exception);
2820 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2821  #pragma omp parallel for schedule(static) shared(progress,status) \
2822  magick_number_threads(image,image,image->rows,1)
2823 #endif
2824  for (y=0; y < (ssize_t) image->rows; y++)
2825  {
2826  CacheView
2827  *image_view;
2828 
2829  const Image
2830  *next;
2831 
2832  const int
2833  id = GetOpenMPThreadId();
2834 
2835  IndexPacket
2836  *magick_restrict polynomial_indexes;
2837 
2839  *polynomial_pixel;
2840 
2841  PixelPacket
2842  *magick_restrict q;
2843 
2844  ssize_t
2845  i,
2846  x;
2847 
2848  size_t
2849  number_images;
2850 
2851  if (status == MagickFalse)
2852  continue;
2853  q=QueueCacheViewAuthenticPixels(polynomial_view,0,y,image->columns,1,
2854  exception);
2855  if (q == (PixelPacket *) NULL)
2856  {
2857  status=MagickFalse;
2858  continue;
2859  }
2860  polynomial_indexes=GetCacheViewAuthenticIndexQueue(polynomial_view);
2861  polynomial_pixel=polynomial_pixels[id];
2862  for (x=0; x < (ssize_t) image->columns; x++)
2863  polynomial_pixel[x]=zero;
2864  next=images;
2865  number_images=GetImageListLength(images);
2866  for (i=0; i < (ssize_t) number_images; i++)
2867  {
2868  const IndexPacket
2869  *indexes;
2870 
2871  const PixelPacket
2872  *p;
2873 
2874  if (i >= (ssize_t) number_terms)
2875  break;
2876  image_view=AcquireVirtualCacheView(next,exception);
2877  p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2878  if (p == (const PixelPacket *) NULL)
2879  {
2880  image_view=DestroyCacheView(image_view);
2881  break;
2882  }
2883  indexes=GetCacheViewVirtualIndexQueue(image_view);
2884  for (x=0; x < (ssize_t) image->columns; x++)
2885  {
2886  double
2887  coefficient,
2888  degree;
2889 
2890  coefficient=terms[i << 1];
2891  degree=terms[(i << 1)+1];
2892  if ((channel & RedChannel) != 0)
2893  polynomial_pixel[x].red+=coefficient*pow(QuantumScale*p->red,degree);
2894  if ((channel & GreenChannel) != 0)
2895  polynomial_pixel[x].green+=coefficient*pow(QuantumScale*p->green,
2896  degree);
2897  if ((channel & BlueChannel) != 0)
2898  polynomial_pixel[x].blue+=coefficient*pow(QuantumScale*p->blue,
2899  degree);
2900  if ((channel & OpacityChannel) != 0)
2901  polynomial_pixel[x].opacity+=coefficient*pow(QuantumScale*
2902  (QuantumRange-p->opacity),degree);
2903  if (((channel & IndexChannel) != 0) &&
2904  (image->colorspace == CMYKColorspace))
2905  polynomial_pixel[x].index+=coefficient*pow(QuantumScale*indexes[x],
2906  degree);
2907  p++;
2908  }
2909  image_view=DestroyCacheView(image_view);
2910  next=GetNextImageInList(next);
2911  }
2912  for (x=0; x < (ssize_t) image->columns; x++)
2913  {
2914  SetPixelRed(q,ClampToQuantum(QuantumRange*polynomial_pixel[x].red));
2915  SetPixelGreen(q,ClampToQuantum(QuantumRange*polynomial_pixel[x].green));
2916  SetPixelBlue(q,ClampToQuantum(QuantumRange*polynomial_pixel[x].blue));
2917  if (image->matte == MagickFalse)
2918  SetPixelOpacity(q,ClampToQuantum(QuantumRange-QuantumRange*
2919  polynomial_pixel[x].opacity));
2920  else
2921  SetPixelAlpha(q,ClampToQuantum(QuantumRange-QuantumRange*
2922  polynomial_pixel[x].opacity));
2923  if (image->colorspace == CMYKColorspace)
2924  SetPixelIndex(polynomial_indexes+x,ClampToQuantum(QuantumRange*
2925  polynomial_pixel[x].index));
2926  q++;
2927  }
2928  if (SyncCacheViewAuthenticPixels(polynomial_view,exception) == MagickFalse)
2929  status=MagickFalse;
2930  if (images->progress_monitor != (MagickProgressMonitor) NULL)
2931  {
2932  MagickBooleanType
2933  proceed;
2934 
2935  proceed=SetImageProgress(images,PolynomialImageTag,progress++,
2936  image->rows);
2937  if (proceed == MagickFalse)
2938  status=MagickFalse;
2939  }
2940  }
2941  polynomial_view=DestroyCacheView(polynomial_view);
2942  polynomial_pixels=DestroyPixelTLS(images,polynomial_pixels);
2943  if (status == MagickFalse)
2944  image=DestroyImage(image);
2945  return(image);
2946 }
2947 
2948 /*
2949 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2950 % %
2951 % %
2952 % %
2953 % S t a t i s t i c I m a g e %
2954 % %
2955 % %
2956 % %
2957 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2958 %
2959 % StatisticImage() makes each pixel the min / max / median / mode / etc. of
2960 % the neighborhood of the specified width and height.
2961 %
2962 % The format of the StatisticImage method is:
2963 %
2964 % Image *StatisticImage(const Image *image,const StatisticType type,
2965 % const size_t width,const size_t height,ExceptionInfo *exception)
2966 % Image *StatisticImageChannel(const Image *image,
2967 % const ChannelType channel,const StatisticType type,
2968 % const size_t width,const size_t height,ExceptionInfo *exception)
2969 %
2970 % A description of each parameter follows:
2971 %
2972 % o image: the image.
2973 %
2974 % o channel: the image channel.
2975 %
2976 % o type: the statistic type (median, mode, etc.).
2977 %
2978 % o width: the width of the pixel neighborhood.
2979 %
2980 % o height: the height of the pixel neighborhood.
2981 %
2982 % o exception: return any errors or warnings in this structure.
2983 %
2984 */
2985 
2986 #define ListChannels 5
2987 
2988 typedef struct _ListNode
2989 {
2990  size_t
2991  next[9],
2992  count,
2993  signature;
2994 } ListNode;
2995 
2996 typedef struct _SkipList
2997 {
2998  ssize_t
2999  level;
3000 
3001  ListNode
3002  *nodes;
3003 } SkipList;
3004 
3005 typedef struct _PixelList
3006 {
3007  size_t
3008  length,
3009  seed,
3010  signature;
3011 
3012  SkipList
3013  lists[ListChannels];
3014 } PixelList;
3015 
3016 static PixelList *DestroyPixelList(PixelList *pixel_list)
3017 {
3018  ssize_t
3019  i;
3020 
3021  if (pixel_list == (PixelList *) NULL)
3022  return((PixelList *) NULL);
3023  for (i=0; i < ListChannels; i++)
3024  if (pixel_list->lists[i].nodes != (ListNode *) NULL)
3025  pixel_list->lists[i].nodes=(ListNode *) RelinquishAlignedMemory(
3026  pixel_list->lists[i].nodes);
3027  pixel_list=(PixelList *) RelinquishMagickMemory(pixel_list);
3028  return(pixel_list);
3029 }
3030 
3031 static PixelList **DestroyPixelListTLS(PixelList **pixel_list)
3032 {
3033  ssize_t
3034  i;
3035 
3036  assert(pixel_list != (PixelList **) NULL);
3037  for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
3038  if (pixel_list[i] != (PixelList *) NULL)
3039  pixel_list[i]=DestroyPixelList(pixel_list[i]);
3040  pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list);
3041  return(pixel_list);
3042 }
3043 
3044 static PixelList *AcquirePixelList(const size_t width,const size_t height)
3045 {
3046  PixelList
3047  *pixel_list;
3048 
3049  ssize_t
3050  i;
3051 
3052  pixel_list=(PixelList *) AcquireMagickMemory(sizeof(*pixel_list));
3053  if (pixel_list == (PixelList *) NULL)
3054  return(pixel_list);
3055  (void) memset((void *) pixel_list,0,sizeof(*pixel_list));
3056  pixel_list->length=width*height;
3057  for (i=0; i < ListChannels; i++)
3058  {
3059  pixel_list->lists[i].nodes=(ListNode *) AcquireAlignedMemory(65537UL,
3060  sizeof(*pixel_list->lists[i].nodes));
3061  if (pixel_list->lists[i].nodes == (ListNode *) NULL)
3062  return(DestroyPixelList(pixel_list));
3063  (void) memset(pixel_list->lists[i].nodes,0,65537UL*
3064  sizeof(*pixel_list->lists[i].nodes));
3065  }
3066  pixel_list->signature=MagickCoreSignature;
3067  return(pixel_list);
3068 }
3069 
3070 static PixelList **AcquirePixelListTLS(const size_t width,
3071  const size_t height)
3072 {
3073  PixelList
3074  **pixel_list;
3075 
3076  ssize_t
3077  i;
3078 
3079  size_t
3080  number_threads;
3081 
3082  number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
3083  pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
3084  sizeof(*pixel_list));
3085  if (pixel_list == (PixelList **) NULL)
3086  return((PixelList **) NULL);
3087  (void) memset(pixel_list,0,number_threads*sizeof(*pixel_list));
3088  for (i=0; i < (ssize_t) number_threads; i++)
3089  {
3090  pixel_list[i]=AcquirePixelList(width,height);
3091  if (pixel_list[i] == (PixelList *) NULL)
3092  return(DestroyPixelListTLS(pixel_list));
3093  }
3094  return(pixel_list);
3095 }
3096 
3097 static void AddNodePixelList(PixelList *pixel_list,const ssize_t channel,
3098  const size_t color)
3099 {
3100  SkipList
3101  *list;
3102 
3103  ssize_t
3104  level;
3105 
3106  size_t
3107  search,
3108  update[9];
3109 
3110  /*
3111  Initialize the node.
3112  */
3113  list=pixel_list->lists+channel;
3114  list->nodes[color].signature=pixel_list->signature;
3115  list->nodes[color].count=1;
3116  /*
3117  Determine where it belongs in the list.
3118  */
3119  search=65536UL;
3120  for (level=list->level; level >= 0; level--)
3121  {
3122  while (list->nodes[search].next[level] < color)
3123  search=list->nodes[search].next[level];
3124  update[level]=search;
3125  }
3126  /*
3127  Generate a pseudo-random level for this node.
3128  */
3129  for (level=0; ; level++)
3130  {
3131  pixel_list->seed=(pixel_list->seed*42893621L)+1L;
3132  if ((pixel_list->seed & 0x300) != 0x300)
3133  break;
3134  }
3135  if (level > 8)
3136  level=8;
3137  if (level > (list->level+2))
3138  level=list->level+2;
3139  /*
3140  If we're raising the list's level, link back to the root node.
3141  */
3142  while (level > list->level)
3143  {
3144  list->level++;
3145  update[list->level]=65536UL;
3146  }
3147  /*
3148  Link the node into the skip-list.
3149  */
3150  do
3151  {
3152  list->nodes[color].next[level]=list->nodes[update[level]].next[level];
3153  list->nodes[update[level]].next[level]=color;
3154  } while (level-- > 0);
3155 }
3156 
3157 static void GetMaximumPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3158 {
3159  SkipList
3160  *list;
3161 
3162  ssize_t
3163  channel;
3164 
3165  size_t
3166  color,
3167  maximum;
3168 
3169  ssize_t
3170  count;
3171 
3172  unsigned short
3173  channels[ListChannels];
3174 
3175  /*
3176  Find the maximum value for each of the color.
3177  */
3178  for (channel=0; channel < 5; channel++)
3179  {
3180  list=pixel_list->lists+channel;
3181  color=65536L;
3182  count=0;
3183  maximum=list->nodes[color].next[0];
3184  do
3185  {
3186  color=list->nodes[color].next[0];
3187  if (color > maximum)
3188  maximum=color;
3189  count+=list->nodes[color].count;
3190  } while (count < (ssize_t) pixel_list->length);
3191  channels[channel]=(unsigned short) maximum;
3192  }
3193  pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3194  pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3195  pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3196  pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3197  pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3198 }
3199 
3200 static void GetMeanPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3201 {
3202  MagickRealType
3203  sum;
3204 
3205  SkipList
3206  *list;
3207 
3208  ssize_t
3209  channel;
3210 
3211  size_t
3212  color;
3213 
3214  ssize_t
3215  count;
3216 
3217  unsigned short
3218  channels[ListChannels];
3219 
3220  /*
3221  Find the mean value for each of the color.
3222  */
3223  for (channel=0; channel < 5; channel++)
3224  {
3225  list=pixel_list->lists+channel;
3226  color=65536L;
3227  count=0;
3228  sum=0.0;
3229  do
3230  {
3231  color=list->nodes[color].next[0];
3232  sum+=(MagickRealType) list->nodes[color].count*color;
3233  count+=list->nodes[color].count;
3234  } while (count < (ssize_t) pixel_list->length);
3235  sum/=pixel_list->length;
3236  channels[channel]=(unsigned short) sum;
3237  }
3238  pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3239  pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3240  pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3241  pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3242  pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3243 }
3244 
3245 static void GetMedianPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3246 {
3247  SkipList
3248  *list;
3249 
3250  ssize_t
3251  channel;
3252 
3253  size_t
3254  color;
3255 
3256  ssize_t
3257  count;
3258 
3259  unsigned short
3260  channels[ListChannels];
3261 
3262  /*
3263  Find the median value for each of the color.
3264  */
3265  for (channel=0; channel < 5; channel++)
3266  {
3267  list=pixel_list->lists+channel;
3268  color=65536L;
3269  count=0;
3270  do
3271  {
3272  color=list->nodes[color].next[0];
3273  count+=list->nodes[color].count;
3274  } while (count <= (ssize_t) (pixel_list->length >> 1));
3275  channels[channel]=(unsigned short) color;
3276  }
3277  GetMagickPixelPacket((const Image *) NULL,pixel);
3278  pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3279  pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3280  pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3281  pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3282  pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3283 }
3284 
3285 static void GetMinimumPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3286 {
3287  SkipList
3288  *list;
3289 
3290  ssize_t
3291  channel;
3292 
3293  size_t
3294  color,
3295  minimum;
3296 
3297  ssize_t
3298  count;
3299 
3300  unsigned short
3301  channels[ListChannels];
3302 
3303  /*
3304  Find the minimum value for each of the color.
3305  */
3306  for (channel=0; channel < 5; channel++)
3307  {
3308  list=pixel_list->lists+channel;
3309  count=0;
3310  color=65536UL;
3311  minimum=list->nodes[color].next[0];
3312  do
3313  {
3314  color=list->nodes[color].next[0];
3315  if (color < minimum)
3316  minimum=color;
3317  count+=list->nodes[color].count;
3318  } while (count < (ssize_t) pixel_list->length);
3319  channels[channel]=(unsigned short) minimum;
3320  }
3321  pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3322  pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3323  pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3324  pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3325  pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3326 }
3327 
3328 static void GetModePixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3329 {
3330  SkipList
3331  *list;
3332 
3333  ssize_t
3334  channel;
3335 
3336  size_t
3337  color,
3338  max_count,
3339  mode;
3340 
3341  ssize_t
3342  count;
3343 
3344  unsigned short
3345  channels[5];
3346 
3347  /*
3348  Make each pixel the 'predominant color' of the specified neighborhood.
3349  */
3350  for (channel=0; channel < 5; channel++)
3351  {
3352  list=pixel_list->lists+channel;
3353  color=65536L;
3354  mode=color;
3355  max_count=list->nodes[mode].count;
3356  count=0;
3357  do
3358  {
3359  color=list->nodes[color].next[0];
3360  if (list->nodes[color].count > max_count)
3361  {
3362  mode=color;
3363  max_count=list->nodes[mode].count;
3364  }
3365  count+=list->nodes[color].count;
3366  } while (count < (ssize_t) pixel_list->length);
3367  channels[channel]=(unsigned short) mode;
3368  }
3369  pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3370  pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3371  pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3372  pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3373  pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3374 }
3375 
3376 static void GetNonpeakPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3377 {
3378  SkipList
3379  *list;
3380 
3381  ssize_t
3382  channel;
3383 
3384  size_t
3385  color,
3386  next,
3387  previous;
3388 
3389  ssize_t
3390  count;
3391 
3392  unsigned short
3393  channels[5];
3394 
3395  /*
3396  Finds the non peak value for each of the colors.
3397  */
3398  for (channel=0; channel < 5; channel++)
3399  {
3400  list=pixel_list->lists+channel;
3401  color=65536L;
3402  next=list->nodes[color].next[0];
3403  count=0;
3404  do
3405  {
3406  previous=color;
3407  color=next;
3408  next=list->nodes[color].next[0];
3409  count+=list->nodes[color].count;
3410  } while (count <= (ssize_t) (pixel_list->length >> 1));
3411  if ((previous == 65536UL) && (next != 65536UL))
3412  color=next;
3413  else
3414  if ((previous != 65536UL) && (next == 65536UL))
3415  color=previous;
3416  channels[channel]=(unsigned short) color;
3417  }
3418  pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3419  pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3420  pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3421  pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3422  pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3423 }
3424 
3425 static void GetRootMeanSquarePixelList(PixelList *pixel_list,
3426  MagickPixelPacket *pixel)
3427 {
3428  MagickRealType
3429  sum;
3430 
3431  SkipList
3432  *list;
3433 
3434  ssize_t
3435  channel;
3436 
3437  size_t
3438  color;
3439 
3440  ssize_t
3441  count;
3442 
3443  unsigned short
3444  channels[ListChannels];
3445 
3446  /*
3447  Find the root mean square value for each of the color.
3448  */
3449  for (channel=0; channel < 5; channel++)
3450  {
3451  list=pixel_list->lists+channel;
3452  color=65536L;
3453  count=0;
3454  sum=0.0;
3455  do
3456  {
3457  color=list->nodes[color].next[0];
3458  sum+=(MagickRealType) (list->nodes[color].count*color*color);
3459  count+=list->nodes[color].count;
3460  } while (count < (ssize_t) pixel_list->length);
3461  sum/=pixel_list->length;
3462  channels[channel]=(unsigned short) sqrt(sum);
3463  }
3464  pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3465  pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3466  pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3467  pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3468  pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3469 }
3470 
3471 static void GetStandardDeviationPixelList(PixelList *pixel_list,
3472  MagickPixelPacket *pixel)
3473 {
3474  MagickRealType
3475  sum,
3476  sum_squared;
3477 
3478  SkipList
3479  *list;
3480 
3481  ssize_t
3482  channel;
3483 
3484  size_t
3485  color;
3486 
3487  ssize_t
3488  count;
3489 
3490  unsigned short
3491  channels[ListChannels];
3492 
3493  /*
3494  Find the standard-deviation value for each of the color.
3495  */
3496  for (channel=0; channel < 5; channel++)
3497  {
3498  list=pixel_list->lists+channel;
3499  color=65536L;
3500  count=0;
3501  sum=0.0;
3502  sum_squared=0.0;
3503  do
3504  {
3505  ssize_t
3506  i;
3507 
3508  color=list->nodes[color].next[0];
3509  sum+=(MagickRealType) list->nodes[color].count*color;
3510  for (i=0; i < (ssize_t) list->nodes[color].count; i++)
3511  sum_squared+=((MagickRealType) color)*((MagickRealType) color);
3512  count+=list->nodes[color].count;
3513  } while (count < (ssize_t) pixel_list->length);
3514  sum/=pixel_list->length;
3515  sum_squared/=pixel_list->length;
3516  channels[channel]=(unsigned short) sqrt(sum_squared-(sum*sum));
3517  }
3518  pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3519  pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3520  pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3521  pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3522  pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3523 }
3524 
3525 static inline void InsertPixelList(const Image *image,const PixelPacket *pixel,
3526  const IndexPacket *indexes,PixelList *pixel_list)
3527 {
3528  size_t
3529  signature;
3530 
3531  unsigned short
3532  index;
3533 
3534  index=ScaleQuantumToShort(GetPixelRed(pixel));
3535  signature=pixel_list->lists[0].nodes[index].signature;
3536  if (signature == pixel_list->signature)
3537  pixel_list->lists[0].nodes[index].count++;
3538  else
3539  AddNodePixelList(pixel_list,0,index);
3540  index=ScaleQuantumToShort(GetPixelGreen(pixel));
3541  signature=pixel_list->lists[1].nodes[index].signature;
3542  if (signature == pixel_list->signature)
3543  pixel_list->lists[1].nodes[index].count++;
3544  else
3545  AddNodePixelList(pixel_list,1,index);
3546  index=ScaleQuantumToShort(GetPixelBlue(pixel));
3547  signature=pixel_list->lists[2].nodes[index].signature;
3548  if (signature == pixel_list->signature)
3549  pixel_list->lists[2].nodes[index].count++;
3550  else
3551  AddNodePixelList(pixel_list,2,index);
3552  index=ScaleQuantumToShort(GetPixelOpacity(pixel));
3553  signature=pixel_list->lists[3].nodes[index].signature;
3554  if (signature == pixel_list->signature)
3555  pixel_list->lists[3].nodes[index].count++;
3556  else
3557  AddNodePixelList(pixel_list,3,index);
3558  if (image->colorspace == CMYKColorspace)
3559  index=ScaleQuantumToShort(GetPixelIndex(indexes));
3560  signature=pixel_list->lists[4].nodes[index].signature;
3561  if (signature == pixel_list->signature)
3562  pixel_list->lists[4].nodes[index].count++;
3563  else
3564  AddNodePixelList(pixel_list,4,index);
3565 }
3566 
3567 static void ResetPixelList(PixelList *pixel_list)
3568 {
3569  int
3570  level;
3571 
3572  ListNode
3573  *root;
3574 
3575  SkipList
3576  *list;
3577 
3578  ssize_t
3579  channel;
3580 
3581  /*
3582  Reset the skip-list.
3583  */
3584  for (channel=0; channel < 5; channel++)
3585  {
3586  list=pixel_list->lists+channel;
3587  root=list->nodes+65536UL;
3588  list->level=0;
3589  for (level=0; level < 9; level++)
3590  root->next[level]=65536UL;
3591  }
3592  pixel_list->seed=pixel_list->signature++;
3593 }
3594 
3595 MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
3596  const size_t width,const size_t height,ExceptionInfo *exception)
3597 {
3598  Image
3599  *statistic_image;
3600 
3601  statistic_image=StatisticImageChannel(image,DefaultChannels,type,width,
3602  height,exception);
3603  return(statistic_image);
3604 }
3605 
3606 MagickExport Image *StatisticImageChannel(const Image *image,
3607  const ChannelType channel,const StatisticType type,const size_t width,
3608  const size_t height,ExceptionInfo *exception)
3609 {
3610 #define StatisticImageTag "Statistic/Image"
3611 
3612  CacheView
3613  *image_view,
3614  *statistic_view;
3615 
3616  Image
3617  *statistic_image;
3618 
3619  MagickBooleanType
3620  status;
3621 
3622  MagickOffsetType
3623  progress;
3624 
3625  PixelList
3626  **magick_restrict pixel_list;
3627 
3628  size_t
3629  neighbor_height,
3630  neighbor_width;
3631 
3632  ssize_t
3633  y;
3634 
3635  /*
3636  Initialize statistics image attributes.
3637  */
3638  assert(image != (Image *) NULL);
3639  assert(image->signature == MagickCoreSignature);
3640  if (IsEventLogging() != MagickFalse)
3641  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3642  assert(exception != (ExceptionInfo *) NULL);
3643  assert(exception->signature == MagickCoreSignature);
3644  statistic_image=CloneImage(image,0,0,MagickTrue,exception);
3645  if (statistic_image == (Image *) NULL)
3646  return((Image *) NULL);
3647  if (SetImageStorageClass(statistic_image,DirectClass) == MagickFalse)
3648  {
3649  InheritException(exception,&statistic_image->exception);
3650  statistic_image=DestroyImage(statistic_image);
3651  return((Image *) NULL);
3652  }
3653  neighbor_width=width == 0 ? GetOptimalKernelWidth2D((double) width,0.5) :
3654  width;
3655  neighbor_height=height == 0 ? GetOptimalKernelWidth2D((double) height,0.5) :
3656  height;
3657  pixel_list=AcquirePixelListTLS(neighbor_width,neighbor_height);
3658  if (pixel_list == (PixelList **) NULL)
3659  {
3660  statistic_image=DestroyImage(statistic_image);
3661  ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3662  }
3663  /*
3664  Make each pixel the min / max / median / mode / etc. of the neighborhood.
3665  */
3666  status=MagickTrue;
3667  progress=0;
3668  image_view=AcquireVirtualCacheView(image,exception);
3669  statistic_view=AcquireAuthenticCacheView(statistic_image,exception);
3670 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3671  #pragma omp parallel for schedule(static) shared(progress,status) \
3672  magick_number_threads(image,statistic_image,statistic_image->rows,1)
3673 #endif
3674  for (y=0; y < (ssize_t) statistic_image->rows; y++)
3675  {
3676  const int
3677  id = GetOpenMPThreadId();
3678 
3679  const IndexPacket
3680  *magick_restrict indexes;
3681 
3682  const PixelPacket
3683  *magick_restrict p;
3684 
3685  IndexPacket
3686  *magick_restrict statistic_indexes;
3687 
3688  PixelPacket
3689  *magick_restrict q;
3690 
3691  ssize_t
3692  x;
3693 
3694  if (status == MagickFalse)
3695  continue;
3696  p=GetCacheViewVirtualPixels(image_view,-((ssize_t) neighbor_width/2L),y-
3697  (ssize_t) (neighbor_height/2L),image->columns+neighbor_width,
3698  neighbor_height,exception);
3699  q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns, 1,exception);
3700  if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
3701  {
3702  status=MagickFalse;
3703  continue;
3704  }
3705  indexes=GetCacheViewVirtualIndexQueue(image_view);
3706  statistic_indexes=GetCacheViewAuthenticIndexQueue(statistic_view);
3707  for (x=0; x < (ssize_t) statistic_image->columns; x++)
3708  {
3710  pixel;
3711 
3712  const IndexPacket
3713  *magick_restrict s;
3714 
3715  const PixelPacket
3716  *magick_restrict r;
3717 
3718  ssize_t
3719  u,
3720  v;
3721 
3722  r=p;
3723  s=indexes+x;
3724  ResetPixelList(pixel_list[id]);
3725  for (v=0; v < (ssize_t) neighbor_height; v++)
3726  {
3727  for (u=0; u < (ssize_t) neighbor_width; u++)
3728  InsertPixelList(image,r+u,s+u,pixel_list[id]);
3729  r+=image->columns+neighbor_width;
3730  s+=image->columns+neighbor_width;
3731  }
3732  GetMagickPixelPacket(image,&pixel);
3733  SetMagickPixelPacket(image,p+neighbor_width*neighbor_height/2,indexes+x+
3734  neighbor_width*neighbor_height/2,&pixel);
3735  switch (type)
3736  {
3737  case GradientStatistic:
3738  {
3740  maximum,
3741  minimum;
3742 
3743  GetMinimumPixelList(pixel_list[id],&pixel);
3744  minimum=pixel;
3745  GetMaximumPixelList(pixel_list[id],&pixel);
3746  maximum=pixel;
3747  pixel.red=MagickAbsoluteValue(maximum.red-minimum.red);
3748  pixel.green=MagickAbsoluteValue(maximum.green-minimum.green);
3749  pixel.blue=MagickAbsoluteValue(maximum.blue-minimum.blue);
3750  pixel.opacity=MagickAbsoluteValue(maximum.opacity-minimum.opacity);
3751  if (image->colorspace == CMYKColorspace)
3752  pixel.index=MagickAbsoluteValue(maximum.index-minimum.index);
3753  break;
3754  }
3755  case MaximumStatistic:
3756  {
3757  GetMaximumPixelList(pixel_list[id],&pixel);
3758  break;
3759  }
3760  case MeanStatistic:
3761  {
3762  GetMeanPixelList(pixel_list[id],&pixel);
3763  break;
3764  }
3765  case MedianStatistic:
3766  default:
3767  {
3768  GetMedianPixelList(pixel_list[id],&pixel);
3769  break;
3770  }
3771  case MinimumStatistic:
3772  {
3773  GetMinimumPixelList(pixel_list[id],&pixel);
3774  break;
3775  }
3776  case ModeStatistic:
3777  {
3778  GetModePixelList(pixel_list[id],&pixel);
3779  break;
3780  }
3781  case NonpeakStatistic:
3782  {
3783  GetNonpeakPixelList(pixel_list[id],&pixel);
3784  break;
3785  }
3786  case RootMeanSquareStatistic:
3787  {
3788  GetRootMeanSquarePixelList(pixel_list[id],&pixel);
3789  break;
3790  }
3791  case StandardDeviationStatistic:
3792  {
3793  GetStandardDeviationPixelList(pixel_list[id],&pixel);
3794  break;
3795  }
3796  }
3797  if ((channel & RedChannel) != 0)
3798  SetPixelRed(q,ClampToQuantum(pixel.red));
3799  if ((channel & GreenChannel) != 0)
3800  SetPixelGreen(q,ClampToQuantum(pixel.green));
3801  if ((channel & BlueChannel) != 0)
3802  SetPixelBlue(q,ClampToQuantum(pixel.blue));
3803  if ((channel & OpacityChannel) != 0)
3804  SetPixelOpacity(q,ClampToQuantum(pixel.opacity));
3805  if (((channel & IndexChannel) != 0) &&
3806  (image->colorspace == CMYKColorspace))
3807  SetPixelIndex(statistic_indexes+x,ClampToQuantum(pixel.index));
3808  p++;
3809  q++;
3810  }
3811  if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
3812  status=MagickFalse;
3813  if (image->progress_monitor != (MagickProgressMonitor) NULL)
3814  {
3815  MagickBooleanType
3816  proceed;
3817 
3818  proceed=SetImageProgress(image,StatisticImageTag,progress++,
3819  image->rows);
3820  if (proceed == MagickFalse)
3821  status=MagickFalse;
3822  }
3823  }
3824  statistic_view=DestroyCacheView(statistic_view);
3825  image_view=DestroyCacheView(image_view);
3826  pixel_list=DestroyPixelListTLS(pixel_list);
3827  if (status == MagickFalse)
3828  statistic_image=DestroyImage(statistic_image);
3829  return(statistic_image);
3830 }
Definition: image.h:152