MagickCore 7.1.2-27
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/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
136typedef struct _PixelChannels
137{
138 double
139 channel[MaxPixelChannels];
140} PixelChannels;
141
142static 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
161static 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
206static 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)
214extern "C" {
215#endif
216
217static 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
241static 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
453static 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
478MagickExport 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
828MagickExport 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
971static 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
1069MagickExport 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*/
1186MagickExport 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*/
1237MagickExport 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*/
1287MagickExport 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*/
1341MagickExport 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*/
1393MagickExport 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*/
1441MagickExport 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*/
1728MagickExport 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*/
1841MagickExport 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
1961static 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
2012static 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
2026MagickExport 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*/
2421MagickExport 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
2624typedef struct _SkipNode
2625{
2626 size_t
2627 next[9],
2628 count,
2629 signature;
2630} SkipNode;
2631
2632typedef struct _SkipList
2633{
2634 ssize_t
2635 level;
2636
2637 SkipNode
2638 *nodes;
2639} SkipList;
2640
2641typedef 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
2654static 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
2665static 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
2678static 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
2698static 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
2725static 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
2785static 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
2810static 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
2844static 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
2879static 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
2897static 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
2919MagickExport 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