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