MagickCore  6.9.13-52
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/attribute.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  channel,
2143  i;
2144 
2145  /*
2146  Blur then transform to xyY colorspace.
2147  */
2148  hash_image=BlurImage(image,0.0,1.0,exception);
2149  if (hash_image == (Image *) NULL)
2150  return((ChannelPerceptualHash *) NULL);
2151  status=TransformImageColorspace(hash_image,xyYColorspace);
2152  if (status == MagickFalse)
2153  {
2154  hash_image=DestroyImage(hash_image);
2155  return((ChannelPerceptualHash *) NULL);
2156  }
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  (void) memset(perceptual_hash,0,(CompositeChannels+1UL)*
2166  sizeof(*perceptual_hash));
2167  for (channel=0; channel <= CompositeChannels; channel++)
2168  for (i=0; i < MaximumNumberOfPerceptualHashes; i++)
2169  perceptual_hash[channel].P[i]=(-MagickSafeLog10(fabs(
2170  moments[channel].I[i])));
2171  moments=(ChannelMoments *) RelinquishMagickMemory(moments);
2172  /*
2173  Blur then transform to HSB colorspace.
2174  */
2175  hash_image=BlurImage(image,0.0,1.0,exception);
2176  if (hash_image == (Image *) NULL)
2177  {
2178  perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
2179  perceptual_hash);
2180  return((ChannelPerceptualHash *) NULL);
2181  }
2182  status=TransformImageColorspace(hash_image,HSBColorspace);
2183  if (status == MagickFalse)
2184  {
2185  hash_image=DestroyImage(hash_image);
2186  perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
2187  perceptual_hash);
2188  return((ChannelPerceptualHash *) NULL);
2189  }
2190  moments=GetImageChannelMoments(hash_image,exception);
2191  hash_image=DestroyImage(hash_image);
2192  if (moments == (ChannelMoments *) NULL)
2193  {
2194  perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
2195  perceptual_hash);
2196  return((ChannelPerceptualHash *) NULL);
2197  }
2198  for (channel=0; channel <= CompositeChannels; channel++)
2199  for (i=0; i < MaximumNumberOfPerceptualHashes; i++)
2200  perceptual_hash[channel].Q[i]=(-MagickSafeLog10(fabs(
2201  moments[channel].I[i])));
2202  moments=(ChannelMoments *) RelinquishMagickMemory(moments);
2203  return(perceptual_hash);
2204 }
2205 
2206 /*
2207 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2208 % %
2209 % %
2210 % %
2211 % G e t I m a g e C h a n n e l R a n g e %
2212 % %
2213 % %
2214 % %
2215 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2216 %
2217 % GetImageChannelRange() returns the range of one or more image channels.
2218 %
2219 % The format of the GetImageChannelRange method is:
2220 %
2221 % MagickBooleanType GetImageChannelRange(const Image *image,
2222 % const ChannelType channel,double *minima,double *maxima,
2223 % ExceptionInfo *exception)
2224 %
2225 % A description of each parameter follows:
2226 %
2227 % o image: the image.
2228 %
2229 % o channel: the channel.
2230 %
2231 % o minima: the minimum value in the channel.
2232 %
2233 % o maxima: the maximum value in the channel.
2234 %
2235 % o exception: return any errors or warnings in this structure.
2236 %
2237 */
2238 
2239 MagickExport MagickBooleanType GetImageRange(const Image *image,
2240  double *minima,double *maxima,ExceptionInfo *exception)
2241 {
2242  return(GetImageChannelRange(image,CompositeChannels,minima,maxima,exception));
2243 }
2244 
2245 MagickExport MagickBooleanType GetImageChannelRange(const Image *image,
2246  const ChannelType channel,double *minima,double *maxima,
2247  ExceptionInfo *exception)
2248 {
2249  MagickPixelPacket
2250  pixel;
2251 
2252  ssize_t
2253  y;
2254 
2255  assert(image != (Image *) NULL);
2256  assert(image->signature == MagickCoreSignature);
2257  if (IsEventLogging() != MagickFalse)
2258  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2259  *maxima=(-MagickMaximumValue);
2260  *minima=MagickMaximumValue;
2261  GetMagickPixelPacket(image,&pixel);
2262  for (y=0; y < (ssize_t) image->rows; y++)
2263  {
2264  const IndexPacket
2265  *magick_restrict indexes;
2266 
2267  const PixelPacket
2268  *magick_restrict p;
2269 
2270  ssize_t
2271  x;
2272 
2273  p=GetVirtualPixels(image,0,y,image->columns,1,exception);
2274  if (p == (const PixelPacket *) NULL)
2275  break;
2276  indexes=GetVirtualIndexQueue(image);
2277  for (x=0; x < (ssize_t) image->columns; x++)
2278  {
2279  SetMagickPixelPacket(image,p,indexes+x,&pixel);
2280  if ((channel & RedChannel) != 0)
2281  {
2282  if (pixel.red < *minima)
2283  *minima=(double) pixel.red;
2284  if (pixel.red > *maxima)
2285  *maxima=(double) pixel.red;
2286  }
2287  if ((channel & GreenChannel) != 0)
2288  {
2289  if (pixel.green < *minima)
2290  *minima=(double) pixel.green;
2291  if (pixel.green > *maxima)
2292  *maxima=(double) pixel.green;
2293  }
2294  if ((channel & BlueChannel) != 0)
2295  {
2296  if (pixel.blue < *minima)
2297  *minima=(double) pixel.blue;
2298  if (pixel.blue > *maxima)
2299  *maxima=(double) pixel.blue;
2300  }
2301  if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
2302  {
2303  if (((MagickRealType) QuantumRange-(MagickRealType) pixel.opacity) < *minima)
2304  *minima=(double) ((MagickRealType) QuantumRange-(MagickRealType)
2305  pixel.opacity);
2306  if (((MagickRealType) QuantumRange-(MagickRealType) pixel.opacity) > *maxima)
2307  *maxima=(double) ((MagickRealType) QuantumRange-(MagickRealType)
2308  pixel.opacity);
2309  }
2310  if (((channel & IndexChannel) != 0) &&
2311  (image->colorspace == CMYKColorspace))
2312  {
2313  if ((double) pixel.index < *minima)
2314  *minima=(double) pixel.index;
2315  if ((double) pixel.index > *maxima)
2316  *maxima=(double) pixel.index;
2317  }
2318  p++;
2319  }
2320  }
2321  return(y == (ssize_t) image->rows ? MagickTrue : MagickFalse);
2322 }
2323 
2324 /*
2325 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2326 % %
2327 % %
2328 % %
2329 % 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 %
2330 % %
2331 % %
2332 % %
2333 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2334 %
2335 % GetImageChannelStatistics() returns statistics for each channel in the
2336 % image. The statistics include the channel depth, its minima, maxima, mean,
2337 % standard deviation, kurtosis and skewness. You can access the red channel
2338 % mean, for example, like this:
2339 %
2340 % channel_statistics=GetImageChannelStatistics(image,exception);
2341 % red_mean=channel_statistics[RedChannel].mean;
2342 %
2343 % Use MagickRelinquishMemory() to free the statistics buffer.
2344 %
2345 % The format of the GetImageChannelStatistics method is:
2346 %
2347 % ChannelStatistics *GetImageChannelStatistics(const Image *image,
2348 % ExceptionInfo *exception)
2349 %
2350 % A description of each parameter follows:
2351 %
2352 % o image: the image.
2353 %
2354 % o exception: return any errors or warnings in this structure.
2355 %
2356 */
2357 MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
2358  ExceptionInfo *exception)
2359 {
2360  ChannelStatistics
2361  *channel_statistics;
2362 
2363  double
2364  area,
2365  standard_deviation;
2366 
2367  MagickPixelPacket
2368  number_bins,
2369  *histogram;
2370 
2371  QuantumAny
2372  range;
2373 
2374  size_t
2375  channels,
2376  depth,
2377  length;
2378 
2379  ssize_t
2380  i,
2381  y;
2382 
2383  assert(image != (Image *) NULL);
2384  assert(image->signature == MagickCoreSignature);
2385  if (IsEventLogging() != MagickFalse)
2386  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2387  length=CompositeChannels+1UL;
2388  channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(length,
2389  sizeof(*channel_statistics));
2390  histogram=(MagickPixelPacket *) AcquireQuantumMemory(MaxMap+1U,
2391  sizeof(*histogram));
2392  if ((channel_statistics == (ChannelStatistics *) NULL) ||
2393  (histogram == (MagickPixelPacket *) NULL))
2394  {
2395  if (histogram != (MagickPixelPacket *) NULL)
2396  histogram=(MagickPixelPacket *) RelinquishMagickMemory(histogram);
2397  if (channel_statistics != (ChannelStatistics *) NULL)
2398  channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
2399  channel_statistics);
2400  return(channel_statistics);
2401  }
2402  (void) memset(channel_statistics,0,length*
2403  sizeof(*channel_statistics));
2404  for (i=0; i <= (ssize_t) CompositeChannels; i++)
2405  {
2406  ChannelStatistics *cs = channel_statistics+i;
2407  cs->depth=1;
2408  cs->maxima=(-MagickMaximumValue);
2409  cs->minima=MagickMaximumValue;
2410  cs->sum=0.0;
2411  cs->mean=0.0;
2412  cs->standard_deviation=0.0;
2413  cs->variance=0.0;
2414  cs->skewness=0.0;
2415  cs->kurtosis=0.0;
2416  cs->entropy=0.0;
2417  }
2418  (void) memset(histogram,0,(MaxMap+1U)*sizeof(*histogram));
2419  (void) memset(&number_bins,0,sizeof(number_bins));
2420  for (y=0; y < (ssize_t) image->rows; y++)
2421  {
2422  const IndexPacket
2423  *magick_restrict indexes;
2424 
2425  const PixelPacket
2426  *magick_restrict p;
2427 
2428  ssize_t
2429  x;
2430 
2431  /*
2432  Compute pixel statistics.
2433  */
2434  p=GetVirtualPixels(image,0,y,image->columns,1,exception);
2435  if (p == (const PixelPacket *) NULL)
2436  break;
2437  indexes=GetVirtualIndexQueue(image);
2438  for (x=0; x < (ssize_t) image->columns; )
2439  {
2440  if (channel_statistics[RedChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2441  {
2442  depth=channel_statistics[RedChannel].depth;
2443  range=GetQuantumRange(depth);
2444  if (IsPixelAtDepth(GetPixelRed(p),range) == MagickFalse)
2445  {
2446  channel_statistics[RedChannel].depth++;
2447  continue;
2448  }
2449  }
2450  if (channel_statistics[GreenChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2451  {
2452  depth=channel_statistics[GreenChannel].depth;
2453  range=GetQuantumRange(depth);
2454  if (IsPixelAtDepth(GetPixelGreen(p),range) == MagickFalse)
2455  {
2456  channel_statistics[GreenChannel].depth++;
2457  continue;
2458  }
2459  }
2460  if (channel_statistics[BlueChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2461  {
2462  depth=channel_statistics[BlueChannel].depth;
2463  range=GetQuantumRange(depth);
2464  if (IsPixelAtDepth(GetPixelBlue(p),range) == MagickFalse)
2465  {
2466  channel_statistics[BlueChannel].depth++;
2467  continue;
2468  }
2469  }
2470  if (image->matte != MagickFalse)
2471  {
2472  if (channel_statistics[OpacityChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2473  {
2474  depth=channel_statistics[OpacityChannel].depth;
2475  range=GetQuantumRange(depth);
2476  if (IsPixelAtDepth(GetPixelAlpha(p),range) == MagickFalse)
2477  {
2478  channel_statistics[OpacityChannel].depth++;
2479  continue;
2480  }
2481  }
2482  }
2483  if (image->colorspace == CMYKColorspace)
2484  {
2485  if (channel_statistics[BlackChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2486  {
2487  depth=channel_statistics[BlackChannel].depth;
2488  range=GetQuantumRange(depth);
2489  if (IsPixelAtDepth(GetPixelIndex(indexes+x),range) == MagickFalse)
2490  {
2491  channel_statistics[BlackChannel].depth++;
2492  continue;
2493  }
2494  }
2495  }
2496  if ((double) GetPixelRed(p) < channel_statistics[RedChannel].minima)
2497  channel_statistics[RedChannel].minima=(double) GetPixelRed(p);
2498  if ((double) GetPixelRed(p) > channel_statistics[RedChannel].maxima)
2499  channel_statistics[RedChannel].maxima=(double) GetPixelRed(p);
2500  channel_statistics[RedChannel].sum+=QuantumScale*GetPixelRed(p);
2501  channel_statistics[RedChannel].sum_squared+=QuantumScale*GetPixelRed(p)*
2502  QuantumScale*GetPixelRed(p);
2503  channel_statistics[RedChannel].sum_cubed+=QuantumScale*GetPixelRed(p)*
2504  QuantumScale*GetPixelRed(p)*QuantumScale*GetPixelRed(p);
2505  channel_statistics[RedChannel].sum_fourth_power+=QuantumScale*
2506  GetPixelRed(p)*QuantumScale*GetPixelRed(p)*QuantumScale*GetPixelRed(p)*
2507  QuantumScale*GetPixelRed(p);
2508  if ((double) GetPixelGreen(p) < channel_statistics[GreenChannel].minima)
2509  channel_statistics[GreenChannel].minima=(double) GetPixelGreen(p);
2510  if ((double) GetPixelGreen(p) > channel_statistics[GreenChannel].maxima)
2511  channel_statistics[GreenChannel].maxima=(double) GetPixelGreen(p);
2512  channel_statistics[GreenChannel].sum+=QuantumScale*GetPixelGreen(p);
2513  channel_statistics[GreenChannel].sum_squared+=QuantumScale*GetPixelGreen(p)*
2514  QuantumScale*GetPixelGreen(p);
2515  channel_statistics[GreenChannel].sum_cubed+=QuantumScale*GetPixelGreen(p)*
2516  QuantumScale*GetPixelGreen(p)*QuantumScale*GetPixelGreen(p);
2517  channel_statistics[GreenChannel].sum_fourth_power+=QuantumScale*
2518  GetPixelGreen(p)*QuantumScale*GetPixelGreen(p)*QuantumScale*
2519  GetPixelGreen(p)*QuantumScale*GetPixelGreen(p);
2520  if ((double) GetPixelBlue(p) < channel_statistics[BlueChannel].minima)
2521  channel_statistics[BlueChannel].minima=(double) GetPixelBlue(p);
2522  if ((double) GetPixelBlue(p) > channel_statistics[BlueChannel].maxima)
2523  channel_statistics[BlueChannel].maxima=(double) GetPixelBlue(p);
2524  channel_statistics[BlueChannel].sum+=QuantumScale*GetPixelBlue(p);
2525  channel_statistics[BlueChannel].sum_squared+=QuantumScale*GetPixelBlue(p)*
2526  QuantumScale*GetPixelBlue(p);
2527  channel_statistics[BlueChannel].sum_cubed+=QuantumScale*GetPixelBlue(p)*
2528  QuantumScale*GetPixelBlue(p)*QuantumScale*GetPixelBlue(p);
2529  channel_statistics[BlueChannel].sum_fourth_power+=QuantumScale*
2530  GetPixelBlue(p)*QuantumScale*GetPixelBlue(p)*QuantumScale*
2531  GetPixelBlue(p)*QuantumScale*GetPixelBlue(p);
2532  histogram[ScaleQuantumToMap(GetPixelRed(p))].red++;
2533  histogram[ScaleQuantumToMap(GetPixelGreen(p))].green++;
2534  histogram[ScaleQuantumToMap(GetPixelBlue(p))].blue++;
2535  if (image->matte != MagickFalse)
2536  {
2537  if ((double) GetPixelAlpha(p) < channel_statistics[OpacityChannel].minima)
2538  channel_statistics[OpacityChannel].minima=(double) GetPixelAlpha(p);
2539  if ((double) GetPixelAlpha(p) > channel_statistics[OpacityChannel].maxima)
2540  channel_statistics[OpacityChannel].maxima=(double) GetPixelAlpha(p);
2541  channel_statistics[OpacityChannel].sum+=QuantumScale*GetPixelAlpha(p);
2542  channel_statistics[OpacityChannel].sum_squared+=QuantumScale*
2543  GetPixelAlpha(p)*QuantumScale*GetPixelAlpha(p);
2544  channel_statistics[OpacityChannel].sum_cubed+=QuantumScale*
2545  GetPixelAlpha(p)*QuantumScale*GetPixelAlpha(p)*QuantumScale*
2546  GetPixelAlpha(p);
2547  channel_statistics[OpacityChannel].sum_fourth_power+=QuantumScale*
2548  GetPixelAlpha(p)*QuantumScale*GetPixelAlpha(p)*QuantumScale*
2549  GetPixelAlpha(p)*QuantumScale*GetPixelAlpha(p);
2550  histogram[ScaleQuantumToMap(GetPixelAlpha(p))].opacity++;
2551  }
2552  if (image->colorspace == CMYKColorspace)
2553  {
2554  if ((double) GetPixelIndex(indexes+x) < channel_statistics[BlackChannel].minima)
2555  channel_statistics[BlackChannel].minima=(double)
2556  GetPixelIndex(indexes+x);
2557  if ((double) GetPixelIndex(indexes+x) > channel_statistics[BlackChannel].maxima)
2558  channel_statistics[BlackChannel].maxima=(double)
2559  GetPixelIndex(indexes+x);
2560  channel_statistics[BlackChannel].sum+=QuantumScale*
2561  GetPixelIndex(indexes+x);
2562  channel_statistics[BlackChannel].sum_squared+=QuantumScale*
2563  GetPixelIndex(indexes+x)*QuantumScale*GetPixelIndex(indexes+x);
2564  channel_statistics[BlackChannel].sum_cubed+=QuantumScale*
2565  GetPixelIndex(indexes+x)*QuantumScale*GetPixelIndex(indexes+x)*
2566  QuantumScale*GetPixelIndex(indexes+x);
2567  channel_statistics[BlackChannel].sum_fourth_power+=QuantumScale*
2568  GetPixelIndex(indexes+x)*QuantumScale*GetPixelIndex(indexes+x)*
2569  QuantumScale*GetPixelIndex(indexes+x)*QuantumScale*
2570  GetPixelIndex(indexes+x);
2571  histogram[ScaleQuantumToMap(GetPixelIndex(indexes+x))].index++;
2572  }
2573  x++;
2574  p++;
2575  }
2576  }
2577  for (i=0; i < (ssize_t) CompositeChannels; i++)
2578  {
2579  double
2580  area,
2581  mean,
2582  standard_deviation;
2583 
2584  /*
2585  Normalize pixel statistics.
2586  */
2587  area=MagickSafeReciprocal((double) image->columns*image->rows);
2588  mean=channel_statistics[i].sum*area;
2589  channel_statistics[i].sum=mean;
2590  channel_statistics[i].sum_squared*=area;
2591  channel_statistics[i].sum_cubed*=area;
2592  channel_statistics[i].sum_fourth_power*=area;
2593  channel_statistics[i].mean=mean;
2594  channel_statistics[i].variance=channel_statistics[i].sum_squared;
2595  standard_deviation=sqrt(channel_statistics[i].variance-(mean*mean));
2596  area=MagickSafeReciprocal((double) image->columns*image->rows-1.0)*
2597  ((double) image->columns*image->rows);
2598  standard_deviation=sqrt(area*standard_deviation*standard_deviation);
2599  channel_statistics[i].standard_deviation=standard_deviation;
2600  }
2601  for (i=0; i < (ssize_t) (MaxMap+1U); i++)
2602  {
2603  if (histogram[i].red > 0.0)
2604  number_bins.red++;
2605  if (histogram[i].green > 0.0)
2606  number_bins.green++;
2607  if (histogram[i].blue > 0.0)
2608  number_bins.blue++;
2609  if ((image->matte != MagickFalse) && (histogram[i].opacity > 0.0))
2610  number_bins.opacity++;
2611  if ((image->colorspace == CMYKColorspace) && (histogram[i].index > 0.0))
2612  number_bins.index++;
2613  }
2614  area=MagickSafeReciprocal((double) image->columns*image->rows);
2615  for (i=0; i < (ssize_t) (MaxMap+1U); i++)
2616  {
2617  double
2618  entropy;
2619 
2620  /*
2621  Compute pixel entropy.
2622  */
2623  histogram[i].red*=area;
2624  entropy=-histogram[i].red*log2(histogram[i].red)*
2625  MagickSafeReciprocal(log2((double) number_bins.red));
2626  if (IsNaN(entropy) == 0)
2627  channel_statistics[RedChannel].entropy+=entropy;
2628  histogram[i].green*=area;
2629  entropy=-histogram[i].green*log2(histogram[i].green)*
2630  MagickSafeReciprocal(log2((double) number_bins.green));
2631  if (IsNaN(entropy) == 0)
2632  channel_statistics[GreenChannel].entropy+=entropy;
2633  histogram[i].blue*=area;
2634  entropy=-histogram[i].blue*log2(histogram[i].blue)*
2635  MagickSafeReciprocal(log2((double) number_bins.blue));
2636  if (IsNaN(entropy) == 0)
2637  channel_statistics[BlueChannel].entropy+=entropy;
2638  if (image->matte != MagickFalse)
2639  {
2640  histogram[i].opacity*=area;
2641  entropy=-histogram[i].opacity*log2(histogram[i].opacity)*
2642  MagickSafeReciprocal(log2((double) number_bins.opacity));
2643  if (IsNaN(entropy) == 0)
2644  channel_statistics[OpacityChannel].entropy+=entropy;
2645  }
2646  if (image->colorspace == CMYKColorspace)
2647  {
2648  histogram[i].index*=area;
2649  entropy=-histogram[i].index*log2(histogram[i].index)*
2650  MagickSafeReciprocal(log2((double) number_bins.index));
2651  if (IsNaN(entropy) == 0)
2652  channel_statistics[IndexChannel].entropy+=entropy;
2653  }
2654  }
2655  /*
2656  Compute overall statistics.
2657  */
2658  for (i=0; i < (ssize_t) CompositeChannels; i++)
2659  {
2660  channel_statistics[CompositeChannels].depth=(size_t) EvaluateMax((double)
2661  channel_statistics[CompositeChannels].depth,(double)
2662  channel_statistics[i].depth);
2663  channel_statistics[CompositeChannels].minima=MagickMin(
2664  channel_statistics[CompositeChannels].minima,
2665  channel_statistics[i].minima);
2666  channel_statistics[CompositeChannels].maxima=EvaluateMax(
2667  channel_statistics[CompositeChannels].maxima,
2668  channel_statistics[i].maxima);
2669  channel_statistics[CompositeChannels].sum+=channel_statistics[i].sum;
2670  channel_statistics[CompositeChannels].sum_squared+=
2671  channel_statistics[i].sum_squared;
2672  channel_statistics[CompositeChannels].sum_cubed+=
2673  channel_statistics[i].sum_cubed;
2674  channel_statistics[CompositeChannels].sum_fourth_power+=
2675  channel_statistics[i].sum_fourth_power;
2676  channel_statistics[CompositeChannels].mean+=channel_statistics[i].mean;
2677  channel_statistics[CompositeChannels].variance+=
2678  channel_statistics[i].variance-channel_statistics[i].mean*
2679  channel_statistics[i].mean;
2680  standard_deviation=sqrt(channel_statistics[i].variance-
2681  (channel_statistics[i].mean*channel_statistics[i].mean));
2682  area=MagickSafeReciprocal((double) image->columns*image->rows-1.0)*
2683  ((double) image->columns*image->rows);
2684  standard_deviation=sqrt(area*standard_deviation*standard_deviation);
2685  channel_statistics[CompositeChannels].standard_deviation=standard_deviation;
2686  channel_statistics[CompositeChannels].entropy+=
2687  channel_statistics[i].entropy;
2688  }
2689  channels=3;
2690  if (image->matte != MagickFalse)
2691  channels++;
2692  if (image->colorspace == CMYKColorspace)
2693  channels++;
2694  channel_statistics[CompositeChannels].sum/=channels;
2695  channel_statistics[CompositeChannels].sum_squared/=channels;
2696  channel_statistics[CompositeChannels].sum_cubed/=channels;
2697  channel_statistics[CompositeChannels].sum_fourth_power/=channels;
2698  channel_statistics[CompositeChannels].mean/=channels;
2699  channel_statistics[CompositeChannels].kurtosis/=channels;
2700  channel_statistics[CompositeChannels].skewness/=channels;
2701  channel_statistics[CompositeChannels].entropy/=channels;
2702  i=CompositeChannels;
2703  area=MagickSafeReciprocal((double) channels*image->columns*image->rows);
2704  channel_statistics[i].variance=channel_statistics[i].sum_squared;
2705  channel_statistics[i].mean=channel_statistics[i].sum;
2706  standard_deviation=sqrt(channel_statistics[i].variance-
2707  (channel_statistics[i].mean*channel_statistics[i].mean));
2708  standard_deviation=sqrt(MagickSafeReciprocal((double) channels*
2709  image->columns*image->rows-1.0)*channels*image->columns*image->rows*
2710  standard_deviation*standard_deviation);
2711  channel_statistics[i].standard_deviation=standard_deviation;
2712  for (i=0; i <= (ssize_t) CompositeChannels; i++)
2713  {
2714  /*
2715  Compute kurtosis & skewness statistics.
2716  */
2717  standard_deviation=MagickSafeReciprocal(
2718  channel_statistics[i].standard_deviation);
2719  channel_statistics[i].skewness=(channel_statistics[i].sum_cubed-3.0*
2720  channel_statistics[i].mean*channel_statistics[i].sum_squared+2.0*
2721  channel_statistics[i].mean*channel_statistics[i].mean*
2722  channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2723  standard_deviation);
2724  channel_statistics[i].kurtosis=(channel_statistics[i].sum_fourth_power-4.0*
2725  channel_statistics[i].mean*channel_statistics[i].sum_cubed+6.0*
2726  channel_statistics[i].mean*channel_statistics[i].mean*
2727  channel_statistics[i].sum_squared-3.0*channel_statistics[i].mean*
2728  channel_statistics[i].mean*1.0*channel_statistics[i].mean*
2729  channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2730  standard_deviation*standard_deviation)-3.0;
2731  }
2732  for (i=0; i <= (ssize_t) CompositeChannels; i++)
2733  {
2734  channel_statistics[i].mean*=QuantumRange;
2735  channel_statistics[i].variance*=QuantumRange;
2736  channel_statistics[i].standard_deviation*=QuantumRange;
2737  channel_statistics[i].sum*=QuantumRange;
2738  channel_statistics[i].sum_squared*=QuantumRange;
2739  channel_statistics[i].sum_cubed*=QuantumRange;
2740  channel_statistics[i].sum_fourth_power*=QuantumRange;
2741  }
2742  channel_statistics[CompositeChannels].mean=0.0;
2743  channel_statistics[CompositeChannels].standard_deviation=0.0;
2744  for (i=0; i < (ssize_t) CompositeChannels; i++)
2745  {
2746  channel_statistics[CompositeChannels].mean+=
2747  channel_statistics[i].mean;
2748  channel_statistics[CompositeChannels].standard_deviation+=
2749  channel_statistics[i].standard_deviation;
2750  }
2751  channel_statistics[CompositeChannels].mean/=(double) channels;
2752  channel_statistics[CompositeChannels].standard_deviation/=(double) channels;
2753  histogram=(MagickPixelPacket *) RelinquishMagickMemory(histogram);
2754  if (y < (ssize_t) image->rows)
2755  channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
2756  channel_statistics);
2757  return(channel_statistics);
2758 }
2759 
2760 /*
2761 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2762 % %
2763 % %
2764 % %
2765 % P o l y n o m i a l I m a g e %
2766 % %
2767 % %
2768 % %
2769 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2770 %
2771 % PolynomialImage() returns a new image where each pixel is the sum of the
2772 % pixels in the image sequence after applying its corresponding terms
2773 % (coefficient and degree pairs).
2774 %
2775 % The format of the PolynomialImage method is:
2776 %
2777 % Image *PolynomialImage(const Image *images,const size_t number_terms,
2778 % const double *terms,ExceptionInfo *exception)
2779 % Image *PolynomialImageChannel(const Image *images,
2780 % const size_t number_terms,const ChannelType channel,
2781 % const double *terms,ExceptionInfo *exception)
2782 %
2783 % A description of each parameter follows:
2784 %
2785 % o images: the image sequence.
2786 %
2787 % o channel: the channel.
2788 %
2789 % o number_terms: the number of terms in the list. The actual list length
2790 % is 2 x number_terms + 1 (the constant).
2791 %
2792 % o terms: the list of polynomial coefficients and degree pairs and a
2793 % constant.
2794 %
2795 % o exception: return any errors or warnings in this structure.
2796 %
2797 */
2798 MagickExport Image *PolynomialImage(const Image *images,
2799  const size_t number_terms,const double *terms,ExceptionInfo *exception)
2800 {
2801  Image
2802  *polynomial_image;
2803 
2804  polynomial_image=PolynomialImageChannel(images,DefaultChannels,number_terms,
2805  terms,exception);
2806  return(polynomial_image);
2807 }
2808 
2809 MagickExport Image *PolynomialImageChannel(const Image *images,
2810  const ChannelType channel,const size_t number_terms,const double *terms,
2811  ExceptionInfo *exception)
2812 {
2813 #define PolynomialImageTag "Polynomial/Image"
2814 
2815  CacheView
2816  *polynomial_view;
2817 
2818  Image
2819  *image;
2820 
2821  MagickBooleanType
2822  status;
2823 
2824  MagickOffsetType
2825  progress;
2826 
2827  MagickPixelPacket
2828  **magick_restrict polynomial_pixels,
2829  zero;
2830 
2831  ssize_t
2832  y;
2833 
2834  assert(images != (Image *) NULL);
2835  assert(images->signature == MagickCoreSignature);
2836  if (IsEventLogging() != MagickFalse)
2837  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
2838  assert(exception != (ExceptionInfo *) NULL);
2839  assert(exception->signature == MagickCoreSignature);
2840  image=AcquireImageCanvas(images,exception);
2841  if (image == (Image *) NULL)
2842  return((Image *) NULL);
2843  if (SetImageStorageClass(image,DirectClass) == MagickFalse)
2844  {
2845  InheritException(exception,&image->exception);
2846  image=DestroyImage(image);
2847  return((Image *) NULL);
2848  }
2849  polynomial_pixels=AcquirePixelTLS(images);
2850  if (polynomial_pixels == (MagickPixelPacket **) NULL)
2851  {
2852  image=DestroyImage(image);
2853  (void) ThrowMagickException(exception,GetMagickModule(),
2854  ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
2855  return((Image *) NULL);
2856  }
2857  /*
2858  Polynomial image pixels.
2859  */
2860  status=MagickTrue;
2861  progress=0;
2862  GetMagickPixelPacket(images,&zero);
2863  polynomial_view=AcquireAuthenticCacheView(image,exception);
2864 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2865  #pragma omp parallel for schedule(static) shared(progress,status) \
2866  magick_number_threads(image,image,image->rows,1)
2867 #endif
2868  for (y=0; y < (ssize_t) image->rows; y++)
2869  {
2870  CacheView
2871  *image_view;
2872 
2873  const Image
2874  *next;
2875 
2876  const int
2877  id = GetOpenMPThreadId();
2878 
2879  IndexPacket
2880  *magick_restrict polynomial_indexes;
2881 
2882  MagickPixelPacket
2883  *polynomial_pixel;
2884 
2885  PixelPacket
2886  *magick_restrict q;
2887 
2888  ssize_t
2889  i,
2890  x;
2891 
2892  size_t
2893  number_images;
2894 
2895  if (status == MagickFalse)
2896  continue;
2897  q=QueueCacheViewAuthenticPixels(polynomial_view,0,y,image->columns,1,
2898  exception);
2899  if (q == (PixelPacket *) NULL)
2900  {
2901  status=MagickFalse;
2902  continue;
2903  }
2904  polynomial_indexes=GetCacheViewAuthenticIndexQueue(polynomial_view);
2905  polynomial_pixel=polynomial_pixels[id];
2906  for (x=0; x < (ssize_t) image->columns; x++)
2907  polynomial_pixel[x]=zero;
2908  next=images;
2909  number_images=GetImageListLength(images);
2910  for (i=0; i < (ssize_t) number_images; i++)
2911  {
2912  const IndexPacket
2913  *indexes;
2914 
2915  const PixelPacket
2916  *p;
2917 
2918  if (i >= (ssize_t) number_terms)
2919  break;
2920  image_view=AcquireVirtualCacheView(next,exception);
2921  p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2922  if (p == (const PixelPacket *) NULL)
2923  {
2924  image_view=DestroyCacheView(image_view);
2925  break;
2926  }
2927  indexes=GetCacheViewVirtualIndexQueue(image_view);
2928  for (x=0; x < (ssize_t) image->columns; x++)
2929  {
2930  double
2931  coefficient,
2932  degree;
2933 
2934  coefficient=terms[i << 1];
2935  degree=terms[(i << 1)+1];
2936  if ((channel & RedChannel) != 0)
2937  polynomial_pixel[x].red+=coefficient*pow(QuantumScale*(double)
2938  p->red,degree);
2939  if ((channel & GreenChannel) != 0)
2940  polynomial_pixel[x].green+=coefficient*pow(QuantumScale*(double)
2941  p->green,
2942  degree);
2943  if ((channel & BlueChannel) != 0)
2944  polynomial_pixel[x].blue+=coefficient*pow(QuantumScale*(double)
2945  p->blue,degree);
2946  if ((channel & OpacityChannel) != 0)
2947  polynomial_pixel[x].opacity+=coefficient*pow(QuantumScale*
2948  ((double) QuantumRange-(double) p->opacity),degree);
2949  if (((channel & IndexChannel) != 0) &&
2950  (image->colorspace == CMYKColorspace))
2951  polynomial_pixel[x].index+=coefficient*pow(QuantumScale*(double)
2952  indexes[x],degree);
2953  p++;
2954  }
2955  image_view=DestroyCacheView(image_view);
2956  next=GetNextImageInList(next);
2957  }
2958  for (x=0; x < (ssize_t) image->columns; x++)
2959  {
2960  SetPixelRed(q,ClampToQuantum((MagickRealType) QuantumRange*
2961  polynomial_pixel[x].red));
2962  SetPixelGreen(q,ClampToQuantum((MagickRealType) QuantumRange*
2963  polynomial_pixel[x].green));
2964  SetPixelBlue(q,ClampToQuantum((MagickRealType) QuantumRange*
2965  polynomial_pixel[x].blue));
2966  if (image->matte == MagickFalse)
2967  SetPixelOpacity(q,ClampToQuantum((MagickRealType) QuantumRange-
2968  (MagickRealType) QuantumRange*polynomial_pixel[x].opacity));
2969  else
2970  SetPixelAlpha(q,ClampToQuantum((MagickRealType) QuantumRange-
2971  (MagickRealType) QuantumRange*polynomial_pixel[x].opacity));
2972  if (image->colorspace == CMYKColorspace)
2973  SetPixelIndex(polynomial_indexes+x,ClampToQuantum((MagickRealType)
2974  QuantumRange*polynomial_pixel[x].index));
2975  q++;
2976  }
2977  if (SyncCacheViewAuthenticPixels(polynomial_view,exception) == MagickFalse)
2978  status=MagickFalse;
2979  if (images->progress_monitor != (MagickProgressMonitor) NULL)
2980  {
2981  MagickBooleanType
2982  proceed;
2983 
2984  proceed=SetImageProgress(images,PolynomialImageTag,progress++,
2985  image->rows);
2986  if (proceed == MagickFalse)
2987  status=MagickFalse;
2988  }
2989  }
2990  polynomial_view=DestroyCacheView(polynomial_view);
2991  polynomial_pixels=DestroyPixelTLS(images,polynomial_pixels);
2992  if (status == MagickFalse)
2993  image=DestroyImage(image);
2994  return(image);
2995 }
2996 
2997 /*
2998 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2999 % %
3000 % %
3001 % %
3002 % S t a t i s t i c I m a g e %
3003 % %
3004 % %
3005 % %
3006 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3007 %
3008 % StatisticImage() makes each pixel the min / max / median / mode / etc. of
3009 % the neighborhood of the specified width and height.
3010 %
3011 % The format of the StatisticImage method is:
3012 %
3013 % Image *StatisticImage(const Image *image,const StatisticType type,
3014 % const size_t width,const size_t height,ExceptionInfo *exception)
3015 % Image *StatisticImageChannel(const Image *image,
3016 % const ChannelType channel,const StatisticType type,
3017 % const size_t width,const size_t height,ExceptionInfo *exception)
3018 %
3019 % A description of each parameter follows:
3020 %
3021 % o image: the image.
3022 %
3023 % o channel: the image channel.
3024 %
3025 % o type: the statistic type (median, mode, etc.).
3026 %
3027 % o width: the width of the pixel neighborhood.
3028 %
3029 % o height: the height of the pixel neighborhood.
3030 %
3031 % o exception: return any errors or warnings in this structure.
3032 %
3033 */
3034 
3035 #define ListChannels 5
3036 
3037 typedef struct _ListNode
3038 {
3039  size_t
3040  next[9],
3041  count,
3042  signature;
3043 } ListNode;
3044 
3045 typedef struct _SkipList
3046 {
3047  ssize_t
3048  level;
3049 
3050  ListNode
3051  *nodes;
3052 } SkipList;
3053 
3054 typedef struct _PixelList
3055 {
3056  size_t
3057  length,
3058  seed,
3059  signature;
3060 
3061  SkipList
3062  lists[ListChannels];
3063 } PixelList;
3064 
3065 static PixelList *DestroyPixelList(PixelList *pixel_list)
3066 {
3067  ssize_t
3068  i;
3069 
3070  if (pixel_list == (PixelList *) NULL)
3071  return((PixelList *) NULL);
3072  for (i=0; i < ListChannels; i++)
3073  if (pixel_list->lists[i].nodes != (ListNode *) NULL)
3074  pixel_list->lists[i].nodes=(ListNode *) RelinquishAlignedMemory(
3075  pixel_list->lists[i].nodes);
3076  pixel_list=(PixelList *) RelinquishMagickMemory(pixel_list);
3077  return(pixel_list);
3078 }
3079 
3080 static PixelList **DestroyPixelListTLS(PixelList **pixel_list)
3081 {
3082  ssize_t
3083  i;
3084 
3085  assert(pixel_list != (PixelList **) NULL);
3086  for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
3087  if (pixel_list[i] != (PixelList *) NULL)
3088  pixel_list[i]=DestroyPixelList(pixel_list[i]);
3089  pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list);
3090  return(pixel_list);
3091 }
3092 
3093 static PixelList *AcquirePixelList(const size_t width,const size_t height)
3094 {
3095  PixelList
3096  *pixel_list;
3097 
3098  ssize_t
3099  i;
3100 
3101  pixel_list=(PixelList *) AcquireMagickMemory(sizeof(*pixel_list));
3102  if (pixel_list == (PixelList *) NULL)
3103  return(pixel_list);
3104  (void) memset((void *) pixel_list,0,sizeof(*pixel_list));
3105  pixel_list->length=width*height;
3106  for (i=0; i < ListChannels; i++)
3107  {
3108  pixel_list->lists[i].nodes=(ListNode *) AcquireAlignedMemory(65537UL,
3109  sizeof(*pixel_list->lists[i].nodes));
3110  if (pixel_list->lists[i].nodes == (ListNode *) NULL)
3111  return(DestroyPixelList(pixel_list));
3112  (void) memset(pixel_list->lists[i].nodes,0,65537UL*
3113  sizeof(*pixel_list->lists[i].nodes));
3114  }
3115  pixel_list->signature=MagickCoreSignature;
3116  return(pixel_list);
3117 }
3118 
3119 static PixelList **AcquirePixelListTLS(const size_t width,const size_t height)
3120 {
3121  PixelList
3122  **pixel_list;
3123 
3124  ssize_t
3125  i;
3126 
3127  size_t
3128  number_threads;
3129 
3130  number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
3131  pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
3132  sizeof(*pixel_list));
3133  if (pixel_list == (PixelList **) NULL)
3134  return((PixelList **) NULL);
3135  (void) memset(pixel_list,0,number_threads*sizeof(*pixel_list));
3136  for (i=0; i < (ssize_t) number_threads; i++)
3137  {
3138  pixel_list[i]=AcquirePixelList(width,height);
3139  if (pixel_list[i] == (PixelList *) NULL)
3140  return(DestroyPixelListTLS(pixel_list));
3141  }
3142  return(pixel_list);
3143 }
3144 
3145 static void AddNodePixelList(PixelList *pixel_list,const ssize_t channel,
3146  const size_t color)
3147 {
3148  SkipList
3149  *list;
3150 
3151  ssize_t
3152  level;
3153 
3154  size_t
3155  search,
3156  update[9];
3157 
3158  /*
3159  Initialize the node.
3160  */
3161  list=pixel_list->lists+channel;
3162  list->nodes[color].signature=pixel_list->signature;
3163  list->nodes[color].count=1;
3164  /*
3165  Determine where it belongs in the list.
3166  */
3167  search=65536UL;
3168  (void) memset(update,0,sizeof(update));
3169  for (level=list->level; level >= 0; level--)
3170  {
3171  while (list->nodes[search].next[level] < color)
3172  search=list->nodes[search].next[level];
3173  update[level]=search;
3174  }
3175  /*
3176  Generate a pseudo-random level for this node.
3177  */
3178  for (level=0; ; level++)
3179  {
3180  pixel_list->seed=(pixel_list->seed*42893621L)+1L;
3181  if ((pixel_list->seed & 0x300) != 0x300)
3182  break;
3183  }
3184  if (level > 8)
3185  level=8;
3186  if (level > (list->level+2))
3187  level=list->level+2;
3188  /*
3189  If we're raising the list's level, link back to the root node.
3190  */
3191  while (level > list->level)
3192  {
3193  list->level++;
3194  update[list->level]=65536UL;
3195  }
3196  /*
3197  Link the node into the skip-list.
3198  */
3199  do
3200  {
3201  list->nodes[color].next[level]=list->nodes[update[level]].next[level];
3202  list->nodes[update[level]].next[level]=color;
3203  } while (level-- > 0);
3204 }
3205 
3206 static void GetMaximumPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3207 {
3208  SkipList
3209  *list;
3210 
3211  ssize_t
3212  channel;
3213 
3214  size_t
3215  color,
3216  maximum;
3217 
3218  ssize_t
3219  count;
3220 
3221  unsigned short
3222  channels[ListChannels];
3223 
3224  /*
3225  Find the maximum value for each of the color.
3226  */
3227  for (channel=0; channel < 5; channel++)
3228  {
3229  list=pixel_list->lists+channel;
3230  color=65536L;
3231  count=0;
3232  maximum=list->nodes[color].next[0];
3233  do
3234  {
3235  color=list->nodes[color].next[0];
3236  if (color > maximum)
3237  maximum=color;
3238  count+=list->nodes[color].count;
3239  } while (count < (ssize_t) pixel_list->length);
3240  channels[channel]=(unsigned short) maximum;
3241  }
3242  pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3243  pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3244  pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3245  pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3246  pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3247 }
3248 
3249 static void GetMeanPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3250 {
3251  MagickRealType
3252  sum;
3253 
3254  SkipList
3255  *list;
3256 
3257  ssize_t
3258  channel;
3259 
3260  size_t
3261  color;
3262 
3263  ssize_t
3264  count;
3265 
3266  unsigned short
3267  channels[ListChannels];
3268 
3269  /*
3270  Find the mean value for each of the color.
3271  */
3272  for (channel=0; channel < 5; channel++)
3273  {
3274  list=pixel_list->lists+channel;
3275  color=65536L;
3276  count=0;
3277  sum=0.0;
3278  do
3279  {
3280  color=list->nodes[color].next[0];
3281  sum+=(MagickRealType) list->nodes[color].count*color;
3282  count+=list->nodes[color].count;
3283  } while (count < (ssize_t) pixel_list->length);
3284  sum/=pixel_list->length;
3285  channels[channel]=(unsigned short) sum;
3286  }
3287  pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3288  pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3289  pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3290  pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3291  pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3292 }
3293 
3294 static void GetMedianPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3295 {
3296  SkipList
3297  *list;
3298 
3299  ssize_t
3300  channel;
3301 
3302  size_t
3303  color;
3304 
3305  ssize_t
3306  count;
3307 
3308  unsigned short
3309  channels[ListChannels];
3310 
3311  /*
3312  Find the median value for each of the color.
3313  */
3314  for (channel=0; channel < 5; channel++)
3315  {
3316  list=pixel_list->lists+channel;
3317  color=65536L;
3318  count=0;
3319  do
3320  {
3321  color=list->nodes[color].next[0];
3322  count+=list->nodes[color].count;
3323  } while (count <= (ssize_t) (pixel_list->length >> 1));
3324  channels[channel]=(unsigned short) color;
3325  }
3326  GetMagickPixelPacket((const Image *) NULL,pixel);
3327  pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3328  pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3329  pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3330  pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3331  pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3332 }
3333 
3334 static void GetMinimumPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3335 {
3336  SkipList
3337  *list;
3338 
3339  ssize_t
3340  channel;
3341 
3342  size_t
3343  color,
3344  minimum;
3345 
3346  ssize_t
3347  count;
3348 
3349  unsigned short
3350  channels[ListChannels];
3351 
3352  /*
3353  Find the minimum value for each of the color.
3354  */
3355  for (channel=0; channel < 5; channel++)
3356  {
3357  list=pixel_list->lists+channel;
3358  count=0;
3359  color=65536UL;
3360  minimum=list->nodes[color].next[0];
3361  do
3362  {
3363  color=list->nodes[color].next[0];
3364  if (color < minimum)
3365  minimum=color;
3366  count+=list->nodes[color].count;
3367  } while (count < (ssize_t) pixel_list->length);
3368  channels[channel]=(unsigned short) minimum;
3369  }
3370  pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3371  pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3372  pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3373  pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3374  pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3375 }
3376 
3377 static void GetModePixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3378 {
3379  SkipList
3380  *list;
3381 
3382  ssize_t
3383  channel;
3384 
3385  size_t
3386  color,
3387  max_count,
3388  mode;
3389 
3390  ssize_t
3391  count;
3392 
3393  unsigned short
3394  channels[5];
3395 
3396  /*
3397  Make each pixel the 'predominant color' of the specified neighborhood.
3398  */
3399  for (channel=0; channel < 5; channel++)
3400  {
3401  list=pixel_list->lists+channel;
3402  color=65536L;
3403  mode=color;
3404  max_count=list->nodes[mode].count;
3405  count=0;
3406  do
3407  {
3408  color=list->nodes[color].next[0];
3409  if (list->nodes[color].count > max_count)
3410  {
3411  mode=color;
3412  max_count=list->nodes[mode].count;
3413  }
3414  count+=list->nodes[color].count;
3415  } while (count < (ssize_t) pixel_list->length);
3416  channels[channel]=(unsigned short) mode;
3417  }
3418  pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3419  pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3420  pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3421  pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3422  pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3423 }
3424 
3425 static void GetNonpeakPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3426 {
3427  SkipList
3428  *list;
3429 
3430  ssize_t
3431  channel;
3432 
3433  size_t
3434  color,
3435  next,
3436  previous;
3437 
3438  ssize_t
3439  count;
3440 
3441  unsigned short
3442  channels[5];
3443 
3444  /*
3445  Finds the non peak value for each of the colors.
3446  */
3447  for (channel=0; channel < 5; channel++)
3448  {
3449  list=pixel_list->lists+channel;
3450  color=65536L;
3451  next=list->nodes[color].next[0];
3452  count=0;
3453  do
3454  {
3455  previous=color;
3456  color=next;
3457  next=list->nodes[color].next[0];
3458  count+=list->nodes[color].count;
3459  } while (count <= (ssize_t) (pixel_list->length >> 1));
3460  if ((previous == 65536UL) && (next != 65536UL))
3461  color=next;
3462  else
3463  if ((previous != 65536UL) && (next == 65536UL))
3464  color=previous;
3465  channels[channel]=(unsigned short) color;
3466  }
3467  pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3468  pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3469  pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3470  pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3471  pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3472 }
3473 
3474 static void GetRootMeanSquarePixelList(PixelList *pixel_list,
3475  MagickPixelPacket *pixel)
3476 {
3477  MagickRealType
3478  sum;
3479 
3480  SkipList
3481  *list;
3482 
3483  ssize_t
3484  channel;
3485 
3486  size_t
3487  color;
3488 
3489  ssize_t
3490  count;
3491 
3492  unsigned short
3493  channels[ListChannels];
3494 
3495  /*
3496  Find the root mean square value for each of the color.
3497  */
3498  for (channel=0; channel < 5; channel++)
3499  {
3500  list=pixel_list->lists+channel;
3501  color=65536L;
3502  count=0;
3503  sum=0.0;
3504  do
3505  {
3506  color=list->nodes[color].next[0];
3507  sum+=(MagickRealType) (list->nodes[color].count*color*color);
3508  count+=list->nodes[color].count;
3509  } while (count < (ssize_t) pixel_list->length);
3510  sum/=pixel_list->length;
3511  channels[channel]=(unsigned short) sqrt(sum);
3512  }
3513  pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3514  pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3515  pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3516  pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3517  pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3518 }
3519 
3520 static void GetStandardDeviationPixelList(PixelList *pixel_list,
3521  MagickPixelPacket *pixel)
3522 {
3523  MagickRealType
3524  sum,
3525  sum_squared;
3526 
3527  SkipList
3528  *list;
3529 
3530  size_t
3531  color;
3532 
3533  ssize_t
3534  channel,
3535  count;
3536 
3537  unsigned short
3538  channels[ListChannels];
3539 
3540  /*
3541  Find the standard-deviation value for each of the color.
3542  */
3543  for (channel=0; channel < 5; channel++)
3544  {
3545  list=pixel_list->lists+channel;
3546  color=65536L;
3547  count=0;
3548  sum=0.0;
3549  sum_squared=0.0;
3550  do
3551  {
3552  ssize_t
3553  i;
3554 
3555  color=list->nodes[color].next[0];
3556  sum+=(MagickRealType) list->nodes[color].count*color;
3557  for (i=0; i < (ssize_t) list->nodes[color].count; i++)
3558  sum_squared+=((MagickRealType) color)*((MagickRealType) color);
3559  count+=list->nodes[color].count;
3560  } while (count < (ssize_t) pixel_list->length);
3561  sum/=pixel_list->length;
3562  sum_squared/=pixel_list->length;
3563  channels[channel]=(unsigned short) sqrt(sum_squared-(sum*sum));
3564  }
3565  pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3566  pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3567  pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3568  pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3569  pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3570 }
3571 
3572 static inline void InsertPixelList(const Image *image,const PixelPacket *pixel,
3573  const IndexPacket *indexes,PixelList *pixel_list)
3574 {
3575  size_t
3576  signature;
3577 
3578  unsigned short
3579  index;
3580 
3581  index=ScaleQuantumToShort(GetPixelRed(pixel));
3582  signature=pixel_list->lists[0].nodes[index].signature;
3583  if (signature == pixel_list->signature)
3584  pixel_list->lists[0].nodes[index].count++;
3585  else
3586  AddNodePixelList(pixel_list,0,index);
3587  index=ScaleQuantumToShort(GetPixelGreen(pixel));
3588  signature=pixel_list->lists[1].nodes[index].signature;
3589  if (signature == pixel_list->signature)
3590  pixel_list->lists[1].nodes[index].count++;
3591  else
3592  AddNodePixelList(pixel_list,1,index);
3593  index=ScaleQuantumToShort(GetPixelBlue(pixel));
3594  signature=pixel_list->lists[2].nodes[index].signature;
3595  if (signature == pixel_list->signature)
3596  pixel_list->lists[2].nodes[index].count++;
3597  else
3598  AddNodePixelList(pixel_list,2,index);
3599  index=ScaleQuantumToShort(GetPixelOpacity(pixel));
3600  signature=pixel_list->lists[3].nodes[index].signature;
3601  if (signature == pixel_list->signature)
3602  pixel_list->lists[3].nodes[index].count++;
3603  else
3604  AddNodePixelList(pixel_list,3,index);
3605  if (image->colorspace == CMYKColorspace)
3606  index=ScaleQuantumToShort(GetPixelIndex(indexes));
3607  signature=pixel_list->lists[4].nodes[index].signature;
3608  if (signature == pixel_list->signature)
3609  pixel_list->lists[4].nodes[index].count++;
3610  else
3611  AddNodePixelList(pixel_list,4,index);
3612 }
3613 
3614 static void ResetPixelList(PixelList *pixel_list)
3615 {
3616  int
3617  level;
3618 
3619  ListNode
3620  *root;
3621 
3622  SkipList
3623  *list;
3624 
3625  ssize_t
3626  channel;
3627 
3628  /*
3629  Reset the skip-list.
3630  */
3631  for (channel=0; channel < 5; channel++)
3632  {
3633  list=pixel_list->lists+channel;
3634  root=list->nodes+65536UL;
3635  list->level=0;
3636  for (level=0; level < 9; level++)
3637  root->next[level]=65536UL;
3638  }
3639  pixel_list->seed=pixel_list->signature++;
3640 }
3641 
3642 MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
3643  const size_t width,const size_t height,ExceptionInfo *exception)
3644 {
3645  Image
3646  *statistic_image;
3647 
3648  statistic_image=StatisticImageChannel(image,DefaultChannels,type,width,
3649  height,exception);
3650  return(statistic_image);
3651 }
3652 
3653 MagickExport Image *StatisticImageChannel(const Image *image,
3654  const ChannelType channel,const StatisticType type,const size_t width,
3655  const size_t height,ExceptionInfo *exception)
3656 {
3657 #define StatisticImageTag "Statistic/Image"
3658 
3659  CacheView
3660  *image_view,
3661  *statistic_view;
3662 
3663  Image
3664  *statistic_image;
3665 
3666  MagickBooleanType
3667  status;
3668 
3669  MagickOffsetType
3670  progress;
3671 
3672  PixelList
3673  **magick_restrict pixel_list;
3674 
3675  size_t
3676  neighbor_height,
3677  neighbor_width;
3678 
3679  ssize_t
3680  y;
3681 
3682  /*
3683  Initialize statistics image attributes.
3684  */
3685  assert(image != (Image *) NULL);
3686  assert(image->signature == MagickCoreSignature);
3687  if (IsEventLogging() != MagickFalse)
3688  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3689  assert(exception != (ExceptionInfo *) NULL);
3690  assert(exception->signature == MagickCoreSignature);
3691  statistic_image=CloneImage(image,0,0,MagickTrue,exception);
3692  if (statistic_image == (Image *) NULL)
3693  return((Image *) NULL);
3694  if (SetImageStorageClass(statistic_image,DirectClass) == MagickFalse)
3695  {
3696  InheritException(exception,&statistic_image->exception);
3697  statistic_image=DestroyImage(statistic_image);
3698  return((Image *) NULL);
3699  }
3700  neighbor_width=width == 0 ? GetOptimalKernelWidth2D((double) width,0.5) :
3701  width;
3702  neighbor_height=height == 0 ? GetOptimalKernelWidth2D((double) height,0.5) :
3703  height;
3704  pixel_list=AcquirePixelListTLS(neighbor_width,neighbor_height);
3705  if (pixel_list == (PixelList **) NULL)
3706  {
3707  statistic_image=DestroyImage(statistic_image);
3708  ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3709  }
3710  /*
3711  Make each pixel the min / max / median / mode / etc. of the neighborhood.
3712  */
3713  status=MagickTrue;
3714  progress=0;
3715  image_view=AcquireVirtualCacheView(image,exception);
3716  statistic_view=AcquireAuthenticCacheView(statistic_image,exception);
3717 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3718  #pragma omp parallel for schedule(static) shared(progress,status) \
3719  magick_number_threads(image,statistic_image,statistic_image->rows,1)
3720 #endif
3721  for (y=0; y < (ssize_t) statistic_image->rows; y++)
3722  {
3723  const int
3724  id = GetOpenMPThreadId();
3725 
3726  const IndexPacket
3727  *magick_restrict indexes;
3728 
3729  const PixelPacket
3730  *magick_restrict p;
3731 
3732  IndexPacket
3733  *magick_restrict statistic_indexes;
3734 
3735  PixelPacket
3736  *magick_restrict q;
3737 
3738  ssize_t
3739  x;
3740 
3741  if (status == MagickFalse)
3742  continue;
3743  p=GetCacheViewVirtualPixels(image_view,-((ssize_t) neighbor_width/2L),y-
3744  (ssize_t) (neighbor_height/2L),image->columns+neighbor_width,
3745  neighbor_height,exception);
3746  q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns, 1,exception);
3747  if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
3748  {
3749  status=MagickFalse;
3750  continue;
3751  }
3752  indexes=GetCacheViewVirtualIndexQueue(image_view);
3753  statistic_indexes=GetCacheViewAuthenticIndexQueue(statistic_view);
3754  for (x=0; x < (ssize_t) statistic_image->columns; x++)
3755  {
3756  MagickPixelPacket
3757  pixel;
3758 
3759  const IndexPacket
3760  *magick_restrict s;
3761 
3762  const PixelPacket
3763  *magick_restrict r;
3764 
3765  ssize_t
3766  u,
3767  v;
3768 
3769  r=p;
3770  s=indexes+x;
3771  ResetPixelList(pixel_list[id]);
3772  for (v=0; v < (ssize_t) neighbor_height; v++)
3773  {
3774  for (u=0; u < (ssize_t) neighbor_width; u++)
3775  InsertPixelList(image,r+u,s+u,pixel_list[id]);
3776  r+=(ptrdiff_t) image->columns+neighbor_width;
3777  s+=(ptrdiff_t) image->columns+neighbor_width;
3778  }
3779  GetMagickPixelPacket(image,&pixel);
3780  SetMagickPixelPacket(image,p+neighbor_width*neighbor_height/2,indexes+x+
3781  neighbor_width*neighbor_height/2,&pixel);
3782  switch (type)
3783  {
3784  case GradientStatistic:
3785  {
3786  MagickPixelPacket
3787  maximum,
3788  minimum;
3789 
3790  GetMinimumPixelList(pixel_list[id],&pixel);
3791  minimum=pixel;
3792  GetMaximumPixelList(pixel_list[id],&pixel);
3793  maximum=pixel;
3794  pixel.red=MagickAbsoluteValue(maximum.red-minimum.red);
3795  pixel.green=MagickAbsoluteValue(maximum.green-minimum.green);
3796  pixel.blue=MagickAbsoluteValue(maximum.blue-minimum.blue);
3797  pixel.opacity=MagickAbsoluteValue(maximum.opacity-minimum.opacity);
3798  if (image->colorspace == CMYKColorspace)
3799  pixel.index=MagickAbsoluteValue(maximum.index-minimum.index);
3800  break;
3801  }
3802  case MaximumStatistic:
3803  {
3804  GetMaximumPixelList(pixel_list[id],&pixel);
3805  break;
3806  }
3807  case MeanStatistic:
3808  {
3809  GetMeanPixelList(pixel_list[id],&pixel);
3810  break;
3811  }
3812  case MedianStatistic:
3813  default:
3814  {
3815  GetMedianPixelList(pixel_list[id],&pixel);
3816  break;
3817  }
3818  case MinimumStatistic:
3819  {
3820  GetMinimumPixelList(pixel_list[id],&pixel);
3821  break;
3822  }
3823  case ModeStatistic:
3824  {
3825  GetModePixelList(pixel_list[id],&pixel);
3826  break;
3827  }
3828  case NonpeakStatistic:
3829  {
3830  GetNonpeakPixelList(pixel_list[id],&pixel);
3831  break;
3832  }
3833  case RootMeanSquareStatistic:
3834  {
3835  GetRootMeanSquarePixelList(pixel_list[id],&pixel);
3836  break;
3837  }
3838  case StandardDeviationStatistic:
3839  {
3840  GetStandardDeviationPixelList(pixel_list[id],&pixel);
3841  break;
3842  }
3843  }
3844  if ((channel & RedChannel) != 0)
3845  SetPixelRed(q,ClampToQuantum(pixel.red));
3846  if ((channel & GreenChannel) != 0)
3847  SetPixelGreen(q,ClampToQuantum(pixel.green));
3848  if ((channel & BlueChannel) != 0)
3849  SetPixelBlue(q,ClampToQuantum(pixel.blue));
3850  if ((channel & OpacityChannel) != 0)
3851  SetPixelOpacity(q,ClampToQuantum(pixel.opacity));
3852  if (((channel & IndexChannel) != 0) &&
3853  (image->colorspace == CMYKColorspace))
3854  SetPixelIndex(statistic_indexes+x,ClampToQuantum(pixel.index));
3855  p++;
3856  q++;
3857  }
3858  if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
3859  status=MagickFalse;
3860  if (image->progress_monitor != (MagickProgressMonitor) NULL)
3861  {
3862  MagickBooleanType
3863  proceed;
3864 
3865  proceed=SetImageProgress(image,StatisticImageTag,progress++,
3866  image->rows);
3867  if (proceed == MagickFalse)
3868  status=MagickFalse;
3869  }
3870  }
3871  statistic_view=DestroyCacheView(statistic_view);
3872  image_view=DestroyCacheView(image_view);
3873  pixel_list=DestroyPixelListTLS(pixel_list);
3874  if (status == MagickFalse)
3875  statistic_image=DestroyImage(statistic_image);
3876  return(statistic_image);
3877 }
_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:3037
_PixelList
Definition: statistic.c:3054
_ChannelPerceptualHash
Definition: statistic.h:67
_ChannelStatistics
Definition: statistic.h:30
_PointInfo
Definition: draw.h:139
_SkipList
Definition: statistic.c:3045