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