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