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