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