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