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