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