MagickCore  7.0.8
Convert, Edit, Or Compose Bitmap Images
statistic.c
Go to the documentation of this file.
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-2018 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://www.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 "MagickCore/studio.h"
45 #include "MagickCore/animate.h"
46 #include "MagickCore/artifact.h"
47 #include "MagickCore/blob.h"
49 #include "MagickCore/cache.h"
51 #include "MagickCore/cache-view.h"
52 #include "MagickCore/client.h"
53 #include "MagickCore/color.h"
55 #include "MagickCore/colorspace.h"
57 #include "MagickCore/composite.h"
59 #include "MagickCore/compress.h"
60 #include "MagickCore/constitute.h"
61 #include "MagickCore/display.h"
62 #include "MagickCore/draw.h"
63 #include "MagickCore/enhance.h"
64 #include "MagickCore/exception.h"
66 #include "MagickCore/gem.h"
67 #include "MagickCore/gem-private.h"
68 #include "MagickCore/geometry.h"
69 #include "MagickCore/list.h"
71 #include "MagickCore/magic.h"
72 #include "MagickCore/magick.h"
73 #include "MagickCore/memory_.h"
74 #include "MagickCore/module.h"
75 #include "MagickCore/monitor.h"
77 #include "MagickCore/option.h"
78 #include "MagickCore/paint.h"
80 #include "MagickCore/profile.h"
81 #include "MagickCore/property.h"
82 #include "MagickCore/quantize.h"
84 #include "MagickCore/random_.h"
86 #include "MagickCore/resource_.h"
87 #include "MagickCore/segment.h"
88 #include "MagickCore/semaphore.h"
90 #include "MagickCore/statistic.h"
91 #include "MagickCore/string_.h"
93 #include "MagickCore/timer.h"
94 #include "MagickCore/utility.h"
95 #include "MagickCore/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 EvaluateImage 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 %
122 % A description of each parameter follows:
123 %
124 % o image: the image.
125 %
126 % o op: A channel op.
127 %
128 % o value: A value value.
129 %
130 % o exception: return any errors or warnings in this structure.
131 %
132 */
133 
134 typedef struct _PixelChannels
135 {
136  double
138 } PixelChannels;
139 
141 {
142  register ssize_t
143  i;
144 
145  assert(pixels != (PixelChannels **) NULL);
146  for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
147  if (pixels[i] != (PixelChannels *) NULL)
148  pixels[i]=(PixelChannels *) RelinquishMagickMemory(pixels[i]);
149  pixels=(PixelChannels **) RelinquishMagickMemory(pixels);
150  return(pixels);
151 }
152 
154 {
156  **pixels;
157 
158  register ssize_t
159  i;
160 
161  size_t
162  number_threads;
163 
164  number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
165  pixels=(PixelChannels **) AcquireQuantumMemory(number_threads,
166  sizeof(*pixels));
167  if (pixels == (PixelChannels **) NULL)
168  return((PixelChannels **) NULL);
169  (void) memset(pixels,0,number_threads*sizeof(*pixels));
170  for (i=0; i < (ssize_t) number_threads; i++)
171  {
172  register ssize_t
173  j;
174 
175  pixels[i]=(PixelChannels *) AcquireQuantumMemory(image->columns,
176  sizeof(**pixels));
177  if (pixels[i] == (PixelChannels *) NULL)
178  return(DestroyPixelThreadSet(pixels));
179  for (j=0; j < (ssize_t) image->columns; j++)
180  {
181  register ssize_t
182  k;
183 
184  for (k=0; k < MaxPixelChannels; k++)
185  pixels[i][j].channel[k]=0.0;
186  }
187  }
188  return(pixels);
189 }
190 
191 static inline double EvaluateMax(const double x,const double y)
192 {
193  if (x > y)
194  return(x);
195  return(y);
196 }
197 
198 #if defined(__cplusplus) || defined(c_plusplus)
199 extern "C" {
200 #endif
201 
202 static int IntensityCompare(const void *x,const void *y)
203 {
204  const PixelChannels
205  *color_1,
206  *color_2;
207 
208  double
209  distance;
210 
211  register ssize_t
212  i;
213 
214  color_1=(const PixelChannels *) x;
215  color_2=(const PixelChannels *) y;
216  distance=0.0;
217  for (i=0; i < MaxPixelChannels; i++)
218  distance+=color_1->channel[i]-(double) color_2->channel[i];
219  return(distance < 0 ? -1 : distance > 0 ? 1 : 0);
220 }
221 
222 #if defined(__cplusplus) || defined(c_plusplus)
223 }
224 #endif
225 
227  const MagickEvaluateOperator op,const double value)
228 {
229  double
230  result;
231 
232  result=0.0;
233  switch (op)
234  {
236  break;
237  case AbsEvaluateOperator:
238  {
239  result=(double) fabs((double) (pixel+value));
240  break;
241  }
242  case AddEvaluateOperator:
243  {
244  result=(double) (pixel+value);
245  break;
246  }
248  {
249  /*
250  This returns a 'floored modulus' of the addition which is a positive
251  result. It differs from % or fmod() that returns a 'truncated modulus'
252  result, where floor() is replaced by trunc() and could return a
253  negative result (which is clipped).
254  */
255  result=pixel+value;
256  result-=(QuantumRange+1.0)*floor((double) result/(QuantumRange+1.0));
257  break;
258  }
259  case AndEvaluateOperator:
260  {
261  result=(double) ((size_t) pixel & (size_t) (value+0.5));
262  break;
263  }
265  {
266  result=(double) (QuantumRange*(0.5*cos((double) (2.0*MagickPI*
267  QuantumScale*pixel*value))+0.5));
268  break;
269  }
271  {
272  result=pixel/(value == 0.0 ? 1.0 : value);
273  break;
274  }
276  {
277  result=(double) (QuantumRange*exp((double) (value*QuantumScale*pixel)));
278  break;
279  }
281  {
282  result=(double) GenerateDifferentialNoise(random_info,pixel,
283  GaussianNoise,value);
284  break;
285  }
287  {
288  result=(double) GenerateDifferentialNoise(random_info,pixel,ImpulseNoise,
289  value);
290  break;
291  }
293  {
294  result=(double) GenerateDifferentialNoise(random_info,pixel,
295  LaplacianNoise,value);
296  break;
297  }
299  {
300  result=(double) ((size_t) pixel << (size_t) (value+0.5));
301  break;
302  }
303  case LogEvaluateOperator:
304  {
305  if ((QuantumScale*pixel) >= MagickEpsilon)
306  result=(double) (QuantumRange*log((double) (QuantumScale*value*pixel+
307  1.0))/log((double) (value+1.0)));
308  break;
309  }
310  case MaxEvaluateOperator:
311  {
312  result=(double) EvaluateMax((double) pixel,value);
313  break;
314  }
316  {
317  result=(double) (pixel+value);
318  break;
319  }
321  {
322  result=(double) (pixel+value);
323  break;
324  }
325  case MinEvaluateOperator:
326  {
327  result=(double) MagickMin((double) pixel,value);
328  break;
329  }
331  {
332  result=(double) GenerateDifferentialNoise(random_info,pixel,
334  break;
335  }
337  {
338  result=(double) (value*pixel);
339  break;
340  }
341  case OrEvaluateOperator:
342  {
343  result=(double) ((size_t) pixel | (size_t) (value+0.5));
344  break;
345  }
347  {
348  result=(double) GenerateDifferentialNoise(random_info,pixel,PoissonNoise,
349  value);
350  break;
351  }
352  case PowEvaluateOperator:
353  {
354  result=(double) (QuantumRange*pow((double) (QuantumScale*pixel),(double)
355  value));
356  break;
357  }
359  {
360  result=(double) ((size_t) pixel >> (size_t) (value+0.5));
361  break;
362  }
364  {
365  result=(double) (pixel*pixel+value);
366  break;
367  }
368  case SetEvaluateOperator:
369  {
370  result=value;
371  break;
372  }
374  {
375  result=(double) (QuantumRange*(0.5*sin((double) (2.0*MagickPI*
376  QuantumScale*pixel*value))+0.5));
377  break;
378  }
380  {
381  result=(double) (pixel-value);
382  break;
383  }
384  case SumEvaluateOperator:
385  {
386  result=(double) (pixel+value);
387  break;
388  }
390  {
391  result=(double) (((double) pixel <= value) ? 0 : QuantumRange);
392  break;
393  }
395  {
396  result=(double) (((double) pixel <= value) ? 0 : pixel);
397  break;
398  }
400  {
401  result=(double) (((double) pixel > value) ? QuantumRange : pixel);
402  break;
403  }
405  {
406  result=(double) GenerateDifferentialNoise(random_info,pixel,UniformNoise,
407  value);
408  break;
409  }
410  case XorEvaluateOperator:
411  {
412  result=(double) ((size_t) pixel ^ (size_t) (value+0.5));
413  break;
414  }
415  }
416  return(result);
417 }
418 
419 static Image *AcquireImageCanvas(const Image *images,ExceptionInfo *exception)
420 {
421  const Image
422  *p,
423  *q;
424 
425  size_t
426  columns,
427  rows;
428 
429  q=images;
430  columns=images->columns;
431  rows=images->rows;
432  for (p=images; p != (Image *) NULL; p=p->next)
433  {
434  if (p->number_channels > q->number_channels)
435  q=p;
436  if (p->columns > columns)
437  columns=p->columns;
438  if (p->rows > rows)
439  rows=p->rows;
440  }
441  return(CloneImage(q,columns,rows,MagickTrue,exception));
442 }
443 
445  const MagickEvaluateOperator op,ExceptionInfo *exception)
446 {
447 #define EvaluateImageTag "Evaluate/Image"
448 
449  CacheView
450  *evaluate_view;
451 
452  Image
453  *image;
454 
456  status;
457 
459  progress;
460 
462  **magick_restrict evaluate_pixels;
463 
464  RandomInfo
466 
467  size_t
468  number_images;
469 
470  ssize_t
471  y;
472 
473 #if defined(MAGICKCORE_OPENMP_SUPPORT)
474  unsigned long
475  key;
476 #endif
477 
478  assert(images != (Image *) NULL);
479  assert(images->signature == MagickCoreSignature);
480  if (images->debug != MagickFalse)
481  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
482  assert(exception != (ExceptionInfo *) NULL);
483  assert(exception->signature == MagickCoreSignature);
484  image=AcquireImageCanvas(images,exception);
485  if (image == (Image *) NULL)
486  return((Image *) NULL);
487  if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
488  {
489  image=DestroyImage(image);
490  return((Image *) NULL);
491  }
492  number_images=GetImageListLength(images);
493  evaluate_pixels=AcquirePixelThreadSet(images);
494  if (evaluate_pixels == (PixelChannels **) NULL)
495  {
496  image=DestroyImage(image);
497  (void) ThrowMagickException(exception,GetMagickModule(),
498  ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
499  return((Image *) NULL);
500  }
501  /*
502  Evaluate image pixels.
503  */
504  status=MagickTrue;
505  progress=0;
507  evaluate_view=AcquireAuthenticCacheView(image,exception);
508  if (op == MedianEvaluateOperator)
509  {
510 #if defined(MAGICKCORE_OPENMP_SUPPORT)
512  #pragma omp parallel for schedule(static) shared(progress,status) \
513  magick_number_threads(image,images,image->rows,key == ~0UL)
514 #endif
515  for (y=0; y < (ssize_t) image->rows; y++)
516  {
517  CacheView
518  *image_view;
519 
520  const Image
521  *next;
522 
523  const int
524  id = GetOpenMPThreadId();
525 
526  register PixelChannels
527  *evaluate_pixel;
528 
529  register Quantum
530  *magick_restrict q;
531 
532  register ssize_t
533  x;
534 
535  if (status == MagickFalse)
536  continue;
537  q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,image->columns,1,
538  exception);
539  if (q == (Quantum *) NULL)
540  {
541  status=MagickFalse;
542  continue;
543  }
544  evaluate_pixel=evaluate_pixels[id];
545  for (x=0; x < (ssize_t) image->columns; x++)
546  {
547  register ssize_t
548  j,
549  k;
550 
551  for (j=0; j < (ssize_t) number_images; j++)
552  for (k=0; k < MaxPixelChannels; k++)
553  evaluate_pixel[j].channel[k]=0.0;
554  next=images;
555  for (j=0; j < (ssize_t) number_images; j++)
556  {
557  register const Quantum
558  *p;
559 
560  register ssize_t
561  i;
562 
563  image_view=AcquireVirtualCacheView(next,exception);
564  p=GetCacheViewVirtualPixels(image_view,x,y,1,1,exception);
565  if (p == (const Quantum *) NULL)
566  {
567  image_view=DestroyCacheView(image_view);
568  break;
569  }
570  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
571  {
572  PixelChannel channel = GetPixelChannelChannel(image,i);
573  PixelTrait evaluate_traits=GetPixelChannelTraits(image,channel);
574  PixelTrait traits = GetPixelChannelTraits(next,channel);
575  if ((traits == UndefinedPixelTrait) ||
576  (evaluate_traits == UndefinedPixelTrait))
577  continue;
578  if ((traits & UpdatePixelTrait) == 0)
579  continue;
580  evaluate_pixel[j].channel[i]=ApplyEvaluateOperator(
581  random_info[id],GetPixelChannel(image,channel,p),op,
582  evaluate_pixel[j].channel[i]);
583  }
584  image_view=DestroyCacheView(image_view);
585  next=GetNextImageInList(next);
586  }
587  qsort((void *) evaluate_pixel,number_images,sizeof(*evaluate_pixel),
589  for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
590  q[k]=ClampToQuantum(evaluate_pixel[j/2].channel[k]);
591  q+=GetPixelChannels(image);
592  }
593  if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
594  status=MagickFalse;
595  if (images->progress_monitor != (MagickProgressMonitor) NULL)
596  {
598  proceed;
599 
600 #if defined(MAGICKCORE_OPENMP_SUPPORT)
601  #pragma omp critical (MagickCore_EvaluateImages)
602 #endif
603  proceed=SetImageProgress(images,EvaluateImageTag,progress++,
604  image->rows);
605  if (proceed == MagickFalse)
606  status=MagickFalse;
607  }
608  }
609  }
610  else
611  {
612 #if defined(MAGICKCORE_OPENMP_SUPPORT)
614  #pragma omp parallel for schedule(static) shared(progress,status) \
615  magick_number_threads(image,images,image->rows,key == ~0UL)
616 #endif
617  for (y=0; y < (ssize_t) image->rows; y++)
618  {
619  CacheView
620  *image_view;
621 
622  const Image
623  *next;
624 
625  const int
626  id = GetOpenMPThreadId();
627 
628  register ssize_t
629  i,
630  x;
631 
632  register PixelChannels
633  *evaluate_pixel;
634 
635  register Quantum
636  *magick_restrict q;
637 
638  ssize_t
639  j;
640 
641  if (status == MagickFalse)
642  continue;
643  q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,image->columns,1,
644  exception);
645  if (q == (Quantum *) NULL)
646  {
647  status=MagickFalse;
648  continue;
649  }
650  evaluate_pixel=evaluate_pixels[id];
651  for (j=0; j < (ssize_t) image->columns; j++)
652  for (i=0; i < MaxPixelChannels; i++)
653  evaluate_pixel[j].channel[i]=0.0;
654  next=images;
655  for (j=0; j < (ssize_t) number_images; j++)
656  {
657  register const Quantum
658  *p;
659 
660  image_view=AcquireVirtualCacheView(next,exception);
661  p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,
662  exception);
663  if (p == (const Quantum *) NULL)
664  {
665  image_view=DestroyCacheView(image_view);
666  break;
667  }
668  for (x=0; x < (ssize_t) image->columns; x++)
669  {
670  register ssize_t
671  i;
672 
673  for (i=0; i < (ssize_t) GetPixelChannels(next); i++)
674  {
675  PixelChannel channel = GetPixelChannelChannel(image,i);
676  PixelTrait traits = GetPixelChannelTraits(next,channel);
677  PixelTrait evaluate_traits=GetPixelChannelTraits(image,channel);
678  if ((traits == UndefinedPixelTrait) ||
679  (evaluate_traits == UndefinedPixelTrait))
680  continue;
681  if ((traits & UpdatePixelTrait) == 0)
682  continue;
683  evaluate_pixel[x].channel[i]=ApplyEvaluateOperator(
684  random_info[id],GetPixelChannel(image,channel,p),j == 0 ?
685  AddEvaluateOperator : op,evaluate_pixel[x].channel[i]);
686  }
687  p+=GetPixelChannels(next);
688  }
689  image_view=DestroyCacheView(image_view);
690  next=GetNextImageInList(next);
691  }
692  for (x=0; x < (ssize_t) image->columns; x++)
693  {
694  register ssize_t
695  i;
696 
697  switch (op)
698  {
700  {
701  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
702  evaluate_pixel[x].channel[i]/=(double) number_images;
703  break;
704  }
706  {
707  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
708  {
709  register ssize_t
710  j;
711 
712  for (j=0; j < (ssize_t) (number_images-1); j++)
713  evaluate_pixel[x].channel[i]*=QuantumScale;
714  }
715  break;
716  }
718  {
719  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
720  evaluate_pixel[x].channel[i]=sqrt(evaluate_pixel[x].channel[i]/
721  number_images);
722  break;
723  }
724  default:
725  break;
726  }
727  }
728  for (x=0; x < (ssize_t) image->columns; x++)
729  {
730  register ssize_t
731  i;
732 
733  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
734  {
735  PixelChannel channel = GetPixelChannelChannel(image,i);
736  PixelTrait traits = GetPixelChannelTraits(image,channel);
737  if (traits == UndefinedPixelTrait)
738  continue;
739  if ((traits & UpdatePixelTrait) == 0)
740  continue;
741  q[i]=ClampToQuantum(evaluate_pixel[x].channel[i]);
742  }
743  q+=GetPixelChannels(image);
744  }
745  if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
746  status=MagickFalse;
747  if (images->progress_monitor != (MagickProgressMonitor) NULL)
748  {
750  proceed;
751 
752 #if defined(MAGICKCORE_OPENMP_SUPPORT)
753  #pragma omp critical (MagickCore_EvaluateImages)
754 #endif
755  proceed=SetImageProgress(images,EvaluateImageTag,progress++,
756  image->rows);
757  if (proceed == MagickFalse)
758  status=MagickFalse;
759  }
760  }
761  }
762  evaluate_view=DestroyCacheView(evaluate_view);
763  evaluate_pixels=DestroyPixelThreadSet(evaluate_pixels);
765  if (status == MagickFalse)
766  image=DestroyImage(image);
767  return(image);
768 }
769 
771  const MagickEvaluateOperator op,const double value,ExceptionInfo *exception)
772 {
773  CacheView
774  *image_view;
775 
777  status;
778 
780  progress;
781 
782  RandomInfo
784 
785  ssize_t
786  y;
787 
788 #if defined(MAGICKCORE_OPENMP_SUPPORT)
789  unsigned long
790  key;
791 #endif
792 
793  assert(image != (Image *) NULL);
794  assert(image->signature == MagickCoreSignature);
795  if (image->debug != MagickFalse)
796  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
797  assert(exception != (ExceptionInfo *) NULL);
798  assert(exception->signature == MagickCoreSignature);
799  if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
800  return(MagickFalse);
801  status=MagickTrue;
802  progress=0;
804  image_view=AcquireAuthenticCacheView(image,exception);
805 #if defined(MAGICKCORE_OPENMP_SUPPORT)
807  #pragma omp parallel for schedule(static) shared(progress,status) \
808  magick_number_threads(image,image,image->rows,key == ~0UL)
809 #endif
810  for (y=0; y < (ssize_t) image->rows; y++)
811  {
812  const int
813  id = GetOpenMPThreadId();
814 
815  register Quantum
816  *magick_restrict q;
817 
818  register ssize_t
819  x;
820 
821  if (status == MagickFalse)
822  continue;
823  q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
824  if (q == (Quantum *) NULL)
825  {
826  status=MagickFalse;
827  continue;
828  }
829  for (x=0; x < (ssize_t) image->columns; x++)
830  {
831  double
832  result;
833 
834  register ssize_t
835  i;
836 
837  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
838  {
839  PixelChannel channel = GetPixelChannelChannel(image,i);
840  PixelTrait traits = GetPixelChannelTraits(image,channel);
841  if (traits == UndefinedPixelTrait)
842  continue;
843  if ((traits & CopyPixelTrait) != 0)
844  continue;
845  if ((traits & UpdatePixelTrait) == 0)
846  continue;
847  result=ApplyEvaluateOperator(random_info[id],q[i],op,value);
848  if (op == MeanEvaluateOperator)
849  result/=2.0;
850  q[i]=ClampToQuantum(result);
851  }
852  q+=GetPixelChannels(image);
853  }
854  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
855  status=MagickFalse;
856  if (image->progress_monitor != (MagickProgressMonitor) NULL)
857  {
859  proceed;
860 
861 #if defined(MAGICKCORE_OPENMP_SUPPORT)
862  #pragma omp critical (MagickCore_EvaluateImage)
863 #endif
864  proceed=SetImageProgress(image,EvaluateImageTag,progress++,image->rows);
865  if (proceed == MagickFalse)
866  status=MagickFalse;
867  }
868  }
869  image_view=DestroyCacheView(image_view);
871  return(status);
872 }
873 
874 /*
875 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
876 % %
877 % %
878 % %
879 % F u n c t i o n I m a g e %
880 % %
881 % %
882 % %
883 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
884 %
885 % FunctionImage() applies a value to the image with an arithmetic, relational,
886 % or logical operator to an image. Use these operations to lighten or darken
887 % an image, to increase or decrease contrast in an image, or to produce the
888 % "negative" of an image.
889 %
890 % The format of the FunctionImage method is:
891 %
892 % MagickBooleanType FunctionImage(Image *image,
893 % const MagickFunction function,const ssize_t number_parameters,
894 % const double *parameters,ExceptionInfo *exception)
895 %
896 % A description of each parameter follows:
897 %
898 % o image: the image.
899 %
900 % o function: A channel function.
901 %
902 % o parameters: one or more parameters.
903 %
904 % o exception: return any errors or warnings in this structure.
905 %
906 */
907 
908 static Quantum ApplyFunction(Quantum pixel,const MagickFunction function,
909  const size_t number_parameters,const double *parameters,
910  ExceptionInfo *exception)
911 {
912  double
913  result;
914 
915  register ssize_t
916  i;
917 
918  (void) exception;
919  result=0.0;
920  switch (function)
921  {
922  case PolynomialFunction:
923  {
924  /*
925  Polynomial: polynomial constants, highest to lowest order (e.g. c0*x^3+
926  c1*x^2+c2*x+c3).
927  */
928  result=0.0;
929  for (i=0; i < (ssize_t) number_parameters; i++)
930  result=result*QuantumScale*pixel+parameters[i];
931  result*=QuantumRange;
932  break;
933  }
934  case SinusoidFunction:
935  {
936  double
937  amplitude,
938  bias,
939  frequency,
940  phase;
941 
942  /*
943  Sinusoid: frequency, phase, amplitude, bias.
944  */
945  frequency=(number_parameters >= 1) ? parameters[0] : 1.0;
946  phase=(number_parameters >= 2) ? parameters[1] : 0.0;
947  amplitude=(number_parameters >= 3) ? parameters[2] : 0.5;
948  bias=(number_parameters >= 4) ? parameters[3] : 0.5;
949  result=(double) (QuantumRange*(amplitude*sin((double) (2.0*
950  MagickPI*(frequency*QuantumScale*pixel+phase/360.0)))+bias));
951  break;
952  }
953  case ArcsinFunction:
954  {
955  double
956  bias,
957  center,
958  range,
959  width;
960 
961  /*
962  Arcsin (peged at range limits for invalid results): width, center,
963  range, and bias.
964  */
965  width=(number_parameters >= 1) ? parameters[0] : 1.0;
966  center=(number_parameters >= 2) ? parameters[1] : 0.5;
967  range=(number_parameters >= 3) ? parameters[2] : 1.0;
968  bias=(number_parameters >= 4) ? parameters[3] : 0.5;
969  result=2.0/width*(QuantumScale*pixel-center);
970  if ( result <= -1.0 )
971  result=bias-range/2.0;
972  else
973  if (result >= 1.0)
974  result=bias+range/2.0;
975  else
976  result=(double) (range/MagickPI*asin((double) result)+bias);
977  result*=QuantumRange;
978  break;
979  }
980  case ArctanFunction:
981  {
982  double
983  center,
984  bias,
985  range,
986  slope;
987 
988  /*
989  Arctan: slope, center, range, and bias.
990  */
991  slope=(number_parameters >= 1) ? parameters[0] : 1.0;
992  center=(number_parameters >= 2) ? parameters[1] : 0.5;
993  range=(number_parameters >= 3) ? parameters[2] : 1.0;
994  bias=(number_parameters >= 4) ? parameters[3] : 0.5;
995  result=(double) (MagickPI*slope*(QuantumScale*pixel-center));
996  result=(double) (QuantumRange*(range/MagickPI*atan((double)
997  result)+bias));
998  break;
999  }
1000  case UndefinedFunction:
1001  break;
1002  }
1003  return(ClampToQuantum(result));
1004 }
1005 
1007  const MagickFunction function,const size_t number_parameters,
1008  const double *parameters,ExceptionInfo *exception)
1009 {
1010 #define FunctionImageTag "Function/Image "
1011 
1012  CacheView
1013  *image_view;
1014 
1016  status;
1017 
1019  progress;
1020 
1021  ssize_t
1022  y;
1023 
1024  assert(image != (Image *) NULL);
1025  assert(image->signature == MagickCoreSignature);
1026  if (image->debug != MagickFalse)
1027  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1028  assert(exception != (ExceptionInfo *) NULL);
1029  assert(exception->signature == MagickCoreSignature);
1030 #if defined(MAGICKCORE_OPENCL_SUPPORT)
1031  if (AccelerateFunctionImage(image,function,number_parameters,parameters,
1032  exception) != MagickFalse)
1033  return(MagickTrue);
1034 #endif
1035  if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
1036  return(MagickFalse);
1037  status=MagickTrue;
1038  progress=0;
1039  image_view=AcquireAuthenticCacheView(image,exception);
1040 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1041  #pragma omp parallel for schedule(static) shared(progress,status) \
1042  magick_number_threads(image,image,image->rows,1)
1043 #endif
1044  for (y=0; y < (ssize_t) image->rows; y++)
1045  {
1046  register Quantum
1047  *magick_restrict q;
1048 
1049  register ssize_t
1050  x;
1051 
1052  if (status == MagickFalse)
1053  continue;
1054  q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
1055  if (q == (Quantum *) NULL)
1056  {
1057  status=MagickFalse;
1058  continue;
1059  }
1060  for (x=0; x < (ssize_t) image->columns; x++)
1061  {
1062  register ssize_t
1063  i;
1064 
1065  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1066  {
1067  PixelChannel channel = GetPixelChannelChannel(image,i);
1068  PixelTrait traits = GetPixelChannelTraits(image,channel);
1069  if (traits == UndefinedPixelTrait)
1070  continue;
1071  if ((traits & UpdatePixelTrait) == 0)
1072  continue;
1073  q[i]=ApplyFunction(q[i],function,number_parameters,parameters,
1074  exception);
1075  }
1076  q+=GetPixelChannels(image);
1077  }
1078  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1079  status=MagickFalse;
1080  if (image->progress_monitor != (MagickProgressMonitor) NULL)
1081  {
1083  proceed;
1084 
1085 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1086  #pragma omp critical (MagickCore_FunctionImage)
1087 #endif
1088  proceed=SetImageProgress(image,FunctionImageTag,progress++,image->rows);
1089  if (proceed == MagickFalse)
1090  status=MagickFalse;
1091  }
1092  }
1093  image_view=DestroyCacheView(image_view);
1094  return(status);
1095 }
1096 
1097 /*
1098 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1099 % %
1100 % %
1101 % %
1102 % G e t I m a g e E n t r o p y %
1103 % %
1104 % %
1105 % %
1106 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1107 %
1108 % GetImageEntropy() returns the entropy of one or more image channels.
1109 %
1110 % The format of the GetImageEntropy method is:
1111 %
1112 % MagickBooleanType GetImageEntropy(const Image *image,double *entropy,
1113 % ExceptionInfo *exception)
1114 %
1115 % A description of each parameter follows:
1116 %
1117 % o image: the image.
1118 %
1119 % o entropy: the average entropy of the selected channels.
1120 %
1121 % o exception: return any errors or warnings in this structure.
1122 %
1123 */
1125  double *entropy,ExceptionInfo *exception)
1126 {
1128  *channel_statistics;
1129 
1130  assert(image != (Image *) NULL);
1131  assert(image->signature == MagickCoreSignature);
1132  if (image->debug != MagickFalse)
1133  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1134  channel_statistics=GetImageStatistics(image,exception);
1135  if (channel_statistics == (ChannelStatistics *) NULL)
1136  return(MagickFalse);
1137  *entropy=channel_statistics[CompositePixelChannel].entropy;
1138  channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1139  channel_statistics);
1140  return(MagickTrue);
1141 }
1142 
1143 /*
1144 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1145 % %
1146 % %
1147 % %
1148 % G e t I m a g e E x t r e m a %
1149 % %
1150 % %
1151 % %
1152 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1153 %
1154 % GetImageExtrema() returns the extrema of one or more image channels.
1155 %
1156 % The format of the GetImageExtrema method is:
1157 %
1158 % MagickBooleanType GetImageExtrema(const Image *image,size_t *minima,
1159 % size_t *maxima,ExceptionInfo *exception)
1160 %
1161 % A description of each parameter follows:
1162 %
1163 % o image: the image.
1164 %
1165 % o minima: the minimum value in the channel.
1166 %
1167 % o maxima: the maximum value in the channel.
1168 %
1169 % o exception: return any errors or warnings in this structure.
1170 %
1171 */
1173  size_t *minima,size_t *maxima,ExceptionInfo *exception)
1174 {
1175  double
1176  max,
1177  min;
1178 
1180  status;
1181 
1182  assert(image != (Image *) NULL);
1183  assert(image->signature == MagickCoreSignature);
1184  if (image->debug != MagickFalse)
1185  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1186  status=GetImageRange(image,&min,&max,exception);
1187  *minima=(size_t) ceil(min-0.5);
1188  *maxima=(size_t) floor(max+0.5);
1189  return(status);
1190 }
1191 
1192 /*
1193 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1194 % %
1195 % %
1196 % %
1197 % G e t I m a g e K u r t o s i s %
1198 % %
1199 % %
1200 % %
1201 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1202 %
1203 % GetImageKurtosis() returns the kurtosis and skewness of one or more image
1204 % channels.
1205 %
1206 % The format of the GetImageKurtosis method is:
1207 %
1208 % MagickBooleanType GetImageKurtosis(const Image *image,double *kurtosis,
1209 % double *skewness,ExceptionInfo *exception)
1210 %
1211 % A description of each parameter follows:
1212 %
1213 % o image: the image.
1214 %
1215 % o kurtosis: the kurtosis of the channel.
1216 %
1217 % o skewness: the skewness of the channel.
1218 %
1219 % o exception: return any errors or warnings in this structure.
1220 %
1221 */
1223  double *kurtosis,double *skewness,ExceptionInfo *exception)
1224 {
1226  *channel_statistics;
1227 
1228  assert(image != (Image *) NULL);
1229  assert(image->signature == MagickCoreSignature);
1230  if (image->debug != MagickFalse)
1231  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1232  channel_statistics=GetImageStatistics(image,exception);
1233  if (channel_statistics == (ChannelStatistics *) NULL)
1234  return(MagickFalse);
1235  *kurtosis=channel_statistics[CompositePixelChannel].kurtosis;
1236  *skewness=channel_statistics[CompositePixelChannel].skewness;
1237  channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1238  channel_statistics);
1239  return(MagickTrue);
1240 }
1241 
1242 /*
1243 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1244 % %
1245 % %
1246 % %
1247 % G e t I m a g e M e a n %
1248 % %
1249 % %
1250 % %
1251 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1252 %
1253 % GetImageMean() returns the mean and standard deviation of one or more image
1254 % channels.
1255 %
1256 % The format of the GetImageMean method is:
1257 %
1258 % MagickBooleanType GetImageMean(const Image *image,double *mean,
1259 % double *standard_deviation,ExceptionInfo *exception)
1260 %
1261 % A description of each parameter follows:
1262 %
1263 % o image: the image.
1264 %
1265 % o mean: the average value in the channel.
1266 %
1267 % o standard_deviation: the standard deviation of the channel.
1268 %
1269 % o exception: return any errors or warnings in this structure.
1270 %
1271 */
1273  double *standard_deviation,ExceptionInfo *exception)
1274 {
1276  *channel_statistics;
1277 
1278  assert(image != (Image *) NULL);
1279  assert(image->signature == MagickCoreSignature);
1280  if (image->debug != MagickFalse)
1281  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1282  channel_statistics=GetImageStatistics(image,exception);
1283  if (channel_statistics == (ChannelStatistics *) NULL)
1284  return(MagickFalse);
1285  *mean=channel_statistics[CompositePixelChannel].mean;
1286  *standard_deviation=
1287  channel_statistics[CompositePixelChannel].standard_deviation;
1288  channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1289  channel_statistics);
1290  return(MagickTrue);
1291 }
1292 
1293 /*
1294 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1295 % %
1296 % %
1297 % %
1298 % G e t I m a g e M o m e n t s %
1299 % %
1300 % %
1301 % %
1302 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1303 %
1304 % GetImageMoments() returns the normalized moments of one or more image
1305 % channels.
1306 %
1307 % The format of the GetImageMoments method is:
1308 %
1309 % ChannelMoments *GetImageMoments(const Image *image,
1310 % ExceptionInfo *exception)
1311 %
1312 % A description of each parameter follows:
1313 %
1314 % o image: the image.
1315 %
1316 % o exception: return any errors or warnings in this structure.
1317 %
1318 */
1319 
1320 static size_t GetImageChannels(const Image *image)
1321 {
1322  register ssize_t
1323  i;
1324 
1325  size_t
1326  channels;
1327 
1328  channels=0;
1329  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1330  {
1331  PixelChannel channel = GetPixelChannelChannel(image,i);
1332  PixelTrait traits = GetPixelChannelTraits(image,channel);
1333  if (traits == UndefinedPixelTrait)
1334  continue;
1335  if ((traits & UpdatePixelTrait) == 0)
1336  continue;
1337  channels++;
1338  }
1339  return((size_t) (channels == 0 ? 1 : channels));
1340 }
1341 
1343  ExceptionInfo *exception)
1344 {
1345 #define MaxNumberImageMoments 8
1346 
1347  CacheView
1348  *image_view;
1349 
1351  *channel_moments;
1352 
1353  double
1354  M00[MaxPixelChannels+1],
1355  M01[MaxPixelChannels+1],
1356  M02[MaxPixelChannels+1],
1357  M03[MaxPixelChannels+1],
1358  M10[MaxPixelChannels+1],
1359  M11[MaxPixelChannels+1],
1360  M12[MaxPixelChannels+1],
1361  M20[MaxPixelChannels+1],
1362  M21[MaxPixelChannels+1],
1363  M22[MaxPixelChannels+1],
1364  M30[MaxPixelChannels+1];
1365 
1366  PointInfo
1367  centroid[MaxPixelChannels+1];
1368 
1369  ssize_t
1370  channel,
1371  y;
1372 
1373  assert(image != (Image *) NULL);
1374  assert(image->signature == MagickCoreSignature);
1375  if (image->debug != MagickFalse)
1376  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1378  sizeof(*channel_moments));
1379  if (channel_moments == (ChannelMoments *) NULL)
1380  return(channel_moments);
1381  (void) memset(channel_moments,0,(MaxPixelChannels+1)*
1382  sizeof(*channel_moments));
1383  (void) memset(centroid,0,sizeof(centroid));
1384  (void) memset(M00,0,sizeof(M00));
1385  (void) memset(M01,0,sizeof(M01));
1386  (void) memset(M02,0,sizeof(M02));
1387  (void) memset(M03,0,sizeof(M03));
1388  (void) memset(M10,0,sizeof(M10));
1389  (void) memset(M11,0,sizeof(M11));
1390  (void) memset(M12,0,sizeof(M12));
1391  (void) memset(M20,0,sizeof(M20));
1392  (void) memset(M21,0,sizeof(M21));
1393  (void) memset(M22,0,sizeof(M22));
1394  (void) memset(M30,0,sizeof(M30));
1395  image_view=AcquireVirtualCacheView(image,exception);
1396  for (y=0; y < (ssize_t) image->rows; y++)
1397  {
1398  register const Quantum
1399  *magick_restrict p;
1400 
1401  register ssize_t
1402  x;
1403 
1404  /*
1405  Compute center of mass (centroid).
1406  */
1407  p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1408  if (p == (const Quantum *) NULL)
1409  break;
1410  for (x=0; x < (ssize_t) image->columns; x++)
1411  {
1412  register ssize_t
1413  i;
1414 
1415  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1416  {
1417  PixelChannel channel = GetPixelChannelChannel(image,i);
1418  PixelTrait traits = GetPixelChannelTraits(image,channel);
1419  if (traits == UndefinedPixelTrait)
1420  continue;
1421  if ((traits & UpdatePixelTrait) == 0)
1422  continue;
1423  M00[channel]+=QuantumScale*p[i];
1424  M00[MaxPixelChannels]+=QuantumScale*p[i];
1425  M10[channel]+=x*QuantumScale*p[i];
1426  M10[MaxPixelChannels]+=x*QuantumScale*p[i];
1427  M01[channel]+=y*QuantumScale*p[i];
1428  M01[MaxPixelChannels]+=y*QuantumScale*p[i];
1429  }
1430  p+=GetPixelChannels(image);
1431  }
1432  }
1433  for (channel=0; channel <= MaxPixelChannels; channel++)
1434  {
1435  /*
1436  Compute center of mass (centroid).
1437  */
1438  if (M00[channel] < MagickEpsilon)
1439  {
1440  M00[channel]+=MagickEpsilon;
1441  centroid[channel].x=(double) image->columns/2.0;
1442  centroid[channel].y=(double) image->rows/2.0;
1443  continue;
1444  }
1445  M00[channel]+=MagickEpsilon;
1446  centroid[channel].x=M10[channel]/M00[channel];
1447  centroid[channel].y=M01[channel]/M00[channel];
1448  }
1449  for (y=0; y < (ssize_t) image->rows; y++)
1450  {
1451  register const Quantum
1452  *magick_restrict p;
1453 
1454  register ssize_t
1455  x;
1456 
1457  /*
1458  Compute the image moments.
1459  */
1460  p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1461  if (p == (const Quantum *) NULL)
1462  break;
1463  for (x=0; x < (ssize_t) image->columns; x++)
1464  {
1465  register ssize_t
1466  i;
1467 
1468  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1469  {
1470  PixelChannel channel = GetPixelChannelChannel(image,i);
1471  PixelTrait traits = GetPixelChannelTraits(image,channel);
1472  if (traits == UndefinedPixelTrait)
1473  continue;
1474  if ((traits & UpdatePixelTrait) == 0)
1475  continue;
1476  M11[channel]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
1477  QuantumScale*p[i];
1478  M11[MaxPixelChannels]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
1479  QuantumScale*p[i];
1480  M20[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1481  QuantumScale*p[i];
1482  M20[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1483  QuantumScale*p[i];
1484  M02[channel]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
1485  QuantumScale*p[i];
1486  M02[MaxPixelChannels]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
1487  QuantumScale*p[i];
1488  M21[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1489  (y-centroid[channel].y)*QuantumScale*p[i];
1490  M21[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1491  (y-centroid[channel].y)*QuantumScale*p[i];
1492  M12[channel]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
1493  (y-centroid[channel].y)*QuantumScale*p[i];
1494  M12[MaxPixelChannels]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
1495  (y-centroid[channel].y)*QuantumScale*p[i];
1496  M22[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1497  (y-centroid[channel].y)*(y-centroid[channel].y)*QuantumScale*p[i];
1498  M22[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1499  (y-centroid[channel].y)*(y-centroid[channel].y)*QuantumScale*p[i];
1500  M30[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1501  (x-centroid[channel].x)*QuantumScale*p[i];
1502  M30[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1503  (x-centroid[channel].x)*QuantumScale*p[i];
1504  M03[channel]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
1505  (y-centroid[channel].y)*QuantumScale*p[i];
1506  M03[MaxPixelChannels]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
1507  (y-centroid[channel].y)*QuantumScale*p[i];
1508  }
1509  p+=GetPixelChannels(image);
1510  }
1511  }
1512  M00[MaxPixelChannels]/=GetImageChannels(image);
1513  M01[MaxPixelChannels]/=GetImageChannels(image);
1514  M02[MaxPixelChannels]/=GetImageChannels(image);
1515  M03[MaxPixelChannels]/=GetImageChannels(image);
1516  M10[MaxPixelChannels]/=GetImageChannels(image);
1517  M11[MaxPixelChannels]/=GetImageChannels(image);
1518  M12[MaxPixelChannels]/=GetImageChannels(image);
1519  M20[MaxPixelChannels]/=GetImageChannels(image);
1520  M21[MaxPixelChannels]/=GetImageChannels(image);
1521  M22[MaxPixelChannels]/=GetImageChannels(image);
1522  M30[MaxPixelChannels]/=GetImageChannels(image);
1523  for (channel=0; channel <= MaxPixelChannels; channel++)
1524  {
1525  /*
1526  Compute elliptical angle, major and minor axes, eccentricity, & intensity.
1527  */
1528  channel_moments[channel].centroid=centroid[channel];
1529  channel_moments[channel].ellipse_axis.x=sqrt((2.0/M00[channel])*
1530  ((M20[channel]+M02[channel])+sqrt(4.0*M11[channel]*M11[channel]+
1531  (M20[channel]-M02[channel])*(M20[channel]-M02[channel]))));
1532  channel_moments[channel].ellipse_axis.y=sqrt((2.0/M00[channel])*
1533  ((M20[channel]+M02[channel])-sqrt(4.0*M11[channel]*M11[channel]+
1534  (M20[channel]-M02[channel])*(M20[channel]-M02[channel]))));
1535  channel_moments[channel].ellipse_angle=RadiansToDegrees(0.5*atan(2.0*
1536  M11[channel]/(M20[channel]-M02[channel]+MagickEpsilon)));
1537  if (fabs(M11[channel]) < MagickEpsilon)
1538  {
1539  if (fabs(M20[channel]-M02[channel]) < MagickEpsilon)
1540  channel_moments[channel].ellipse_angle+=0.0;
1541  else
1542  if ((M20[channel]-M02[channel]) < 0.0)
1543  channel_moments[channel].ellipse_angle+=90.0;
1544  else
1545  channel_moments[channel].ellipse_angle+=0.0;
1546  }
1547  else
1548  if (M11[channel] < 0.0)
1549  {
1550  if (fabs(M20[channel]-M02[channel]) < MagickEpsilon)
1551  channel_moments[channel].ellipse_angle+=0.0;
1552  else
1553  if ((M20[channel]-M02[channel]) < 0.0)
1554  channel_moments[channel].ellipse_angle+=90.0;
1555  else
1556  channel_moments[channel].ellipse_angle+=180.0;
1557  }
1558  else
1559  {
1560  if (fabs(M20[channel]-M02[channel]) < MagickEpsilon)
1561  channel_moments[channel].ellipse_angle+=0.0;
1562  else
1563  if ((M20[channel]-M02[channel]) < 0.0)
1564  channel_moments[channel].ellipse_angle+=90.0;
1565  else
1566  channel_moments[channel].ellipse_angle+=0.0;
1567  }
1568  channel_moments[channel].ellipse_eccentricity=sqrt(1.0-(
1569  channel_moments[channel].ellipse_axis.y/
1570  (channel_moments[channel].ellipse_axis.x+MagickEpsilon)));
1571  channel_moments[channel].ellipse_intensity=M00[channel]/
1572  (MagickPI*channel_moments[channel].ellipse_axis.x*
1573  channel_moments[channel].ellipse_axis.y+MagickEpsilon);
1574  }
1575  for (channel=0; channel <= MaxPixelChannels; channel++)
1576  {
1577  /*
1578  Normalize image moments.
1579  */
1580  M10[channel]=0.0;
1581  M01[channel]=0.0;
1582  M11[channel]/=pow(M00[channel],1.0+(1.0+1.0)/2.0);
1583  M20[channel]/=pow(M00[channel],1.0+(2.0+0.0)/2.0);
1584  M02[channel]/=pow(M00[channel],1.0+(0.0+2.0)/2.0);
1585  M21[channel]/=pow(M00[channel],1.0+(2.0+1.0)/2.0);
1586  M12[channel]/=pow(M00[channel],1.0+(1.0+2.0)/2.0);
1587  M22[channel]/=pow(M00[channel],1.0+(2.0+2.0)/2.0);
1588  M30[channel]/=pow(M00[channel],1.0+(3.0+0.0)/2.0);
1589  M03[channel]/=pow(M00[channel],1.0+(0.0+3.0)/2.0);
1590  M00[channel]=1.0;
1591  }
1592  image_view=DestroyCacheView(image_view);
1593  for (channel=0; channel <= MaxPixelChannels; channel++)
1594  {
1595  /*
1596  Compute Hu invariant moments.
1597  */
1598  channel_moments[channel].invariant[0]=M20[channel]+M02[channel];
1599  channel_moments[channel].invariant[1]=(M20[channel]-M02[channel])*
1600  (M20[channel]-M02[channel])+4.0*M11[channel]*M11[channel];
1601  channel_moments[channel].invariant[2]=(M30[channel]-3.0*M12[channel])*
1602  (M30[channel]-3.0*M12[channel])+(3.0*M21[channel]-M03[channel])*
1603  (3.0*M21[channel]-M03[channel]);
1604  channel_moments[channel].invariant[3]=(M30[channel]+M12[channel])*
1605  (M30[channel]+M12[channel])+(M21[channel]+M03[channel])*
1606  (M21[channel]+M03[channel]);
1607  channel_moments[channel].invariant[4]=(M30[channel]-3.0*M12[channel])*
1608  (M30[channel]+M12[channel])*((M30[channel]+M12[channel])*
1609  (M30[channel]+M12[channel])-3.0*(M21[channel]+M03[channel])*
1610  (M21[channel]+M03[channel]))+(3.0*M21[channel]-M03[channel])*
1611  (M21[channel]+M03[channel])*(3.0*(M30[channel]+M12[channel])*
1612  (M30[channel]+M12[channel])-(M21[channel]+M03[channel])*
1613  (M21[channel]+M03[channel]));
1614  channel_moments[channel].invariant[5]=(M20[channel]-M02[channel])*
1615  ((M30[channel]+M12[channel])*(M30[channel]+M12[channel])-
1616  (M21[channel]+M03[channel])*(M21[channel]+M03[channel]))+
1617  4.0*M11[channel]*(M30[channel]+M12[channel])*(M21[channel]+M03[channel]);
1618  channel_moments[channel].invariant[6]=(3.0*M21[channel]-M03[channel])*
1619  (M30[channel]+M12[channel])*((M30[channel]+M12[channel])*
1620  (M30[channel]+M12[channel])-3.0*(M21[channel]+M03[channel])*
1621  (M21[channel]+M03[channel]))-(M30[channel]-3*M12[channel])*
1622  (M21[channel]+M03[channel])*(3.0*(M30[channel]+M12[channel])*
1623  (M30[channel]+M12[channel])-(M21[channel]+M03[channel])*
1624  (M21[channel]+M03[channel]));
1625  channel_moments[channel].invariant[7]=M11[channel]*((M30[channel]+
1626  M12[channel])*(M30[channel]+M12[channel])-(M03[channel]+M21[channel])*
1627  (M03[channel]+M21[channel]))-(M20[channel]-M02[channel])*
1628  (M30[channel]+M12[channel])*(M03[channel]+M21[channel]);
1629  }
1630  if (y < (ssize_t) image->rows)
1631  channel_moments=(ChannelMoments *) RelinquishMagickMemory(channel_moments);
1632  return(channel_moments);
1633 }
1634 
1635 /*
1636 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1637 % %
1638 % %
1639 % %
1640 % 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 %
1641 % %
1642 % %
1643 % %
1644 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1645 %
1646 % GetImagePerceptualHash() returns the perceptual hash of one or more
1647 % image channels.
1648 %
1649 % The format of the GetImagePerceptualHash method is:
1650 %
1651 % ChannelPerceptualHash *GetImagePerceptualHash(const Image *image,
1652 % ExceptionInfo *exception)
1653 %
1654 % A description of each parameter follows:
1655 %
1656 % o image: the image.
1657 %
1658 % o exception: return any errors or warnings in this structure.
1659 %
1660 */
1661 
1662 static inline double MagickLog10(const double x)
1663 {
1664 #define Log10Epsilon (1.0e-11)
1665 
1666  if (fabs(x) < Log10Epsilon)
1667  return(log10(Log10Epsilon));
1668  return(log10(fabs(x)));
1669 }
1670 
1672  ExceptionInfo *exception)
1673 {
1675  *perceptual_hash;
1676 
1677  char
1678  *colorspaces,
1679  *q;
1680 
1681  const char
1682  *artifact;
1683 
1685  status;
1686 
1687  register char
1688  *p;
1689 
1690  register ssize_t
1691  i;
1692 
1693  perceptual_hash=(ChannelPerceptualHash *) AcquireQuantumMemory(
1694  MaxPixelChannels+1UL,sizeof(*perceptual_hash));
1695  if (perceptual_hash == (ChannelPerceptualHash *) NULL)
1696  return((ChannelPerceptualHash *) NULL);
1697  artifact=GetImageArtifact(image,"phash:colorspaces");
1698  if (artifact != NULL)
1699  colorspaces=AcquireString(artifact);
1700  else
1701  colorspaces=AcquireString("sRGB,HCLp");
1702  perceptual_hash[0].number_colorspaces=0;
1703  perceptual_hash[0].number_channels=0;
1704  q=colorspaces;
1705  for (i=0; (p=StringToken(",",&q)) != (char *) NULL; i++)
1706  {
1708  *moments;
1709 
1710  Image
1711  *hash_image;
1712 
1713  size_t
1714  j;
1715 
1716  ssize_t
1717  channel,
1718  colorspace;
1719 
1721  break;
1723  if (colorspace < 0)
1724  break;
1725  perceptual_hash[0].colorspace[i]=(ColorspaceType) colorspace;
1726  hash_image=BlurImage(image,0.0,1.0,exception);
1727  if (hash_image == (Image *) NULL)
1728  break;
1729  hash_image->depth=8;
1730  status=TransformImageColorspace(hash_image,(ColorspaceType) colorspace,
1731  exception);
1732  if (status == MagickFalse)
1733  break;
1734  moments=GetImageMoments(hash_image,exception);
1735  perceptual_hash[0].number_colorspaces++;
1736  perceptual_hash[0].number_channels+=GetImageChannels(hash_image);
1737  hash_image=DestroyImage(hash_image);
1738  if (moments == (ChannelMoments *) NULL)
1739  break;
1740  for (channel=0; channel <= MaxPixelChannels; channel++)
1741  for (j=0; j < MaximumNumberOfImageMoments; j++)
1742  perceptual_hash[channel].phash[i][j]=
1743  (-MagickLog10(moments[channel].invariant[j]));
1744  moments=(ChannelMoments *) RelinquishMagickMemory(moments);
1745  }
1746  colorspaces=DestroyString(colorspaces);
1747  return(perceptual_hash);
1748 }
1749 
1750 /*
1751 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1752 % %
1753 % %
1754 % %
1755 % G e t I m a g e R a n g e %
1756 % %
1757 % %
1758 % %
1759 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1760 %
1761 % GetImageRange() returns the range of one or more image channels.
1762 %
1763 % The format of the GetImageRange method is:
1764 %
1765 % MagickBooleanType GetImageRange(const Image *image,double *minima,
1766 % double *maxima,ExceptionInfo *exception)
1767 %
1768 % A description of each parameter follows:
1769 %
1770 % o image: the image.
1771 %
1772 % o minima: the minimum value in the channel.
1773 %
1774 % o maxima: the maximum value in the channel.
1775 %
1776 % o exception: return any errors or warnings in this structure.
1777 %
1778 */
1780  double *maxima,ExceptionInfo *exception)
1781 {
1782  CacheView
1783  *image_view;
1784 
1786  initialize,
1787  status;
1788 
1789  ssize_t
1790  y;
1791 
1792  assert(image != (Image *) NULL);
1793  assert(image->signature == MagickCoreSignature);
1794  if (image->debug != MagickFalse)
1795  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1796  status=MagickTrue;
1797  initialize=MagickTrue;
1798  *maxima=0.0;
1799  *minima=0.0;
1800  image_view=AcquireVirtualCacheView(image,exception);
1801 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1802  #pragma omp parallel for schedule(static) shared(status,initialize) \
1803  magick_number_threads(image,image,image->rows,1)
1804 #endif
1805  for (y=0; y < (ssize_t) image->rows; y++)
1806  {
1807  double
1808  row_maxima = 0.0,
1809  row_minima = 0.0;
1810 
1812  row_initialize;
1813 
1814  register const Quantum
1815  *magick_restrict p;
1816 
1817  register ssize_t
1818  x;
1819 
1820  if (status == MagickFalse)
1821  continue;
1822  p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1823  if (p == (const Quantum *) NULL)
1824  {
1825  status=MagickFalse;
1826  continue;
1827  }
1828  row_initialize=MagickTrue;
1829  for (x=0; x < (ssize_t) image->columns; x++)
1830  {
1831  register ssize_t
1832  i;
1833 
1834  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1835  {
1836  PixelChannel channel = GetPixelChannelChannel(image,i);
1837  PixelTrait traits = GetPixelChannelTraits(image,channel);
1838  if (traits == UndefinedPixelTrait)
1839  continue;
1840  if ((traits & UpdatePixelTrait) == 0)
1841  continue;
1842  if (row_initialize != MagickFalse)
1843  {
1844  row_minima=(double) p[i];
1845  row_maxima=(double) p[i];
1846  row_initialize=MagickFalse;
1847  }
1848  else
1849  {
1850  if ((double) p[i] < row_minima)
1851  row_minima=(double) p[i];
1852  if ((double) p[i] > row_maxima)
1853  row_maxima=(double) p[i];
1854  }
1855  }
1856  p+=GetPixelChannels(image);
1857  }
1858 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1859 #pragma omp critical (MagickCore_GetImageRange)
1860 #endif
1861  {
1862  if (initialize != MagickFalse)
1863  {
1864  *minima=row_minima;
1865  *maxima=row_maxima;
1866  initialize=MagickFalse;
1867  }
1868  else
1869  {
1870  if (row_minima < *minima)
1871  *minima=row_minima;
1872  if (row_maxima > *maxima)
1873  *maxima=row_maxima;
1874  }
1875  }
1876  }
1877  image_view=DestroyCacheView(image_view);
1878  return(status);
1879 }
1880 
1881 /*
1882 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1883 % %
1884 % %
1885 % %
1886 % G e t I m a g e S t a t i s t i c s %
1887 % %
1888 % %
1889 % %
1890 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1891 %
1892 % GetImageStatistics() returns statistics for each channel in the image. The
1893 % statistics include the channel depth, its minima, maxima, mean, standard
1894 % deviation, kurtosis and skewness. You can access the red channel mean, for
1895 % example, like this:
1896 %
1897 % channel_statistics=GetImageStatistics(image,exception);
1898 % red_mean=channel_statistics[RedPixelChannel].mean;
1899 %
1900 % Use MagickRelinquishMemory() to free the statistics buffer.
1901 %
1902 % The format of the GetImageStatistics method is:
1903 %
1904 % ChannelStatistics *GetImageStatistics(const Image *image,
1905 % ExceptionInfo *exception)
1906 %
1907 % A description of each parameter follows:
1908 %
1909 % o image: the image.
1910 %
1911 % o exception: return any errors or warnings in this structure.
1912 %
1913 */
1915  ExceptionInfo *exception)
1916 {
1918  *channel_statistics;
1919 
1920  double
1921  area,
1922  *histogram,
1923  standard_deviation;
1924 
1926  status;
1927 
1928  QuantumAny
1929  range;
1930 
1931  register ssize_t
1932  i;
1933 
1934  size_t
1935  depth;
1936 
1937  ssize_t
1938  y;
1939 
1940  assert(image != (Image *) NULL);
1941  assert(image->signature == MagickCoreSignature);
1942  if (image->debug != MagickFalse)
1943  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1944  histogram=(double *) AcquireQuantumMemory(MaxMap+1UL,GetPixelChannels(image)*
1945  sizeof(*histogram));
1946  channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(
1947  MaxPixelChannels+1,sizeof(*channel_statistics));
1948  if ((channel_statistics == (ChannelStatistics *) NULL) ||
1949  (histogram == (double *) NULL))
1950  {
1951  if (histogram != (double *) NULL)
1952  histogram=(double *) RelinquishMagickMemory(histogram);
1953  if (channel_statistics != (ChannelStatistics *) NULL)
1954  channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1955  channel_statistics);
1956  return(channel_statistics);
1957  }
1958  (void) memset(channel_statistics,0,(MaxPixelChannels+1)*
1959  sizeof(*channel_statistics));
1960  for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
1961  {
1962  channel_statistics[i].depth=1;
1963  channel_statistics[i].maxima=(-MagickMaximumValue);
1964  channel_statistics[i].minima=MagickMaximumValue;
1965  }
1966  (void) memset(histogram,0,(MaxMap+1)*GetPixelChannels(image)*
1967  sizeof(*histogram));
1968  for (y=0; y < (ssize_t) image->rows; y++)
1969  {
1970  register const Quantum
1971  *magick_restrict p;
1972 
1973  register ssize_t
1974  x;
1975 
1976  /*
1977  Compute pixel statistics.
1978  */
1979  p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1980  if (p == (const Quantum *) NULL)
1981  break;
1982  for (x=0; x < (ssize_t) image->columns; x++)
1983  {
1984  register ssize_t
1985  i;
1986 
1987  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1988  {
1989  PixelChannel channel = GetPixelChannelChannel(image,i);
1990  PixelTrait traits = GetPixelChannelTraits(image,channel);
1991  if (traits == UndefinedPixelTrait)
1992  continue;
1993  if ((traits & UpdatePixelTrait) == 0)
1994  continue;
1995  if (channel_statistics[channel].depth != MAGICKCORE_QUANTUM_DEPTH)
1996  {
1997  depth=channel_statistics[channel].depth;
1998  range=GetQuantumRange(depth);
1999  status=p[i] != ScaleAnyToQuantum(ScaleQuantumToAny(p[i],range),
2000  range) ? MagickTrue : MagickFalse;
2001  if (status != MagickFalse)
2002  {
2003  channel_statistics[channel].depth++;
2004  i--;
2005  continue;
2006  }
2007  }
2008  if ((double) p[i] < channel_statistics[channel].minima)
2009  channel_statistics[channel].minima=(double) p[i];
2010  if ((double) p[i] > channel_statistics[channel].maxima)
2011  channel_statistics[channel].maxima=(double) p[i];
2012  channel_statistics[channel].sum+=p[i];
2013  channel_statistics[channel].sum_squared+=(double) p[i]*p[i];
2014  channel_statistics[channel].sum_cubed+=(double) p[i]*p[i]*p[i];
2015  channel_statistics[channel].sum_fourth_power+=(double) p[i]*p[i]*p[i]*
2016  p[i];
2017  channel_statistics[channel].area++;
2018  if ((double) p[i] < channel_statistics[CompositePixelChannel].minima)
2019  channel_statistics[CompositePixelChannel].minima=(double) p[i];
2020  if ((double) p[i] > channel_statistics[CompositePixelChannel].maxima)
2021  channel_statistics[CompositePixelChannel].maxima=(double) p[i];
2022  histogram[GetPixelChannels(image)*ScaleQuantumToMap(
2023  ClampToQuantum((double) p[i]))+i]++;
2024  channel_statistics[CompositePixelChannel].sum+=(double) p[i];
2025  channel_statistics[CompositePixelChannel].sum_squared+=(double)
2026  p[i]*p[i];
2027  channel_statistics[CompositePixelChannel].sum_cubed+=(double)
2028  p[i]*p[i]*p[i];
2029  channel_statistics[CompositePixelChannel].sum_fourth_power+=(double)
2030  p[i]*p[i]*p[i]*p[i];
2031  channel_statistics[CompositePixelChannel].area++;
2032  }
2033  p+=GetPixelChannels(image);
2034  }
2035  }
2036  for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
2037  {
2038  /*
2039  Normalize pixel statistics.
2040  */
2041  area=PerceptibleReciprocal(channel_statistics[i].area);
2042  channel_statistics[i].sum*=area;
2043  channel_statistics[i].sum_squared*=area;
2044  channel_statistics[i].sum_cubed*=area;
2045  channel_statistics[i].sum_fourth_power*=area;
2046  channel_statistics[i].mean=channel_statistics[i].sum;
2047  channel_statistics[i].variance=channel_statistics[i].sum_squared;
2048  standard_deviation=sqrt(channel_statistics[i].variance-
2049  (channel_statistics[i].mean*channel_statistics[i].mean));
2050  standard_deviation=sqrt(PerceptibleReciprocal(channel_statistics[i].area-
2051  1.0)*channel_statistics[i].area*standard_deviation*standard_deviation);
2052  channel_statistics[i].standard_deviation=standard_deviation;
2053  }
2054  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2055  {
2056  double
2057  number_bins;
2058 
2059  register ssize_t
2060  j;
2061 
2062  /*
2063  Compute pixel entropy.
2064  */
2065  PixelChannel channel = GetPixelChannelChannel(image,i);
2066  number_bins=0.0;
2067  for (j=0; j <= (ssize_t) MaxMap; j++)
2068  if (histogram[GetPixelChannels(image)*j+i] > 0.0)
2069  number_bins++;
2070  area=PerceptibleReciprocal(channel_statistics[channel].area);
2071  for (j=0; j <= (ssize_t) MaxMap; j++)
2072  {
2073  double
2074  count;
2075 
2076  count=area*histogram[GetPixelChannels(image)*j+i];
2077  channel_statistics[channel].entropy+=-count*MagickLog10(count)*
2078  PerceptibleReciprocal(MagickLog10(number_bins));
2079  channel_statistics[CompositePixelChannel].entropy+=-count*
2080  MagickLog10(count)*PerceptibleReciprocal(MagickLog10(number_bins))/
2081  GetPixelChannels(image);
2082  }
2083  }
2084  histogram=(double *) RelinquishMagickMemory(histogram);
2085  for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
2086  {
2087  /*
2088  Compute kurtosis & skewness statistics.
2089  */
2090  standard_deviation=PerceptibleReciprocal(
2091  channel_statistics[i].standard_deviation);
2092  channel_statistics[i].skewness=(channel_statistics[i].sum_cubed-3.0*
2093  channel_statistics[i].mean*channel_statistics[i].sum_squared+2.0*
2094  channel_statistics[i].mean*channel_statistics[i].mean*
2095  channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2096  standard_deviation);
2097  channel_statistics[i].kurtosis=(channel_statistics[i].sum_fourth_power-4.0*
2098  channel_statistics[i].mean*channel_statistics[i].sum_cubed+6.0*
2099  channel_statistics[i].mean*channel_statistics[i].mean*
2100  channel_statistics[i].sum_squared-3.0*channel_statistics[i].mean*
2101  channel_statistics[i].mean*1.0*channel_statistics[i].mean*
2102  channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2103  standard_deviation*standard_deviation)-3.0;
2104  }
2105  channel_statistics[CompositePixelChannel].mean=0.0;
2106  channel_statistics[CompositePixelChannel].standard_deviation=0.0;
2107  channel_statistics[CompositePixelChannel].entropy=0.0;
2108  for (i=0; i < (ssize_t) MaxPixelChannels; i++)
2109  {
2110  channel_statistics[CompositePixelChannel].mean+=
2111  channel_statistics[i].mean;
2112  channel_statistics[CompositePixelChannel].standard_deviation+=
2113  channel_statistics[i].standard_deviation;
2114  channel_statistics[CompositePixelChannel].entropy+=
2115  channel_statistics[i].entropy;
2116  }
2117  channel_statistics[CompositePixelChannel].mean/=(double)
2118  GetImageChannels(image);
2119  channel_statistics[CompositePixelChannel].standard_deviation/=(double)
2120  GetImageChannels(image);
2121  channel_statistics[CompositePixelChannel].entropy/=(double)
2122  GetImageChannels(image);
2123  if (y < (ssize_t) image->rows)
2124  channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
2125  channel_statistics);
2126  return(channel_statistics);
2127 }
2128 
2129 /*
2130 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2131 % %
2132 % %
2133 % %
2134 % P o l y n o m i a l I m a g e %
2135 % %
2136 % %
2137 % %
2138 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2139 %
2140 % PolynomialImage() returns a new image where each pixel is the sum of the
2141 % pixels in the image sequence after applying its corresponding terms
2142 % (coefficient and degree pairs).
2143 %
2144 % The format of the PolynomialImage method is:
2145 %
2146 % Image *PolynomialImage(const Image *images,const size_t number_terms,
2147 % const double *terms,ExceptionInfo *exception)
2148 %
2149 % A description of each parameter follows:
2150 %
2151 % o images: the image sequence.
2152 %
2153 % o number_terms: the number of terms in the list. The actual list length
2154 % is 2 x number_terms + 1 (the constant).
2155 %
2156 % o terms: the list of polynomial coefficients and degree pairs and a
2157 % constant.
2158 %
2159 % o exception: return any errors or warnings in this structure.
2160 %
2161 */
2162 
2164  const size_t number_terms,const double *terms,ExceptionInfo *exception)
2165 {
2166 #define PolynomialImageTag "Polynomial/Image"
2167 
2168  CacheView
2169  *polynomial_view;
2170 
2171  Image
2172  *image;
2173 
2175  status;
2176 
2178  progress;
2179 
2181  **magick_restrict polynomial_pixels;
2182 
2183  size_t
2184  number_images;
2185 
2186  ssize_t
2187  y;
2188 
2189  assert(images != (Image *) NULL);
2190  assert(images->signature == MagickCoreSignature);
2191  if (images->debug != MagickFalse)
2192  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
2193  assert(exception != (ExceptionInfo *) NULL);
2194  assert(exception->signature == MagickCoreSignature);
2195  image=AcquireImageCanvas(images,exception);
2196  if (image == (Image *) NULL)
2197  return((Image *) NULL);
2198  if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
2199  {
2200  image=DestroyImage(image);
2201  return((Image *) NULL);
2202  }
2203  number_images=GetImageListLength(images);
2204  polynomial_pixels=AcquirePixelThreadSet(images);
2205  if (polynomial_pixels == (PixelChannels **) NULL)
2206  {
2207  image=DestroyImage(image);
2208  (void) ThrowMagickException(exception,GetMagickModule(),
2209  ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
2210  return((Image *) NULL);
2211  }
2212  /*
2213  Polynomial image pixels.
2214  */
2215  status=MagickTrue;
2216  progress=0;
2217  polynomial_view=AcquireAuthenticCacheView(image,exception);
2218 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2219  #pragma omp parallel for schedule(static) shared(progress,status) \
2220  magick_number_threads(image,image,image->rows,1)
2221 #endif
2222  for (y=0; y < (ssize_t) image->rows; y++)
2223  {
2224  CacheView
2225  *image_view;
2226 
2227  const Image
2228  *next;
2229 
2230  const int
2231  id = GetOpenMPThreadId();
2232 
2233  register ssize_t
2234  i,
2235  x;
2236 
2237  register PixelChannels
2238  *polynomial_pixel;
2239 
2240  register Quantum
2241  *magick_restrict q;
2242 
2243  ssize_t
2244  j;
2245 
2246  if (status == MagickFalse)
2247  continue;
2248  q=QueueCacheViewAuthenticPixels(polynomial_view,0,y,image->columns,1,
2249  exception);
2250  if (q == (Quantum *) NULL)
2251  {
2252  status=MagickFalse;
2253  continue;
2254  }
2255  polynomial_pixel=polynomial_pixels[id];
2256  for (j=0; j < (ssize_t) image->columns; j++)
2257  for (i=0; i < MaxPixelChannels; i++)
2258  polynomial_pixel[j].channel[i]=0.0;
2259  next=images;
2260  for (j=0; j < (ssize_t) number_images; j++)
2261  {
2262  register const Quantum
2263  *p;
2264 
2265  if (j >= (ssize_t) number_terms)
2266  continue;
2267  image_view=AcquireVirtualCacheView(next,exception);
2268  p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2269  if (p == (const Quantum *) NULL)
2270  {
2271  image_view=DestroyCacheView(image_view);
2272  break;
2273  }
2274  for (x=0; x < (ssize_t) image->columns; x++)
2275  {
2276  register ssize_t
2277  i;
2278 
2279  for (i=0; i < (ssize_t) GetPixelChannels(next); i++)
2280  {
2282  coefficient,
2283  degree;
2284 
2285  PixelChannel channel = GetPixelChannelChannel(image,i);
2286  PixelTrait traits = GetPixelChannelTraits(next,channel);
2287  PixelTrait polynomial_traits=GetPixelChannelTraits(image,channel);
2288  if ((traits == UndefinedPixelTrait) ||
2289  (polynomial_traits == UndefinedPixelTrait))
2290  continue;
2291  if ((traits & UpdatePixelTrait) == 0)
2292  continue;
2293  coefficient=(MagickRealType) terms[2*j];
2294  degree=(MagickRealType) terms[(j << 1)+1];
2295  polynomial_pixel[x].channel[i]+=coefficient*
2296  pow(QuantumScale*GetPixelChannel(image,channel,p),degree);
2297  }
2298  p+=GetPixelChannels(next);
2299  }
2300  image_view=DestroyCacheView(image_view);
2301  next=GetNextImageInList(next);
2302  }
2303  for (x=0; x < (ssize_t) image->columns; x++)
2304  {
2305  register ssize_t
2306  i;
2307 
2308  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2309  {
2310  PixelChannel channel = GetPixelChannelChannel(image,i);
2311  PixelTrait traits = GetPixelChannelTraits(image,channel);
2312  if (traits == UndefinedPixelTrait)
2313  continue;
2314  if ((traits & UpdatePixelTrait) == 0)
2315  continue;
2316  q[i]=ClampToQuantum(QuantumRange*polynomial_pixel[x].channel[i]);
2317  }
2318  q+=GetPixelChannels(image);
2319  }
2320  if (SyncCacheViewAuthenticPixels(polynomial_view,exception) == MagickFalse)
2321  status=MagickFalse;
2322  if (images->progress_monitor != (MagickProgressMonitor) NULL)
2323  {
2325  proceed;
2326 
2327 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2328  #pragma omp critical (MagickCore_PolynomialImages)
2329 #endif
2330  proceed=SetImageProgress(images,PolynomialImageTag,progress++,
2331  image->rows);
2332  if (proceed == MagickFalse)
2333  status=MagickFalse;
2334  }
2335  }
2336  polynomial_view=DestroyCacheView(polynomial_view);
2337  polynomial_pixels=DestroyPixelThreadSet(polynomial_pixels);
2338  if (status == MagickFalse)
2339  image=DestroyImage(image);
2340  return(image);
2341 }
2342 
2343 /*
2344 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2345 % %
2346 % %
2347 % %
2348 % S t a t i s t i c I m a g e %
2349 % %
2350 % %
2351 % %
2352 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2353 %
2354 % StatisticImage() makes each pixel the min / max / median / mode / etc. of
2355 % the neighborhood of the specified width and height.
2356 %
2357 % The format of the StatisticImage method is:
2358 %
2359 % Image *StatisticImage(const Image *image,const StatisticType type,
2360 % const size_t width,const size_t height,ExceptionInfo *exception)
2361 %
2362 % A description of each parameter follows:
2363 %
2364 % o image: the image.
2365 %
2366 % o type: the statistic type (median, mode, etc.).
2367 %
2368 % o width: the width of the pixel neighborhood.
2369 %
2370 % o height: the height of the pixel neighborhood.
2371 %
2372 % o exception: return any errors or warnings in this structure.
2373 %
2374 */
2375 
2376 typedef struct _SkipNode
2377 {
2378  size_t
2379  next[9],
2380  count,
2381  signature;
2382 } SkipNode;
2383 
2384 typedef struct _SkipList
2385 {
2386  ssize_t
2388 
2389  SkipNode
2391 } SkipList;
2392 
2393 typedef struct _PixelList
2394 {
2395  size_t
2397  seed;
2398 
2399  SkipList
2401 
2402  size_t
2404 } PixelList;
2405 
2407 {
2408  if (pixel_list == (PixelList *) NULL)
2409  return((PixelList *) NULL);
2410  if (pixel_list->skip_list.nodes != (SkipNode *) NULL)
2412  pixel_list->skip_list.nodes);
2413  pixel_list=(PixelList *) RelinquishMagickMemory(pixel_list);
2414  return(pixel_list);
2415 }
2416 
2418 {
2419  register ssize_t
2420  i;
2421 
2422  assert(pixel_list != (PixelList **) NULL);
2423  for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
2424  if (pixel_list[i] != (PixelList *) NULL)
2425  pixel_list[i]=DestroyPixelList(pixel_list[i]);
2426  pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list);
2427  return(pixel_list);
2428 }
2429 
2430 static PixelList *AcquirePixelList(const size_t width,const size_t height)
2431 {
2432  PixelList
2433  *pixel_list;
2434 
2435  pixel_list=(PixelList *) AcquireMagickMemory(sizeof(*pixel_list));
2436  if (pixel_list == (PixelList *) NULL)
2437  return(pixel_list);
2438  (void) memset((void *) pixel_list,0,sizeof(*pixel_list));
2439  pixel_list->length=width*height;
2440  pixel_list->skip_list.nodes=(SkipNode *) AcquireAlignedMemory(65537UL,
2441  sizeof(*pixel_list->skip_list.nodes));
2442  if (pixel_list->skip_list.nodes == (SkipNode *) NULL)
2443  return(DestroyPixelList(pixel_list));
2444  (void) memset(pixel_list->skip_list.nodes,0,65537UL*
2445  sizeof(*pixel_list->skip_list.nodes));
2446  pixel_list->signature=MagickCoreSignature;
2447  return(pixel_list);
2448 }
2449 
2450 static PixelList **AcquirePixelListThreadSet(const size_t width,
2451  const size_t height)
2452 {
2453  PixelList
2454  **pixel_list;
2455 
2456  register ssize_t
2457  i;
2458 
2459  size_t
2460  number_threads;
2461 
2462  number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
2463  pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
2464  sizeof(*pixel_list));
2465  if (pixel_list == (PixelList **) NULL)
2466  return((PixelList **) NULL);
2467  (void) memset(pixel_list,0,number_threads*sizeof(*pixel_list));
2468  for (i=0; i < (ssize_t) number_threads; i++)
2469  {
2470  pixel_list[i]=AcquirePixelList(width,height);
2471  if (pixel_list[i] == (PixelList *) NULL)
2472  return(DestroyPixelListThreadSet(pixel_list));
2473  }
2474  return(pixel_list);
2475 }
2476 
2477 static void AddNodePixelList(PixelList *pixel_list,const size_t color)
2478 {
2479  register SkipList
2480  *p;
2481 
2482  register ssize_t
2483  level;
2484 
2485  size_t
2486  search,
2487  update[9];
2488 
2489  /*
2490  Initialize the node.
2491  */
2492  p=(&pixel_list->skip_list);
2493  p->nodes[color].signature=pixel_list->signature;
2494  p->nodes[color].count=1;
2495  /*
2496  Determine where it belongs in the list.
2497  */
2498  search=65536UL;
2499  for (level=p->level; level >= 0; level--)
2500  {
2501  while (p->nodes[search].next[level] < color)
2502  search=p->nodes[search].next[level];
2503  update[level]=search;
2504  }
2505  /*
2506  Generate a pseudo-random level for this node.
2507  */
2508  for (level=0; ; level++)
2509  {
2510  pixel_list->seed=(pixel_list->seed*42893621L)+1L;
2511  if ((pixel_list->seed & 0x300) != 0x300)
2512  break;
2513  }
2514  if (level > 8)
2515  level=8;
2516  if (level > (p->level+2))
2517  level=p->level+2;
2518  /*
2519  If we're raising the list's level, link back to the root node.
2520  */
2521  while (level > p->level)
2522  {
2523  p->level++;
2524  update[p->level]=65536UL;
2525  }
2526  /*
2527  Link the node into the skip-list.
2528  */
2529  do
2530  {
2531  p->nodes[color].next[level]=p->nodes[update[level]].next[level];
2532  p->nodes[update[level]].next[level]=color;
2533  } while (level-- > 0);
2534 }
2535 
2536 static inline void GetMaximumPixelList(PixelList *pixel_list,Quantum *pixel)
2537 {
2538  register SkipList
2539  *p;
2540 
2541  size_t
2542  color,
2543  maximum;
2544 
2545  ssize_t
2546  count;
2547 
2548  /*
2549  Find the maximum value for each of the color.
2550  */
2551  p=(&pixel_list->skip_list);
2552  color=65536L;
2553  count=0;
2554  maximum=p->nodes[color].next[0];
2555  do
2556  {
2557  color=p->nodes[color].next[0];
2558  if (color > maximum)
2559  maximum=color;
2560  count+=p->nodes[color].count;
2561  } while (count < (ssize_t) pixel_list->length);
2562  *pixel=ScaleShortToQuantum((unsigned short) maximum);
2563 }
2564 
2565 static inline void GetMeanPixelList(PixelList *pixel_list,Quantum *pixel)
2566 {
2567  double
2568  sum;
2569 
2570  register SkipList
2571  *p;
2572 
2573  size_t
2574  color;
2575 
2576  ssize_t
2577  count;
2578 
2579  /*
2580  Find the mean value for each of the color.
2581  */
2582  p=(&pixel_list->skip_list);
2583  color=65536L;
2584  count=0;
2585  sum=0.0;
2586  do
2587  {
2588  color=p->nodes[color].next[0];
2589  sum+=(double) p->nodes[color].count*color;
2590  count+=p->nodes[color].count;
2591  } while (count < (ssize_t) pixel_list->length);
2592  sum/=pixel_list->length;
2593  *pixel=ScaleShortToQuantum((unsigned short) sum);
2594 }
2595 
2596 static inline void GetMedianPixelList(PixelList *pixel_list,Quantum *pixel)
2597 {
2598  register SkipList
2599  *p;
2600 
2601  size_t
2602  color;
2603 
2604  ssize_t
2605  count;
2606 
2607  /*
2608  Find the median value for each of the color.
2609  */
2610  p=(&pixel_list->skip_list);
2611  color=65536L;
2612  count=0;
2613  do
2614  {
2615  color=p->nodes[color].next[0];
2616  count+=p->nodes[color].count;
2617  } while (count <= (ssize_t) (pixel_list->length >> 1));
2618  *pixel=ScaleShortToQuantum((unsigned short) color);
2619 }
2620 
2621 static inline void GetMinimumPixelList(PixelList *pixel_list,Quantum *pixel)
2622 {
2623  register SkipList
2624  *p;
2625 
2626  size_t
2627  color,
2628  minimum;
2629 
2630  ssize_t
2631  count;
2632 
2633  /*
2634  Find the minimum value for each of the color.
2635  */
2636  p=(&pixel_list->skip_list);
2637  count=0;
2638  color=65536UL;
2639  minimum=p->nodes[color].next[0];
2640  do
2641  {
2642  color=p->nodes[color].next[0];
2643  if (color < minimum)
2644  minimum=color;
2645  count+=p->nodes[color].count;
2646  } while (count < (ssize_t) pixel_list->length);
2647  *pixel=ScaleShortToQuantum((unsigned short) minimum);
2648 }
2649 
2650 static inline void GetModePixelList(PixelList *pixel_list,Quantum *pixel)
2651 {
2652  register SkipList
2653  *p;
2654 
2655  size_t
2656  color,
2657  max_count,
2658  mode;
2659 
2660  ssize_t
2661  count;
2662 
2663  /*
2664  Make each pixel the 'predominant color' of the specified neighborhood.
2665  */
2666  p=(&pixel_list->skip_list);
2667  color=65536L;
2668  mode=color;
2669  max_count=p->nodes[mode].count;
2670  count=0;
2671  do
2672  {
2673  color=p->nodes[color].next[0];
2674  if (p->nodes[color].count > max_count)
2675  {
2676  mode=color;
2677  max_count=p->nodes[mode].count;
2678  }
2679  count+=p->nodes[color].count;
2680  } while (count < (ssize_t) pixel_list->length);
2681  *pixel=ScaleShortToQuantum((unsigned short) mode);
2682 }
2683 
2684 static inline void GetNonpeakPixelList(PixelList *pixel_list,Quantum *pixel)
2685 {
2686  register SkipList
2687  *p;
2688 
2689  size_t
2690  color,
2691  next,
2692  previous;
2693 
2694  ssize_t
2695  count;
2696 
2697  /*
2698  Finds the non peak value for each of the colors.
2699  */
2700  p=(&pixel_list->skip_list);
2701  color=65536L;
2702  next=p->nodes[color].next[0];
2703  count=0;
2704  do
2705  {
2706  previous=color;
2707  color=next;
2708  next=p->nodes[color].next[0];
2709  count+=p->nodes[color].count;
2710  } while (count <= (ssize_t) (pixel_list->length >> 1));
2711  if ((previous == 65536UL) && (next != 65536UL))
2712  color=next;
2713  else
2714  if ((previous != 65536UL) && (next == 65536UL))
2715  color=previous;
2716  *pixel=ScaleShortToQuantum((unsigned short) color);
2717 }
2718 
2719 static inline void GetRootMeanSquarePixelList(PixelList *pixel_list,
2720  Quantum *pixel)
2721 {
2722  double
2723  sum;
2724 
2725  register SkipList
2726  *p;
2727 
2728  size_t
2729  color;
2730 
2731  ssize_t
2732  count;
2733 
2734  /*
2735  Find the root mean square value for each of the color.
2736  */
2737  p=(&pixel_list->skip_list);
2738  color=65536L;
2739  count=0;
2740  sum=0.0;
2741  do
2742  {
2743  color=p->nodes[color].next[0];
2744  sum+=(double) (p->nodes[color].count*color*color);
2745  count+=p->nodes[color].count;
2746  } while (count < (ssize_t) pixel_list->length);
2747  sum/=pixel_list->length;
2748  *pixel=ScaleShortToQuantum((unsigned short) sqrt(sum));
2749 }
2750 
2751 static inline void GetStandardDeviationPixelList(PixelList *pixel_list,
2752  Quantum *pixel)
2753 {
2754  double
2755  sum,
2756  sum_squared;
2757 
2758  register SkipList
2759  *p;
2760 
2761  size_t
2762  color;
2763 
2764  ssize_t
2765  count;
2766 
2767  /*
2768  Find the standard-deviation value for each of the color.
2769  */
2770  p=(&pixel_list->skip_list);
2771  color=65536L;
2772  count=0;
2773  sum=0.0;
2774  sum_squared=0.0;
2775  do
2776  {
2777  register ssize_t
2778  i;
2779 
2780  color=p->nodes[color].next[0];
2781  sum+=(double) p->nodes[color].count*color;
2782  for (i=0; i < (ssize_t) p->nodes[color].count; i++)
2783  sum_squared+=((double) color)*((double) color);
2784  count+=p->nodes[color].count;
2785  } while (count < (ssize_t) pixel_list->length);
2786  sum/=pixel_list->length;
2787  sum_squared/=pixel_list->length;
2788  *pixel=ScaleShortToQuantum((unsigned short) sqrt(sum_squared-(sum*sum)));
2789 }
2790 
2791 static inline void InsertPixelList(const Quantum pixel,PixelList *pixel_list)
2792 {
2793  size_t
2794  signature;
2795 
2796  unsigned short
2797  index;
2798 
2799  index=ScaleQuantumToShort(pixel);
2800  signature=pixel_list->skip_list.nodes[index].signature;
2801  if (signature == pixel_list->signature)
2802  {
2803  pixel_list->skip_list.nodes[index].count++;
2804  return;
2805  }
2806  AddNodePixelList(pixel_list,index);
2807 }
2808 
2809 static void ResetPixelList(PixelList *pixel_list)
2810 {
2811  int
2812  level;
2813 
2814  register SkipNode
2815  *root;
2816 
2817  register SkipList
2818  *p;
2819 
2820  /*
2821  Reset the skip-list.
2822  */
2823  p=(&pixel_list->skip_list);
2824  root=p->nodes+65536UL;
2825  p->level=0;
2826  for (level=0; level < 9; level++)
2827  root->next[level]=65536UL;
2828  pixel_list->seed=pixel_list->signature++;
2829 }
2830 
2832  const size_t width,const size_t height,ExceptionInfo *exception)
2833 {
2834 #define StatisticImageTag "Statistic/Image"
2835 
2836  CacheView
2837  *image_view,
2838  *statistic_view;
2839 
2840  Image
2841  *statistic_image;
2842 
2844  status;
2845 
2847  progress;
2848 
2849  PixelList
2850  **magick_restrict pixel_list;
2851 
2852  ssize_t
2853  center,
2854  y;
2855 
2856  /*
2857  Initialize statistics image attributes.
2858  */
2859  assert(image != (Image *) NULL);
2860  assert(image->signature == MagickCoreSignature);
2861  if (image->debug != MagickFalse)
2862  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2863  assert(exception != (ExceptionInfo *) NULL);
2864  assert(exception->signature == MagickCoreSignature);
2865  statistic_image=CloneImage(image,0,0,MagickTrue,
2866  exception);
2867  if (statistic_image == (Image *) NULL)
2868  return((Image *) NULL);
2869  status=SetImageStorageClass(statistic_image,DirectClass,exception);
2870  if (status == MagickFalse)
2871  {
2872  statistic_image=DestroyImage(statistic_image);
2873  return((Image *) NULL);
2874  }
2875  pixel_list=AcquirePixelListThreadSet(MagickMax(width,1),MagickMax(height,1));
2876  if (pixel_list == (PixelList **) NULL)
2877  {
2878  statistic_image=DestroyImage(statistic_image);
2879  ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2880  }
2881  /*
2882  Make each pixel the min / max / median / mode / etc. of the neighborhood.
2883  */
2884  center=(ssize_t) GetPixelChannels(image)*(image->columns+MagickMax(width,1))*
2885  (MagickMax(height,1)/2L)+GetPixelChannels(image)*(MagickMax(width,1)/2L);
2886  status=MagickTrue;
2887  progress=0;
2888  image_view=AcquireVirtualCacheView(image,exception);
2889  statistic_view=AcquireAuthenticCacheView(statistic_image,exception);
2890 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2891  #pragma omp parallel for schedule(static) shared(progress,status) \
2892  magick_number_threads(image,statistic_image,statistic_image->rows,1)
2893 #endif
2894  for (y=0; y < (ssize_t) statistic_image->rows; y++)
2895  {
2896  const int
2897  id = GetOpenMPThreadId();
2898 
2899  register const Quantum
2900  *magick_restrict p;
2901 
2902  register Quantum
2903  *magick_restrict q;
2904 
2905  register ssize_t
2906  x;
2907 
2908  if (status == MagickFalse)
2909  continue;
2910  p=GetCacheViewVirtualPixels(image_view,-((ssize_t) MagickMax(width,1)/2L),y-
2911  (ssize_t) (MagickMax(height,1)/2L),image->columns+MagickMax(width,1),
2912  MagickMax(height,1),exception);
2913  q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns, 1,exception);
2914  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2915  {
2916  status=MagickFalse;
2917  continue;
2918  }
2919  for (x=0; x < (ssize_t) statistic_image->columns; x++)
2920  {
2921  register ssize_t
2922  i;
2923 
2924  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2925  {
2926  Quantum
2927  pixel;
2928 
2929  register const Quantum
2930  *magick_restrict pixels;
2931 
2932  register ssize_t
2933  u;
2934 
2935  ssize_t
2936  v;
2937 
2938  PixelChannel channel = GetPixelChannelChannel(image,i);
2939  PixelTrait traits = GetPixelChannelTraits(image,channel);
2940  PixelTrait statistic_traits=GetPixelChannelTraits(statistic_image,
2941  channel);
2942  if ((traits == UndefinedPixelTrait) ||
2943  (statistic_traits == UndefinedPixelTrait))
2944  continue;
2945  if (((statistic_traits & CopyPixelTrait) != 0) ||
2946  (GetPixelWriteMask(image,p) <= (QuantumRange/2)))
2947  {
2948  SetPixelChannel(statistic_image,channel,p[center+i],q);
2949  continue;
2950  }
2951  if ((statistic_traits & UpdatePixelTrait) == 0)
2952  continue;
2953  pixels=p;
2954  ResetPixelList(pixel_list[id]);
2955  for (v=0; v < (ssize_t) MagickMax(height,1); v++)
2956  {
2957  for (u=0; u < (ssize_t) MagickMax(width,1); u++)
2958  {
2959  InsertPixelList(pixels[i],pixel_list[id]);
2960  pixels+=GetPixelChannels(image);
2961  }
2962  pixels+=GetPixelChannels(image)*image->columns;
2963  }
2964  switch (type)
2965  {
2966  case GradientStatistic:
2967  {
2968  double
2969  maximum,
2970  minimum;
2971 
2972  GetMinimumPixelList(pixel_list[id],&pixel);
2973  minimum=(double) pixel;
2974  GetMaximumPixelList(pixel_list[id],&pixel);
2975  maximum=(double) pixel;
2976  pixel=ClampToQuantum(MagickAbsoluteValue(maximum-minimum));
2977  break;
2978  }
2979  case MaximumStatistic:
2980  {
2981  GetMaximumPixelList(pixel_list[id],&pixel);
2982  break;
2983  }
2984  case MeanStatistic:
2985  {
2986  GetMeanPixelList(pixel_list[id],&pixel);
2987  break;
2988  }
2989  case MedianStatistic:
2990  default:
2991  {
2992  GetMedianPixelList(pixel_list[id],&pixel);
2993  break;
2994  }
2995  case MinimumStatistic:
2996  {
2997  GetMinimumPixelList(pixel_list[id],&pixel);
2998  break;
2999  }
3000  case ModeStatistic:
3001  {
3002  GetModePixelList(pixel_list[id],&pixel);
3003  break;
3004  }
3005  case NonpeakStatistic:
3006  {
3007  GetNonpeakPixelList(pixel_list[id],&pixel);
3008  break;
3009  }
3011  {
3012  GetRootMeanSquarePixelList(pixel_list[id],&pixel);
3013  break;
3014  }
3016  {
3017  GetStandardDeviationPixelList(pixel_list[id],&pixel);
3018  break;
3019  }
3020  }
3021  SetPixelChannel(statistic_image,channel,pixel,q);
3022  }
3023  p+=GetPixelChannels(image);
3024  q+=GetPixelChannels(statistic_image);
3025  }
3026  if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
3027  status=MagickFalse;
3028  if (image->progress_monitor != (MagickProgressMonitor) NULL)
3029  {
3031  proceed;
3032 
3033 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3034  #pragma omp critical (MagickCore_StatisticImage)
3035 #endif
3036  proceed=SetImageProgress(image,StatisticImageTag,progress++,
3037  image->rows);
3038  if (proceed == MagickFalse)
3039  status=MagickFalse;
3040  }
3041  }
3042  statistic_view=DestroyCacheView(statistic_view);
3043  image_view=DestroyCacheView(image_view);
3044  pixel_list=DestroyPixelListThreadSet(pixel_list);
3045  if (status == MagickFalse)
3046  statistic_image=DestroyImage(statistic_image);
3047  return(statistic_image);
3048 }
size_t rows
Definition: image.h:172
#define magick_restrict
Definition: MagickCore.h:41
MagickExport Image * BlurImage(const Image *image, const double radius, const double sigma, ExceptionInfo *exception)
Definition: effect.c:768
MagickDoubleType MagickRealType
Definition: magick-type.h:118
MagickExport CacheView * DestroyCacheView(CacheView *cache_view)
Definition: cache-view.c:252
double channel[CompositePixelChannel]
Definition: statistic.c:137
StatisticType
Definition: statistic.h:130
static MagickBooleanType SetImageProgress(const Image *image, const char *tag, const MagickOffsetType offset, const MagickSizeType extent)
double standard_deviation
Definition: statistic.h:35
static MagickSizeType GetQuantumRange(const size_t depth)
MagickProgressMonitor progress_monitor
Definition: image.h:303
MagickExport MagickBooleanType TransformImageColorspace(Image *image, const ColorspaceType colorspace, ExceptionInfo *exception)
Definition: colorspace.c:1325
double ellipse_angle
Definition: statistic.h:60
#define Log10Epsilon
MagickExport ssize_t ParseCommandOption(const CommandOption option, const MagickBooleanType list, const char *options)
Definition: option.c:2957
MagickExport MagickBooleanType FunctionImage(Image *image, const MagickFunction function, const size_t number_parameters, const double *parameters, ExceptionInfo *exception)
Definition: statistic.c:1006
MagickExport ChannelStatistics * GetImageStatistics(const Image *image, ExceptionInfo *exception)
Definition: statistic.c:1914
MagickExport MagickBooleanType EvaluateImage(Image *image, const MagickEvaluateOperator op, const double value, ExceptionInfo *exception)
Definition: statistic.c:770
static PixelChannels ** DestroyPixelThreadSet(PixelChannels **pixels)
Definition: statistic.c:140
size_t signature
Definition: exception.h:123
struct _SkipNode SkipNode
static double ApplyEvaluateOperator(RandomInfo *random_info, const Quantum pixel, const MagickEvaluateOperator op, const double value)
Definition: statistic.c:226
MagickPrivate double GenerateDifferentialNoise(RandomInfo *, const Quantum, const NoiseType, const double)
Definition: gem.c:1455
static Quantum GetPixelChannel(const Image *magick_restrict image, const PixelChannel channel, const Quantum *magick_restrict pixel)
static PixelList * DestroyPixelList(PixelList *pixel_list)
Definition: statistic.c:2406
static void GetMaximumPixelList(PixelList *pixel_list, Quantum *pixel)
Definition: statistic.c:2536
MagickExport const char * GetImageArtifact(const Image *image, const char *artifact)
Definition: artifact.c:273
static int IntensityCompare(const void *x, const void *y)
Definition: statistic.c:202
struct _PixelChannels PixelChannels
static RandomInfo ** DestroyRandomInfoThreadSet(RandomInfo **random_info)
static PixelTrait GetPixelChannelTraits(const Image *magick_restrict image, const PixelChannel channel)
#define MagickPI
Definition: image-private.h:30
#define MagickAbsoluteValue(x)
Definition: image-private.h:25
#define MAGICKCORE_QUANTUM_DEPTH
Definition: magick-type.h:28
MagickExport const Quantum * GetCacheViewVirtualPixels(const CacheView *cache_view, const ssize_t x, const ssize_t y, const size_t columns, const size_t rows, ExceptionInfo *exception)
Definition: cache-view.c:651
static RandomInfo ** AcquireRandomInfoThreadSet(void)
#define PolynomialImageTag
size_t count
Definition: statistic.c:2379
#define MaximumNumberOfImageMoments
Definition: statistic.h:25
static void GetMeanPixelList(PixelList *pixel_list, Quantum *pixel)
Definition: statistic.c:2565
#define MagickEpsilon
Definition: magick-type.h:110
MagickExport MagickBooleanType GetImageEntropy(const Image *image, double *entropy, ExceptionInfo *exception)
Definition: statistic.c:1124
Definition: log.h:52
ssize_t MagickOffsetType
Definition: magick-type.h:127
MagickExport unsigned long GetRandomSecretKey(const RandomInfo *random_info)
Definition: random.c:743
Definition: image.h:151
static double EvaluateMax(const double x, const double y)
Definition: statistic.c:191
size_t length
Definition: statistic.c:2396
#define EvaluateImageTag
double x
Definition: geometry.h:123
#define MagickCoreSignature
static double MagickLog10(const double x)
Definition: statistic.c:1662
MagickExport Quantum * GetCacheViewAuthenticPixels(CacheView *cache_view, const ssize_t x, const ssize_t y, const size_t columns, const size_t rows, ExceptionInfo *exception)
Definition: cache-view.c:299
SkipNode * nodes
Definition: statistic.c:2390
static Image * AcquireImageCanvas(const Image *images, ExceptionInfo *exception)
Definition: statistic.c:419
double invariant[MaximumNumberOfImageMoments+1]
Definition: statistic.h:53
static void GetNonpeakPixelList(PixelList *pixel_list, Quantum *pixel)
Definition: statistic.c:2684
MagickBooleanType
Definition: magick-type.h:156
unsigned int MagickStatusType
Definition: magick-type.h:119
MagickExport char * AcquireString(const char *source)
Definition: string.c:129
double ellipse_intensity
Definition: statistic.h:60
static double PerceptibleReciprocal(const double x)
static Quantum ScaleAnyToQuantum(const QuantumAny quantum, const QuantumAny range)
size_t signature
Definition: statistic.c:2379
size_t signature
Definition: statistic.c:2403
MagickEvaluateOperator
Definition: statistic.h:84
#define MaximumNumberOfPerceptualColorspaces
Definition: statistic.h:26
static Quantum GetPixelWriteMask(const Image *magick_restrict image, const Quantum *magick_restrict pixel)
MagickExport const Quantum * GetVirtualPixels(const Image *image, const ssize_t x, const ssize_t y, const size_t columns, const size_t rows, ExceptionInfo *exception)
Definition: cache.c:3252
static Quantum ApplyFunction(Quantum pixel, const MagickFunction function, const size_t number_parameters, const double *parameters, ExceptionInfo *exception)
Definition: statistic.c:908
MagickExport void * AcquireQuantumMemory(const size_t count, const size_t quantum)
Definition: memory.c:533
double y
Definition: geometry.h:123
static PixelList ** DestroyPixelListThreadSet(PixelList **pixel_list)
Definition: statistic.c:2417
static int GetOpenMPThreadId(void)
MagickExport void * RelinquishAlignedMemory(void *memory)
Definition: memory.c:1016
MagickExport ChannelMoments * GetImageMoments(const Image *image, ExceptionInfo *exception)
Definition: statistic.c:1342
#define MagickMaximumValue
Definition: magick-type.h:111
struct _PixelList PixelList
MagickExport Quantum * QueueCacheViewAuthenticPixels(CacheView *cache_view, const ssize_t x, const ssize_t y, const size_t columns, const size_t rows, ExceptionInfo *exception)
Definition: cache-view.c:977
MagickExport MagickBooleanType ThrowMagickException(ExceptionInfo *exception, const char *module, const char *function, const size_t line, const ExceptionType severity, const char *tag, const char *format,...)
Definition: exception.c:1064
MagickExport MagickBooleanType LogMagickEvent(const LogEventType type, const char *module, const char *function, const size_t line, const char *format,...)
Definition: log.c:1397
size_t signature
Definition: image.h:354
#define QuantumScale
Definition: magick-type.h:113
size_t columns
Definition: image.h:172
MagickExport MagickSizeType GetMagickResourceLimit(const ResourceType type)
Definition: resource.c:745
static void GetMinimumPixelList(PixelList *pixel_list, Quantum *pixel)
Definition: statistic.c:2621
MagickBooleanType(* MagickProgressMonitor)(const char *, const MagickOffsetType, const MagickSizeType, void *)
Definition: monitor.h:26
static void InsertPixelList(const Quantum pixel, PixelList *pixel_list)
Definition: statistic.c:2791
struct _Image * next
Definition: image.h:348
static PixelList * AcquirePixelList(const size_t width, const size_t height)
Definition: statistic.c:2430
MagickExport MagickBooleanType GetImageRange(const Image *image, double *minima, double *maxima, ExceptionInfo *exception)
Definition: statistic.c:1779
MagickExport MagickBooleanType SetImageStorageClass(Image *image, const ClassType storage_class, ExceptionInfo *exception)
Definition: image.c:2613
MagickExport MagickBooleanType GetImageKurtosis(const Image *image, double *kurtosis, double *skewness, ExceptionInfo *exception)
Definition: statistic.c:1222
PixelChannel
Definition: pixel.h:67
MagickExport void * AcquireAlignedMemory(const size_t count, const size_t quantum)
Definition: memory.c:242
MagickExport MagickBooleanType GetImageMean(const Image *image, double *mean, double *standard_deviation, ExceptionInfo *exception)
Definition: statistic.c:1272
#define MaxMap
Definition: magick-type.h:75
#define MagickMax(x, y)
Definition: image-private.h:26
static size_t GetPixelChannels(const Image *magick_restrict image)
static PixelChannels ** AcquirePixelThreadSet(const Image *image)
Definition: statistic.c:153
char filename[MagickPathExtent]
Definition: image.h:319
MagickExport MagickBooleanType GetImageExtrema(const Image *image, size_t *minima, size_t *maxima, ExceptionInfo *exception)
Definition: statistic.c:1172
#define GetMagickModule()
Definition: log.h:28
size_t seed
Definition: statistic.c:2396
ssize_t level
Definition: statistic.c:2387
#define ThrowImageException(severity, tag)
static Quantum ClampToQuantum(const MagickRealType value)
Definition: quantum.h:84
static PixelChannel GetPixelChannelChannel(const Image *magick_restrict image, const ssize_t offset)
MagickExport CacheView * AcquireVirtualCacheView(const Image *image, ExceptionInfo *exception)
Definition: cache-view.c:149
static double RadiansToDegrees(const double radians)
Definition: image-private.h:61
static void GetMedianPixelList(PixelList *pixel_list, Quantum *pixel)
Definition: statistic.c:2596
PointInfo centroid
Definition: statistic.h:56
MagickExport char * StringToken(const char *delimiters, char **string)
Definition: string.c:2314
static void GetModePixelList(PixelList *pixel_list, Quantum *pixel)
Definition: statistic.c:2650
unsigned short Quantum
Definition: magick-type.h:82
MagickExport Image * GetNextImageInList(const Image *images)
Definition: list.c:755
MagickExport char * DestroyString(char *string)
Definition: string.c:823
MagickExport void * AcquireMagickMemory(const size_t size)
Definition: memory.c:462
static size_t GetImageChannels(const Image *image)
Definition: statistic.c:1320
size_t number_channels
Definition: image.h:283
static void SetPixelChannel(const Image *magick_restrict image, const PixelChannel channel, const Quantum quantum, Quantum *magick_restrict pixel)
SkipList skip_list
Definition: statistic.c:2400
static void GetStandardDeviationPixelList(PixelList *pixel_list, Quantum *pixel)
Definition: statistic.c:2751
ColorspaceType colorspace[MaximumNumberOfPerceptualColorspaces+1]
Definition: statistic.h:75
#define StatisticImageTag
#define FunctionImageTag
#define MagickMin(x, y)
Definition: image-private.h:27
ColorspaceType
Definition: colorspace.h:25
static RandomInfo * random_info
Definition: resource.c:111
size_t next[9]
Definition: statistic.c:2379
PointInfo ellipse_axis
Definition: statistic.h:56
MagickExport void * RelinquishMagickMemory(void *memory)
Definition: memory.c:1053
#define MaxPixelChannels
Definition: pixel.h:27
double sum_squared
Definition: statistic.h:35
double ellipse_eccentricity
Definition: statistic.h:60
#define MagickExport
static void AddNodePixelList(PixelList *pixel_list, const size_t color)
Definition: statistic.c:2477
MagickExport MagickBooleanType SyncCacheViewAuthenticPixels(CacheView *magick_restrict cache_view, ExceptionInfo *exception)
Definition: cache-view.c:1100
MagickExport CacheView * AcquireAuthenticCacheView(const Image *image, ExceptionInfo *exception)
Definition: cache-view.c:112
static void ResetPixelList(PixelList *pixel_list)
Definition: statistic.c:2809
MagickExport Image * PolynomialImage(const Image *images, const size_t number_terms, const double *terms, ExceptionInfo *exception)
Definition: statistic.c:2163
static void GetRootMeanSquarePixelList(PixelList *pixel_list, Quantum *pixel)
Definition: statistic.c:2719
PixelTrait
Definition: pixel.h:135
MagickExport ChannelPerceptualHash * GetImagePerceptualHash(const Image *image, ExceptionInfo *exception)
Definition: statistic.c:1671
MagickFunction
Definition: statistic.h:121
static QuantumAny ScaleQuantumToAny(const Quantum quantum, const QuantumAny range)
MagickExport size_t GetImageListLength(const Image *images)
Definition: list.c:680
MagickSizeType QuantumAny
Definition: magick-type.h:142
static PixelList ** AcquirePixelListThreadSet(const size_t width, const size_t height)
Definition: statistic.c:2450
MagickExport Image * DestroyImage(Image *image)
Definition: image.c:1177
MagickExport Image * CloneImage(const Image *image, const size_t columns, const size_t rows, const MagickBooleanType detach, ExceptionInfo *exception)
Definition: image.c:794
MagickExport Image * EvaluateImages(const Image *images, const MagickEvaluateOperator op, ExceptionInfo *exception)
Definition: statistic.c:444
MagickExport Image * StatisticImage(const Image *image, const StatisticType type, const size_t width, const size_t height, ExceptionInfo *exception)
Definition: statistic.c:2831
struct _SkipList SkipList
#define QuantumRange
Definition: magick-type.h:83
MagickBooleanType debug
Definition: image.h:334
double sum_fourth_power
Definition: statistic.h:35
size_t depth
Definition: image.h:172