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