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