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