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