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) memset(pixels,0,number_threads*sizeof(*pixels));
170  for (i=0; i < (ssize_t) number_threads; i++)
171  {
172  register ssize_t
173  j;
174 
175  pixels[i]=(PixelChannels *) AcquireQuantumMemory(image->columns,
176  sizeof(**pixels));
177  if (pixels[i] == (PixelChannels *) NULL)
178  return(DestroyPixelThreadSet(pixels));
179  for (j=0; j < (ssize_t) image->columns; j++)
180  {
181  register ssize_t
182  k;
183 
184  for (k=0; k < MaxPixelChannels; k++)
185  pixels[i][j].channel[k]=0.0;
186  }
187  }
188  return(pixels);
189 }
190 
191 static inline double EvaluateMax(const double x,const double y)
192 {
193  if (x > y)
194  return(x);
195  return(y);
196 }
197 
198 #if defined(__cplusplus) || defined(c_plusplus)
199 extern "C" {
200 #endif
201 
202 static int IntensityCompare(const void *x,const void *y)
203 {
204  const PixelChannels
205  *color_1,
206  *color_2;
207 
208  double
209  distance;
210 
211  register ssize_t
212  i;
213 
214  color_1=(const PixelChannels *) x;
215  color_2=(const PixelChannels *) y;
216  distance=0.0;
217  for (i=0; i < MaxPixelChannels; i++)
218  distance+=color_1->channel[i]-(double) color_2->channel[i];
219  return(distance < 0 ? -1 : distance > 0 ? 1 : 0);
220 }
221 
222 #if defined(__cplusplus) || defined(c_plusplus)
223 }
224 #endif
225 
227  const MagickEvaluateOperator op,const double value)
228 {
229  double
230  result;
231 
232  result=0.0;
233  switch (op)
234  {
236  break;
237  case AbsEvaluateOperator:
238  {
239  result=(double) fabs((double) (pixel+value));
240  break;
241  }
242  case AddEvaluateOperator:
243  {
244  result=(double) (pixel+value);
245  break;
246  }
248  {
249  /*
250  This returns a 'floored modulus' of the addition which is a positive
251  result. It differs from % or fmod() that returns a 'truncated modulus'
252  result, where floor() is replaced by trunc() and could return a
253  negative result (which is clipped).
254  */
255  result=pixel+value;
256  result-=(QuantumRange+1.0)*floor((double) result/(QuantumRange+1.0));
257  break;
258  }
259  case AndEvaluateOperator:
260  {
261  result=(double) ((size_t) pixel & (size_t) (value+0.5));
262  break;
263  }
265  {
266  result=(double) (QuantumRange*(0.5*cos((double) (2.0*MagickPI*
267  QuantumScale*pixel*value))+0.5));
268  break;
269  }
271  {
272  result=pixel/(value == 0.0 ? 1.0 : value);
273  break;
274  }
276  {
277  result=(double) (QuantumRange*exp((double) (value*QuantumScale*pixel)));
278  break;
279  }
281  {
282  result=(double) GenerateDifferentialNoise(random_info,pixel,
283  GaussianNoise,value);
284  break;
285  }
287  {
288  result=(double) GenerateDifferentialNoise(random_info,pixel,ImpulseNoise,
289  value);
290  break;
291  }
293  {
294  result=(double) GenerateDifferentialNoise(random_info,pixel,
295  LaplacianNoise,value);
296  break;
297  }
299  {
300  result=(double) ((size_t) pixel << (size_t) (value+0.5));
301  break;
302  }
303  case LogEvaluateOperator:
304  {
305  if ((QuantumScale*pixel) >= MagickEpsilon)
306  result=(double) (QuantumRange*log((double) (QuantumScale*value*pixel+
307  1.0))/log((double) (value+1.0)));
308  break;
309  }
310  case MaxEvaluateOperator:
311  {
312  result=(double) EvaluateMax((double) pixel,value);
313  break;
314  }
316  {
317  result=(double) (pixel+value);
318  break;
319  }
321  {
322  result=(double) (pixel+value);
323  break;
324  }
325  case MinEvaluateOperator:
326  {
327  result=(double) MagickMin((double) pixel,value);
328  break;
329  }
331  {
332  result=(double) GenerateDifferentialNoise(random_info,pixel,
334  break;
335  }
337  {
338  result=(double) (value*pixel);
339  break;
340  }
341  case OrEvaluateOperator:
342  {
343  result=(double) ((size_t) pixel | (size_t) (value+0.5));
344  break;
345  }
347  {
348  result=(double) GenerateDifferentialNoise(random_info,pixel,PoissonNoise,
349  value);
350  break;
351  }
352  case PowEvaluateOperator:
353  {
354  result=(double) (QuantumRange*pow((double) (QuantumScale*pixel),(double)
355  value));
356  break;
357  }
359  {
360  result=(double) ((size_t) pixel >> (size_t) (value+0.5));
361  break;
362  }
364  {
365  result=(double) (pixel*pixel+value);
366  break;
367  }
368  case SetEvaluateOperator:
369  {
370  result=value;
371  break;
372  }
374  {
375  result=(double) (QuantumRange*(0.5*sin((double) (2.0*MagickPI*
376  QuantumScale*pixel*value))+0.5));
377  break;
378  }
380  {
381  result=(double) (pixel-value);
382  break;
383  }
384  case SumEvaluateOperator:
385  {
386  result=(double) (pixel+value);
387  break;
388  }
390  {
391  result=(double) (((double) pixel <= value) ? 0 : QuantumRange);
392  break;
393  }
395  {
396  result=(double) (((double) pixel <= value) ? 0 : pixel);
397  break;
398  }
400  {
401  result=(double) (((double) pixel > value) ? QuantumRange : pixel);
402  break;
403  }
405  {
406  result=(double) GenerateDifferentialNoise(random_info,pixel,UniformNoise,
407  value);
408  break;
409  }
410  case XorEvaluateOperator:
411  {
412  result=(double) ((size_t) pixel ^ (size_t) (value+0.5));
413  break;
414  }
415  }
416  return(result);
417 }
418 
419 static Image *AcquireImageCanvas(const Image *images,ExceptionInfo *exception)
420 {
421  const Image
422  *p,
423  *q;
424 
425  size_t
426  columns,
427  rows;
428 
429  q=images;
430  columns=images->columns;
431  rows=images->rows;
432  for (p=images; p != (Image *) NULL; p=p->next)
433  {
434  if (p->number_channels > q->number_channels)
435  q=p;
436  if (p->columns > columns)
437  columns=p->columns;
438  if (p->rows > rows)
439  rows=p->rows;
440  }
441  return(CloneImage(q,columns,rows,MagickTrue,exception));
442 }
443 
445  const MagickEvaluateOperator op,ExceptionInfo *exception)
446 {
447 #define EvaluateImageTag "Evaluate/Image"
448 
449  CacheView
450  *evaluate_view;
451 
452  Image
453  *image;
454 
456  status;
457 
459  progress;
460 
462  **magick_restrict evaluate_pixels;
463 
464  RandomInfo
466 
467  size_t
468  number_images;
469 
470  ssize_t
471  y;
472 
473 #if defined(MAGICKCORE_OPENMP_SUPPORT)
474  unsigned long
475  key;
476 #endif
477 
478  assert(images != (Image *) NULL);
479  assert(images->signature == MagickCoreSignature);
480  if (images->debug != MagickFalse)
481  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
482  assert(exception != (ExceptionInfo *) NULL);
483  assert(exception->signature == MagickCoreSignature);
484  image=AcquireImageCanvas(images,exception);
485  if (image == (Image *) NULL)
486  return((Image *) NULL);
487  if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
488  {
489  image=DestroyImage(image);
490  return((Image *) NULL);
491  }
492  number_images=GetImageListLength(images);
493  evaluate_pixels=AcquirePixelThreadSet(images);
494  if (evaluate_pixels == (PixelChannels **) NULL)
495  {
496  image=DestroyImage(image);
497  (void) ThrowMagickException(exception,GetMagickModule(),
498  ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
499  return((Image *) NULL);
500  }
501  /*
502  Evaluate image pixels.
503  */
504  status=MagickTrue;
505  progress=0;
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) 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) 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) 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) 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) memset(channel_moments,0,(MaxPixelChannels+1)*
1398  sizeof(*channel_moments));
1399  (void) memset(centroid,0,sizeof(centroid));
1400  (void) memset(M00,0,sizeof(M00));
1401  (void) memset(M01,0,sizeof(M01));
1402  (void) memset(M02,0,sizeof(M02));
1403  (void) memset(M03,0,sizeof(M03));
1404  (void) memset(M10,0,sizeof(M10));
1405  (void) memset(M11,0,sizeof(M11));
1406  (void) memset(M12,0,sizeof(M12));
1407  (void) memset(M20,0,sizeof(M20));
1408  (void) memset(M21,0,sizeof(M21));
1409  (void) memset(M22,0,sizeof(M22));
1410  (void) memset(M30,0,sizeof(M30));
1411  image_view=AcquireVirtualCacheView(image,exception);
1412  for (y=0; y < (ssize_t) image->rows; y++)
1413  {
1414  register const Quantum
1415  *magick_restrict p;
1416 
1417  register ssize_t
1418  x;
1419 
1420  /*
1421  Compute center of mass (centroid).
1422  */
1423  p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1424  if (p == (const Quantum *) NULL)
1425  break;
1426  for (x=0; x < (ssize_t) image->columns; x++)
1427  {
1428  register ssize_t
1429  i;
1430 
1431  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) 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) memset(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) memset(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  channel_statistics[channel].entropy+=-count*MagickLog10(count)*
2114  PerceptibleReciprocal(MagickLog10(number_bins));
2115  channel_statistics[CompositePixelChannel].entropy+=-count*
2116  MagickLog10(count)*PerceptibleReciprocal(MagickLog10(number_bins))/
2117  GetPixelChannels(image);
2118  }
2119  }
2120  histogram=(double *) RelinquishMagickMemory(histogram);
2121  for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
2122  {
2123  /*
2124  Compute kurtosis & skewness statistics.
2125  */
2126  standard_deviation=PerceptibleReciprocal(
2127  channel_statistics[i].standard_deviation);
2128  channel_statistics[i].skewness=(channel_statistics[i].sum_cubed-3.0*
2129  channel_statistics[i].mean*channel_statistics[i].sum_squared+2.0*
2130  channel_statistics[i].mean*channel_statistics[i].mean*
2131  channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2132  standard_deviation);
2133  channel_statistics[i].kurtosis=(channel_statistics[i].sum_fourth_power-4.0*
2134  channel_statistics[i].mean*channel_statistics[i].sum_cubed+6.0*
2135  channel_statistics[i].mean*channel_statistics[i].mean*
2136  channel_statistics[i].sum_squared-3.0*channel_statistics[i].mean*
2137  channel_statistics[i].mean*1.0*channel_statistics[i].mean*
2138  channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2139  standard_deviation*standard_deviation)-3.0;
2140  }
2141  channel_statistics[CompositePixelChannel].mean=0.0;
2142  channel_statistics[CompositePixelChannel].standard_deviation=0.0;
2143  channel_statistics[CompositePixelChannel].entropy=0.0;
2144  for (i=0; i < (ssize_t) MaxPixelChannels; i++)
2145  {
2146  channel_statistics[CompositePixelChannel].mean+=
2147  channel_statistics[i].mean;
2148  channel_statistics[CompositePixelChannel].standard_deviation+=
2149  channel_statistics[i].standard_deviation;
2150  channel_statistics[CompositePixelChannel].entropy+=
2151  channel_statistics[i].entropy;
2152  }
2153  channel_statistics[CompositePixelChannel].mean/=(double)
2154  GetImageChannels(image);
2155  channel_statistics[CompositePixelChannel].standard_deviation/=(double)
2156  GetImageChannels(image);
2157  channel_statistics[CompositePixelChannel].entropy/=(double)
2158  GetImageChannels(image);
2159  if (y < (ssize_t) image->rows)
2160  channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
2161  channel_statistics);
2162  return(channel_statistics);
2163 }
2164 
2165 /*
2166 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2167 % %
2168 % %
2169 % %
2170 % P o l y n o m i a l I m a g e %
2171 % %
2172 % %
2173 % %
2174 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2175 %
2176 % PolynomialImage() returns a new image where each pixel is the sum of the
2177 % pixels in the image sequence after applying its corresponding terms
2178 % (coefficient and degree pairs).
2179 %
2180 % The format of the PolynomialImage method is:
2181 %
2182 % Image *PolynomialImage(const Image *images,const size_t number_terms,
2183 % const double *terms,ExceptionInfo *exception)
2184 %
2185 % A description of each parameter follows:
2186 %
2187 % o images: the image sequence.
2188 %
2189 % o number_terms: the number of terms in the list. The actual list length
2190 % is 2 x number_terms + 1 (the constant).
2191 %
2192 % o terms: the list of polynomial coefficients and degree pairs and a
2193 % constant.
2194 %
2195 % o exception: return any errors or warnings in this structure.
2196 %
2197 */
2198 
2200  const size_t number_terms,const double *terms,ExceptionInfo *exception)
2201 {
2202 #define PolynomialImageTag "Polynomial/Image"
2203 
2204  CacheView
2205  *polynomial_view;
2206 
2207  Image
2208  *image;
2209 
2211  status;
2212 
2214  progress;
2215 
2217  **magick_restrict polynomial_pixels;
2218 
2219  size_t
2220  number_images;
2221 
2222  ssize_t
2223  y;
2224 
2225  assert(images != (Image *) NULL);
2226  assert(images->signature == MagickCoreSignature);
2227  if (images->debug != MagickFalse)
2228  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
2229  assert(exception != (ExceptionInfo *) NULL);
2230  assert(exception->signature == MagickCoreSignature);
2231  image=AcquireImageCanvas(images,exception);
2232  if (image == (Image *) NULL)
2233  return((Image *) NULL);
2234  if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
2235  {
2236  image=DestroyImage(image);
2237  return((Image *) NULL);
2238  }
2239  number_images=GetImageListLength(images);
2240  polynomial_pixels=AcquirePixelThreadSet(images);
2241  if (polynomial_pixels == (PixelChannels **) NULL)
2242  {
2243  image=DestroyImage(image);
2244  (void) ThrowMagickException(exception,GetMagickModule(),
2245  ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
2246  return((Image *) NULL);
2247  }
2248  /*
2249  Polynomial image pixels.
2250  */
2251  status=MagickTrue;
2252  progress=0;
2253  polynomial_view=AcquireAuthenticCacheView(image,exception);
2254 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2255  #pragma omp parallel for schedule(static) shared(progress,status) \
2256  magick_number_threads(image,image,image->rows,1)
2257 #endif
2258  for (y=0; y < (ssize_t) image->rows; y++)
2259  {
2260  CacheView
2261  *image_view;
2262 
2263  const Image
2264  *next;
2265 
2266  const int
2267  id = GetOpenMPThreadId();
2268 
2269  register ssize_t
2270  i,
2271  x;
2272 
2273  register PixelChannels
2274  *polynomial_pixel;
2275 
2276  register Quantum
2277  *magick_restrict q;
2278 
2279  ssize_t
2280  j;
2281 
2282  if (status == MagickFalse)
2283  continue;
2284  q=QueueCacheViewAuthenticPixels(polynomial_view,0,y,image->columns,1,
2285  exception);
2286  if (q == (Quantum *) NULL)
2287  {
2288  status=MagickFalse;
2289  continue;
2290  }
2291  polynomial_pixel=polynomial_pixels[id];
2292  for (j=0; j < (ssize_t) image->columns; j++)
2293  for (i=0; i < MaxPixelChannels; i++)
2294  polynomial_pixel[j].channel[i]=0.0;
2295  next=images;
2296  for (j=0; j < (ssize_t) number_images; j++)
2297  {
2298  register const Quantum
2299  *p;
2300 
2301  if (j >= (ssize_t) number_terms)
2302  continue;
2303  image_view=AcquireVirtualCacheView(next,exception);
2304  p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2305  if (p == (const Quantum *) NULL)
2306  {
2307  image_view=DestroyCacheView(image_view);
2308  break;
2309  }
2310  for (x=0; x < (ssize_t) image->columns; x++)
2311  {
2312  register ssize_t
2313  i;
2314 
2315  if (GetPixelWriteMask(next,p) <= (QuantumRange/2))
2316  {
2317  p+=GetPixelChannels(next);
2318  continue;
2319  }
2320  for (i=0; i < (ssize_t) GetPixelChannels(next); i++)
2321  {
2323  coefficient,
2324  degree;
2325 
2327  PixelTrait traits = GetPixelChannelTraits(next,channel);
2328  PixelTrait polynomial_traits=GetPixelChannelTraits(image,channel);
2329  if ((traits == UndefinedPixelTrait) ||
2330  (polynomial_traits == UndefinedPixelTrait))
2331  continue;
2332  if ((traits & UpdatePixelTrait) == 0)
2333  continue;
2334  coefficient=(MagickRealType) terms[2*j];
2335  degree=(MagickRealType) terms[(j << 1)+1];
2336  polynomial_pixel[x].channel[i]+=coefficient*
2337  pow(QuantumScale*GetPixelChannel(image,channel,p),degree);
2338  }
2339  p+=GetPixelChannels(next);
2340  }
2341  image_view=DestroyCacheView(image_view);
2342  next=GetNextImageInList(next);
2343  }
2344  for (x=0; x < (ssize_t) image->columns; x++)
2345  {
2346  register ssize_t
2347  i;
2348 
2349  if (GetPixelWriteMask(image,q) <= (QuantumRange/2))
2350  {
2351  q+=GetPixelChannels(image);
2352  continue;
2353  }
2354  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2355  {
2357  PixelTrait traits = GetPixelChannelTraits(image,channel);
2358  if (traits == UndefinedPixelTrait)
2359  continue;
2360  if ((traits & UpdatePixelTrait) == 0)
2361  continue;
2362  q[i]=ClampToQuantum(QuantumRange*polynomial_pixel[x].channel[i]);
2363  }
2364  q+=GetPixelChannels(image);
2365  }
2366  if (SyncCacheViewAuthenticPixels(polynomial_view,exception) == MagickFalse)
2367  status=MagickFalse;
2368  if (images->progress_monitor != (MagickProgressMonitor) NULL)
2369  {
2371  proceed;
2372 
2373 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2374  #pragma omp critical (MagickCore_PolynomialImages)
2375 #endif
2376  proceed=SetImageProgress(images,PolynomialImageTag,progress++,
2377  image->rows);
2378  if (proceed == MagickFalse)
2379  status=MagickFalse;
2380  }
2381  }
2382  polynomial_view=DestroyCacheView(polynomial_view);
2383  polynomial_pixels=DestroyPixelThreadSet(polynomial_pixels);
2384  if (status == MagickFalse)
2385  image=DestroyImage(image);
2386  return(image);
2387 }
2388 
2389 /*
2390 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2391 % %
2392 % %
2393 % %
2394 % S t a t i s t i c I m a g e %
2395 % %
2396 % %
2397 % %
2398 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2399 %
2400 % StatisticImage() makes each pixel the min / max / median / mode / etc. of
2401 % the neighborhood of the specified width and height.
2402 %
2403 % The format of the StatisticImage method is:
2404 %
2405 % Image *StatisticImage(const Image *image,const StatisticType type,
2406 % const size_t width,const size_t height,ExceptionInfo *exception)
2407 %
2408 % A description of each parameter follows:
2409 %
2410 % o image: the image.
2411 %
2412 % o type: the statistic type (median, mode, etc.).
2413 %
2414 % o width: the width of the pixel neighborhood.
2415 %
2416 % o height: the height of the pixel neighborhood.
2417 %
2418 % o exception: return any errors or warnings in this structure.
2419 %
2420 */
2421 
2422 typedef struct _SkipNode
2423 {
2424  size_t
2425  next[9],
2426  count,
2427  signature;
2428 } SkipNode;
2429 
2430 typedef struct _SkipList
2431 {
2432  ssize_t
2434 
2435  SkipNode
2437 } SkipList;
2438 
2439 typedef struct _PixelList
2440 {
2441  size_t
2442  length,
2443  seed;
2444 
2445  SkipList
2447 
2448  size_t
2450 } PixelList;
2451 
2453 {
2454  if (pixel_list == (PixelList *) NULL)
2455  return((PixelList *) NULL);
2456  if (pixel_list->skip_list.nodes != (SkipNode *) NULL)
2458  pixel_list->skip_list.nodes);
2459  pixel_list=(PixelList *) RelinquishMagickMemory(pixel_list);
2460  return(pixel_list);
2461 }
2462 
2464 {
2465  register ssize_t
2466  i;
2467 
2468  assert(pixel_list != (PixelList **) NULL);
2469  for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
2470  if (pixel_list[i] != (PixelList *) NULL)
2471  pixel_list[i]=DestroyPixelList(pixel_list[i]);
2472  pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list);
2473  return(pixel_list);
2474 }
2475 
2476 static PixelList *AcquirePixelList(const size_t width,const size_t height)
2477 {
2478  PixelList
2479  *pixel_list;
2480 
2481  pixel_list=(PixelList *) AcquireMagickMemory(sizeof(*pixel_list));
2482  if (pixel_list == (PixelList *) NULL)
2483  return(pixel_list);
2484  (void) memset((void *) pixel_list,0,sizeof(*pixel_list));
2485  pixel_list->length=width*height;
2486  pixel_list->skip_list.nodes=(SkipNode *) AcquireAlignedMemory(65537UL,
2487  sizeof(*pixel_list->skip_list.nodes));
2488  if (pixel_list->skip_list.nodes == (SkipNode *) NULL)
2489  return(DestroyPixelList(pixel_list));
2490  (void) memset(pixel_list->skip_list.nodes,0,65537UL*
2491  sizeof(*pixel_list->skip_list.nodes));
2492  pixel_list->signature=MagickCoreSignature;
2493  return(pixel_list);
2494 }
2495 
2496 static PixelList **AcquirePixelListThreadSet(const size_t width,
2497  const size_t height)
2498 {
2499  PixelList
2500  **pixel_list;
2501 
2502  register ssize_t
2503  i;
2504 
2505  size_t
2506  number_threads;
2507 
2508  number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
2509  pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
2510  sizeof(*pixel_list));
2511  if (pixel_list == (PixelList **) NULL)
2512  return((PixelList **) NULL);
2513  (void) memset(pixel_list,0,number_threads*sizeof(*pixel_list));
2514  for (i=0; i < (ssize_t) number_threads; i++)
2515  {
2516  pixel_list[i]=AcquirePixelList(width,height);
2517  if (pixel_list[i] == (PixelList *) NULL)
2518  return(DestroyPixelListThreadSet(pixel_list));
2519  }
2520  return(pixel_list);
2521 }
2522 
2523 static void AddNodePixelList(PixelList *pixel_list,const size_t color)
2524 {
2525  register SkipList
2526  *p;
2527 
2528  register ssize_t
2529  level;
2530 
2531  size_t
2532  search,
2533  update[9];
2534 
2535  /*
2536  Initialize the node.
2537  */
2538  p=(&pixel_list->skip_list);
2539  p->nodes[color].signature=pixel_list->signature;
2540  p->nodes[color].count=1;
2541  /*
2542  Determine where it belongs in the list.
2543  */
2544  search=65536UL;
2545  for (level=p->level; level >= 0; level--)
2546  {
2547  while (p->nodes[search].next[level] < color)
2548  search=p->nodes[search].next[level];
2549  update[level]=search;
2550  }
2551  /*
2552  Generate a pseudo-random level for this node.
2553  */
2554  for (level=0; ; level++)
2555  {
2556  pixel_list->seed=(pixel_list->seed*42893621L)+1L;
2557  if ((pixel_list->seed & 0x300) != 0x300)
2558  break;
2559  }
2560  if (level > 8)
2561  level=8;
2562  if (level > (p->level+2))
2563  level=p->level+2;
2564  /*
2565  If we're raising the list's level, link back to the root node.
2566  */
2567  while (level > p->level)
2568  {
2569  p->level++;
2570  update[p->level]=65536UL;
2571  }
2572  /*
2573  Link the node into the skip-list.
2574  */
2575  do
2576  {
2577  p->nodes[color].next[level]=p->nodes[update[level]].next[level];
2578  p->nodes[update[level]].next[level]=color;
2579  } while (level-- > 0);
2580 }
2581 
2582 static inline void GetMaximumPixelList(PixelList *pixel_list,Quantum *pixel)
2583 {
2584  register SkipList
2585  *p;
2586 
2587  size_t
2588  color,
2589  maximum;
2590 
2591  ssize_t
2592  count;
2593 
2594  /*
2595  Find the maximum value for each of the color.
2596  */
2597  p=(&pixel_list->skip_list);
2598  color=65536L;
2599  count=0;
2600  maximum=p->nodes[color].next[0];
2601  do
2602  {
2603  color=p->nodes[color].next[0];
2604  if (color > maximum)
2605  maximum=color;
2606  count+=p->nodes[color].count;
2607  } while (count < (ssize_t) pixel_list->length);
2608  *pixel=ScaleShortToQuantum((unsigned short) maximum);
2609 }
2610 
2611 static inline void GetMeanPixelList(PixelList *pixel_list,Quantum *pixel)
2612 {
2613  double
2614  sum;
2615 
2616  register SkipList
2617  *p;
2618 
2619  size_t
2620  color;
2621 
2622  ssize_t
2623  count;
2624 
2625  /*
2626  Find the mean value for each of the color.
2627  */
2628  p=(&pixel_list->skip_list);
2629  color=65536L;
2630  count=0;
2631  sum=0.0;
2632  do
2633  {
2634  color=p->nodes[color].next[0];
2635  sum+=(double) p->nodes[color].count*color;
2636  count+=p->nodes[color].count;
2637  } while (count < (ssize_t) pixel_list->length);
2638  sum/=pixel_list->length;
2639  *pixel=ScaleShortToQuantum((unsigned short) sum);
2640 }
2641 
2642 static inline void GetMedianPixelList(PixelList *pixel_list,Quantum *pixel)
2643 {
2644  register SkipList
2645  *p;
2646 
2647  size_t
2648  color;
2649 
2650  ssize_t
2651  count;
2652 
2653  /*
2654  Find the median value for each of the color.
2655  */
2656  p=(&pixel_list->skip_list);
2657  color=65536L;
2658  count=0;
2659  do
2660  {
2661  color=p->nodes[color].next[0];
2662  count+=p->nodes[color].count;
2663  } while (count <= (ssize_t) (pixel_list->length >> 1));
2664  *pixel=ScaleShortToQuantum((unsigned short) color);
2665 }
2666 
2667 static inline void GetMinimumPixelList(PixelList *pixel_list,Quantum *pixel)
2668 {
2669  register SkipList
2670  *p;
2671 
2672  size_t
2673  color,
2674  minimum;
2675 
2676  ssize_t
2677  count;
2678 
2679  /*
2680  Find the minimum value for each of the color.
2681  */
2682  p=(&pixel_list->skip_list);
2683  count=0;
2684  color=65536UL;
2685  minimum=p->nodes[color].next[0];
2686  do
2687  {
2688  color=p->nodes[color].next[0];
2689  if (color < minimum)
2690  minimum=color;
2691  count+=p->nodes[color].count;
2692  } while (count < (ssize_t) pixel_list->length);
2693  *pixel=ScaleShortToQuantum((unsigned short) minimum);
2694 }
2695 
2696 static inline void GetModePixelList(PixelList *pixel_list,Quantum *pixel)
2697 {
2698  register SkipList
2699  *p;
2700 
2701  size_t
2702  color,
2703  max_count,
2704  mode;
2705 
2706  ssize_t
2707  count;
2708 
2709  /*
2710  Make each pixel the 'predominant color' of the specified neighborhood.
2711  */
2712  p=(&pixel_list->skip_list);
2713  color=65536L;
2714  mode=color;
2715  max_count=p->nodes[mode].count;
2716  count=0;
2717  do
2718  {
2719  color=p->nodes[color].next[0];
2720  if (p->nodes[color].count > max_count)
2721  {
2722  mode=color;
2723  max_count=p->nodes[mode].count;
2724  }
2725  count+=p->nodes[color].count;
2726  } while (count < (ssize_t) pixel_list->length);
2727  *pixel=ScaleShortToQuantum((unsigned short) mode);
2728 }
2729 
2730 static inline void GetNonpeakPixelList(PixelList *pixel_list,Quantum *pixel)
2731 {
2732  register SkipList
2733  *p;
2734 
2735  size_t
2736  color,
2737  next,
2738  previous;
2739 
2740  ssize_t
2741  count;
2742 
2743  /*
2744  Finds the non peak value for each of the colors.
2745  */
2746  p=(&pixel_list->skip_list);
2747  color=65536L;
2748  next=p->nodes[color].next[0];
2749  count=0;
2750  do
2751  {
2752  previous=color;
2753  color=next;
2754  next=p->nodes[color].next[0];
2755  count+=p->nodes[color].count;
2756  } while (count <= (ssize_t) (pixel_list->length >> 1));
2757  if ((previous == 65536UL) && (next != 65536UL))
2758  color=next;
2759  else
2760  if ((previous != 65536UL) && (next == 65536UL))
2761  color=previous;
2762  *pixel=ScaleShortToQuantum((unsigned short) color);
2763 }
2764 
2765 static inline void GetRootMeanSquarePixelList(PixelList *pixel_list,
2766  Quantum *pixel)
2767 {
2768  double
2769  sum;
2770 
2771  register SkipList
2772  *p;
2773 
2774  size_t
2775  color;
2776 
2777  ssize_t
2778  count;
2779 
2780  /*
2781  Find the root mean square value for each of the color.
2782  */
2783  p=(&pixel_list->skip_list);
2784  color=65536L;
2785  count=0;
2786  sum=0.0;
2787  do
2788  {
2789  color=p->nodes[color].next[0];
2790  sum+=(double) (p->nodes[color].count*color*color);
2791  count+=p->nodes[color].count;
2792  } while (count < (ssize_t) pixel_list->length);
2793  sum/=pixel_list->length;
2794  *pixel=ScaleShortToQuantum((unsigned short) sqrt(sum));
2795 }
2796 
2797 static inline void GetStandardDeviationPixelList(PixelList *pixel_list,
2798  Quantum *pixel)
2799 {
2800  double
2801  sum,
2802  sum_squared;
2803 
2804  register SkipList
2805  *p;
2806 
2807  size_t
2808  color;
2809 
2810  ssize_t
2811  count;
2812 
2813  /*
2814  Find the standard-deviation value for each of the color.
2815  */
2816  p=(&pixel_list->skip_list);
2817  color=65536L;
2818  count=0;
2819  sum=0.0;
2820  sum_squared=0.0;
2821  do
2822  {
2823  register ssize_t
2824  i;
2825 
2826  color=p->nodes[color].next[0];
2827  sum+=(double) p->nodes[color].count*color;
2828  for (i=0; i < (ssize_t) p->nodes[color].count; i++)
2829  sum_squared+=((double) color)*((double) color);
2830  count+=p->nodes[color].count;
2831  } while (count < (ssize_t) pixel_list->length);
2832  sum/=pixel_list->length;
2833  sum_squared/=pixel_list->length;
2834  *pixel=ScaleShortToQuantum((unsigned short) sqrt(sum_squared-(sum*sum)));
2835 }
2836 
2837 static inline void InsertPixelList(const Quantum pixel,PixelList *pixel_list)
2838 {
2839  size_t
2840  signature;
2841 
2842  unsigned short
2843  index;
2844 
2845  index=ScaleQuantumToShort(pixel);
2846  signature=pixel_list->skip_list.nodes[index].signature;
2847  if (signature == pixel_list->signature)
2848  {
2849  pixel_list->skip_list.nodes[index].count++;
2850  return;
2851  }
2852  AddNodePixelList(pixel_list,index);
2853 }
2854 
2855 static void ResetPixelList(PixelList *pixel_list)
2856 {
2857  int
2858  level;
2859 
2860  register SkipNode
2861  *root;
2862 
2863  register SkipList
2864  *p;
2865 
2866  /*
2867  Reset the skip-list.
2868  */
2869  p=(&pixel_list->skip_list);
2870  root=p->nodes+65536UL;
2871  p->level=0;
2872  for (level=0; level < 9; level++)
2873  root->next[level]=65536UL;
2874  pixel_list->seed=pixel_list->signature++;
2875 }
2876 
2878  const size_t width,const size_t height,ExceptionInfo *exception)
2879 {
2880 #define StatisticImageTag "Statistic/Image"
2881 
2882  CacheView
2883  *image_view,
2884  *statistic_view;
2885 
2886  Image
2887  *statistic_image;
2888 
2890  status;
2891 
2893  progress;
2894 
2895  PixelList
2896  **magick_restrict pixel_list;
2897 
2898  ssize_t
2899  center,
2900  y;
2901 
2902  /*
2903  Initialize statistics image attributes.
2904  */
2905  assert(image != (Image *) NULL);
2906  assert(image->signature == MagickCoreSignature);
2907  if (image->debug != MagickFalse)
2908  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2909  assert(exception != (ExceptionInfo *) NULL);
2910  assert(exception->signature == MagickCoreSignature);
2911  statistic_image=CloneImage(image,image->columns,image->rows,MagickTrue,
2912  exception);
2913  if (statistic_image == (Image *) NULL)
2914  return((Image *) NULL);
2915  status=SetImageStorageClass(statistic_image,DirectClass,exception);
2916  if (status == MagickFalse)
2917  {
2918  statistic_image=DestroyImage(statistic_image);
2919  return((Image *) NULL);
2920  }
2921  pixel_list=AcquirePixelListThreadSet(MagickMax(width,1),MagickMax(height,1));
2922  if (pixel_list == (PixelList **) NULL)
2923  {
2924  statistic_image=DestroyImage(statistic_image);
2925  ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2926  }
2927  /*
2928  Make each pixel the min / max / median / mode / etc. of the neighborhood.
2929  */
2930  center=(ssize_t) GetPixelChannels(image)*(image->columns+MagickMax(width,1))*
2931  (MagickMax(height,1)/2L)+GetPixelChannels(image)*(MagickMax(width,1)/2L);
2932  status=MagickTrue;
2933  progress=0;
2934  image_view=AcquireVirtualCacheView(image,exception);
2935  statistic_view=AcquireAuthenticCacheView(statistic_image,exception);
2936 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2937  #pragma omp parallel for schedule(static) shared(progress,status) \
2938  magick_number_threads(image,statistic_image,statistic_image->rows,1)
2939 #endif
2940  for (y=0; y < (ssize_t) statistic_image->rows; y++)
2941  {
2942  const int
2943  id = GetOpenMPThreadId();
2944 
2945  register const Quantum
2946  *magick_restrict p;
2947 
2948  register Quantum
2949  *magick_restrict q;
2950 
2951  register ssize_t
2952  x;
2953 
2954  if (status == MagickFalse)
2955  continue;
2956  p=GetCacheViewVirtualPixels(image_view,-((ssize_t) MagickMax(width,1)/2L),y-
2957  (ssize_t) (MagickMax(height,1)/2L),image->columns+MagickMax(width,1),
2958  MagickMax(height,1),exception);
2959  q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns, 1,exception);
2960  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2961  {
2962  status=MagickFalse;
2963  continue;
2964  }
2965  for (x=0; x < (ssize_t) statistic_image->columns; x++)
2966  {
2967  register ssize_t
2968  i;
2969 
2970  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2971  {
2972  Quantum
2973  pixel;
2974 
2975  register const Quantum
2976  *magick_restrict pixels;
2977 
2978  register ssize_t
2979  u;
2980 
2981  ssize_t
2982  v;
2983 
2985  PixelTrait traits = GetPixelChannelTraits(image,channel);
2986  PixelTrait statistic_traits=GetPixelChannelTraits(statistic_image,
2987  channel);
2988  if ((traits == UndefinedPixelTrait) ||
2989  (statistic_traits == UndefinedPixelTrait))
2990  continue;
2991  if (((statistic_traits & CopyPixelTrait) != 0) ||
2992  (GetPixelWriteMask(image,p) <= (QuantumRange/2)))
2993  {
2994  SetPixelChannel(statistic_image,channel,p[center+i],q);
2995  continue;
2996  }
2997  if ((statistic_traits & UpdatePixelTrait) == 0)
2998  continue;
2999  pixels=p;
3000  ResetPixelList(pixel_list[id]);
3001  for (v=0; v < (ssize_t) MagickMax(height,1); v++)
3002  {
3003  for (u=0; u < (ssize_t) MagickMax(width,1); u++)
3004  {
3005  InsertPixelList(pixels[i],pixel_list[id]);
3006  pixels+=GetPixelChannels(image);
3007  }
3008  pixels+=GetPixelChannels(image)*image->columns;
3009  }
3010  switch (type)
3011  {
3012  case GradientStatistic:
3013  {
3014  double
3015  maximum,
3016  minimum;
3017 
3018  GetMinimumPixelList(pixel_list[id],&pixel);
3019  minimum=(double) pixel;
3020  GetMaximumPixelList(pixel_list[id],&pixel);
3021  maximum=(double) pixel;
3022  pixel=ClampToQuantum(MagickAbsoluteValue(maximum-minimum));
3023  break;
3024  }
3025  case MaximumStatistic:
3026  {
3027  GetMaximumPixelList(pixel_list[id],&pixel);
3028  break;
3029  }
3030  case MeanStatistic:
3031  {
3032  GetMeanPixelList(pixel_list[id],&pixel);
3033  break;
3034  }
3035  case MedianStatistic:
3036  default:
3037  {
3038  GetMedianPixelList(pixel_list[id],&pixel);
3039  break;
3040  }
3041  case MinimumStatistic:
3042  {
3043  GetMinimumPixelList(pixel_list[id],&pixel);
3044  break;
3045  }
3046  case ModeStatistic:
3047  {
3048  GetModePixelList(pixel_list[id],&pixel);
3049  break;
3050  }
3051  case NonpeakStatistic:
3052  {
3053  GetNonpeakPixelList(pixel_list[id],&pixel);
3054  break;
3055  }
3057  {
3058  GetRootMeanSquarePixelList(pixel_list[id],&pixel);
3059  break;
3060  }
3062  {
3063  GetStandardDeviationPixelList(pixel_list[id],&pixel);
3064  break;
3065  }
3066  }
3067  SetPixelChannel(statistic_image,channel,pixel,q);
3068  }
3069  p+=GetPixelChannels(image);
3070  q+=GetPixelChannels(statistic_image);
3071  }
3072  if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
3073  status=MagickFalse;
3074  if (image->progress_monitor != (MagickProgressMonitor) NULL)
3075  {
3077  proceed;
3078 
3079 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3080  #pragma omp critical (MagickCore_StatisticImage)
3081 #endif
3082  proceed=SetImageProgress(image,StatisticImageTag,progress++,
3083  image->rows);
3084  if (proceed == MagickFalse)
3085  status=MagickFalse;
3086  }
3087  }
3088  statistic_view=DestroyCacheView(statistic_view);
3089  image_view=DestroyCacheView(image_view);
3090  pixel_list=DestroyPixelListThreadSet(pixel_list);
3091  if (status == MagickFalse)
3092  statistic_image=DestroyImage(statistic_image);
3093  return(statistic_image);
3094 }
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:1325
double ellipse_angle
Definition: statistic.h:60
#define Log10Epsilon
MagickExport ssize_t ParseCommandOption(const CommandOption option, const MagickBooleanType list, const char *options)
Definition: option.c:2956
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:2452
static void GetMaximumPixelList(PixelList *pixel_list, Quantum *pixel)
Definition: statistic.c:2582
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:2425
#define MaximumNumberOfImageMoments
Definition: statistic.h:25
static void GetMeanPixelList(PixelList *pixel_list, Quantum *pixel)
Definition: statistic.c:2611
#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:2442
#define EvaluateImageTag
double x
Definition: geometry.h:123
#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:2436
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:2730
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:2425
size_t signature
Definition: statistic.c:2449
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:3153
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:533
double y
Definition: geometry.h:123
static PixelList ** DestroyPixelListThreadSet(PixelList **pixel_list)
Definition: statistic.c:2463
static int GetOpenMPThreadId(void)
MagickExport void * RelinquishAlignedMemory(void *memory)
Definition: memory.c:990
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:745
static void GetMinimumPixelList(PixelList *pixel_list, Quantum *pixel)
Definition: statistic.c:2667
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:2837
struct _Image * next
Definition: image.h:348
static PixelList * AcquirePixelList(const size_t width, const size_t height)
Definition: statistic.c:2476
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:2602
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:242
MagickExport MagickBooleanType GetImageMean(const Image *image, double *mean, double *standard_deviation, ExceptionInfo *exception)
Definition: statistic.c:1288
#define MaxMap
Definition: magick-type.h:75
#define MagickMax(x, y)
Definition: image-private.h:26
static size_t GetPixelChannels(const Image *magick_restrict image)
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:2442
ssize_t level
Definition: statistic.c:2433
#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:2642
PointInfo centroid
Definition: statistic.h:56
MagickExport char * StringToken(const char *delimiters, char **string)
Definition: string.c:2307
static void GetModePixelList(PixelList *pixel_list, Quantum *pixel)
Definition: statistic.c:2696
unsigned short Quantum
Definition: magick-type.h:82
MagickExport Image * GetNextImageInList(const Image *images)
Definition: list.c:762
MagickExport char * DestroyString(char *string)
Definition: string.c:816
MagickExport void * AcquireMagickMemory(const size_t size)
Definition: memory.c:462
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:2446
static void GetStandardDeviationPixelList(PixelList *pixel_list, Quantum *pixel)
Definition: statistic.c:2797
ColorspaceType colorspace[MaximumNumberOfPerceptualColorspaces+1]
Definition: statistic.h:75
#define StatisticImageTag
#define FunctionImageTag
#define MagickMin(x, y)
Definition: image-private.h:27
ColorspaceType
Definition: colorspace.h:25
static RandomInfo * random_info
Definition: resource.c:111
size_t next[9]
Definition: statistic.c:2425
PointInfo ellipse_axis
Definition: statistic.h:56
MagickExport void * RelinquishMagickMemory(void *memory)
Definition: memory.c:1027
#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:2523
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:2855
MagickExport Image * PolynomialImage(const Image *images, const size_t number_terms, const double *terms, ExceptionInfo *exception)
Definition: statistic.c:2199
static void GetRootMeanSquarePixelList(PixelList *pixel_list, Quantum *pixel)
Definition: statistic.c:2765
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:687
MagickSizeType QuantumAny
Definition: magick-type.h:142
static PixelList ** AcquirePixelListThreadSet(const size_t width, const size_t height)
Definition: statistic.c:2496
MagickExport Image * DestroyImage(Image *image)
Definition: image.c:1183
MagickExport Image * CloneImage(const Image *image, const size_t columns, const size_t rows, const MagickBooleanType detach, ExceptionInfo *exception)
Definition: image.c:800
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:2877
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