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