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