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 MagickBooleanType
1847 status;
1848
1849 RangeInfo
1850 range_info = { -MagickMaximumValue, MagickMaximumValue };
1851
1852 ssize_t
1853 y;
1854
1855 assert(image != (Image *) NULL);
1856 assert(image->signature == MagickCoreSignature);
1857 if (IsEventLogging() != MagickFalse)
1858 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1859 status=MagickTrue;
1860 image_view=AcquireVirtualCacheView(image,exception);
1861#if defined(MAGICKCORE_OPENMP_SUPPORT)
1862 #pragma omp parallel for schedule(static) shared(range_info,status) \
1863 magick_number_threads(image,image,image->rows,1)
1864#endif
1865 for (y=0; y < (ssize_t) image->rows; y++)
1866 {
1867 const Quantum
1868 *magick_restrict p;
1869
1870 RangeInfo
1871 channel_range;
1872
1873 ssize_t
1874 x;
1875
1876 if (status == MagickFalse)
1877 continue;
1878 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1879 if (p == (const Quantum *) NULL)
1880 {
1881 status=MagickFalse;
1882 continue;
1883 }
1884 channel_range=range_info;
1885 for (x=0; x < (ssize_t) image->columns; x++)
1886 {
1887 ssize_t
1888 i;
1889
1890 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1891 {
1892 PixelChannel channel = GetPixelChannelChannel(image,i);
1893 PixelTrait traits = GetPixelChannelTraits(image,channel);
1894 if ((traits & UpdatePixelTrait) == 0)
1895 continue;
1896 if ((double) p[i] > channel_range.maxima)
1897 channel_range.maxima=(double) p[i];
1898 if ((double) p[i] < channel_range.minima)
1899 channel_range.minima=(double) p[i];
1900 }
1901 p+=(ptrdiff_t) GetPixelChannels(image);
1902 }
1903#if defined(MAGICKCORE_OPENMP_SUPPORT)
1904 #pragma omp critical (MagickCore_GetImageRange)
1905#endif
1906 {
1907 if (channel_range.maxima > range_info.maxima)
1908 range_info.maxima=channel_range.maxima;
1909 if (channel_range.minima < range_info.minima)
1910 range_info.minima=channel_range.minima;
1911 }
1912 }
1913 image_view=DestroyCacheView(image_view);
1914 *maxima=range_info.maxima;
1915 *minima=range_info.minima;
1916 return(status);
1917}
1918␌
1919/*
1920%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1921% %
1922% %
1923% %
1924% G e t I m a g e S t a t i s t i c s %
1925% %
1926% %
1927% %
1928%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1929%
1930% GetImageStatistics() returns statistics for each channel in the image. The
1931% statistics include the channel depth, its minima, maxima, mean, standard
1932% deviation, kurtosis and skewness. You can access the red channel mean, for
1933% example, like this:
1934%
1935% channel_statistics=GetImageStatistics(image,exception);
1936% red_mean=channel_statistics[RedPixelChannel].mean;
1937%
1938% Use MagickRelinquishMemory() to free the statistics buffer.
1939%
1940% The format of the GetImageStatistics method is:
1941%
1942% ChannelStatistics *GetImageStatistics(const Image *image,
1943% ExceptionInfo *exception)
1944%
1945% A description of each parameter follows:
1946%
1947% o image: the image.
1948%
1949% o exception: return any errors or warnings in this structure.
1950%
1951*/
1952
1953static ssize_t GetMedianPixel(Quantum *pixels,const size_t n)
1954{
1955#define SwapPixels(alpha,beta) \
1956{ \
1957 Quantum gamma=(alpha); \
1958 (alpha)=(beta);(beta)=gamma; \
1959}
1960
1961 ssize_t
1962 low = 0,
1963 high = (ssize_t) n-1,
1964 median = (low+high)/2;
1965
1966 for ( ; ; )
1967 {
1968 ssize_t
1969 l = low+1,
1970 h = high,
1971 mid = (low+high)/2;
1972
1973 if (high <= low)
1974 return(median);
1975 if (high == (low+1))
1976 {
1977 if (pixels[low] > pixels[high])
1978 SwapPixels(pixels[low],pixels[high]);
1979 return(median);
1980 }
1981 if (pixels[mid] > pixels[high])
1982 SwapPixels(pixels[mid],pixels[high]);
1983 if (pixels[low] > pixels[high])
1984 SwapPixels(pixels[low], pixels[high]);
1985 if (pixels[mid] > pixels[low])
1986 SwapPixels(pixels[mid],pixels[low]);
1987 SwapPixels(pixels[mid],pixels[low+1]);
1988 for ( ; ; )
1989 {
1990 do l++; while (pixels[low] > pixels[l]);
1991 do h--; while (pixels[h] > pixels[low]);
1992 if (h < l)
1993 break;
1994 SwapPixels(pixels[l],pixels[h]);
1995 }
1996 SwapPixels(pixels[low],pixels[h]);
1997 if (h <= median)
1998 low=l;
1999 if (h >= median)
2000 high=h-1;
2001 }
2002}
2003
2004static inline long double MagickSafeReciprocalLD(const long double x)
2005{
2006 long double
2007 sign;
2008
2009 /*
2010 Return 1/x where x is perceptible (not unlimited or infinitesimal).
2011 */
2012 sign=x < 0.0 ? -1.0 : 1.0;
2013 if ((sign*x) >= MagickEpsilon)
2014 return(1.0/x);
2015 return(sign/MagickEpsilon);
2016}
2017
2018MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
2019 ExceptionInfo *exception)
2020{
2021 ChannelStatistics
2022 *channel_statistics;
2023
2024 double
2025 channels,
2026 *histogram;
2027
2028 long double
2029 area;
2030
2031 MagickStatusType
2032 status;
2033
2034 MemoryInfo
2035 *median_info;
2036
2037 Quantum
2038 *median;
2039
2040 QuantumAny
2041 range;
2042
2043 size_t
2044 depth;
2045
2046 ssize_t
2047 i,
2048 y;
2049
2050 assert(image != (Image *) NULL);
2051 assert(image->signature == MagickCoreSignature);
2052 if (IsEventLogging() != MagickFalse)
2053 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2054 histogram=(double *) AcquireQuantumMemory(MaxMap+1UL,
2055 MagickMax(GetPixelChannels(image),1)*sizeof(*histogram));
2056 channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(
2057 MaxPixelChannels+1,sizeof(*channel_statistics));
2058 if ((channel_statistics == (ChannelStatistics *) NULL) ||
2059 (histogram == (double *) NULL))
2060 {
2061 if (histogram != (double *) NULL)
2062 histogram=(double *) RelinquishMagickMemory(histogram);
2063 if (channel_statistics != (ChannelStatistics *) NULL)
2064 channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
2065 channel_statistics);
2066 (void) ThrowMagickException(exception,GetMagickModule(),
2067 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
2068 return(channel_statistics);
2069 }
2070 (void) memset(channel_statistics,0,(MaxPixelChannels+1)*
2071 sizeof(*channel_statistics));
2072 for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
2073 {
2074 ChannelStatistics *cs = channel_statistics+i;
2075 cs->area=0.0;
2076 cs->depth=1;
2077 cs->maxima=(-MagickMaximumValue);
2078 cs->minima=MagickMaximumValue;
2079 cs->sum=0.0;
2080 cs->sumLD=0.0;
2081 cs->mean=0.0;
2082 cs->standard_deviation=0.0;
2083 cs->variance=0.0;
2084 cs->skewness=0.0;
2085 cs->kurtosis=0.0;
2086 cs->entropy=0.0;
2087 }
2088 (void) memset(histogram,0,(MaxMap+1)*GetPixelChannels(image)*
2089 sizeof(*histogram));
2090 for (y=0; y < (ssize_t) image->rows; y++)
2091 {
2092 const Quantum
2093 *magick_restrict p;
2094
2095 ssize_t
2096 x;
2097
2098 /*
2099 Compute pixel statistics.
2100 */
2101 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
2102 if (p == (const Quantum *) NULL)
2103 break;
2104 for (x=0; x < (ssize_t) image->columns; x++)
2105 {
2106 if (GetPixelReadMask(image,p) <= (QuantumRange/2))
2107 {
2108 p+=(ptrdiff_t) GetPixelChannels(image);
2109 continue;
2110 }
2111 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2112 {
2113 ChannelStatistics
2114 *cs;
2115
2116 PixelChannel channel = GetPixelChannelChannel(image,i);
2117 PixelTrait traits = GetPixelChannelTraits(image,channel);
2118 if ((traits & UpdatePixelTrait) == 0)
2119 continue;
2120 cs=channel_statistics+channel;
2121 if (cs->depth != MAGICKCORE_QUANTUM_DEPTH)
2122 {
2123 depth=cs->depth;
2124 range=GetQuantumRange(depth);
2125 status=p[i] != ScaleAnyToQuantum(ScaleQuantumToAny(p[i],range),
2126 range) ? MagickTrue : MagickFalse;
2127 if (status != MagickFalse)
2128 {
2129 cs->depth++;
2130 if (cs->depth > channel_statistics[CompositePixelChannel].depth)
2131 channel_statistics[CompositePixelChannel].depth=cs->depth;
2132 i--;
2133 continue;
2134 }
2135 }
2136 cs->area++;
2137 if ((double) p[i] < cs->minima)
2138 cs->minima=(double) p[i];
2139 if ((double) p[i] > cs->maxima)
2140 cs->maxima=(double) p[i];
2141 histogram[(ssize_t) GetPixelChannels(image)*ScaleQuantumToMap(
2142 ClampToQuantum((double) p[i]))+i]++;
2143 cs->sumLD+=(long double) p[i];
2144 /*
2145 sum_squared, sum_cubed and sum_fourth_power are not used in
2146 MagickCore or MagickWand, but are made available in
2147 Magick++/lib/Statistic.cpp, so we need to calculate these.
2148 */
2149 cs->sum_squared+=(double) p[i]*(double) p[i];
2150 cs->sum_cubed+=(double) p[i]*(double) p[i]*(double) p[i];
2151 cs->sum_fourth_power+=(double) p[i]*(double) p[i]*(double) p[i]*
2152 (double) p[i];
2153 {
2154 /* Calculate running totals for Welford's method.
2155 */
2156 double
2157 n = cs->area,
2158 n1 = cs->area-1;
2159
2160 long double
2161 delta,
2162 delta_n,
2163 delta_n2,
2164 term1;
2165
2166 delta=(double) p[i]-cs->M1;
2167 delta_n=delta/n;
2168 delta_n2=delta_n*delta_n;
2169 term1=delta*delta_n*n1;
2170 cs->M4+=term1*delta_n2*(n*n-3.0*n+3.0)+6.0*delta_n2*cs->M2-4.0*
2171 delta_n*cs->M3;
2172 cs->M3+=term1*delta_n*(n-2.0)-3.0*delta_n*cs->M2;
2173 cs->M2+=term1;
2174 cs->M1+=delta_n;
2175 }
2176 }
2177 p+=(ptrdiff_t) GetPixelChannels(image);
2178 }
2179 }
2180 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2181 {
2182 ChannelStatistics
2183 *cs;
2184
2185 double
2186 adj_area = 1.0;
2187
2188 PixelChannel channel = GetPixelChannelChannel(image,i);
2189 PixelTrait traits = GetPixelChannelTraits(image,channel);
2190 if ((traits & UpdatePixelTrait) == 0)
2191 continue;
2192 cs=channel_statistics+(ssize_t) channel;
2193 cs->mean=0.0;
2194 if (cs->area > 0)
2195 {
2196 cs->mean=(double) (cs->sumLD/(long double) cs->area);
2197 if (cs->area > 1.0)
2198 adj_area=cs->area/(cs->area-1.0);
2199 }
2200 cs->sum=(double) cs->sum;
2201 if (cs->M2 == 0.0)
2202 {
2203 cs->standard_deviation=0.0;
2204 cs->variance=0.0;
2205 cs->skewness=0.0;
2206 cs->kurtosis=0.0;
2207 }
2208 else
2209 {
2210 if (cs->area > 1.0)
2211 cs->standard_deviation=(double) sqrtl(cs->M2/((long double)
2212 cs->area-1.0));
2213 else
2214 cs->standard_deviation=(double) sqrtl(cs->M2/((long double)
2215 cs->area));
2216 cs->variance=cs->standard_deviation*cs->standard_deviation;
2217 cs->skewness=(double) (sqrtl(cs->area)*cs->M3/powl(cs->M2*adj_area,
2218 1.5));
2219 cs->kurtosis=(double) (cs->area*cs->M4/(cs->M2*cs->M2*adj_area*
2220 adj_area)-3.0);
2221 }
2222 }
2223 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2224 {
2225 ChannelStatistics
2226 *cs;
2227
2228 double
2229 number_bins;
2230
2231 ssize_t
2232 j;
2233
2234 PixelChannel channel = GetPixelChannelChannel(image,i);
2235 PixelTrait traits = GetPixelChannelTraits(image,channel);
2236 if ((traits & UpdatePixelTrait) == 0)
2237 continue;
2238 cs=channel_statistics+(ssize_t) channel;
2239 if (cs->area > 0.0)
2240 {
2241 cs->sum/=cs->area;
2242 cs->sum_squared/=cs->area;
2243 cs->sum_cubed/=cs->area;
2244 cs->sum_fourth_power/=cs->area;
2245 }
2246 /*
2247 Compute pixel entropy.
2248 */
2249 number_bins=0.0;
2250 for (j=0; j <= (ssize_t) MaxMap; j++)
2251 if (histogram[(ssize_t) GetPixelChannels(image)*j+i] > 0.0)
2252 number_bins++;
2253 area=MagickSafeReciprocalLD(channel_statistics[channel].area);
2254 number_bins=(double) MagickSafeReciprocalLD((long double) log2(number_bins));
2255 for (j=0; j <= (ssize_t) MaxMap; j++)
2256 {
2257 double
2258 entropy,
2259 count;
2260
2261 count=(double) (area*histogram[(ssize_t) GetPixelChannels(image)*j+i]);
2262 entropy=-count*log2(count)*number_bins;
2263 if (IsNaN(entropy) != 0)
2264 continue;
2265 channel_statistics[channel].entropy+=(double) entropy;
2266 channel_statistics[CompositePixelChannel].entropy+=((double) entropy/
2267 GetPixelChannels(image));
2268 }
2269 }
2270 histogram=(double *) RelinquishMagickMemory(histogram);
2271 median_info=AcquireVirtualMemory(image->columns,image->rows*sizeof(*median));
2272 if (median_info == (MemoryInfo *) NULL)
2273 (void) ThrowMagickException(exception,GetMagickModule(),
2274 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
2275 else
2276 {
2277 median=(Quantum *) GetVirtualMemoryBlob(median_info);
2278 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2279 {
2280 size_t
2281 n = 0;
2282
2283 /*
2284 Compute median statistics for each channel.
2285 */
2286 PixelChannel channel = GetPixelChannelChannel(image,i);
2287 PixelTrait traits = GetPixelChannelTraits(image,channel);
2288 if ((traits & UpdatePixelTrait) == 0)
2289 continue;
2290 for (y=0; y < (ssize_t) image->rows; y++)
2291 {
2292 const Quantum
2293 *magick_restrict p;
2294
2295 ssize_t
2296 x;
2297
2298 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
2299 if (p == (const Quantum *) NULL)
2300 break;
2301 for (x=0; x < (ssize_t) image->columns; x++)
2302 {
2303 if (GetPixelReadMask(image,p) <= (QuantumRange/2))
2304 {
2305 p+=(ptrdiff_t) GetPixelChannels(image);
2306 continue;
2307 }
2308 median[n++]=p[i];
2309 p+=(ptrdiff_t) GetPixelChannels(image);
2310 }
2311 }
2312 channel_statistics[channel].median=(double) median[
2313 GetMedianPixel(median,n)];
2314 }
2315 median_info=RelinquishVirtualMemory(median_info);
2316 }
2317 {
2318 ChannelStatistics *cs_comp = channel_statistics+CompositePixelChannel;
2319 cs_comp->sum=0.0;
2320 cs_comp->sum_squared=0.0;
2321 cs_comp->sum_cubed=0.0;
2322 cs_comp->sum_fourth_power=0.0;
2323 cs_comp->maxima=(-MagickMaximumValue);
2324 cs_comp->minima=MagickMaximumValue;
2325 cs_comp->area=0.0;
2326 cs_comp->mean=0.0;
2327 cs_comp->median=0.0;
2328 cs_comp->variance=0.0;
2329 cs_comp->standard_deviation=0.0;
2330 cs_comp->entropy=0.0;
2331 cs_comp->skewness=0.0;
2332 cs_comp->kurtosis=0.0;
2333 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2334 {
2335 ChannelStatistics
2336 *cs;
2337
2338 PixelChannel channel = GetPixelChannelChannel(image,i);
2339 PixelTrait traits = GetPixelChannelTraits(image,channel);
2340 if ((traits & UpdatePixelTrait) == 0)
2341 continue;
2342 cs=channel_statistics+channel;
2343 if (cs_comp->maxima < cs->maxima)
2344 cs_comp->maxima=cs->maxima;
2345 if (cs_comp->minima > cs->minima)
2346 cs_comp->minima=cs->minima;
2347 cs_comp->sum+=cs->sum;
2348 cs_comp->sum_squared+=cs->sum_squared;
2349 cs_comp->sum_cubed+=cs->sum_cubed;
2350 cs_comp->sum_fourth_power+=cs->sum_fourth_power;
2351 cs_comp->median+=cs->median;
2352 cs_comp->area+=cs->area;
2353 cs_comp->mean+=cs->mean;
2354 cs_comp->variance+=cs->variance;
2355 cs_comp->standard_deviation+=cs->standard_deviation;
2356 cs_comp->skewness+=cs->skewness;
2357 cs_comp->kurtosis+=cs->kurtosis;
2358 cs_comp->entropy+=cs->entropy;
2359 }
2360 channels=(double) GetImageChannels(image);
2361 cs_comp->sum/=channels;
2362 cs_comp->sum_squared/=channels;
2363 cs_comp->sum_cubed/=channels;
2364 cs_comp->sum_fourth_power/=channels;
2365 cs_comp->median/=channels;
2366 cs_comp->area/=channels;
2367 cs_comp->mean/=channels;
2368 cs_comp->variance/=channels;
2369 cs_comp->standard_deviation/=channels;
2370 cs_comp->skewness/=channels;
2371 cs_comp->kurtosis/=channels;
2372 cs_comp->entropy/=channels;
2373 }
2374 if (y < (ssize_t) image->rows)
2375 channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
2376 channel_statistics);
2377 return(channel_statistics);
2378}
2379␌
2380/*
2381%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2382% %
2383% %
2384% %
2385% P o l y n o m i a l I m a g e %
2386% %
2387% %
2388% %
2389%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2390%
2391% PolynomialImage() returns a new image where each pixel is the sum of the
2392% pixels in the image sequence after applying its corresponding terms
2393% (coefficient and degree pairs).
2394%
2395% The format of the PolynomialImage method is:
2396%
2397% Image *PolynomialImage(const Image *images,const size_t number_terms,
2398% const double *terms,ExceptionInfo *exception)
2399%
2400% A description of each parameter follows:
2401%
2402% o images: the image sequence.
2403%
2404% o number_terms: the number of terms in the list. The actual list length
2405% is 2 x number_terms + 1 (the constant).
2406%
2407% o terms: the list of polynomial coefficients and degree pairs and a
2408% constant.
2409%
2410% o exception: return any errors or warnings in this structure.
2411%
2412*/
2413MagickExport Image *PolynomialImage(const Image *images,
2414 const size_t number_terms,const double *terms,ExceptionInfo *exception)
2415{
2416#define PolynomialImageTag "Polynomial/Image"
2417
2418 CacheView
2419 *polynomial_view;
2420
2421 Image
2422 *image;
2423
2424 MagickBooleanType
2425 status;
2426
2427 MagickOffsetType
2428 progress;
2429
2430 PixelChannels
2431 **magick_restrict polynomial_pixels;
2432
2433 size_t
2434 number_images;
2435
2436 ssize_t
2437 y;
2438
2439 assert(images != (Image *) NULL);
2440 assert(images->signature == MagickCoreSignature);
2441 if (IsEventLogging() != MagickFalse)
2442 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
2443 assert(exception != (ExceptionInfo *) NULL);
2444 assert(exception->signature == MagickCoreSignature);
2445 image=AcquireImageCanvas(images,exception);
2446 if (image == (Image *) NULL)
2447 return((Image *) NULL);
2448 if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
2449 {
2450 image=DestroyImage(image);
2451 return((Image *) NULL);
2452 }
2453 number_images=GetImageListLength(images);
2454 polynomial_pixels=AcquirePixelTLS(images);
2455 if (polynomial_pixels == (PixelChannels **) NULL)
2456 {
2457 image=DestroyImage(image);
2458 (void) ThrowMagickException(exception,GetMagickModule(),
2459 ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
2460 return((Image *) NULL);
2461 }
2462 /*
2463 Polynomial image pixels.
2464 */
2465 status=MagickTrue;
2466 progress=0;
2467 polynomial_view=AcquireAuthenticCacheView(image,exception);
2468#if defined(MAGICKCORE_OPENMP_SUPPORT)
2469 #pragma omp parallel for schedule(static) shared(progress,status) \
2470 magick_number_threads(image,image,image->rows,1)
2471#endif
2472 for (y=0; y < (ssize_t) image->rows; y++)
2473 {
2474 CacheView
2475 *image_view;
2476
2477 const Image
2478 *next;
2479
2480 const int
2481 id = GetOpenMPThreadId();
2482
2483 PixelChannels
2484 *polynomial_pixel;
2485
2486 Quantum
2487 *magick_restrict q;
2488
2489 ssize_t
2490 i,
2491 j,
2492 x;
2493
2494 if (status == MagickFalse)
2495 continue;
2496 q=QueueCacheViewAuthenticPixels(polynomial_view,0,y,image->columns,1,
2497 exception);
2498 if (q == (Quantum *) NULL)
2499 {
2500 status=MagickFalse;
2501 continue;
2502 }
2503 polynomial_pixel=polynomial_pixels[id];
2504 for (j=0; j < (ssize_t) image->columns; j++)
2505 for (i=0; i < MaxPixelChannels; i++)
2506 polynomial_pixel[j].channel[i]=0.0;
2507 next=images;
2508 for (j=0; j < (ssize_t) number_images; j++)
2509 {
2510 const Quantum
2511 *p;
2512
2513 if (j >= (ssize_t) number_terms)
2514 continue;
2515 image_view=AcquireVirtualCacheView(next,exception);
2516 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2517 if (p == (const Quantum *) NULL)
2518 {
2519 image_view=DestroyCacheView(image_view);
2520 break;
2521 }
2522 for (x=0; x < (ssize_t) image->columns; x++)
2523 {
2524 for (i=0; i < (ssize_t) GetPixelChannels(next); i++)
2525 {
2526 MagickRealType
2527 coefficient,
2528 degree;
2529
2530 PixelChannel channel = GetPixelChannelChannel(image,i);
2531 PixelTrait traits = GetPixelChannelTraits(next,channel);
2532 PixelTrait polynomial_traits = GetPixelChannelTraits(image,channel);
2533 if (((traits & UpdatePixelTrait) == 0) ||
2534 ((polynomial_traits & UpdatePixelTrait) == 0))
2535 continue;
2536 coefficient=(MagickRealType) terms[2*j];
2537 degree=(MagickRealType) terms[(j << 1)+1];
2538 polynomial_pixel[x].channel[i]+=coefficient*
2539 pow(QuantumScale*(double) GetPixelChannel(image,channel,p),degree);
2540 }
2541 p+=(ptrdiff_t) GetPixelChannels(next);
2542 }
2543 image_view=DestroyCacheView(image_view);
2544 next=GetNextImageInList(next);
2545 }
2546 for (x=0; x < (ssize_t) image->columns; x++)
2547 {
2548 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2549 {
2550 PixelChannel channel = GetPixelChannelChannel(image,i);
2551 PixelTrait traits = GetPixelChannelTraits(image,channel);
2552 if ((traits & UpdatePixelTrait) == 0)
2553 continue;
2554 q[i]=ClampToQuantum((double) QuantumRange*
2555 polynomial_pixel[x].channel[i]);
2556 }
2557 q+=(ptrdiff_t) GetPixelChannels(image);
2558 }
2559 if (SyncCacheViewAuthenticPixels(polynomial_view,exception) == MagickFalse)
2560 status=MagickFalse;
2561 if (images->progress_monitor != (MagickProgressMonitor) NULL)
2562 {
2563 MagickBooleanType
2564 proceed;
2565
2566#if defined(MAGICKCORE_OPENMP_SUPPORT)
2567 #pragma omp atomic
2568#endif
2569 progress++;
2570 proceed=SetImageProgress(images,PolynomialImageTag,progress,
2571 image->rows);
2572 if (proceed == MagickFalse)
2573 status=MagickFalse;
2574 }
2575 }
2576 polynomial_view=DestroyCacheView(polynomial_view);
2577 polynomial_pixels=DestroyPixelTLS(images,polynomial_pixels);
2578 if (status == MagickFalse)
2579 image=DestroyImage(image);
2580 return(image);
2581}
2582␌
2583/*
2584%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2585% %
2586% %
2587% %
2588% S t a t i s t i c I m a g e %
2589% %
2590% %
2591% %
2592%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2593%
2594% StatisticImage() makes each pixel the min / max / median / mode / etc. of
2595% the neighborhood of the specified width and height.
2596%
2597% The format of the StatisticImage method is:
2598%
2599% Image *StatisticImage(const Image *image,const StatisticType type,
2600% const size_t width,const size_t height,ExceptionInfo *exception)
2601%
2602% A description of each parameter follows:
2603%
2604% o image: the image.
2605%
2606% o type: the statistic type (median, mode, etc.).
2607%
2608% o width: the width of the pixel neighborhood.
2609%
2610% o height: the height of the pixel neighborhood.
2611%
2612% o exception: return any errors or warnings in this structure.
2613%
2614*/
2615
2616typedef struct _SkipNode
2617{
2618 size_t
2619 next[9],
2620 count,
2621 signature;
2622} SkipNode;
2623
2624typedef struct _SkipList
2625{
2626 ssize_t
2627 level;
2628
2629 SkipNode
2630 *nodes;
2631} SkipList;
2632
2633typedef struct _PixelList
2634{
2635 size_t
2636 length,
2637 seed;
2638
2639 SkipList
2640 skip_list;
2641
2642 size_t
2643 signature;
2644} PixelList;
2645
2646static PixelList *DestroyPixelList(PixelList *pixel_list)
2647{
2648 if (pixel_list == (PixelList *) NULL)
2649 return((PixelList *) NULL);
2650 if (pixel_list->skip_list.nodes != (SkipNode *) NULL)
2651 pixel_list->skip_list.nodes=(SkipNode *) RelinquishAlignedMemory(
2652 pixel_list->skip_list.nodes);
2653 pixel_list=(PixelList *) RelinquishMagickMemory(pixel_list);
2654 return(pixel_list);
2655}
2656
2657static PixelList **DestroyPixelListTLS(PixelList **pixel_list)
2658{
2659 ssize_t
2660 i;
2661
2662 assert(pixel_list != (PixelList **) NULL);
2663 for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
2664 if (pixel_list[i] != (PixelList *) NULL)
2665 pixel_list[i]=DestroyPixelList(pixel_list[i]);
2666 pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list);
2667 return(pixel_list);
2668}
2669
2670static PixelList *AcquirePixelList(const size_t width,const size_t height)
2671{
2672 PixelList
2673 *pixel_list;
2674
2675 pixel_list=(PixelList *) AcquireMagickMemory(sizeof(*pixel_list));
2676 if (pixel_list == (PixelList *) NULL)
2677 return(pixel_list);
2678 (void) memset((void *) pixel_list,0,sizeof(*pixel_list));
2679 pixel_list->length=width*height;
2680 pixel_list->skip_list.nodes=(SkipNode *) AcquireAlignedMemory(65537UL,
2681 sizeof(*pixel_list->skip_list.nodes));
2682 if (pixel_list->skip_list.nodes == (SkipNode *) NULL)
2683 return(DestroyPixelList(pixel_list));
2684 (void) memset(pixel_list->skip_list.nodes,0,65537UL*
2685 sizeof(*pixel_list->skip_list.nodes));
2686 pixel_list->signature=MagickCoreSignature;
2687 return(pixel_list);
2688}
2689
2690static PixelList **AcquirePixelListTLS(const size_t width,
2691 const size_t height)
2692{
2693 PixelList
2694 **pixel_list;
2695
2696 ssize_t
2697 i;
2698
2699 size_t
2700 number_threads;
2701
2702 number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
2703 pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
2704 sizeof(*pixel_list));
2705 if (pixel_list == (PixelList **) NULL)
2706 return((PixelList **) NULL);
2707 (void) memset(pixel_list,0,number_threads*sizeof(*pixel_list));
2708 for (i=0; i < (ssize_t) number_threads; i++)
2709 {
2710 pixel_list[i]=AcquirePixelList(width,height);
2711 if (pixel_list[i] == (PixelList *) NULL)
2712 return(DestroyPixelListTLS(pixel_list));
2713 }
2714 return(pixel_list);
2715}
2716
2717static void AddNodePixelList(PixelList *pixel_list,const size_t color)
2718{
2719 SkipList
2720 *p;
2721
2722 ssize_t
2723 level;
2724
2725 size_t
2726 search,
2727 update[9];
2728
2729 /*
2730 Initialize the node.
2731 */
2732 p=(&pixel_list->skip_list);
2733 p->nodes[color].signature=pixel_list->signature;
2734 p->nodes[color].count=1;
2735 /*
2736 Determine where it belongs in the list.
2737 */
2738 search=65536UL;
2739 (void) memset(update,0,sizeof(update));
2740 for (level=p->level; level >= 0; level--)
2741 {
2742 while (p->nodes[search].next[level] < color)
2743 search=p->nodes[search].next[level];
2744 update[level]=search;
2745 }
2746 /*
2747 Generate a pseudo-random level for this node.
2748 */
2749 for (level=0; ; level++)
2750 {
2751 pixel_list->seed=(pixel_list->seed*42893621L)+1L;
2752 if ((pixel_list->seed & 0x300) != 0x300)
2753 break;
2754 }
2755 if (level > 8)
2756 level=8;
2757 if (level > (p->level+2))
2758 level=p->level+2;
2759 /*
2760 If we're raising the list's level, link back to the root node.
2761 */
2762 while (level > p->level)
2763 {
2764 p->level++;
2765 update[p->level]=65536UL;
2766 }
2767 /*
2768 Link the node into the skip-list.
2769 */
2770 do
2771 {
2772 p->nodes[color].next[level]=p->nodes[update[level]].next[level];
2773 p->nodes[update[level]].next[level]=color;
2774 } while (level-- > 0);
2775}
2776
2777static inline void GetMedianPixelList(PixelList *pixel_list,Quantum *pixel)
2778{
2779 SkipList
2780 *p;
2781
2782 size_t
2783 color;
2784
2785 ssize_t
2786 count;
2787
2788 /*
2789 Find the median value for each of the color.
2790 */
2791 p=(&pixel_list->skip_list);
2792 color=65536L;
2793 count=0;
2794 do
2795 {
2796 color=p->nodes[color].next[0];
2797 count+=(ssize_t) p->nodes[color].count;
2798 } while (count <= (ssize_t) (pixel_list->length >> 1));
2799 *pixel=ScaleShortToQuantum((unsigned short) color);
2800}
2801
2802static inline void GetModePixelList(PixelList *pixel_list,Quantum *pixel)
2803{
2804 SkipList
2805 *p;
2806
2807 size_t
2808 color,
2809 max_count,
2810 mode;
2811
2812 ssize_t
2813 count;
2814
2815 /*
2816 Make each pixel the 'predominant color' of the specified neighborhood.
2817 */
2818 p=(&pixel_list->skip_list);
2819 color=65536L;
2820 mode=color;
2821 max_count=p->nodes[mode].count;
2822 count=0;
2823 do
2824 {
2825 color=p->nodes[color].next[0];
2826 if (p->nodes[color].count > max_count)
2827 {
2828 mode=color;
2829 max_count=p->nodes[mode].count;
2830 }
2831 count+=(ssize_t) p->nodes[color].count;
2832 } while (count < (ssize_t) pixel_list->length);
2833 *pixel=ScaleShortToQuantum((unsigned short) mode);
2834}
2835
2836static inline void GetNonpeakPixelList(PixelList *pixel_list,Quantum *pixel)
2837{
2838 SkipList
2839 *p;
2840
2841 size_t
2842 color,
2843 next,
2844 previous;
2845
2846 ssize_t
2847 count;
2848
2849 /*
2850 Finds the non peak value for each of the colors.
2851 */
2852 p=(&pixel_list->skip_list);
2853 color=65536L;
2854 next=p->nodes[color].next[0];
2855 count=0;
2856 do
2857 {
2858 previous=color;
2859 color=next;
2860 next=p->nodes[color].next[0];
2861 count+=(ssize_t) p->nodes[color].count;
2862 } while (count <= (ssize_t) (pixel_list->length >> 1));
2863 if ((previous == 65536UL) && (next != 65536UL))
2864 color=next;
2865 else
2866 if ((previous != 65536UL) && (next == 65536UL))
2867 color=previous;
2868 *pixel=ScaleShortToQuantum((unsigned short) color);
2869}
2870
2871static inline void InsertPixelList(const Quantum pixel,PixelList *pixel_list)
2872{
2873 size_t
2874 signature;
2875
2876 unsigned short
2877 index;
2878
2879 index=ScaleQuantumToShort(pixel);
2880 signature=pixel_list->skip_list.nodes[index].signature;
2881 if (signature == pixel_list->signature)
2882 {
2883 pixel_list->skip_list.nodes[index].count++;
2884 return;
2885 }
2886 AddNodePixelList(pixel_list,index);
2887}
2888
2889static void ResetPixelList(PixelList *pixel_list)
2890{
2891 int
2892 level;
2893
2894 SkipNode
2895 *root;
2896
2897 SkipList
2898 *p;
2899
2900 /*
2901 Reset the skip-list.
2902 */
2903 p=(&pixel_list->skip_list);
2904 root=p->nodes+65536UL;
2905 p->level=0;
2906 for (level=0; level < 9; level++)
2907 root->next[level]=65536UL;
2908 pixel_list->seed=pixel_list->signature++;
2909}
2910
2911MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
2912 const size_t width,const size_t height,ExceptionInfo *exception)
2913{
2914#define StatisticImageTag "Statistic/Image"
2915
2916 CacheView
2917 *image_view,
2918 *statistic_view;
2919
2920 Image
2921 *statistic_image;
2922
2923 MagickBooleanType
2924 status;
2925
2926 MagickOffsetType
2927 progress;
2928
2929 PixelList
2930 **magick_restrict pixel_list;
2931
2932 ssize_t
2933 center,
2934 y;
2935
2936 /*
2937 Initialize statistics image attributes.
2938 */
2939 assert(image != (Image *) NULL);
2940 assert(image->signature == MagickCoreSignature);
2941 if (IsEventLogging() != MagickFalse)
2942 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2943 assert(exception != (ExceptionInfo *) NULL);
2944 assert(exception->signature == MagickCoreSignature);
2945 statistic_image=CloneImage(image,0,0,MagickTrue,
2946 exception);
2947 if (statistic_image == (Image *) NULL)
2948 return((Image *) NULL);
2949 status=SetImageStorageClass(statistic_image,DirectClass,exception);
2950 if (status == MagickFalse)
2951 {
2952 statistic_image=DestroyImage(statistic_image);
2953 return((Image *) NULL);
2954 }
2955 pixel_list=AcquirePixelListTLS(MagickMax(width,1),MagickMax(height,1));
2956 if (pixel_list == (PixelList **) NULL)
2957 {
2958 statistic_image=DestroyImage(statistic_image);
2959 ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2960 }
2961 /*
2962 Make each pixel the min / max / median / mode / etc. of the neighborhood.
2963 */
2964 center=(ssize_t) GetPixelChannels(image)*((ssize_t) image->columns+
2965 MagickMax((ssize_t) width,1L))*(MagickMax((ssize_t) height,1)/2L)+(ssize_t)
2966 GetPixelChannels(image)*(MagickMax((ssize_t) width,1L)/2L);
2967 status=MagickTrue;
2968 progress=0;
2969 image_view=AcquireVirtualCacheView(image,exception);
2970 statistic_view=AcquireAuthenticCacheView(statistic_image,exception);
2971#if defined(MAGICKCORE_OPENMP_SUPPORT)
2972 #pragma omp parallel for schedule(static) shared(progress,status) \
2973 magick_number_threads(image,statistic_image,statistic_image->rows,1)
2974#endif
2975 for (y=0; y < (ssize_t) statistic_image->rows; y++)
2976 {
2977 const int
2978 id = GetOpenMPThreadId();
2979
2980 const Quantum
2981 *magick_restrict p;
2982
2983 Quantum
2984 *magick_restrict q;
2985
2986 ssize_t
2987 x;
2988
2989 if (status == MagickFalse)
2990 continue;
2991 p=GetCacheViewVirtualPixels(image_view,-((ssize_t) MagickMax(width,1)/2L),y-
2992 (ssize_t) (MagickMax(height,1)/2L),image->columns+MagickMax(width,1),
2993 MagickMax(height,1),exception);
2994 q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns, 1,exception);
2995 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2996 {
2997 status=MagickFalse;
2998 continue;
2999 }
3000 for (x=0; x < (ssize_t) statistic_image->columns; x++)
3001 {
3002 ssize_t
3003 i;
3004
3005 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3006 {
3007 double
3008 area,
3009 maximum,
3010 minimum,
3011 sum,
3012 sum_squared;
3013
3014 Quantum
3015 pixel;
3016
3017 const Quantum
3018 *magick_restrict pixels;
3019
3020 ssize_t
3021 u;
3022
3023 ssize_t
3024 v;
3025
3026 PixelChannel channel = GetPixelChannelChannel(image,i);
3027 PixelTrait traits = GetPixelChannelTraits(image,channel);
3028 PixelTrait statistic_traits=GetPixelChannelTraits(statistic_image,
3029 channel);
3030 if (((traits & UpdatePixelTrait) == 0) ||
3031 ((statistic_traits & UpdatePixelTrait) == 0))
3032 continue;
3033 if (((statistic_traits & CopyPixelTrait) != 0) ||
3034 (GetPixelWriteMask(image,p) <= (QuantumRange/2)))
3035 {
3036 SetPixelChannel(statistic_image,channel,p[center+i],q);
3037 continue;
3038 }
3039 pixels=p;
3040 area=0.0;
3041 minimum=pixels[i];
3042 maximum=pixels[i];
3043 sum=0.0;
3044 sum_squared=0.0;
3045 ResetPixelList(pixel_list[id]);
3046 for (v=0; v < (ssize_t) MagickMax(height,1); v++)
3047 {
3048 for (u=0; u < (ssize_t) MagickMax(width,1); u++)
3049 {
3050 if ((type == MedianStatistic) || (type == ModeStatistic) ||
3051 (type == NonpeakStatistic))
3052 {
3053 InsertPixelList(pixels[i],pixel_list[id]);
3054 pixels+=(ptrdiff_t) GetPixelChannels(image);
3055 continue;
3056 }
3057 area++;
3058 if ((double) pixels[i] < minimum)
3059 minimum=(double) pixels[i];
3060 if ((double) pixels[i] > maximum)
3061 maximum=(double) pixels[i];
3062 sum+=(double) pixels[i];
3063 sum_squared+=(double) pixels[i]*(double) pixels[i];
3064 pixels+=(ptrdiff_t) GetPixelChannels(image);
3065 }
3066 pixels+=(ptrdiff_t) GetPixelChannels(image)*image->columns;
3067 }
3068 switch (type)
3069 {
3070 case ContrastStatistic:
3071 {
3072 pixel=ClampToQuantum(MagickAbsoluteValue((maximum-minimum)*
3073 MagickSafeReciprocal(maximum+minimum)));
3074 break;
3075 }
3076 case GradientStatistic:
3077 {
3078 pixel=ClampToQuantum(MagickAbsoluteValue(maximum-minimum));
3079 break;
3080 }
3081 case MaximumStatistic:
3082 {
3083 pixel=ClampToQuantum(maximum);
3084 break;
3085 }
3086 case MeanStatistic:
3087 default:
3088 {
3089 pixel=ClampToQuantum(sum/area);
3090 break;
3091 }
3092 case MedianStatistic:
3093 {
3094 GetMedianPixelList(pixel_list[id],&pixel);
3095 break;
3096 }
3097 case MinimumStatistic:
3098 {
3099 pixel=ClampToQuantum(minimum);
3100 break;
3101 }
3102 case ModeStatistic:
3103 {
3104 GetModePixelList(pixel_list[id],&pixel);
3105 break;
3106 }
3107 case NonpeakStatistic:
3108 {
3109 GetNonpeakPixelList(pixel_list[id],&pixel);
3110 break;
3111 }
3112 case RootMeanSquareStatistic:
3113 {
3114 pixel=ClampToQuantum(sqrt(sum_squared/area));
3115 break;
3116 }
3117 case StandardDeviationStatistic:
3118 {
3119 pixel=ClampToQuantum(sqrt(sum_squared/area-(sum/area*sum/area)));
3120 break;
3121 }
3122 }
3123 SetPixelChannel(statistic_image,channel,pixel,q);
3124 }
3125 p+=(ptrdiff_t) GetPixelChannels(image);
3126 q+=(ptrdiff_t) GetPixelChannels(statistic_image);
3127 }
3128 if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
3129 status=MagickFalse;
3130 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3131 {
3132 MagickBooleanType
3133 proceed;
3134
3135#if defined(MAGICKCORE_OPENMP_SUPPORT)
3136 #pragma omp atomic
3137#endif
3138 progress++;
3139 proceed=SetImageProgress(image,StatisticImageTag,progress,image->rows);
3140 if (proceed == MagickFalse)
3141 status=MagickFalse;
3142 }
3143 }
3144 statistic_view=DestroyCacheView(statistic_view);
3145 image_view=DestroyCacheView(image_view);
3146 pixel_list=DestroyPixelListTLS(pixel_list);
3147 if (status == MagickFalse)
3148 statistic_image=DestroyImage(statistic_image);
3149 return(statistic_image);
3150}
_PixelChannels
Definition statistic.c:136
_PixelList
Definition statistic.c:2634
_SkipList
Definition statistic.c:2625
_SkipNode
Definition statistic.c:2617