MagickCore  6.9.13-17
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  {
1604  *mean=NAN;
1605  *standard_deviation=NAN;
1606  return(MagickFalse);
1607  }
1608  channels=0;
1609  channel_statistics[CompositeChannels].mean=0.0;
1610  channel_statistics[CompositeChannels].standard_deviation=0.0;
1611  if ((channel & RedChannel) != 0)
1612  {
1613  channel_statistics[CompositeChannels].mean+=
1614  channel_statistics[RedChannel].mean;
1615  channel_statistics[CompositeChannels].standard_deviation+=
1616  channel_statistics[RedChannel].standard_deviation;
1617  channels++;
1618  }
1619  if ((channel & GreenChannel) != 0)
1620  {
1621  channel_statistics[CompositeChannels].mean+=
1622  channel_statistics[GreenChannel].mean;
1623  channel_statistics[CompositeChannels].standard_deviation+=
1624  channel_statistics[GreenChannel].standard_deviation;
1625  channels++;
1626  }
1627  if ((channel & BlueChannel) != 0)
1628  {
1629  channel_statistics[CompositeChannels].mean+=
1630  channel_statistics[BlueChannel].mean;
1631  channel_statistics[CompositeChannels].standard_deviation+=
1632  channel_statistics[BlueChannel].standard_deviation;
1633  channels++;
1634  }
1635  if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
1636  {
1637  channel_statistics[CompositeChannels].mean+=
1638  channel_statistics[OpacityChannel].mean;
1639  channel_statistics[CompositeChannels].standard_deviation+=
1640  channel_statistics[OpacityChannel].standard_deviation;
1641  channels++;
1642  }
1643  if (((channel & IndexChannel) != 0) && (image->colorspace == CMYKColorspace))
1644  {
1645  channel_statistics[CompositeChannels].mean+=
1646  channel_statistics[BlackChannel].mean;
1647  channel_statistics[CompositeChannels].standard_deviation+=
1648  channel_statistics[CompositeChannels].standard_deviation;
1649  channels++;
1650  }
1651  channel_statistics[CompositeChannels].mean/=channels;
1652  channel_statistics[CompositeChannels].standard_deviation/=channels;
1653  *mean=channel_statistics[CompositeChannels].mean;
1654  *standard_deviation=channel_statistics[CompositeChannels].standard_deviation;
1655  channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1656  channel_statistics);
1657  return(MagickTrue);
1658 }
1659 
1660 /*
1661 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1662 % %
1663 % %
1664 % %
1665 % G e t I m a g e C h a n n e l M o m e n t s %
1666 % %
1667 % %
1668 % %
1669 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1670 %
1671 % GetImageChannelMoments() returns the normalized moments of one or more image
1672 % channels.
1673 %
1674 % The format of the GetImageChannelMoments method is:
1675 %
1676 % ChannelMoments *GetImageChannelMoments(const Image *image,
1677 % ExceptionInfo *exception)
1678 %
1679 % A description of each parameter follows:
1680 %
1681 % o image: the image.
1682 %
1683 % o exception: return any errors or warnings in this structure.
1684 %
1685 */
1686 MagickExport ChannelMoments *GetImageChannelMoments(const Image *image,
1687  ExceptionInfo *exception)
1688 {
1689 #define MaxNumberImageMoments 8
1690 
1692  *channel_moments;
1693 
1694  double
1695  M00[CompositeChannels+1],
1696  M01[CompositeChannels+1],
1697  M02[CompositeChannels+1],
1698  M03[CompositeChannels+1],
1699  M10[CompositeChannels+1],
1700  M11[CompositeChannels+1],
1701  M12[CompositeChannels+1],
1702  M20[CompositeChannels+1],
1703  M21[CompositeChannels+1],
1704  M22[CompositeChannels+1],
1705  M30[CompositeChannels+1];
1706 
1708  pixel;
1709 
1710  PointInfo
1711  centroid[CompositeChannels+1];
1712 
1713  ssize_t
1714  channel,
1715  channels,
1716  y;
1717 
1718  size_t
1719  length;
1720 
1721  assert(image != (Image *) NULL);
1722  assert(image->signature == MagickCoreSignature);
1723  if (IsEventLogging() != MagickFalse)
1724  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1725  length=CompositeChannels+1UL;
1726  channel_moments=(ChannelMoments *) AcquireQuantumMemory(length,
1727  sizeof(*channel_moments));
1728  if (channel_moments == (ChannelMoments *) NULL)
1729  return(channel_moments);
1730  (void) memset(channel_moments,0,length*sizeof(*channel_moments));
1731  (void) memset(centroid,0,sizeof(centroid));
1732  (void) memset(M00,0,sizeof(M00));
1733  (void) memset(M01,0,sizeof(M01));
1734  (void) memset(M02,0,sizeof(M02));
1735  (void) memset(M03,0,sizeof(M03));
1736  (void) memset(M10,0,sizeof(M10));
1737  (void) memset(M11,0,sizeof(M11));
1738  (void) memset(M12,0,sizeof(M12));
1739  (void) memset(M20,0,sizeof(M20));
1740  (void) memset(M21,0,sizeof(M21));
1741  (void) memset(M22,0,sizeof(M22));
1742  (void) memset(M30,0,sizeof(M30));
1743  GetMagickPixelPacket(image,&pixel);
1744  for (y=0; y < (ssize_t) image->rows; y++)
1745  {
1746  const IndexPacket
1747  *magick_restrict indexes;
1748 
1749  const PixelPacket
1750  *magick_restrict p;
1751 
1752  ssize_t
1753  x;
1754 
1755  /*
1756  Compute center of mass (centroid).
1757  */
1758  p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1759  if (p == (const PixelPacket *) NULL)
1760  break;
1761  indexes=GetVirtualIndexQueue(image);
1762  for (x=0; x < (ssize_t) image->columns; x++)
1763  {
1764  SetMagickPixelPacket(image,p,indexes+x,&pixel);
1765  M00[RedChannel]+=QuantumScale*pixel.red;
1766  M10[RedChannel]+=x*QuantumScale*pixel.red;
1767  M01[RedChannel]+=y*QuantumScale*pixel.red;
1768  M00[GreenChannel]+=QuantumScale*pixel.green;
1769  M10[GreenChannel]+=x*QuantumScale*pixel.green;
1770  M01[GreenChannel]+=y*QuantumScale*pixel.green;
1771  M00[BlueChannel]+=QuantumScale*pixel.blue;
1772  M10[BlueChannel]+=x*QuantumScale*pixel.blue;
1773  M01[BlueChannel]+=y*QuantumScale*pixel.blue;
1774  if (image->matte != MagickFalse)
1775  {
1776  M00[OpacityChannel]+=QuantumScale*pixel.opacity;
1777  M10[OpacityChannel]+=x*QuantumScale*pixel.opacity;
1778  M01[OpacityChannel]+=y*QuantumScale*pixel.opacity;
1779  }
1780  if (image->colorspace == CMYKColorspace)
1781  {
1782  M00[IndexChannel]+=QuantumScale*pixel.index;
1783  M10[IndexChannel]+=x*QuantumScale*pixel.index;
1784  M01[IndexChannel]+=y*QuantumScale*pixel.index;
1785  }
1786  p++;
1787  }
1788  }
1789  for (channel=0; channel <= CompositeChannels; channel++)
1790  {
1791  /*
1792  Compute center of mass (centroid).
1793  */
1794  if (M00[channel] < MagickEpsilon)
1795  {
1796  M00[channel]+=MagickEpsilon;
1797  centroid[channel].x=(double) image->columns/2.0;
1798  centroid[channel].y=(double) image->rows/2.0;
1799  continue;
1800  }
1801  M00[channel]+=MagickEpsilon;
1802  centroid[channel].x=M10[channel]/M00[channel];
1803  centroid[channel].y=M01[channel]/M00[channel];
1804  }
1805  for (y=0; y < (ssize_t) image->rows; y++)
1806  {
1807  const IndexPacket
1808  *magick_restrict indexes;
1809 
1810  const PixelPacket
1811  *magick_restrict p;
1812 
1813  ssize_t
1814  x;
1815 
1816  /*
1817  Compute the image moments.
1818  */
1819  p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1820  if (p == (const PixelPacket *) NULL)
1821  break;
1822  indexes=GetVirtualIndexQueue(image);
1823  for (x=0; x < (ssize_t) image->columns; x++)
1824  {
1825  SetMagickPixelPacket(image,p,indexes+x,&pixel);
1826  M11[RedChannel]+=(x-centroid[RedChannel].x)*(y-
1827  centroid[RedChannel].y)*QuantumScale*pixel.red;
1828  M20[RedChannel]+=(x-centroid[RedChannel].x)*(x-
1829  centroid[RedChannel].x)*QuantumScale*pixel.red;
1830  M02[RedChannel]+=(y-centroid[RedChannel].y)*(y-
1831  centroid[RedChannel].y)*QuantumScale*pixel.red;
1832  M21[RedChannel]+=(x-centroid[RedChannel].x)*(x-
1833  centroid[RedChannel].x)*(y-centroid[RedChannel].y)*QuantumScale*
1834  pixel.red;
1835  M12[RedChannel]+=(x-centroid[RedChannel].x)*(y-
1836  centroid[RedChannel].y)*(y-centroid[RedChannel].y)*QuantumScale*
1837  pixel.red;
1838  M22[RedChannel]+=(x-centroid[RedChannel].x)*(x-
1839  centroid[RedChannel].x)*(y-centroid[RedChannel].y)*(y-
1840  centroid[RedChannel].y)*QuantumScale*pixel.red;
1841  M30[RedChannel]+=(x-centroid[RedChannel].x)*(x-
1842  centroid[RedChannel].x)*(x-centroid[RedChannel].x)*QuantumScale*
1843  pixel.red;
1844  M03[RedChannel]+=(y-centroid[RedChannel].y)*(y-
1845  centroid[RedChannel].y)*(y-centroid[RedChannel].y)*QuantumScale*
1846  pixel.red;
1847  M11[GreenChannel]+=(x-centroid[GreenChannel].x)*(y-
1848  centroid[GreenChannel].y)*QuantumScale*pixel.green;
1849  M20[GreenChannel]+=(x-centroid[GreenChannel].x)*(x-
1850  centroid[GreenChannel].x)*QuantumScale*pixel.green;
1851  M02[GreenChannel]+=(y-centroid[GreenChannel].y)*(y-
1852  centroid[GreenChannel].y)*QuantumScale*pixel.green;
1853  M21[GreenChannel]+=(x-centroid[GreenChannel].x)*(x-
1854  centroid[GreenChannel].x)*(y-centroid[GreenChannel].y)*QuantumScale*
1855  pixel.green;
1856  M12[GreenChannel]+=(x-centroid[GreenChannel].x)*(y-
1857  centroid[GreenChannel].y)*(y-centroid[GreenChannel].y)*QuantumScale*
1858  pixel.green;
1859  M22[GreenChannel]+=(x-centroid[GreenChannel].x)*(x-
1860  centroid[GreenChannel].x)*(y-centroid[GreenChannel].y)*(y-
1861  centroid[GreenChannel].y)*QuantumScale*pixel.green;
1862  M30[GreenChannel]+=(x-centroid[GreenChannel].x)*(x-
1863  centroid[GreenChannel].x)*(x-centroid[GreenChannel].x)*QuantumScale*
1864  pixel.green;
1865  M03[GreenChannel]+=(y-centroid[GreenChannel].y)*(y-
1866  centroid[GreenChannel].y)*(y-centroid[GreenChannel].y)*QuantumScale*
1867  pixel.green;
1868  M11[BlueChannel]+=(x-centroid[BlueChannel].x)*(y-
1869  centroid[BlueChannel].y)*QuantumScale*pixel.blue;
1870  M20[BlueChannel]+=(x-centroid[BlueChannel].x)*(x-
1871  centroid[BlueChannel].x)*QuantumScale*pixel.blue;
1872  M02[BlueChannel]+=(y-centroid[BlueChannel].y)*(y-
1873  centroid[BlueChannel].y)*QuantumScale*pixel.blue;
1874  M21[BlueChannel]+=(x-centroid[BlueChannel].x)*(x-
1875  centroid[BlueChannel].x)*(y-centroid[BlueChannel].y)*QuantumScale*
1876  pixel.blue;
1877  M12[BlueChannel]+=(x-centroid[BlueChannel].x)*(y-
1878  centroid[BlueChannel].y)*(y-centroid[BlueChannel].y)*QuantumScale*
1879  pixel.blue;
1880  M22[BlueChannel]+=(x-centroid[BlueChannel].x)*(x-
1881  centroid[BlueChannel].x)*(y-centroid[BlueChannel].y)*(y-
1882  centroid[BlueChannel].y)*QuantumScale*pixel.blue;
1883  M30[BlueChannel]+=(x-centroid[BlueChannel].x)*(x-
1884  centroid[BlueChannel].x)*(x-centroid[BlueChannel].x)*QuantumScale*
1885  pixel.blue;
1886  M03[BlueChannel]+=(y-centroid[BlueChannel].y)*(y-
1887  centroid[BlueChannel].y)*(y-centroid[BlueChannel].y)*QuantumScale*
1888  pixel.blue;
1889  if (image->matte != MagickFalse)
1890  {
1891  M11[OpacityChannel]+=(x-centroid[OpacityChannel].x)*(y-
1892  centroid[OpacityChannel].y)*QuantumScale*pixel.opacity;
1893  M20[OpacityChannel]+=(x-centroid[OpacityChannel].x)*(x-
1894  centroid[OpacityChannel].x)*QuantumScale*pixel.opacity;
1895  M02[OpacityChannel]+=(y-centroid[OpacityChannel].y)*(y-
1896  centroid[OpacityChannel].y)*QuantumScale*pixel.opacity;
1897  M21[OpacityChannel]+=(x-centroid[OpacityChannel].x)*(x-
1898  centroid[OpacityChannel].x)*(y-centroid[OpacityChannel].y)*
1899  QuantumScale*pixel.opacity;
1900  M12[OpacityChannel]+=(x-centroid[OpacityChannel].x)*(y-
1901  centroid[OpacityChannel].y)*(y-centroid[OpacityChannel].y)*
1902  QuantumScale*pixel.opacity;
1903  M22[OpacityChannel]+=(x-centroid[OpacityChannel].x)*(x-
1904  centroid[OpacityChannel].x)*(y-centroid[OpacityChannel].y)*(y-
1905  centroid[OpacityChannel].y)*QuantumScale*pixel.opacity;
1906  M30[OpacityChannel]+=(x-centroid[OpacityChannel].x)*(x-
1907  centroid[OpacityChannel].x)*(x-centroid[OpacityChannel].x)*
1908  QuantumScale*pixel.opacity;
1909  M03[OpacityChannel]+=(y-centroid[OpacityChannel].y)*(y-
1910  centroid[OpacityChannel].y)*(y-centroid[OpacityChannel].y)*
1911  QuantumScale*pixel.opacity;
1912  }
1913  if (image->colorspace == CMYKColorspace)
1914  {
1915  M11[IndexChannel]+=(x-centroid[IndexChannel].x)*(y-
1916  centroid[IndexChannel].y)*QuantumScale*pixel.index;
1917  M20[IndexChannel]+=(x-centroid[IndexChannel].x)*(x-
1918  centroid[IndexChannel].x)*QuantumScale*pixel.index;
1919  M02[IndexChannel]+=(y-centroid[IndexChannel].y)*(y-
1920  centroid[IndexChannel].y)*QuantumScale*pixel.index;
1921  M21[IndexChannel]+=(x-centroid[IndexChannel].x)*(x-
1922  centroid[IndexChannel].x)*(y-centroid[IndexChannel].y)*
1923  QuantumScale*pixel.index;
1924  M12[IndexChannel]+=(x-centroid[IndexChannel].x)*(y-
1925  centroid[IndexChannel].y)*(y-centroid[IndexChannel].y)*
1926  QuantumScale*pixel.index;
1927  M22[IndexChannel]+=(x-centroid[IndexChannel].x)*(x-
1928  centroid[IndexChannel].x)*(y-centroid[IndexChannel].y)*(y-
1929  centroid[IndexChannel].y)*QuantumScale*pixel.index;
1930  M30[IndexChannel]+=(x-centroid[IndexChannel].x)*(x-
1931  centroid[IndexChannel].x)*(x-centroid[IndexChannel].x)*
1932  QuantumScale*pixel.index;
1933  M03[IndexChannel]+=(y-centroid[IndexChannel].y)*(y-
1934  centroid[IndexChannel].y)*(y-centroid[IndexChannel].y)*
1935  QuantumScale*pixel.index;
1936  }
1937  p++;
1938  }
1939  }
1940  channels=3;
1941  M00[CompositeChannels]+=(M00[RedChannel]+M00[GreenChannel]+M00[BlueChannel]);
1942  M01[CompositeChannels]+=(M01[RedChannel]+M01[GreenChannel]+M01[BlueChannel]);
1943  M02[CompositeChannels]+=(M02[RedChannel]+M02[GreenChannel]+M02[BlueChannel]);
1944  M03[CompositeChannels]+=(M03[RedChannel]+M03[GreenChannel]+M03[BlueChannel]);
1945  M10[CompositeChannels]+=(M10[RedChannel]+M10[GreenChannel]+M10[BlueChannel]);
1946  M11[CompositeChannels]+=(M11[RedChannel]+M11[GreenChannel]+M11[BlueChannel]);
1947  M12[CompositeChannels]+=(M12[RedChannel]+M12[GreenChannel]+M12[BlueChannel]);
1948  M20[CompositeChannels]+=(M20[RedChannel]+M20[GreenChannel]+M20[BlueChannel]);
1949  M21[CompositeChannels]+=(M21[RedChannel]+M21[GreenChannel]+M21[BlueChannel]);
1950  M22[CompositeChannels]+=(M22[RedChannel]+M22[GreenChannel]+M22[BlueChannel]);
1951  M30[CompositeChannels]+=(M30[RedChannel]+M30[GreenChannel]+M30[BlueChannel]);
1952  if (image->matte != MagickFalse)
1953  {
1954  channels+=1;
1955  M00[CompositeChannels]+=M00[OpacityChannel];
1956  M01[CompositeChannels]+=M01[OpacityChannel];
1957  M02[CompositeChannels]+=M02[OpacityChannel];
1958  M03[CompositeChannels]+=M03[OpacityChannel];
1959  M10[CompositeChannels]+=M10[OpacityChannel];
1960  M11[CompositeChannels]+=M11[OpacityChannel];
1961  M12[CompositeChannels]+=M12[OpacityChannel];
1962  M20[CompositeChannels]+=M20[OpacityChannel];
1963  M21[CompositeChannels]+=M21[OpacityChannel];
1964  M22[CompositeChannels]+=M22[OpacityChannel];
1965  M30[CompositeChannels]+=M30[OpacityChannel];
1966  }
1967  if (image->colorspace == CMYKColorspace)
1968  {
1969  channels+=1;
1970  M00[CompositeChannels]+=M00[IndexChannel];
1971  M01[CompositeChannels]+=M01[IndexChannel];
1972  M02[CompositeChannels]+=M02[IndexChannel];
1973  M03[CompositeChannels]+=M03[IndexChannel];
1974  M10[CompositeChannels]+=M10[IndexChannel];
1975  M11[CompositeChannels]+=M11[IndexChannel];
1976  M12[CompositeChannels]+=M12[IndexChannel];
1977  M20[CompositeChannels]+=M20[IndexChannel];
1978  M21[CompositeChannels]+=M21[IndexChannel];
1979  M22[CompositeChannels]+=M22[IndexChannel];
1980  M30[CompositeChannels]+=M30[IndexChannel];
1981  }
1982  M00[CompositeChannels]/=(double) channels;
1983  M01[CompositeChannels]/=(double) channels;
1984  M02[CompositeChannels]/=(double) channels;
1985  M03[CompositeChannels]/=(double) channels;
1986  M10[CompositeChannels]/=(double) channels;
1987  M11[CompositeChannels]/=(double) channels;
1988  M12[CompositeChannels]/=(double) channels;
1989  M20[CompositeChannels]/=(double) channels;
1990  M21[CompositeChannels]/=(double) channels;
1991  M22[CompositeChannels]/=(double) channels;
1992  M30[CompositeChannels]/=(double) channels;
1993  for (channel=0; channel <= CompositeChannels; channel++)
1994  {
1995  /*
1996  Compute elliptical angle, major and minor axes, eccentricity, & intensity.
1997  */
1998  channel_moments[channel].centroid=centroid[channel];
1999  channel_moments[channel].ellipse_axis.x=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_axis.y=sqrt((2.0*
2004  PerceptibleReciprocal(M00[channel]))*((M20[channel]+M02[channel])-
2005  sqrt(4.0*M11[channel]*M11[channel]+(M20[channel]-M02[channel])*
2006  (M20[channel]-M02[channel]))));
2007  channel_moments[channel].ellipse_angle=RadiansToDegrees(1.0/2.0*atan(2.0*
2008  M11[channel]*PerceptibleReciprocal(M20[channel]-M02[channel])));
2009  if (fabs(M11[channel]) < 0.0)
2010  {
2011  if ((fabs(M20[channel]-M02[channel]) >= 0.0) &&
2012  ((M20[channel]-M02[channel]) < 0.0))
2013  channel_moments[channel].ellipse_angle+=90.0;
2014  }
2015  else
2016  if (M11[channel] < 0.0)
2017  {
2018  if (fabs(M20[channel]-M02[channel]) >= 0.0)
2019  {
2020  if ((M20[channel]-M02[channel]) < 0.0)
2021  channel_moments[channel].ellipse_angle+=90.0;
2022  else
2023  channel_moments[channel].ellipse_angle+=180.0;
2024  }
2025  }
2026  else
2027  if ((fabs(M20[channel]-M02[channel]) >= 0.0) &&
2028  ((M20[channel]-M02[channel]) < 0.0))
2029  channel_moments[channel].ellipse_angle+=90.0;
2030  channel_moments[channel].ellipse_eccentricity=sqrt(1.0-(
2031  channel_moments[channel].ellipse_axis.y*
2032  channel_moments[channel].ellipse_axis.y*PerceptibleReciprocal(
2033  channel_moments[channel].ellipse_axis.x*
2034  channel_moments[channel].ellipse_axis.x)));
2035  channel_moments[channel].ellipse_intensity=M00[channel]/
2036  (MagickPI*channel_moments[channel].ellipse_axis.x*
2037  channel_moments[channel].ellipse_axis.y+MagickEpsilon);
2038  }
2039  for (channel=0; channel <= CompositeChannels; channel++)
2040  {
2041  /*
2042  Normalize image moments.
2043  */
2044  M10[channel]=0.0;
2045  M01[channel]=0.0;
2046  M11[channel]/=pow(M00[channel],1.0+(1.0+1.0)/2.0);
2047  M20[channel]/=pow(M00[channel],1.0+(2.0+0.0)/2.0);
2048  M02[channel]/=pow(M00[channel],1.0+(0.0+2.0)/2.0);
2049  M21[channel]/=pow(M00[channel],1.0+(2.0+1.0)/2.0);
2050  M12[channel]/=pow(M00[channel],1.0+(1.0+2.0)/2.0);
2051  M22[channel]/=pow(M00[channel],1.0+(2.0+2.0)/2.0);
2052  M30[channel]/=pow(M00[channel],1.0+(3.0+0.0)/2.0);
2053  M03[channel]/=pow(M00[channel],1.0+(0.0+3.0)/2.0);
2054  M00[channel]=1.0;
2055  }
2056  for (channel=0; channel <= CompositeChannels; channel++)
2057  {
2058  /*
2059  Compute Hu invariant moments.
2060  */
2061  channel_moments[channel].I[0]=M20[channel]+M02[channel];
2062  channel_moments[channel].I[1]=(M20[channel]-M02[channel])*
2063  (M20[channel]-M02[channel])+4.0*M11[channel]*M11[channel];
2064  channel_moments[channel].I[2]=(M30[channel]-3.0*M12[channel])*
2065  (M30[channel]-3.0*M12[channel])+(3.0*M21[channel]-M03[channel])*
2066  (3.0*M21[channel]-M03[channel]);
2067  channel_moments[channel].I[3]=(M30[channel]+M12[channel])*
2068  (M30[channel]+M12[channel])+(M21[channel]+M03[channel])*
2069  (M21[channel]+M03[channel]);
2070  channel_moments[channel].I[4]=(M30[channel]-3.0*M12[channel])*
2071  (M30[channel]+M12[channel])*((M30[channel]+M12[channel])*
2072  (M30[channel]+M12[channel])-3.0*(M21[channel]+M03[channel])*
2073  (M21[channel]+M03[channel]))+(3.0*M21[channel]-M03[channel])*
2074  (M21[channel]+M03[channel])*(3.0*(M30[channel]+M12[channel])*
2075  (M30[channel]+M12[channel])-(M21[channel]+M03[channel])*
2076  (M21[channel]+M03[channel]));
2077  channel_moments[channel].I[5]=(M20[channel]-M02[channel])*
2078  ((M30[channel]+M12[channel])*(M30[channel]+M12[channel])-
2079  (M21[channel]+M03[channel])*(M21[channel]+M03[channel]))+
2080  4.0*M11[channel]*(M30[channel]+M12[channel])*(M21[channel]+M03[channel]);
2081  channel_moments[channel].I[6]=(3.0*M21[channel]-M03[channel])*
2082  (M30[channel]+M12[channel])*((M30[channel]+M12[channel])*
2083  (M30[channel]+M12[channel])-3.0*(M21[channel]+M03[channel])*
2084  (M21[channel]+M03[channel]))-(M30[channel]-3*M12[channel])*
2085  (M21[channel]+M03[channel])*(3.0*(M30[channel]+M12[channel])*
2086  (M30[channel]+M12[channel])-(M21[channel]+M03[channel])*
2087  (M21[channel]+M03[channel]));
2088  channel_moments[channel].I[7]=M11[channel]*((M30[channel]+M12[channel])*
2089  (M30[channel]+M12[channel])-(M03[channel]+M21[channel])*
2090  (M03[channel]+M21[channel]))-(M20[channel]-M02[channel])*
2091  (M30[channel]+M12[channel])*(M03[channel]+M21[channel]);
2092  }
2093  if (y < (ssize_t) image->rows)
2094  channel_moments=(ChannelMoments *) RelinquishMagickMemory(channel_moments);
2095  return(channel_moments);
2096 }
2097 
2098 /*
2099 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2100 % %
2101 % %
2102 % %
2103 % 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 %
2104 % %
2105 % %
2106 % %
2107 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2108 %
2109 % GetImageChannelPerceptualHash() returns the perceptual hash of one or more
2110 % image channels.
2111 %
2112 % The format of the GetImageChannelPerceptualHash method is:
2113 %
2114 % ChannelPerceptualHash *GetImageChannelPerceptualHash(const Image *image,
2115 % ExceptionInfo *exception)
2116 %
2117 % A description of each parameter follows:
2118 %
2119 % o image: the image.
2120 %
2121 % o exception: return any errors or warnings in this structure.
2122 %
2123 */
2124 MagickExport ChannelPerceptualHash *GetImageChannelPerceptualHash(
2125  const Image *image,ExceptionInfo *exception)
2126 {
2128  *moments;
2129 
2131  *perceptual_hash;
2132 
2133  Image
2134  *hash_image;
2135 
2136  MagickBooleanType
2137  status;
2138 
2139  ssize_t
2140  i;
2141 
2142  ssize_t
2143  channel;
2144 
2145  /*
2146  Blur then transform to xyY colorspace.
2147  */
2148  hash_image=BlurImage(image,0.0,1.0,exception);
2149  if (hash_image == (Image *) NULL)
2150  return((ChannelPerceptualHash *) NULL);
2151  hash_image->depth=8;
2152  status=TransformImageColorspace(hash_image,xyYColorspace);
2153  if (status == MagickFalse)
2154  return((ChannelPerceptualHash *) NULL);
2155  moments=GetImageChannelMoments(hash_image,exception);
2156  hash_image=DestroyImage(hash_image);
2157  if (moments == (ChannelMoments *) NULL)
2158  return((ChannelPerceptualHash *) NULL);
2159  perceptual_hash=(ChannelPerceptualHash *) AcquireQuantumMemory(
2160  CompositeChannels+1UL,sizeof(*perceptual_hash));
2161  if (perceptual_hash == (ChannelPerceptualHash *) NULL)
2162  return((ChannelPerceptualHash *) NULL);
2163  for (channel=0; channel <= CompositeChannels; channel++)
2164  for (i=0; i < MaximumNumberOfImageMoments; i++)
2165  perceptual_hash[channel].P[i]=(-MagickLog10(moments[channel].I[i]));
2166  moments=(ChannelMoments *) RelinquishMagickMemory(moments);
2167  /*
2168  Blur then transform to HCLp colorspace.
2169  */
2170  hash_image=BlurImage(image,0.0,1.0,exception);
2171  if (hash_image == (Image *) NULL)
2172  {
2173  perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
2174  perceptual_hash);
2175  return((ChannelPerceptualHash *) NULL);
2176  }
2177  hash_image->depth=8;
2178  status=TransformImageColorspace(hash_image,HSBColorspace);
2179  if (status == MagickFalse)
2180  {
2181  perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
2182  perceptual_hash);
2183  return((ChannelPerceptualHash *) NULL);
2184  }
2185  moments=GetImageChannelMoments(hash_image,exception);
2186  hash_image=DestroyImage(hash_image);
2187  if (moments == (ChannelMoments *) NULL)
2188  {
2189  perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
2190  perceptual_hash);
2191  return((ChannelPerceptualHash *) NULL);
2192  }
2193  for (channel=0; channel <= CompositeChannels; channel++)
2194  for (i=0; i < MaximumNumberOfImageMoments; i++)
2195  perceptual_hash[channel].Q[i]=(-MagickLog10(moments[channel].I[i]));
2196  moments=(ChannelMoments *) RelinquishMagickMemory(moments);
2197  return(perceptual_hash);
2198 }
2199 
2200 /*
2201 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2202 % %
2203 % %
2204 % %
2205 % G e t I m a g e C h a n n e l R a n g e %
2206 % %
2207 % %
2208 % %
2209 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2210 %
2211 % GetImageChannelRange() returns the range of one or more image channels.
2212 %
2213 % The format of the GetImageChannelRange method is:
2214 %
2215 % MagickBooleanType GetImageChannelRange(const Image *image,
2216 % const ChannelType channel,double *minima,double *maxima,
2217 % ExceptionInfo *exception)
2218 %
2219 % A description of each parameter follows:
2220 %
2221 % o image: the image.
2222 %
2223 % o channel: the channel.
2224 %
2225 % o minima: the minimum value in the channel.
2226 %
2227 % o maxima: the maximum value in the channel.
2228 %
2229 % o exception: return any errors or warnings in this structure.
2230 %
2231 */
2232 
2233 MagickExport MagickBooleanType GetImageRange(const Image *image,
2234  double *minima,double *maxima,ExceptionInfo *exception)
2235 {
2236  return(GetImageChannelRange(image,CompositeChannels,minima,maxima,exception));
2237 }
2238 
2239 MagickExport MagickBooleanType GetImageChannelRange(const Image *image,
2240  const ChannelType channel,double *minima,double *maxima,
2241  ExceptionInfo *exception)
2242 {
2244  pixel;
2245 
2246  ssize_t
2247  y;
2248 
2249  assert(image != (Image *) NULL);
2250  assert(image->signature == MagickCoreSignature);
2251  if (IsEventLogging() != MagickFalse)
2252  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2253  *maxima=(-MagickMaximumValue);
2254  *minima=MagickMaximumValue;
2255  GetMagickPixelPacket(image,&pixel);
2256  for (y=0; y < (ssize_t) image->rows; y++)
2257  {
2258  const IndexPacket
2259  *magick_restrict indexes;
2260 
2261  const PixelPacket
2262  *magick_restrict p;
2263 
2264  ssize_t
2265  x;
2266 
2267  p=GetVirtualPixels(image,0,y,image->columns,1,exception);
2268  if (p == (const PixelPacket *) NULL)
2269  break;
2270  indexes=GetVirtualIndexQueue(image);
2271  for (x=0; x < (ssize_t) image->columns; x++)
2272  {
2273  SetMagickPixelPacket(image,p,indexes+x,&pixel);
2274  if ((channel & RedChannel) != 0)
2275  {
2276  if (pixel.red < *minima)
2277  *minima=(double) pixel.red;
2278  if (pixel.red > *maxima)
2279  *maxima=(double) pixel.red;
2280  }
2281  if ((channel & GreenChannel) != 0)
2282  {
2283  if (pixel.green < *minima)
2284  *minima=(double) pixel.green;
2285  if (pixel.green > *maxima)
2286  *maxima=(double) pixel.green;
2287  }
2288  if ((channel & BlueChannel) != 0)
2289  {
2290  if (pixel.blue < *minima)
2291  *minima=(double) pixel.blue;
2292  if (pixel.blue > *maxima)
2293  *maxima=(double) pixel.blue;
2294  }
2295  if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
2296  {
2297  if (((MagickRealType) QuantumRange-(MagickRealType) pixel.opacity) < *minima)
2298  *minima=(double) ((MagickRealType) QuantumRange-(MagickRealType)
2299  pixel.opacity);
2300  if (((MagickRealType) QuantumRange-(MagickRealType) pixel.opacity) > *maxima)
2301  *maxima=(double) ((MagickRealType) QuantumRange-(MagickRealType)
2302  pixel.opacity);
2303  }
2304  if (((channel & IndexChannel) != 0) &&
2305  (image->colorspace == CMYKColorspace))
2306  {
2307  if ((double) pixel.index < *minima)
2308  *minima=(double) pixel.index;
2309  if ((double) pixel.index > *maxima)
2310  *maxima=(double) pixel.index;
2311  }
2312  p++;
2313  }
2314  }
2315  return(y == (ssize_t) image->rows ? MagickTrue : MagickFalse);
2316 }
2317 
2318 /*
2319 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2320 % %
2321 % %
2322 % %
2323 % 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 %
2324 % %
2325 % %
2326 % %
2327 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2328 %
2329 % GetImageChannelStatistics() returns statistics for each channel in the
2330 % image. The statistics include the channel depth, its minima, maxima, mean,
2331 % standard deviation, kurtosis and skewness. You can access the red channel
2332 % mean, for example, like this:
2333 %
2334 % channel_statistics=GetImageChannelStatistics(image,exception);
2335 % red_mean=channel_statistics[RedChannel].mean;
2336 %
2337 % Use MagickRelinquishMemory() to free the statistics buffer.
2338 %
2339 % The format of the GetImageChannelStatistics method is:
2340 %
2341 % ChannelStatistics *GetImageChannelStatistics(const Image *image,
2342 % ExceptionInfo *exception)
2343 %
2344 % A description of each parameter follows:
2345 %
2346 % o image: the image.
2347 %
2348 % o exception: return any errors or warnings in this structure.
2349 %
2350 */
2351 MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
2352  ExceptionInfo *exception)
2353 {
2355  *channel_statistics;
2356 
2357  double
2358  area,
2359  standard_deviation;
2360 
2362  number_bins,
2363  *histogram;
2364 
2365  QuantumAny
2366  range;
2367 
2368  size_t
2369  channels,
2370  depth,
2371  length;
2372 
2373  ssize_t
2374  i,
2375  y;
2376 
2377  assert(image != (Image *) NULL);
2378  assert(image->signature == MagickCoreSignature);
2379  if (IsEventLogging() != MagickFalse)
2380  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2381  length=CompositeChannels+1UL;
2382  channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(length,
2383  sizeof(*channel_statistics));
2384  histogram=(MagickPixelPacket *) AcquireQuantumMemory(MaxMap+1U,
2385  sizeof(*histogram));
2386  if ((channel_statistics == (ChannelStatistics *) NULL) ||
2387  (histogram == (MagickPixelPacket *) NULL))
2388  {
2389  if (histogram != (MagickPixelPacket *) NULL)
2390  histogram=(MagickPixelPacket *) RelinquishMagickMemory(histogram);
2391  if (channel_statistics != (ChannelStatistics *) NULL)
2392  channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
2393  channel_statistics);
2394  return(channel_statistics);
2395  }
2396  (void) memset(channel_statistics,0,length*
2397  sizeof(*channel_statistics));
2398  for (i=0; i <= (ssize_t) CompositeChannels; i++)
2399  {
2400  channel_statistics[i].depth=1;
2401  channel_statistics[i].maxima=(-MagickMaximumValue);
2402  channel_statistics[i].minima=MagickMaximumValue;
2403  }
2404  (void) memset(histogram,0,(MaxMap+1U)*sizeof(*histogram));
2405  (void) memset(&number_bins,0,sizeof(number_bins));
2406  for (y=0; y < (ssize_t) image->rows; y++)
2407  {
2408  const IndexPacket
2409  *magick_restrict indexes;
2410 
2411  const PixelPacket
2412  *magick_restrict p;
2413 
2414  ssize_t
2415  x;
2416 
2417  /*
2418  Compute pixel statistics.
2419  */
2420  p=GetVirtualPixels(image,0,y,image->columns,1,exception);
2421  if (p == (const PixelPacket *) NULL)
2422  break;
2423  indexes=GetVirtualIndexQueue(image);
2424  for (x=0; x < (ssize_t) image->columns; )
2425  {
2426  if (channel_statistics[RedChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2427  {
2428  depth=channel_statistics[RedChannel].depth;
2429  range=GetQuantumRange(depth);
2430  if (IsPixelAtDepth(GetPixelRed(p),range) == MagickFalse)
2431  {
2432  channel_statistics[RedChannel].depth++;
2433  continue;
2434  }
2435  }
2436  if (channel_statistics[GreenChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2437  {
2438  depth=channel_statistics[GreenChannel].depth;
2439  range=GetQuantumRange(depth);
2440  if (IsPixelAtDepth(GetPixelGreen(p),range) == MagickFalse)
2441  {
2442  channel_statistics[GreenChannel].depth++;
2443  continue;
2444  }
2445  }
2446  if (channel_statistics[BlueChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2447  {
2448  depth=channel_statistics[BlueChannel].depth;
2449  range=GetQuantumRange(depth);
2450  if (IsPixelAtDepth(GetPixelBlue(p),range) == MagickFalse)
2451  {
2452  channel_statistics[BlueChannel].depth++;
2453  continue;
2454  }
2455  }
2456  if (image->matte != MagickFalse)
2457  {
2458  if (channel_statistics[OpacityChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2459  {
2460  depth=channel_statistics[OpacityChannel].depth;
2461  range=GetQuantumRange(depth);
2462  if (IsPixelAtDepth(GetPixelAlpha(p),range) == MagickFalse)
2463  {
2464  channel_statistics[OpacityChannel].depth++;
2465  continue;
2466  }
2467  }
2468  }
2469  if (image->colorspace == CMYKColorspace)
2470  {
2471  if (channel_statistics[BlackChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2472  {
2473  depth=channel_statistics[BlackChannel].depth;
2474  range=GetQuantumRange(depth);
2475  if (IsPixelAtDepth(GetPixelIndex(indexes+x),range) == MagickFalse)
2476  {
2477  channel_statistics[BlackChannel].depth++;
2478  continue;
2479  }
2480  }
2481  }
2482  if ((double) GetPixelRed(p) < channel_statistics[RedChannel].minima)
2483  channel_statistics[RedChannel].minima=(double) GetPixelRed(p);
2484  if ((double) GetPixelRed(p) > channel_statistics[RedChannel].maxima)
2485  channel_statistics[RedChannel].maxima=(double) GetPixelRed(p);
2486  channel_statistics[RedChannel].sum+=QuantumScale*GetPixelRed(p);
2487  channel_statistics[RedChannel].sum_squared+=QuantumScale*GetPixelRed(p)*
2488  QuantumScale*GetPixelRed(p);
2489  channel_statistics[RedChannel].sum_cubed+=QuantumScale*GetPixelRed(p)*
2490  QuantumScale*GetPixelRed(p)*QuantumScale*GetPixelRed(p);
2491  channel_statistics[RedChannel].sum_fourth_power+=QuantumScale*
2492  GetPixelRed(p)*QuantumScale*GetPixelRed(p)*QuantumScale*GetPixelRed(p)*
2493  QuantumScale*GetPixelRed(p);
2494  if ((double) GetPixelGreen(p) < channel_statistics[GreenChannel].minima)
2495  channel_statistics[GreenChannel].minima=(double) GetPixelGreen(p);
2496  if ((double) GetPixelGreen(p) > channel_statistics[GreenChannel].maxima)
2497  channel_statistics[GreenChannel].maxima=(double) GetPixelGreen(p);
2498  channel_statistics[GreenChannel].sum+=QuantumScale*GetPixelGreen(p);
2499  channel_statistics[GreenChannel].sum_squared+=QuantumScale*GetPixelGreen(p)*
2500  QuantumScale*GetPixelGreen(p);
2501  channel_statistics[GreenChannel].sum_cubed+=QuantumScale*GetPixelGreen(p)*
2502  QuantumScale*GetPixelGreen(p)*QuantumScale*GetPixelGreen(p);
2503  channel_statistics[GreenChannel].sum_fourth_power+=QuantumScale*
2504  GetPixelGreen(p)*QuantumScale*GetPixelGreen(p)*QuantumScale*
2505  GetPixelGreen(p)*QuantumScale*GetPixelGreen(p);
2506  if ((double) GetPixelBlue(p) < channel_statistics[BlueChannel].minima)
2507  channel_statistics[BlueChannel].minima=(double) GetPixelBlue(p);
2508  if ((double) GetPixelBlue(p) > channel_statistics[BlueChannel].maxima)
2509  channel_statistics[BlueChannel].maxima=(double) GetPixelBlue(p);
2510  channel_statistics[BlueChannel].sum+=QuantumScale*GetPixelBlue(p);
2511  channel_statistics[BlueChannel].sum_squared+=QuantumScale*GetPixelBlue(p)*
2512  QuantumScale*GetPixelBlue(p);
2513  channel_statistics[BlueChannel].sum_cubed+=QuantumScale*GetPixelBlue(p)*
2514  QuantumScale*GetPixelBlue(p)*QuantumScale*GetPixelBlue(p);
2515  channel_statistics[BlueChannel].sum_fourth_power+=QuantumScale*
2516  GetPixelBlue(p)*QuantumScale*GetPixelBlue(p)*QuantumScale*
2517  GetPixelBlue(p)*QuantumScale*GetPixelBlue(p);
2518  histogram[ScaleQuantumToMap(GetPixelRed(p))].red++;
2519  histogram[ScaleQuantumToMap(GetPixelGreen(p))].green++;
2520  histogram[ScaleQuantumToMap(GetPixelBlue(p))].blue++;
2521  if (image->matte != MagickFalse)
2522  {
2523  if ((double) GetPixelAlpha(p) < channel_statistics[OpacityChannel].minima)
2524  channel_statistics[OpacityChannel].minima=(double) GetPixelAlpha(p);
2525  if ((double) GetPixelAlpha(p) > channel_statistics[OpacityChannel].maxima)
2526  channel_statistics[OpacityChannel].maxima=(double) GetPixelAlpha(p);
2527  channel_statistics[OpacityChannel].sum+=QuantumScale*GetPixelAlpha(p);
2528  channel_statistics[OpacityChannel].sum_squared+=QuantumScale*
2529  GetPixelAlpha(p)*QuantumScale*GetPixelAlpha(p);
2530  channel_statistics[OpacityChannel].sum_cubed+=QuantumScale*
2531  GetPixelAlpha(p)*QuantumScale*GetPixelAlpha(p)*QuantumScale*
2532  GetPixelAlpha(p);
2533  channel_statistics[OpacityChannel].sum_fourth_power+=QuantumScale*
2534  GetPixelAlpha(p)*QuantumScale*GetPixelAlpha(p)*QuantumScale*
2535  GetPixelAlpha(p)*QuantumScale*GetPixelAlpha(p);
2536  histogram[ScaleQuantumToMap(GetPixelAlpha(p))].opacity++;
2537  }
2538  if (image->colorspace == CMYKColorspace)
2539  {
2540  if ((double) GetPixelIndex(indexes+x) < channel_statistics[BlackChannel].minima)
2541  channel_statistics[BlackChannel].minima=(double)
2542  GetPixelIndex(indexes+x);
2543  if ((double) GetPixelIndex(indexes+x) > channel_statistics[BlackChannel].maxima)
2544  channel_statistics[BlackChannel].maxima=(double)
2545  GetPixelIndex(indexes+x);
2546  channel_statistics[BlackChannel].sum+=QuantumScale*
2547  GetPixelIndex(indexes+x);
2548  channel_statistics[BlackChannel].sum_squared+=QuantumScale*
2549  GetPixelIndex(indexes+x)*QuantumScale*GetPixelIndex(indexes+x);
2550  channel_statistics[BlackChannel].sum_cubed+=QuantumScale*
2551  GetPixelIndex(indexes+x)*QuantumScale*GetPixelIndex(indexes+x)*
2552  QuantumScale*GetPixelIndex(indexes+x);
2553  channel_statistics[BlackChannel].sum_fourth_power+=QuantumScale*
2554  GetPixelIndex(indexes+x)*QuantumScale*GetPixelIndex(indexes+x)*
2555  QuantumScale*GetPixelIndex(indexes+x)*QuantumScale*
2556  GetPixelIndex(indexes+x);
2557  histogram[ScaleQuantumToMap(GetPixelIndex(indexes+x))].index++;
2558  }
2559  x++;
2560  p++;
2561  }
2562  }
2563  for (i=0; i < (ssize_t) CompositeChannels; i++)
2564  {
2565  double
2566  area,
2567  mean,
2568  standard_deviation;
2569 
2570  /*
2571  Normalize pixel statistics.
2572  */
2573  area=PerceptibleReciprocal((double) image->columns*image->rows);
2574  mean=channel_statistics[i].sum*area;
2575  channel_statistics[i].sum=mean;
2576  channel_statistics[i].sum_squared*=area;
2577  channel_statistics[i].sum_cubed*=area;
2578  channel_statistics[i].sum_fourth_power*=area;
2579  channel_statistics[i].mean=mean;
2580  channel_statistics[i].variance=channel_statistics[i].sum_squared;
2581  standard_deviation=sqrt(channel_statistics[i].variance-(mean*mean));
2582  area=PerceptibleReciprocal((double) image->columns*image->rows-1.0)*
2583  ((double) image->columns*image->rows);
2584  standard_deviation=sqrt(area*standard_deviation*standard_deviation);
2585  channel_statistics[i].standard_deviation=standard_deviation;
2586  }
2587  for (i=0; i < (ssize_t) (MaxMap+1U); i++)
2588  {
2589  if (histogram[i].red > 0.0)
2590  number_bins.red++;
2591  if (histogram[i].green > 0.0)
2592  number_bins.green++;
2593  if (histogram[i].blue > 0.0)
2594  number_bins.blue++;
2595  if ((image->matte != MagickFalse) && (histogram[i].opacity > 0.0))
2596  number_bins.opacity++;
2597  if ((image->colorspace == CMYKColorspace) && (histogram[i].index > 0.0))
2598  number_bins.index++;
2599  }
2600  area=PerceptibleReciprocal((double) image->columns*image->rows);
2601  for (i=0; i < (ssize_t) (MaxMap+1U); i++)
2602  {
2603  /*
2604  Compute pixel entropy.
2605  */
2606  histogram[i].red*=area;
2607  channel_statistics[RedChannel].entropy+=-histogram[i].red*
2608  MagickLog10(histogram[i].red)*
2609  PerceptibleReciprocal(MagickLog10((double) number_bins.red));
2610  histogram[i].green*=area;
2611  channel_statistics[GreenChannel].entropy+=-histogram[i].green*
2612  MagickLog10(histogram[i].green)*
2613  PerceptibleReciprocal(MagickLog10((double) number_bins.green));
2614  histogram[i].blue*=area;
2615  channel_statistics[BlueChannel].entropy+=-histogram[i].blue*
2616  MagickLog10(histogram[i].blue)*
2617  PerceptibleReciprocal(MagickLog10((double) number_bins.blue));
2618  if (image->matte != MagickFalse)
2619  {
2620  histogram[i].opacity*=area;
2621  channel_statistics[OpacityChannel].entropy+=-histogram[i].opacity*
2622  MagickLog10(histogram[i].opacity)*
2623  PerceptibleReciprocal(MagickLog10((double) number_bins.opacity));
2624  }
2625  if (image->colorspace == CMYKColorspace)
2626  {
2627  histogram[i].index*=area;
2628  channel_statistics[IndexChannel].entropy+=-histogram[i].index*
2629  MagickLog10(histogram[i].index)*
2630  PerceptibleReciprocal(MagickLog10((double) number_bins.index));
2631  }
2632  }
2633  /*
2634  Compute overall statistics.
2635  */
2636  for (i=0; i < (ssize_t) CompositeChannels; i++)
2637  {
2638  channel_statistics[CompositeChannels].depth=(size_t) EvaluateMax((double)
2639  channel_statistics[CompositeChannels].depth,(double)
2640  channel_statistics[i].depth);
2641  channel_statistics[CompositeChannels].minima=MagickMin(
2642  channel_statistics[CompositeChannels].minima,
2643  channel_statistics[i].minima);
2644  channel_statistics[CompositeChannels].maxima=EvaluateMax(
2645  channel_statistics[CompositeChannels].maxima,
2646  channel_statistics[i].maxima);
2647  channel_statistics[CompositeChannels].sum+=channel_statistics[i].sum;
2648  channel_statistics[CompositeChannels].sum_squared+=
2649  channel_statistics[i].sum_squared;
2650  channel_statistics[CompositeChannels].sum_cubed+=
2651  channel_statistics[i].sum_cubed;
2652  channel_statistics[CompositeChannels].sum_fourth_power+=
2653  channel_statistics[i].sum_fourth_power;
2654  channel_statistics[CompositeChannels].mean+=channel_statistics[i].mean;
2655  channel_statistics[CompositeChannels].variance+=
2656  channel_statistics[i].variance-channel_statistics[i].mean*
2657  channel_statistics[i].mean;
2658  standard_deviation=sqrt(channel_statistics[i].variance-
2659  (channel_statistics[i].mean*channel_statistics[i].mean));
2660  area=PerceptibleReciprocal((double) image->columns*image->rows-1.0)*
2661  ((double) image->columns*image->rows);
2662  standard_deviation=sqrt(area*standard_deviation*standard_deviation);
2663  channel_statistics[CompositeChannels].standard_deviation=standard_deviation;
2664  channel_statistics[CompositeChannels].entropy+=
2665  channel_statistics[i].entropy;
2666  }
2667  channels=3;
2668  if (image->matte != MagickFalse)
2669  channels++;
2670  if (image->colorspace == CMYKColorspace)
2671  channels++;
2672  channel_statistics[CompositeChannels].sum/=channels;
2673  channel_statistics[CompositeChannels].sum_squared/=channels;
2674  channel_statistics[CompositeChannels].sum_cubed/=channels;
2675  channel_statistics[CompositeChannels].sum_fourth_power/=channels;
2676  channel_statistics[CompositeChannels].mean/=channels;
2677  channel_statistics[CompositeChannels].kurtosis/=channels;
2678  channel_statistics[CompositeChannels].skewness/=channels;
2679  channel_statistics[CompositeChannels].entropy/=channels;
2680  i=CompositeChannels;
2681  area=PerceptibleReciprocal((double) channels*image->columns*image->rows);
2682  channel_statistics[i].variance=channel_statistics[i].sum_squared;
2683  channel_statistics[i].mean=channel_statistics[i].sum;
2684  standard_deviation=sqrt(channel_statistics[i].variance-
2685  (channel_statistics[i].mean*channel_statistics[i].mean));
2686  standard_deviation=sqrt(PerceptibleReciprocal((double) channels*
2687  image->columns*image->rows-1.0)*channels*image->columns*image->rows*
2688  standard_deviation*standard_deviation);
2689  channel_statistics[i].standard_deviation=standard_deviation;
2690  for (i=0; i <= (ssize_t) CompositeChannels; i++)
2691  {
2692  /*
2693  Compute kurtosis & skewness statistics.
2694  */
2695  standard_deviation=PerceptibleReciprocal(
2696  channel_statistics[i].standard_deviation);
2697  channel_statistics[i].skewness=(channel_statistics[i].sum_cubed-3.0*
2698  channel_statistics[i].mean*channel_statistics[i].sum_squared+2.0*
2699  channel_statistics[i].mean*channel_statistics[i].mean*
2700  channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2701  standard_deviation);
2702  channel_statistics[i].kurtosis=(channel_statistics[i].sum_fourth_power-4.0*
2703  channel_statistics[i].mean*channel_statistics[i].sum_cubed+6.0*
2704  channel_statistics[i].mean*channel_statistics[i].mean*
2705  channel_statistics[i].sum_squared-3.0*channel_statistics[i].mean*
2706  channel_statistics[i].mean*1.0*channel_statistics[i].mean*
2707  channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2708  standard_deviation*standard_deviation)-3.0;
2709  }
2710  for (i=0; i <= (ssize_t) CompositeChannels; i++)
2711  {
2712  channel_statistics[i].mean*=QuantumRange;
2713  channel_statistics[i].variance*=QuantumRange;
2714  channel_statistics[i].standard_deviation*=QuantumRange;
2715  channel_statistics[i].sum*=QuantumRange;
2716  channel_statistics[i].sum_squared*=QuantumRange;
2717  channel_statistics[i].sum_cubed*=QuantumRange;
2718  channel_statistics[i].sum_fourth_power*=QuantumRange;
2719  }
2720  channel_statistics[CompositeChannels].mean=0.0;
2721  channel_statistics[CompositeChannels].standard_deviation=0.0;
2722  for (i=0; i < (ssize_t) CompositeChannels; i++)
2723  {
2724  channel_statistics[CompositeChannels].mean+=
2725  channel_statistics[i].mean;
2726  channel_statistics[CompositeChannels].standard_deviation+=
2727  channel_statistics[i].standard_deviation;
2728  }
2729  channel_statistics[CompositeChannels].mean/=(double) channels;
2730  channel_statistics[CompositeChannels].standard_deviation/=(double) channels;
2731  histogram=(MagickPixelPacket *) RelinquishMagickMemory(histogram);
2732  if (y < (ssize_t) image->rows)
2733  channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
2734  channel_statistics);
2735  return(channel_statistics);
2736 }
2737 
2738 /*
2739 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2740 % %
2741 % %
2742 % %
2743 % P o l y n o m i a l I m a g e %
2744 % %
2745 % %
2746 % %
2747 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2748 %
2749 % PolynomialImage() returns a new image where each pixel is the sum of the
2750 % pixels in the image sequence after applying its corresponding terms
2751 % (coefficient and degree pairs).
2752 %
2753 % The format of the PolynomialImage method is:
2754 %
2755 % Image *PolynomialImage(const Image *images,const size_t number_terms,
2756 % const double *terms,ExceptionInfo *exception)
2757 % Image *PolynomialImageChannel(const Image *images,
2758 % const size_t number_terms,const ChannelType channel,
2759 % const double *terms,ExceptionInfo *exception)
2760 %
2761 % A description of each parameter follows:
2762 %
2763 % o images: the image sequence.
2764 %
2765 % o channel: the channel.
2766 %
2767 % o number_terms: the number of terms in the list. The actual list length
2768 % is 2 x number_terms + 1 (the constant).
2769 %
2770 % o terms: the list of polynomial coefficients and degree pairs and a
2771 % constant.
2772 %
2773 % o exception: return any errors or warnings in this structure.
2774 %
2775 */
2776 MagickExport Image *PolynomialImage(const Image *images,
2777  const size_t number_terms,const double *terms,ExceptionInfo *exception)
2778 {
2779  Image
2780  *polynomial_image;
2781 
2782  polynomial_image=PolynomialImageChannel(images,DefaultChannels,number_terms,
2783  terms,exception);
2784  return(polynomial_image);
2785 }
2786 
2787 MagickExport Image *PolynomialImageChannel(const Image *images,
2788  const ChannelType channel,const size_t number_terms,const double *terms,
2789  ExceptionInfo *exception)
2790 {
2791 #define PolynomialImageTag "Polynomial/Image"
2792 
2793  CacheView
2794  *polynomial_view;
2795 
2796  Image
2797  *image;
2798 
2799  MagickBooleanType
2800  status;
2801 
2802  MagickOffsetType
2803  progress;
2804 
2806  **magick_restrict polynomial_pixels,
2807  zero;
2808 
2809  ssize_t
2810  y;
2811 
2812  assert(images != (Image *) NULL);
2813  assert(images->signature == MagickCoreSignature);
2814  if (IsEventLogging() != MagickFalse)
2815  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
2816  assert(exception != (ExceptionInfo *) NULL);
2817  assert(exception->signature == MagickCoreSignature);
2818  image=AcquireImageCanvas(images,exception);
2819  if (image == (Image *) NULL)
2820  return((Image *) NULL);
2821  if (SetImageStorageClass(image,DirectClass) == MagickFalse)
2822  {
2823  InheritException(exception,&image->exception);
2824  image=DestroyImage(image);
2825  return((Image *) NULL);
2826  }
2827  polynomial_pixels=AcquirePixelTLS(images);
2828  if (polynomial_pixels == (MagickPixelPacket **) NULL)
2829  {
2830  image=DestroyImage(image);
2831  (void) ThrowMagickException(exception,GetMagickModule(),
2832  ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
2833  return((Image *) NULL);
2834  }
2835  /*
2836  Polynomial image pixels.
2837  */
2838  status=MagickTrue;
2839  progress=0;
2840  GetMagickPixelPacket(images,&zero);
2841  polynomial_view=AcquireAuthenticCacheView(image,exception);
2842 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2843  #pragma omp parallel for schedule(static) shared(progress,status) \
2844  magick_number_threads(image,image,image->rows,1)
2845 #endif
2846  for (y=0; y < (ssize_t) image->rows; y++)
2847  {
2848  CacheView
2849  *image_view;
2850 
2851  const Image
2852  *next;
2853 
2854  const int
2855  id = GetOpenMPThreadId();
2856 
2857  IndexPacket
2858  *magick_restrict polynomial_indexes;
2859 
2861  *polynomial_pixel;
2862 
2863  PixelPacket
2864  *magick_restrict q;
2865 
2866  ssize_t
2867  i,
2868  x;
2869 
2870  size_t
2871  number_images;
2872 
2873  if (status == MagickFalse)
2874  continue;
2875  q=QueueCacheViewAuthenticPixels(polynomial_view,0,y,image->columns,1,
2876  exception);
2877  if (q == (PixelPacket *) NULL)
2878  {
2879  status=MagickFalse;
2880  continue;
2881  }
2882  polynomial_indexes=GetCacheViewAuthenticIndexQueue(polynomial_view);
2883  polynomial_pixel=polynomial_pixels[id];
2884  for (x=0; x < (ssize_t) image->columns; x++)
2885  polynomial_pixel[x]=zero;
2886  next=images;
2887  number_images=GetImageListLength(images);
2888  for (i=0; i < (ssize_t) number_images; i++)
2889  {
2890  const IndexPacket
2891  *indexes;
2892 
2893  const PixelPacket
2894  *p;
2895 
2896  if (i >= (ssize_t) number_terms)
2897  break;
2898  image_view=AcquireVirtualCacheView(next,exception);
2899  p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2900  if (p == (const PixelPacket *) NULL)
2901  {
2902  image_view=DestroyCacheView(image_view);
2903  break;
2904  }
2905  indexes=GetCacheViewVirtualIndexQueue(image_view);
2906  for (x=0; x < (ssize_t) image->columns; x++)
2907  {
2908  double
2909  coefficient,
2910  degree;
2911 
2912  coefficient=terms[i << 1];
2913  degree=terms[(i << 1)+1];
2914  if ((channel & RedChannel) != 0)
2915  polynomial_pixel[x].red+=coefficient*pow(QuantumScale*(double)
2916  p->red,degree);
2917  if ((channel & GreenChannel) != 0)
2918  polynomial_pixel[x].green+=coefficient*pow(QuantumScale*(double)
2919  p->green,
2920  degree);
2921  if ((channel & BlueChannel) != 0)
2922  polynomial_pixel[x].blue+=coefficient*pow(QuantumScale*(double)
2923  p->blue,degree);
2924  if ((channel & OpacityChannel) != 0)
2925  polynomial_pixel[x].opacity+=coefficient*pow(QuantumScale*
2926  ((double) QuantumRange-(double) p->opacity),degree);
2927  if (((channel & IndexChannel) != 0) &&
2928  (image->colorspace == CMYKColorspace))
2929  polynomial_pixel[x].index+=coefficient*pow(QuantumScale*(double)
2930  indexes[x],degree);
2931  p++;
2932  }
2933  image_view=DestroyCacheView(image_view);
2934  next=GetNextImageInList(next);
2935  }
2936  for (x=0; x < (ssize_t) image->columns; x++)
2937  {
2938  SetPixelRed(q,ClampToQuantum((MagickRealType) QuantumRange*
2939  polynomial_pixel[x].red));
2940  SetPixelGreen(q,ClampToQuantum((MagickRealType) QuantumRange*
2941  polynomial_pixel[x].green));
2942  SetPixelBlue(q,ClampToQuantum((MagickRealType) QuantumRange*
2943  polynomial_pixel[x].blue));
2944  if (image->matte == MagickFalse)
2945  SetPixelOpacity(q,ClampToQuantum((MagickRealType) QuantumRange-
2946  (MagickRealType) QuantumRange*polynomial_pixel[x].opacity));
2947  else
2948  SetPixelAlpha(q,ClampToQuantum((MagickRealType) QuantumRange-
2949  (MagickRealType) QuantumRange*polynomial_pixel[x].opacity));
2950  if (image->colorspace == CMYKColorspace)
2951  SetPixelIndex(polynomial_indexes+x,ClampToQuantum((MagickRealType)
2952  QuantumRange*polynomial_pixel[x].index));
2953  q++;
2954  }
2955  if (SyncCacheViewAuthenticPixels(polynomial_view,exception) == MagickFalse)
2956  status=MagickFalse;
2957  if (images->progress_monitor != (MagickProgressMonitor) NULL)
2958  {
2959  MagickBooleanType
2960  proceed;
2961 
2962  proceed=SetImageProgress(images,PolynomialImageTag,progress++,
2963  image->rows);
2964  if (proceed == MagickFalse)
2965  status=MagickFalse;
2966  }
2967  }
2968  polynomial_view=DestroyCacheView(polynomial_view);
2969  polynomial_pixels=DestroyPixelTLS(images,polynomial_pixels);
2970  if (status == MagickFalse)
2971  image=DestroyImage(image);
2972  return(image);
2973 }
2974 
2975 /*
2976 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2977 % %
2978 % %
2979 % %
2980 % S t a t i s t i c I m a g e %
2981 % %
2982 % %
2983 % %
2984 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2985 %
2986 % StatisticImage() makes each pixel the min / max / median / mode / etc. of
2987 % the neighborhood of the specified width and height.
2988 %
2989 % The format of the StatisticImage method is:
2990 %
2991 % Image *StatisticImage(const Image *image,const StatisticType type,
2992 % const size_t width,const size_t height,ExceptionInfo *exception)
2993 % Image *StatisticImageChannel(const Image *image,
2994 % const ChannelType channel,const StatisticType type,
2995 % const size_t width,const size_t height,ExceptionInfo *exception)
2996 %
2997 % A description of each parameter follows:
2998 %
2999 % o image: the image.
3000 %
3001 % o channel: the image channel.
3002 %
3003 % o type: the statistic type (median, mode, etc.).
3004 %
3005 % o width: the width of the pixel neighborhood.
3006 %
3007 % o height: the height of the pixel neighborhood.
3008 %
3009 % o exception: return any errors or warnings in this structure.
3010 %
3011 */
3012 
3013 #define ListChannels 5
3014 
3015 typedef struct _ListNode
3016 {
3017  size_t
3018  next[9],
3019  count,
3020  signature;
3021 } ListNode;
3022 
3023 typedef struct _SkipList
3024 {
3025  ssize_t
3026  level;
3027 
3028  ListNode
3029  *nodes;
3030 } SkipList;
3031 
3032 typedef struct _PixelList
3033 {
3034  size_t
3035  length,
3036  seed,
3037  signature;
3038 
3039  SkipList
3040  lists[ListChannels];
3041 } PixelList;
3042 
3043 static PixelList *DestroyPixelList(PixelList *pixel_list)
3044 {
3045  ssize_t
3046  i;
3047 
3048  if (pixel_list == (PixelList *) NULL)
3049  return((PixelList *) NULL);
3050  for (i=0; i < ListChannels; i++)
3051  if (pixel_list->lists[i].nodes != (ListNode *) NULL)
3052  pixel_list->lists[i].nodes=(ListNode *) RelinquishAlignedMemory(
3053  pixel_list->lists[i].nodes);
3054  pixel_list=(PixelList *) RelinquishMagickMemory(pixel_list);
3055  return(pixel_list);
3056 }
3057 
3058 static PixelList **DestroyPixelListTLS(PixelList **pixel_list)
3059 {
3060  ssize_t
3061  i;
3062 
3063  assert(pixel_list != (PixelList **) NULL);
3064  for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
3065  if (pixel_list[i] != (PixelList *) NULL)
3066  pixel_list[i]=DestroyPixelList(pixel_list[i]);
3067  pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list);
3068  return(pixel_list);
3069 }
3070 
3071 static PixelList *AcquirePixelList(const size_t width,const size_t height)
3072 {
3073  PixelList
3074  *pixel_list;
3075 
3076  ssize_t
3077  i;
3078 
3079  pixel_list=(PixelList *) AcquireMagickMemory(sizeof(*pixel_list));
3080  if (pixel_list == (PixelList *) NULL)
3081  return(pixel_list);
3082  (void) memset((void *) pixel_list,0,sizeof(*pixel_list));
3083  pixel_list->length=width*height;
3084  for (i=0; i < ListChannels; i++)
3085  {
3086  pixel_list->lists[i].nodes=(ListNode *) AcquireAlignedMemory(65537UL,
3087  sizeof(*pixel_list->lists[i].nodes));
3088  if (pixel_list->lists[i].nodes == (ListNode *) NULL)
3089  return(DestroyPixelList(pixel_list));
3090  (void) memset(pixel_list->lists[i].nodes,0,65537UL*
3091  sizeof(*pixel_list->lists[i].nodes));
3092  }
3093  pixel_list->signature=MagickCoreSignature;
3094  return(pixel_list);
3095 }
3096 
3097 static PixelList **AcquirePixelListTLS(const size_t width,
3098  const size_t height)
3099 {
3100  PixelList
3101  **pixel_list;
3102 
3103  ssize_t
3104  i;
3105 
3106  size_t
3107  number_threads;
3108 
3109  number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
3110  pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
3111  sizeof(*pixel_list));
3112  if (pixel_list == (PixelList **) NULL)
3113  return((PixelList **) NULL);
3114  (void) memset(pixel_list,0,number_threads*sizeof(*pixel_list));
3115  for (i=0; i < (ssize_t) number_threads; i++)
3116  {
3117  pixel_list[i]=AcquirePixelList(width,height);
3118  if (pixel_list[i] == (PixelList *) NULL)
3119  return(DestroyPixelListTLS(pixel_list));
3120  }
3121  return(pixel_list);
3122 }
3123 
3124 static void AddNodePixelList(PixelList *pixel_list,const ssize_t channel,
3125  const size_t color)
3126 {
3127  SkipList
3128  *list;
3129 
3130  ssize_t
3131  level;
3132 
3133  size_t
3134  search,
3135  update[9];
3136 
3137  /*
3138  Initialize the node.
3139  */
3140  list=pixel_list->lists+channel;
3141  list->nodes[color].signature=pixel_list->signature;
3142  list->nodes[color].count=1;
3143  /*
3144  Determine where it belongs in the list.
3145  */
3146  search=65536UL;
3147  (void) memset(update,0,sizeof(update));
3148  for (level=list->level; level >= 0; level--)
3149  {
3150  while (list->nodes[search].next[level] < color)
3151  search=list->nodes[search].next[level];
3152  update[level]=search;
3153  }
3154  /*
3155  Generate a pseudo-random level for this node.
3156  */
3157  for (level=0; ; level++)
3158  {
3159  pixel_list->seed=(pixel_list->seed*42893621L)+1L;
3160  if ((pixel_list->seed & 0x300) != 0x300)
3161  break;
3162  }
3163  if (level > 8)
3164  level=8;
3165  if (level > (list->level+2))
3166  level=list->level+2;
3167  /*
3168  If we're raising the list's level, link back to the root node.
3169  */
3170  while (level > list->level)
3171  {
3172  list->level++;
3173  update[list->level]=65536UL;
3174  }
3175  /*
3176  Link the node into the skip-list.
3177  */
3178  do
3179  {
3180  list->nodes[color].next[level]=list->nodes[update[level]].next[level];
3181  list->nodes[update[level]].next[level]=color;
3182  } while (level-- > 0);
3183 }
3184 
3185 static void GetMaximumPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3186 {
3187  SkipList
3188  *list;
3189 
3190  ssize_t
3191  channel;
3192 
3193  size_t
3194  color,
3195  maximum;
3196 
3197  ssize_t
3198  count;
3199 
3200  unsigned short
3201  channels[ListChannels];
3202 
3203  /*
3204  Find the maximum value for each of the color.
3205  */
3206  for (channel=0; channel < 5; channel++)
3207  {
3208  list=pixel_list->lists+channel;
3209  color=65536L;
3210  count=0;
3211  maximum=list->nodes[color].next[0];
3212  do
3213  {
3214  color=list->nodes[color].next[0];
3215  if (color > maximum)
3216  maximum=color;
3217  count+=list->nodes[color].count;
3218  } while (count < (ssize_t) pixel_list->length);
3219  channels[channel]=(unsigned short) maximum;
3220  }
3221  pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3222  pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3223  pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3224  pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3225  pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3226 }
3227 
3228 static void GetMeanPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3229 {
3230  MagickRealType
3231  sum;
3232 
3233  SkipList
3234  *list;
3235 
3236  ssize_t
3237  channel;
3238 
3239  size_t
3240  color;
3241 
3242  ssize_t
3243  count;
3244 
3245  unsigned short
3246  channels[ListChannels];
3247 
3248  /*
3249  Find the mean value for each of the color.
3250  */
3251  for (channel=0; channel < 5; channel++)
3252  {
3253  list=pixel_list->lists+channel;
3254  color=65536L;
3255  count=0;
3256  sum=0.0;
3257  do
3258  {
3259  color=list->nodes[color].next[0];
3260  sum+=(MagickRealType) list->nodes[color].count*color;
3261  count+=list->nodes[color].count;
3262  } while (count < (ssize_t) pixel_list->length);
3263  sum/=pixel_list->length;
3264  channels[channel]=(unsigned short) sum;
3265  }
3266  pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3267  pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3268  pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3269  pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3270  pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3271 }
3272 
3273 static void GetMedianPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3274 {
3275  SkipList
3276  *list;
3277 
3278  ssize_t
3279  channel;
3280 
3281  size_t
3282  color;
3283 
3284  ssize_t
3285  count;
3286 
3287  unsigned short
3288  channels[ListChannels];
3289 
3290  /*
3291  Find the median value for each of the color.
3292  */
3293  for (channel=0; channel < 5; channel++)
3294  {
3295  list=pixel_list->lists+channel;
3296  color=65536L;
3297  count=0;
3298  do
3299  {
3300  color=list->nodes[color].next[0];
3301  count+=list->nodes[color].count;
3302  } while (count <= (ssize_t) (pixel_list->length >> 1));
3303  channels[channel]=(unsigned short) color;
3304  }
3305  GetMagickPixelPacket((const Image *) NULL,pixel);
3306  pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3307  pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3308  pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3309  pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3310  pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3311 }
3312 
3313 static void GetMinimumPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3314 {
3315  SkipList
3316  *list;
3317 
3318  ssize_t
3319  channel;
3320 
3321  size_t
3322  color,
3323  minimum;
3324 
3325  ssize_t
3326  count;
3327 
3328  unsigned short
3329  channels[ListChannels];
3330 
3331  /*
3332  Find the minimum value for each of the color.
3333  */
3334  for (channel=0; channel < 5; channel++)
3335  {
3336  list=pixel_list->lists+channel;
3337  count=0;
3338  color=65536UL;
3339  minimum=list->nodes[color].next[0];
3340  do
3341  {
3342  color=list->nodes[color].next[0];
3343  if (color < minimum)
3344  minimum=color;
3345  count+=list->nodes[color].count;
3346  } while (count < (ssize_t) pixel_list->length);
3347  channels[channel]=(unsigned short) minimum;
3348  }
3349  pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3350  pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3351  pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3352  pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3353  pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3354 }
3355 
3356 static void GetModePixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3357 {
3358  SkipList
3359  *list;
3360 
3361  ssize_t
3362  channel;
3363 
3364  size_t
3365  color,
3366  max_count,
3367  mode;
3368 
3369  ssize_t
3370  count;
3371 
3372  unsigned short
3373  channels[5];
3374 
3375  /*
3376  Make each pixel the 'predominant color' of the specified neighborhood.
3377  */
3378  for (channel=0; channel < 5; channel++)
3379  {
3380  list=pixel_list->lists+channel;
3381  color=65536L;
3382  mode=color;
3383  max_count=list->nodes[mode].count;
3384  count=0;
3385  do
3386  {
3387  color=list->nodes[color].next[0];
3388  if (list->nodes[color].count > max_count)
3389  {
3390  mode=color;
3391  max_count=list->nodes[mode].count;
3392  }
3393  count+=list->nodes[color].count;
3394  } while (count < (ssize_t) pixel_list->length);
3395  channels[channel]=(unsigned short) mode;
3396  }
3397  pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3398  pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3399  pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3400  pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3401  pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3402 }
3403 
3404 static void GetNonpeakPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3405 {
3406  SkipList
3407  *list;
3408 
3409  ssize_t
3410  channel;
3411 
3412  size_t
3413  color,
3414  next,
3415  previous;
3416 
3417  ssize_t
3418  count;
3419 
3420  unsigned short
3421  channels[5];
3422 
3423  /*
3424  Finds the non peak value for each of the colors.
3425  */
3426  for (channel=0; channel < 5; channel++)
3427  {
3428  list=pixel_list->lists+channel;
3429  color=65536L;
3430  next=list->nodes[color].next[0];
3431  count=0;
3432  do
3433  {
3434  previous=color;
3435  color=next;
3436  next=list->nodes[color].next[0];
3437  count+=list->nodes[color].count;
3438  } while (count <= (ssize_t) (pixel_list->length >> 1));
3439  if ((previous == 65536UL) && (next != 65536UL))
3440  color=next;
3441  else
3442  if ((previous != 65536UL) && (next == 65536UL))
3443  color=previous;
3444  channels[channel]=(unsigned short) color;
3445  }
3446  pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3447  pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3448  pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3449  pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3450  pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3451 }
3452 
3453 static void GetRootMeanSquarePixelList(PixelList *pixel_list,
3454  MagickPixelPacket *pixel)
3455 {
3456  MagickRealType
3457  sum;
3458 
3459  SkipList
3460  *list;
3461 
3462  ssize_t
3463  channel;
3464 
3465  size_t
3466  color;
3467 
3468  ssize_t
3469  count;
3470 
3471  unsigned short
3472  channels[ListChannels];
3473 
3474  /*
3475  Find the root mean square value for each of the color.
3476  */
3477  for (channel=0; channel < 5; channel++)
3478  {
3479  list=pixel_list->lists+channel;
3480  color=65536L;
3481  count=0;
3482  sum=0.0;
3483  do
3484  {
3485  color=list->nodes[color].next[0];
3486  sum+=(MagickRealType) (list->nodes[color].count*color*color);
3487  count+=list->nodes[color].count;
3488  } while (count < (ssize_t) pixel_list->length);
3489  sum/=pixel_list->length;
3490  channels[channel]=(unsigned short) sqrt(sum);
3491  }
3492  pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3493  pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3494  pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3495  pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3496  pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3497 }
3498 
3499 static void GetStandardDeviationPixelList(PixelList *pixel_list,
3500  MagickPixelPacket *pixel)
3501 {
3502  MagickRealType
3503  sum,
3504  sum_squared;
3505 
3506  SkipList
3507  *list;
3508 
3509  size_t
3510  color;
3511 
3512  ssize_t
3513  channel,
3514  count;
3515 
3516  unsigned short
3517  channels[ListChannels];
3518 
3519  /*
3520  Find the standard-deviation value for each of the color.
3521  */
3522  for (channel=0; channel < 5; channel++)
3523  {
3524  list=pixel_list->lists+channel;
3525  color=65536L;
3526  count=0;
3527  sum=0.0;
3528  sum_squared=0.0;
3529  do
3530  {
3531  ssize_t
3532  i;
3533 
3534  color=list->nodes[color].next[0];
3535  sum+=(MagickRealType) list->nodes[color].count*color;
3536  for (i=0; i < (ssize_t) list->nodes[color].count; i++)
3537  sum_squared+=((MagickRealType) color)*((MagickRealType) color);
3538  count+=list->nodes[color].count;
3539  } while (count < (ssize_t) pixel_list->length);
3540  sum/=pixel_list->length;
3541  sum_squared/=pixel_list->length;
3542  channels[channel]=(unsigned short) sqrt(sum_squared-(sum*sum));
3543  }
3544  pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3545  pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3546  pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3547  pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3548  pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3549 }
3550 
3551 static inline void InsertPixelList(const Image *image,const PixelPacket *pixel,
3552  const IndexPacket *indexes,PixelList *pixel_list)
3553 {
3554  size_t
3555  signature;
3556 
3557  unsigned short
3558  index;
3559 
3560  index=ScaleQuantumToShort(GetPixelRed(pixel));
3561  signature=pixel_list->lists[0].nodes[index].signature;
3562  if (signature == pixel_list->signature)
3563  pixel_list->lists[0].nodes[index].count++;
3564  else
3565  AddNodePixelList(pixel_list,0,index);
3566  index=ScaleQuantumToShort(GetPixelGreen(pixel));
3567  signature=pixel_list->lists[1].nodes[index].signature;
3568  if (signature == pixel_list->signature)
3569  pixel_list->lists[1].nodes[index].count++;
3570  else
3571  AddNodePixelList(pixel_list,1,index);
3572  index=ScaleQuantumToShort(GetPixelBlue(pixel));
3573  signature=pixel_list->lists[2].nodes[index].signature;
3574  if (signature == pixel_list->signature)
3575  pixel_list->lists[2].nodes[index].count++;
3576  else
3577  AddNodePixelList(pixel_list,2,index);
3578  index=ScaleQuantumToShort(GetPixelOpacity(pixel));
3579  signature=pixel_list->lists[3].nodes[index].signature;
3580  if (signature == pixel_list->signature)
3581  pixel_list->lists[3].nodes[index].count++;
3582  else
3583  AddNodePixelList(pixel_list,3,index);
3584  if (image->colorspace == CMYKColorspace)
3585  index=ScaleQuantumToShort(GetPixelIndex(indexes));
3586  signature=pixel_list->lists[4].nodes[index].signature;
3587  if (signature == pixel_list->signature)
3588  pixel_list->lists[4].nodes[index].count++;
3589  else
3590  AddNodePixelList(pixel_list,4,index);
3591 }
3592 
3593 static void ResetPixelList(PixelList *pixel_list)
3594 {
3595  int
3596  level;
3597 
3598  ListNode
3599  *root;
3600 
3601  SkipList
3602  *list;
3603 
3604  ssize_t
3605  channel;
3606 
3607  /*
3608  Reset the skip-list.
3609  */
3610  for (channel=0; channel < 5; channel++)
3611  {
3612  list=pixel_list->lists+channel;
3613  root=list->nodes+65536UL;
3614  list->level=0;
3615  for (level=0; level < 9; level++)
3616  root->next[level]=65536UL;
3617  }
3618  pixel_list->seed=pixel_list->signature++;
3619 }
3620 
3621 MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
3622  const size_t width,const size_t height,ExceptionInfo *exception)
3623 {
3624  Image
3625  *statistic_image;
3626 
3627  statistic_image=StatisticImageChannel(image,DefaultChannels,type,width,
3628  height,exception);
3629  return(statistic_image);
3630 }
3631 
3632 MagickExport Image *StatisticImageChannel(const Image *image,
3633  const ChannelType channel,const StatisticType type,const size_t width,
3634  const size_t height,ExceptionInfo *exception)
3635 {
3636 #define StatisticImageTag "Statistic/Image"
3637 
3638  CacheView
3639  *image_view,
3640  *statistic_view;
3641 
3642  Image
3643  *statistic_image;
3644 
3645  MagickBooleanType
3646  status;
3647 
3648  MagickOffsetType
3649  progress;
3650 
3651  PixelList
3652  **magick_restrict pixel_list;
3653 
3654  size_t
3655  neighbor_height,
3656  neighbor_width;
3657 
3658  ssize_t
3659  y;
3660 
3661  /*
3662  Initialize statistics image attributes.
3663  */
3664  assert(image != (Image *) NULL);
3665  assert(image->signature == MagickCoreSignature);
3666  if (IsEventLogging() != MagickFalse)
3667  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3668  assert(exception != (ExceptionInfo *) NULL);
3669  assert(exception->signature == MagickCoreSignature);
3670  statistic_image=CloneImage(image,0,0,MagickTrue,exception);
3671  if (statistic_image == (Image *) NULL)
3672  return((Image *) NULL);
3673  if (SetImageStorageClass(statistic_image,DirectClass) == MagickFalse)
3674  {
3675  InheritException(exception,&statistic_image->exception);
3676  statistic_image=DestroyImage(statistic_image);
3677  return((Image *) NULL);
3678  }
3679  neighbor_width=width == 0 ? GetOptimalKernelWidth2D((double) width,0.5) :
3680  width;
3681  neighbor_height=height == 0 ? GetOptimalKernelWidth2D((double) height,0.5) :
3682  height;
3683  pixel_list=AcquirePixelListTLS(neighbor_width,neighbor_height);
3684  if (pixel_list == (PixelList **) NULL)
3685  {
3686  statistic_image=DestroyImage(statistic_image);
3687  ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3688  }
3689  /*
3690  Make each pixel the min / max / median / mode / etc. of the neighborhood.
3691  */
3692  status=MagickTrue;
3693  progress=0;
3694  image_view=AcquireVirtualCacheView(image,exception);
3695  statistic_view=AcquireAuthenticCacheView(statistic_image,exception);
3696 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3697  #pragma omp parallel for schedule(static) shared(progress,status) \
3698  magick_number_threads(image,statistic_image,statistic_image->rows,1)
3699 #endif
3700  for (y=0; y < (ssize_t) statistic_image->rows; y++)
3701  {
3702  const int
3703  id = GetOpenMPThreadId();
3704 
3705  const IndexPacket
3706  *magick_restrict indexes;
3707 
3708  const PixelPacket
3709  *magick_restrict p;
3710 
3711  IndexPacket
3712  *magick_restrict statistic_indexes;
3713 
3714  PixelPacket
3715  *magick_restrict q;
3716 
3717  ssize_t
3718  x;
3719 
3720  if (status == MagickFalse)
3721  continue;
3722  p=GetCacheViewVirtualPixels(image_view,-((ssize_t) neighbor_width/2L),y-
3723  (ssize_t) (neighbor_height/2L),image->columns+neighbor_width,
3724  neighbor_height,exception);
3725  q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns, 1,exception);
3726  if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
3727  {
3728  status=MagickFalse;
3729  continue;
3730  }
3731  indexes=GetCacheViewVirtualIndexQueue(image_view);
3732  statistic_indexes=GetCacheViewAuthenticIndexQueue(statistic_view);
3733  for (x=0; x < (ssize_t) statistic_image->columns; x++)
3734  {
3736  pixel;
3737 
3738  const IndexPacket
3739  *magick_restrict s;
3740 
3741  const PixelPacket
3742  *magick_restrict r;
3743 
3744  ssize_t
3745  u,
3746  v;
3747 
3748  r=p;
3749  s=indexes+x;
3750  ResetPixelList(pixel_list[id]);
3751  for (v=0; v < (ssize_t) neighbor_height; v++)
3752  {
3753  for (u=0; u < (ssize_t) neighbor_width; u++)
3754  InsertPixelList(image,r+u,s+u,pixel_list[id]);
3755  r+=image->columns+neighbor_width;
3756  s+=image->columns+neighbor_width;
3757  }
3758  GetMagickPixelPacket(image,&pixel);
3759  SetMagickPixelPacket(image,p+neighbor_width*neighbor_height/2,indexes+x+
3760  neighbor_width*neighbor_height/2,&pixel);
3761  switch (type)
3762  {
3763  case GradientStatistic:
3764  {
3766  maximum,
3767  minimum;
3768 
3769  GetMinimumPixelList(pixel_list[id],&pixel);
3770  minimum=pixel;
3771  GetMaximumPixelList(pixel_list[id],&pixel);
3772  maximum=pixel;
3773  pixel.red=MagickAbsoluteValue(maximum.red-minimum.red);
3774  pixel.green=MagickAbsoluteValue(maximum.green-minimum.green);
3775  pixel.blue=MagickAbsoluteValue(maximum.blue-minimum.blue);
3776  pixel.opacity=MagickAbsoluteValue(maximum.opacity-minimum.opacity);
3777  if (image->colorspace == CMYKColorspace)
3778  pixel.index=MagickAbsoluteValue(maximum.index-minimum.index);
3779  break;
3780  }
3781  case MaximumStatistic:
3782  {
3783  GetMaximumPixelList(pixel_list[id],&pixel);
3784  break;
3785  }
3786  case MeanStatistic:
3787  {
3788  GetMeanPixelList(pixel_list[id],&pixel);
3789  break;
3790  }
3791  case MedianStatistic:
3792  default:
3793  {
3794  GetMedianPixelList(pixel_list[id],&pixel);
3795  break;
3796  }
3797  case MinimumStatistic:
3798  {
3799  GetMinimumPixelList(pixel_list[id],&pixel);
3800  break;
3801  }
3802  case ModeStatistic:
3803  {
3804  GetModePixelList(pixel_list[id],&pixel);
3805  break;
3806  }
3807  case NonpeakStatistic:
3808  {
3809  GetNonpeakPixelList(pixel_list[id],&pixel);
3810  break;
3811  }
3812  case RootMeanSquareStatistic:
3813  {
3814  GetRootMeanSquarePixelList(pixel_list[id],&pixel);
3815  break;
3816  }
3817  case StandardDeviationStatistic:
3818  {
3819  GetStandardDeviationPixelList(pixel_list[id],&pixel);
3820  break;
3821  }
3822  }
3823  if ((channel & RedChannel) != 0)
3824  SetPixelRed(q,ClampToQuantum(pixel.red));
3825  if ((channel & GreenChannel) != 0)
3826  SetPixelGreen(q,ClampToQuantum(pixel.green));
3827  if ((channel & BlueChannel) != 0)
3828  SetPixelBlue(q,ClampToQuantum(pixel.blue));
3829  if ((channel & OpacityChannel) != 0)
3830  SetPixelOpacity(q,ClampToQuantum(pixel.opacity));
3831  if (((channel & IndexChannel) != 0) &&
3832  (image->colorspace == CMYKColorspace))
3833  SetPixelIndex(statistic_indexes+x,ClampToQuantum(pixel.index));
3834  p++;
3835  q++;
3836  }
3837  if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
3838  status=MagickFalse;
3839  if (image->progress_monitor != (MagickProgressMonitor) NULL)
3840  {
3841  MagickBooleanType
3842  proceed;
3843 
3844  proceed=SetImageProgress(image,StatisticImageTag,progress++,
3845  image->rows);
3846  if (proceed == MagickFalse)
3847  status=MagickFalse;
3848  }
3849  }
3850  statistic_view=DestroyCacheView(statistic_view);
3851  image_view=DestroyCacheView(image_view);
3852  pixel_list=DestroyPixelListTLS(pixel_list);
3853  if (status == MagickFalse)
3854  statistic_image=DestroyImage(statistic_image);
3855  return(statistic_image);
3856 }
Definition: image.h:133