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