MagickCore 7.1.1
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)*PerceptibleReciprocal(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 == UndefinedPixelTrait) ||
634 (evaluate_traits == UndefinedPixelTrait) ||
635 ((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 == UndefinedPixelTrait) ||
651 ((traits & UpdatePixelTrait) == 0))
652 continue;
653 q[i]=ClampToQuantum(evaluate_pixel[number_images/2].channel[i]);
654 }
655 q+=(ptrdiff_t) GetPixelChannels(image);
656 }
657 p=(const Quantum **) RelinquishMagickMemory((void *) p);
658 if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
659 status=MagickFalse;
660 if (images->progress_monitor != (MagickProgressMonitor) NULL)
661 {
662 MagickBooleanType
663 proceed;
664
665#if defined(MAGICKCORE_OPENMP_SUPPORT)
666 #pragma omp atomic
667#endif
668 progress++;
669 proceed=SetImageProgress(images,EvaluateImageTag,progress,
670 image->rows);
671 if (proceed == MagickFalse)
672 status=MagickFalse;
673 }
674 }
675 }
676 else
677 {
678#if defined(MAGICKCORE_OPENMP_SUPPORT)
679 key=GetRandomSecretKey(random_info[0]);
680 #pragma omp parallel for schedule(static) shared(progress,status) \
681 magick_number_threads(image,images,image->rows,key == ~0UL)
682#endif
683 for (y=0; y < (ssize_t) image->rows; y++)
684 {
685 const Image
686 *next;
687
688 const int
689 id = GetOpenMPThreadId();
690
691 const Quantum
692 **p;
693
694 PixelChannels
695 *evaluate_pixel;
696
697 Quantum
698 *magick_restrict q;
699
700 ssize_t
701 i,
702 x;
703
704 ssize_t
705 j;
706
707 if (status == MagickFalse)
708 continue;
709 p=(const Quantum **) AcquireQuantumMemory(number_images,sizeof(*p));
710 if (p == (const Quantum **) NULL)
711 {
712 status=MagickFalse;
713 (void) ThrowMagickException(exception,GetMagickModule(),
714 ResourceLimitError,"MemoryAllocationFailed","`%s'",
715 images->filename);
716 continue;
717 }
718 for (j=0; j < (ssize_t) number_images; j++)
719 {
720 p[j]=GetCacheViewVirtualPixels(image_view[j],0,y,image->columns,1,
721 exception);
722 if (p[j] == (const Quantum *) NULL)
723 break;
724 }
725 q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,image->columns,1,
726 exception);
727 if ((j < (ssize_t) number_images) || (q == (Quantum *) NULL))
728 {
729 status=MagickFalse;
730 continue;
731 }
732 evaluate_pixel=evaluate_pixels[id];
733 for (j=0; j < (ssize_t) image->columns; j++)
734 for (i=0; i < MaxPixelChannels; i++)
735 evaluate_pixel[j].channel[i]=0.0;
736 next=images;
737 for (j=0; j < (ssize_t) number_images; j++)
738 {
739 for (x=0; x < (ssize_t) image->columns; x++)
740 {
741 for (i=0; i < (ssize_t) GetPixelChannels(next); i++)
742 {
743 PixelChannel channel = GetPixelChannelChannel(image,i);
744 PixelTrait traits = GetPixelChannelTraits(next,channel);
745 PixelTrait evaluate_traits = GetPixelChannelTraits(image,channel);
746 if ((traits == UndefinedPixelTrait) ||
747 (evaluate_traits == UndefinedPixelTrait))
748 continue;
749 if ((traits & UpdatePixelTrait) == 0)
750 continue;
751 evaluate_pixel[x].channel[i]=ApplyEvaluateOperator(
752 random_info[id],GetPixelChannel(next,channel,p[j]),j == 0 ?
753 AddEvaluateOperator : op,evaluate_pixel[x].channel[i]);
754 }
755 p[j]+=(ptrdiff_t) GetPixelChannels(next);
756 }
757 next=GetNextImageInList(next);
758 }
759 for (x=0; x < (ssize_t) image->columns; x++)
760 {
761 switch (op)
762 {
763 case MeanEvaluateOperator:
764 {
765 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
766 evaluate_pixel[x].channel[i]/=(double) number_images;
767 break;
768 }
769 case MultiplyEvaluateOperator:
770 {
771 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
772 {
773 for (j=0; j < (ssize_t) (number_images-1); j++)
774 evaluate_pixel[x].channel[i]*=QuantumScale;
775 }
776 break;
777 }
778 case RootMeanSquareEvaluateOperator:
779 {
780 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
781 evaluate_pixel[x].channel[i]=sqrt(evaluate_pixel[x].channel[i]/
782 number_images);
783 break;
784 }
785 default:
786 break;
787 }
788 }
789 for (x=0; x < (ssize_t) image->columns; x++)
790 {
791 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
792 {
793 PixelChannel channel = GetPixelChannelChannel(image,i);
794 PixelTrait traits = GetPixelChannelTraits(image,channel);
795 if ((traits == UndefinedPixelTrait) ||
796 ((traits & UpdatePixelTrait) == 0))
797 continue;
798 q[i]=ClampToQuantum(evaluate_pixel[x].channel[i]);
799 }
800 q+=(ptrdiff_t) GetPixelChannels(image);
801 }
802 p=(const Quantum **) RelinquishMagickMemory((void *) p);
803 if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
804 status=MagickFalse;
805 if (images->progress_monitor != (MagickProgressMonitor) NULL)
806 {
807 MagickBooleanType
808 proceed;
809
810#if defined(MAGICKCORE_OPENMP_SUPPORT)
811 #pragma omp atomic
812#endif
813 progress++;
814 proceed=SetImageProgress(images,EvaluateImageTag,progress,
815 image->rows);
816 if (proceed == MagickFalse)
817 status=MagickFalse;
818 }
819 }
820 }
821 for (n=0; n < (ssize_t) number_images; n++)
822 image_view[n]=DestroyCacheView(image_view[n]);
823 image_view=(CacheView **) RelinquishMagickMemory(image_view);
824 evaluate_view=DestroyCacheView(evaluate_view);
825 evaluate_pixels=DestroyPixelTLS(images,evaluate_pixels);
826 random_info=DestroyRandomInfoTLS(random_info);
827 if (status == MagickFalse)
828 image=DestroyImage(image);
829 return(image);
830}
831
832MagickExport MagickBooleanType EvaluateImage(Image *image,
833 const MagickEvaluateOperator op,const double value,ExceptionInfo *exception)
834{
835 CacheView
836 *image_view;
837
838 const char
839 *artifact;
840
841 MagickBooleanType
842 clamp,
843 status;
844
845 MagickOffsetType
846 progress;
847
848 RandomInfo
849 **magick_restrict random_info;
850
851 ssize_t
852 y;
853
854#if defined(MAGICKCORE_OPENMP_SUPPORT)
855 unsigned long
856 key;
857#endif
858
859 assert(image != (Image *) NULL);
860 assert(image->signature == MagickCoreSignature);
861 assert(exception != (ExceptionInfo *) NULL);
862 assert(exception->signature == MagickCoreSignature);
863 if (IsEventLogging() != MagickFalse)
864 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
865 if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
866 return(MagickFalse);
867 status=MagickTrue;
868 progress=0;
869 clamp=MagickFalse;
870 artifact=GetImageArtifact(image,"evaluate:clamp");
871 if (artifact != (const char *) NULL)
872 clamp=IsStringTrue(artifact);
873 random_info=AcquireRandomInfoTLS();
874 image_view=AcquireAuthenticCacheView(image,exception);
875#if defined(MAGICKCORE_OPENMP_SUPPORT)
876 key=GetRandomSecretKey(random_info[0]);
877 #pragma omp parallel for schedule(static) shared(progress,status) \
878 magick_number_threads(image,image,image->rows,key == ~0UL)
879#endif
880 for (y=0; y < (ssize_t) image->rows; y++)
881 {
882 const int
883 id = GetOpenMPThreadId();
884
885 Quantum
886 *magick_restrict q;
887
888 ssize_t
889 x;
890
891 if (status == MagickFalse)
892 continue;
893 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
894 if (q == (Quantum *) NULL)
895 {
896 status=MagickFalse;
897 continue;
898 }
899 for (x=0; x < (ssize_t) image->columns; x++)
900 {
901 double
902 result;
903
904 ssize_t
905 i;
906
907 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
908 {
909 PixelChannel channel = GetPixelChannelChannel(image,i);
910 PixelTrait traits = GetPixelChannelTraits(image,channel);
911 if (traits == UndefinedPixelTrait)
912 continue;
913 if ((traits & CopyPixelTrait) != 0)
914 continue;
915 if ((traits & UpdatePixelTrait) == 0)
916 continue;
917 result=ApplyEvaluateOperator(random_info[id],q[i],op,value);
918 if (op == MeanEvaluateOperator)
919 result/=2.0;
920 q[i]=clamp != MagickFalse ? ClampPixel(result) : ClampToQuantum(result);
921 }
922 q+=(ptrdiff_t) GetPixelChannels(image);
923 }
924 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
925 status=MagickFalse;
926 if (image->progress_monitor != (MagickProgressMonitor) NULL)
927 {
928 MagickBooleanType
929 proceed;
930
931#if defined(MAGICKCORE_OPENMP_SUPPORT)
932 #pragma omp atomic
933#endif
934 progress++;
935 proceed=SetImageProgress(image,EvaluateImageTag,progress,image->rows);
936 if (proceed == MagickFalse)
937 status=MagickFalse;
938 }
939 }
940 image_view=DestroyCacheView(image_view);
941 random_info=DestroyRandomInfoTLS(random_info);
942 return(status);
943}
944␌
945/*
946%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
947% %
948% %
949% %
950% F u n c t i o n I m a g e %
951% %
952% %
953% %
954%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
955%
956% FunctionImage() applies a value to the image with an arithmetic, relational,
957% or logical operator to an image. Use these operations to lighten or darken
958% an image, to increase or decrease contrast in an image, or to produce the
959% "negative" of an image.
960%
961% The format of the FunctionImage method is:
962%
963% MagickBooleanType FunctionImage(Image *image,
964% const MagickFunction function,const ssize_t number_parameters,
965% const double *parameters,ExceptionInfo *exception)
966%
967% A description of each parameter follows:
968%
969% o image: the image.
970%
971% o function: A channel function.
972%
973% o parameters: one or more parameters.
974%
975% o exception: return any errors or warnings in this structure.
976%
977*/
978
979static Quantum ApplyFunction(Quantum pixel,const MagickFunction function,
980 const size_t number_parameters,const double *parameters,
981 ExceptionInfo *exception)
982{
983 double
984 result;
985
986 ssize_t
987 i;
988
989 (void) exception;
990 result=0.0;
991 switch (function)
992 {
993 case PolynomialFunction:
994 {
995 /*
996 Polynomial: polynomial constants, highest to lowest order (e.g. c0*x^3+
997 c1*x^2+c2*x+c3).
998 */
999 result=0.0;
1000 for (i=0; i < (ssize_t) number_parameters; i++)
1001 result=result*QuantumScale*(double) pixel+parameters[i];
1002 result*=(double) QuantumRange;
1003 break;
1004 }
1005 case SinusoidFunction:
1006 {
1007 double
1008 amplitude,
1009 bias,
1010 frequency,
1011 phase;
1012
1013 /*
1014 Sinusoid: frequency, phase, amplitude, bias.
1015 */
1016 frequency=(number_parameters >= 1) ? parameters[0] : 1.0;
1017 phase=(number_parameters >= 2) ? parameters[1] : 0.0;
1018 amplitude=(number_parameters >= 3) ? parameters[2] : 0.5;
1019 bias=(number_parameters >= 4) ? parameters[3] : 0.5;
1020 result=(double) QuantumRange*(amplitude*sin((double) (2.0*
1021 MagickPI*(frequency*QuantumScale*(double) pixel+phase/360.0)))+bias);
1022 break;
1023 }
1024 case ArcsinFunction:
1025 {
1026 double
1027 bias,
1028 center,
1029 range,
1030 width;
1031
1032 /*
1033 Arcsin (pegged at range limits for invalid results): width, center,
1034 range, and bias.
1035 */
1036 width=(number_parameters >= 1) ? parameters[0] : 1.0;
1037 center=(number_parameters >= 2) ? parameters[1] : 0.5;
1038 range=(number_parameters >= 3) ? parameters[2] : 1.0;
1039 bias=(number_parameters >= 4) ? parameters[3] : 0.5;
1040 result=2.0*PerceptibleReciprocal(width)*(QuantumScale*(double) pixel-
1041 center);
1042 if (result <= -1.0)
1043 result=bias-range/2.0;
1044 else
1045 if (result >= 1.0)
1046 result=bias+range/2.0;
1047 else
1048 result=(double) (range/MagickPI*asin((double) result)+bias);
1049 result*=(double) QuantumRange;
1050 break;
1051 }
1052 case ArctanFunction:
1053 {
1054 double
1055 center,
1056 bias,
1057 range,
1058 slope;
1059
1060 /*
1061 Arctan: slope, center, range, and bias.
1062 */
1063 slope=(number_parameters >= 1) ? parameters[0] : 1.0;
1064 center=(number_parameters >= 2) ? parameters[1] : 0.5;
1065 range=(number_parameters >= 3) ? parameters[2] : 1.0;
1066 bias=(number_parameters >= 4) ? parameters[3] : 0.5;
1067 result=MagickPI*slope*(QuantumScale*(double) pixel-center);
1068 result=(double) QuantumRange*(range/MagickPI*atan((double) result)+bias);
1069 break;
1070 }
1071 case UndefinedFunction:
1072 break;
1073 }
1074 return(ClampToQuantum(result));
1075}
1076
1077MagickExport MagickBooleanType FunctionImage(Image *image,
1078 const MagickFunction function,const size_t number_parameters,
1079 const double *parameters,ExceptionInfo *exception)
1080{
1081#define FunctionImageTag "Function/Image "
1082
1083 CacheView
1084 *image_view;
1085
1086 MagickBooleanType
1087 status;
1088
1089 MagickOffsetType
1090 progress;
1091
1092 ssize_t
1093 y;
1094
1095 assert(image != (Image *) NULL);
1096 assert(image->signature == MagickCoreSignature);
1097 assert(exception != (ExceptionInfo *) NULL);
1098 assert(exception->signature == MagickCoreSignature);
1099 if (IsEventLogging() != MagickFalse)
1100 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1101#if defined(MAGICKCORE_OPENCL_SUPPORT)
1102 if (AccelerateFunctionImage(image,function,number_parameters,parameters,
1103 exception) != MagickFalse)
1104 return(MagickTrue);
1105#endif
1106 if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
1107 return(MagickFalse);
1108 status=MagickTrue;
1109 progress=0;
1110 image_view=AcquireAuthenticCacheView(image,exception);
1111#if defined(MAGICKCORE_OPENMP_SUPPORT)
1112 #pragma omp parallel for schedule(static) shared(progress,status) \
1113 magick_number_threads(image,image,image->rows,1)
1114#endif
1115 for (y=0; y < (ssize_t) image->rows; y++)
1116 {
1117 Quantum
1118 *magick_restrict q;
1119
1120 ssize_t
1121 x;
1122
1123 if (status == MagickFalse)
1124 continue;
1125 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
1126 if (q == (Quantum *) NULL)
1127 {
1128 status=MagickFalse;
1129 continue;
1130 }
1131 for (x=0; x < (ssize_t) image->columns; x++)
1132 {
1133 ssize_t
1134 i;
1135
1136 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1137 {
1138 PixelChannel channel = GetPixelChannelChannel(image,i);
1139 PixelTrait traits = GetPixelChannelTraits(image,channel);
1140 if (traits == UndefinedPixelTrait)
1141 continue;
1142 if ((traits & UpdatePixelTrait) == 0)
1143 continue;
1144 q[i]=ApplyFunction(q[i],function,number_parameters,parameters,
1145 exception);
1146 }
1147 q+=(ptrdiff_t) GetPixelChannels(image);
1148 }
1149 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1150 status=MagickFalse;
1151 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1152 {
1153 MagickBooleanType
1154 proceed;
1155
1156#if defined(MAGICKCORE_OPENMP_SUPPORT)
1157 #pragma omp atomic
1158#endif
1159 progress++;
1160 proceed=SetImageProgress(image,FunctionImageTag,progress,image->rows);
1161 if (proceed == MagickFalse)
1162 status=MagickFalse;
1163 }
1164 }
1165 image_view=DestroyCacheView(image_view);
1166 return(status);
1167}
1168␌
1169/*
1170%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1171% %
1172% %
1173% %
1174% G e t I m a g e E n t r o p y %
1175% %
1176% %
1177% %
1178%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1179%
1180% GetImageEntropy() returns the entropy of one or more image channels.
1181%
1182% The format of the GetImageEntropy method is:
1183%
1184% MagickBooleanType GetImageEntropy(const Image *image,double *entropy,
1185% ExceptionInfo *exception)
1186%
1187% A description of each parameter follows:
1188%
1189% o image: the image.
1190%
1191% o entropy: the average entropy of the selected channels.
1192%
1193% o exception: return any errors or warnings in this structure.
1194%
1195*/
1196MagickExport MagickBooleanType GetImageEntropy(const Image *image,
1197 double *entropy,ExceptionInfo *exception)
1198{
1199 ChannelStatistics
1200 *channel_statistics;
1201
1202 assert(image != (Image *) NULL);
1203 assert(image->signature == MagickCoreSignature);
1204 if (IsEventLogging() != MagickFalse)
1205 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1206 channel_statistics=GetImageStatistics(image,exception);
1207 if (channel_statistics == (ChannelStatistics *) NULL)
1208 {
1209 *entropy=NAN;
1210 return(MagickFalse);
1211 }
1212 *entropy=channel_statistics[CompositePixelChannel].entropy;
1213 channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1214 channel_statistics);
1215 return(MagickTrue);
1216}
1217␌
1218/*
1219%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1220% %
1221% %
1222% %
1223% G e t I m a g e E x t r e m a %
1224% %
1225% %
1226% %
1227%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1228%
1229% GetImageExtrema() returns the extrema of one or more image channels.
1230%
1231% The format of the GetImageExtrema method is:
1232%
1233% MagickBooleanType GetImageExtrema(const Image *image,size_t *minima,
1234% size_t *maxima,ExceptionInfo *exception)
1235%
1236% A description of each parameter follows:
1237%
1238% o image: the image.
1239%
1240% o minima: the minimum value in the channel.
1241%
1242% o maxima: the maximum value in the channel.
1243%
1244% o exception: return any errors or warnings in this structure.
1245%
1246*/
1247MagickExport MagickBooleanType GetImageExtrema(const Image *image,
1248 size_t *minima,size_t *maxima,ExceptionInfo *exception)
1249{
1250 double
1251 max,
1252 min;
1253
1254 MagickBooleanType
1255 status;
1256
1257 assert(image != (Image *) NULL);
1258 assert(image->signature == MagickCoreSignature);
1259 if (IsEventLogging() != MagickFalse)
1260 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1261 status=GetImageRange(image,&min,&max,exception);
1262 *minima=(size_t) ceil(min-0.5);
1263 *maxima=(size_t) floor(max+0.5);
1264 return(status);
1265}
1266␌
1267/*
1268%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1269% %
1270% %
1271% %
1272% G e t I m a g e K u r t o s i s %
1273% %
1274% %
1275% %
1276%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1277%
1278% GetImageKurtosis() returns the kurtosis and skewness of one or more image
1279% channels.
1280%
1281% The format of the GetImageKurtosis method is:
1282%
1283% MagickBooleanType GetImageKurtosis(const Image *image,double *kurtosis,
1284% double *skewness,ExceptionInfo *exception)
1285%
1286% A description of each parameter follows:
1287%
1288% o image: the image.
1289%
1290% o kurtosis: the kurtosis of the channel.
1291%
1292% o skewness: the skewness of the channel.
1293%
1294% o exception: return any errors or warnings in this structure.
1295%
1296*/
1297MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
1298 double *kurtosis,double *skewness,ExceptionInfo *exception)
1299{
1300 ChannelStatistics
1301 *channel_statistics;
1302
1303 assert(image != (Image *) NULL);
1304 assert(image->signature == MagickCoreSignature);
1305 if (IsEventLogging() != MagickFalse)
1306 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1307 channel_statistics=GetImageStatistics(image,exception);
1308 if (channel_statistics == (ChannelStatistics *) NULL)
1309 {
1310 *kurtosis=NAN;
1311 *skewness=NAN;
1312 return(MagickFalse);
1313 }
1314 *kurtosis=channel_statistics[CompositePixelChannel].kurtosis;
1315 *skewness=channel_statistics[CompositePixelChannel].skewness;
1316 channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1317 channel_statistics);
1318 return(MagickTrue);
1319}
1320␌
1321/*
1322%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1323% %
1324% %
1325% %
1326% G e t I m a g e M e a n %
1327% %
1328% %
1329% %
1330%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1331%
1332% GetImageMean() returns the mean and standard deviation of one or more image
1333% channels.
1334%
1335% The format of the GetImageMean method is:
1336%
1337% MagickBooleanType GetImageMean(const Image *image,double *mean,
1338% double *standard_deviation,ExceptionInfo *exception)
1339%
1340% A description of each parameter follows:
1341%
1342% o image: the image.
1343%
1344% o mean: the average value in the channel.
1345%
1346% o standard_deviation: the standard deviation of the channel.
1347%
1348% o exception: return any errors or warnings in this structure.
1349%
1350*/
1351MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean,
1352 double *standard_deviation,ExceptionInfo *exception)
1353{
1354 ChannelStatistics
1355 *channel_statistics;
1356
1357 assert(image != (Image *) NULL);
1358 assert(image->signature == MagickCoreSignature);
1359 if (IsEventLogging() != MagickFalse)
1360 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1361 channel_statistics=GetImageStatistics(image,exception);
1362 if (channel_statistics == (ChannelStatistics *) NULL)
1363 {
1364 *mean=NAN;
1365 *standard_deviation=NAN;
1366 return(MagickFalse);
1367 }
1368 *mean=channel_statistics[CompositePixelChannel].mean;
1369 *standard_deviation=
1370 channel_statistics[CompositePixelChannel].standard_deviation;
1371 channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1372 channel_statistics);
1373 return(MagickTrue);
1374}
1375␌
1376/*
1377%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1378% %
1379% %
1380% %
1381% G e t I m a g e M e d i a n %
1382% %
1383% %
1384% %
1385%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1386%
1387% GetImageMedian() returns the median pixel of one or more image channels.
1388%
1389% The format of the GetImageMedian method is:
1390%
1391% MagickBooleanType GetImageMedian(const Image *image,double *median,
1392% ExceptionInfo *exception)
1393%
1394% A description of each parameter follows:
1395%
1396% o image: the image.
1397%
1398% o median: the average value in the channel.
1399%
1400% o exception: return any errors or warnings in this structure.
1401%
1402*/
1403MagickExport MagickBooleanType GetImageMedian(const Image *image,double *median,
1404 ExceptionInfo *exception)
1405{
1406 ChannelStatistics
1407 *channel_statistics;
1408
1409 assert(image != (Image *) NULL);
1410 assert(image->signature == MagickCoreSignature);
1411 if (IsEventLogging() != MagickFalse)
1412 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1413 channel_statistics=GetImageStatistics(image,exception);
1414 if (channel_statistics == (ChannelStatistics *) NULL)
1415 {
1416 *median=NAN;
1417 return(MagickFalse);
1418 }
1419 *median=channel_statistics[CompositePixelChannel].median;
1420 channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1421 channel_statistics);
1422 return(MagickTrue);
1423}
1424␌
1425/*
1426%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1427% %
1428% %
1429% %
1430% G e t I m a g e M o m e n t s %
1431% %
1432% %
1433% %
1434%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1435%
1436% GetImageMoments() returns the normalized moments of one or more image
1437% channels.
1438%
1439% The format of the GetImageMoments method is:
1440%
1441% ChannelMoments *GetImageMoments(const Image *image,
1442% ExceptionInfo *exception)
1443%
1444% A description of each parameter follows:
1445%
1446% o image: the image.
1447%
1448% o exception: return any errors or warnings in this structure.
1449%
1450*/
1451MagickExport ChannelMoments *GetImageMoments(const Image *image,
1452 ExceptionInfo *exception)
1453{
1454#define MaxNumberImageMoments 8
1455
1456 CacheView
1457 *image_view;
1458
1459 ChannelMoments
1460 *channel_moments;
1461
1462 double
1463 channels,
1464 M00[2*MaxPixelChannels+1],
1465 M01[2*MaxPixelChannels+1],
1466 M02[2*MaxPixelChannels+1],
1467 M03[2*MaxPixelChannels+1],
1468 M10[2*MaxPixelChannels+1],
1469 M11[2*MaxPixelChannels+1],
1470 M12[2*MaxPixelChannels+1],
1471 M20[2*MaxPixelChannels+1],
1472 M21[2*MaxPixelChannels+1],
1473 M22[2*MaxPixelChannels+1],
1474 M30[2*MaxPixelChannels+1];
1475
1476 PointInfo
1477 centroid[2*MaxPixelChannels+1];
1478
1479 ssize_t
1480 c,
1481 y;
1482
1483 assert(image != (Image *) NULL);
1484 assert(image->signature == MagickCoreSignature);
1485 if (IsEventLogging() != MagickFalse)
1486 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1487 channel_moments=(ChannelMoments *) AcquireQuantumMemory(MaxPixelChannels+1,
1488 sizeof(*channel_moments));
1489 if (channel_moments == (ChannelMoments *) NULL)
1490 return(channel_moments);
1491 (void) memset(channel_moments,0,(MaxPixelChannels+1)*
1492 sizeof(*channel_moments));
1493 (void) memset(centroid,0,sizeof(centroid));
1494 (void) memset(M00,0,sizeof(M00));
1495 (void) memset(M01,0,sizeof(M01));
1496 (void) memset(M02,0,sizeof(M02));
1497 (void) memset(M03,0,sizeof(M03));
1498 (void) memset(M10,0,sizeof(M10));
1499 (void) memset(M11,0,sizeof(M11));
1500 (void) memset(M12,0,sizeof(M12));
1501 (void) memset(M20,0,sizeof(M20));
1502 (void) memset(M21,0,sizeof(M21));
1503 (void) memset(M22,0,sizeof(M22));
1504 (void) memset(M30,0,sizeof(M30));
1505 image_view=AcquireVirtualCacheView(image,exception);
1506 for (y=0; y < (ssize_t) image->rows; y++)
1507 {
1508 const Quantum
1509 *magick_restrict p;
1510
1511 ssize_t
1512 x;
1513
1514 /*
1515 Compute center of mass (centroid).
1516 */
1517 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1518 if (p == (const Quantum *) NULL)
1519 break;
1520 for (x=0; x < (ssize_t) image->columns; x++)
1521 {
1522 ssize_t
1523 i;
1524
1525 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1526 {
1527 PixelChannel channel = GetPixelChannelChannel(image,i);
1528 PixelTrait traits = GetPixelChannelTraits(image,channel);
1529 if (traits == UndefinedPixelTrait)
1530 continue;
1531 if ((traits & UpdatePixelTrait) == 0)
1532 continue;
1533 M00[channel]+=QuantumScale*(double) p[i];
1534 M00[MaxPixelChannels]+=QuantumScale*(double) p[i];
1535 M10[channel]+=x*QuantumScale*(double) p[i];
1536 M10[MaxPixelChannels]+=x*QuantumScale*(double) p[i];
1537 M01[channel]+=y*QuantumScale*(double) p[i];
1538 M01[MaxPixelChannels]+=y*QuantumScale*(double) p[i];
1539 }
1540 p+=(ptrdiff_t) GetPixelChannels(image);
1541 }
1542 }
1543 for (c=0; c <= MaxPixelChannels; c++)
1544 {
1545 /*
1546 Compute center of mass (centroid).
1547 */
1548 centroid[c].x=M10[c]*PerceptibleReciprocal(M00[c]);
1549 centroid[c].y=M01[c]*PerceptibleReciprocal(M00[c]);
1550 }
1551 for (y=0; y < (ssize_t) image->rows; y++)
1552 {
1553 const Quantum
1554 *magick_restrict p;
1555
1556 ssize_t
1557 x;
1558
1559 /*
1560 Compute the image moments.
1561 */
1562 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1563 if (p == (const Quantum *) NULL)
1564 break;
1565 for (x=0; x < (ssize_t) image->columns; x++)
1566 {
1567 ssize_t
1568 i;
1569
1570 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1571 {
1572 PixelChannel channel = GetPixelChannelChannel(image,i);
1573 PixelTrait traits = GetPixelChannelTraits(image,channel);
1574 if (traits == UndefinedPixelTrait)
1575 continue;
1576 if ((traits & UpdatePixelTrait) == 0)
1577 continue;
1578 M11[channel]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
1579 QuantumScale*(double) p[i];
1580 M11[MaxPixelChannels]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
1581 QuantumScale*(double) p[i];
1582 M20[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1583 QuantumScale*(double) p[i];
1584 M20[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1585 QuantumScale*(double) p[i];
1586 M02[channel]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
1587 QuantumScale*(double) p[i];
1588 M02[MaxPixelChannels]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
1589 QuantumScale*(double) p[i];
1590 M21[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1591 (y-centroid[channel].y)*QuantumScale*(double) p[i];
1592 M21[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1593 (y-centroid[channel].y)*QuantumScale*(double) p[i];
1594 M12[channel]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
1595 (y-centroid[channel].y)*QuantumScale*(double) p[i];
1596 M12[MaxPixelChannels]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
1597 (y-centroid[channel].y)*QuantumScale*(double) p[i];
1598 M22[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1599 (y-centroid[channel].y)*(y-centroid[channel].y)*QuantumScale*(double)
1600 p[i];
1601 M22[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1602 (y-centroid[channel].y)*(y-centroid[channel].y)*QuantumScale*(double)
1603 p[i];
1604 M30[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1605 (x-centroid[channel].x)*QuantumScale*(double) p[i];
1606 M30[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1607 (x-centroid[channel].x)*QuantumScale*(double) p[i];
1608 M03[channel]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
1609 (y-centroid[channel].y)*QuantumScale*(double) p[i];
1610 M03[MaxPixelChannels]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
1611 (y-centroid[channel].y)*QuantumScale*(double) p[i];
1612 }
1613 p+=(ptrdiff_t) GetPixelChannels(image);
1614 }
1615 }
1616 channels=(double) GetImageChannels(image);
1617 M00[MaxPixelChannels]/=channels;
1618 M01[MaxPixelChannels]/=channels;
1619 M02[MaxPixelChannels]/=channels;
1620 M03[MaxPixelChannels]/=channels;
1621 M10[MaxPixelChannels]/=channels;
1622 M11[MaxPixelChannels]/=channels;
1623 M12[MaxPixelChannels]/=channels;
1624 M20[MaxPixelChannels]/=channels;
1625 M21[MaxPixelChannels]/=channels;
1626 M22[MaxPixelChannels]/=channels;
1627 M30[MaxPixelChannels]/=channels;
1628 for (c=0; c <= MaxPixelChannels; c++)
1629 {
1630 /*
1631 Compute elliptical angle, major and minor axes, eccentricity, & intensity.
1632 */
1633 channel_moments[c].centroid=centroid[c];
1634 channel_moments[c].ellipse_axis.x=sqrt((2.0*PerceptibleReciprocal(M00[c]))*
1635 ((M20[c]+M02[c])+sqrt(4.0*M11[c]*M11[c]+(M20[c]-M02[c])*(M20[c]-M02[c]))));
1636 channel_moments[c].ellipse_axis.y=sqrt((2.0*PerceptibleReciprocal(M00[c]))*
1637 ((M20[c]+M02[c])-sqrt(4.0*M11[c]*M11[c]+(M20[c]-M02[c])*(M20[c]-M02[c]))));
1638 channel_moments[c].ellipse_angle=RadiansToDegrees(1.0/2.0*atan(2.0*
1639 M11[c]*PerceptibleReciprocal(M20[c]-M02[c])));
1640 if (fabs(M11[c]) < 0.0)
1641 {
1642 if ((fabs(M20[c]-M02[c]) >= 0.0) &&
1643 ((M20[c]-M02[c]) < 0.0))
1644 channel_moments[c].ellipse_angle+=90.0;
1645 }
1646 else
1647 if (M11[c] < 0.0)
1648 {
1649 if (fabs(M20[c]-M02[c]) >= 0.0)
1650 {
1651 if ((M20[c]-M02[c]) < 0.0)
1652 channel_moments[c].ellipse_angle+=90.0;
1653 else
1654 channel_moments[c].ellipse_angle+=180.0;
1655 }
1656 }
1657 else
1658 if ((fabs(M20[c]-M02[c]) >= 0.0) && ((M20[c]-M02[c]) < 0.0))
1659 channel_moments[c].ellipse_angle+=90.0;
1660 channel_moments[c].ellipse_eccentricity=sqrt(1.0-(
1661 channel_moments[c].ellipse_axis.y*
1662 channel_moments[c].ellipse_axis.y*PerceptibleReciprocal(
1663 channel_moments[c].ellipse_axis.x*
1664 channel_moments[c].ellipse_axis.x)));
1665 channel_moments[c].ellipse_intensity=M00[c]*
1666 PerceptibleReciprocal(MagickPI*channel_moments[c].ellipse_axis.x*
1667 channel_moments[c].ellipse_axis.y+MagickEpsilon);
1668 }
1669 for (c=0; c <= MaxPixelChannels; c++)
1670 {
1671 /*
1672 Normalize image moments.
1673 */
1674 M10[c]=0.0;
1675 M01[c]=0.0;
1676 M11[c]*=PerceptibleReciprocal(pow(M00[c],1.0+(1.0+1.0)/2.0));
1677 M20[c]*=PerceptibleReciprocal(pow(M00[c],1.0+(2.0+0.0)/2.0));
1678 M02[c]*=PerceptibleReciprocal(pow(M00[c],1.0+(0.0+2.0)/2.0));
1679 M21[c]*=PerceptibleReciprocal(pow(M00[c],1.0+(2.0+1.0)/2.0));
1680 M12[c]*=PerceptibleReciprocal(pow(M00[c],1.0+(1.0+2.0)/2.0));
1681 M22[c]*=PerceptibleReciprocal(pow(M00[c],1.0+(2.0+2.0)/2.0));
1682 M30[c]*=PerceptibleReciprocal(pow(M00[c],1.0+(3.0+0.0)/2.0));
1683 M03[c]*=PerceptibleReciprocal(pow(M00[c],1.0+(0.0+3.0)/2.0));
1684 M00[c]=1.0;
1685 }
1686 image_view=DestroyCacheView(image_view);
1687 for (c=0; c <= MaxPixelChannels; c++)
1688 {
1689 /*
1690 Compute Hu invariant moments.
1691 */
1692 channel_moments[c].invariant[0]=M20[c]+M02[c];
1693 channel_moments[c].invariant[1]=(M20[c]-M02[c])*(M20[c]-M02[c])+4.0*M11[c]*
1694 M11[c];
1695 channel_moments[c].invariant[2]=(M30[c]-3.0*M12[c])*(M30[c]-3.0*M12[c])+
1696 (3.0*M21[c]-M03[c])*(3.0*M21[c]-M03[c]);
1697 channel_moments[c].invariant[3]=(M30[c]+M12[c])*(M30[c]+M12[c])+
1698 (M21[c]+M03[c])*(M21[c]+M03[c]);
1699 channel_moments[c].invariant[4]=(M30[c]-3.0*M12[c])*(M30[c]+M12[c])*
1700 ((M30[c]+M12[c])*(M30[c]+M12[c])-3.0*(M21[c]+M03[c])*(M21[c]+M03[c]))+
1701 (3.0*M21[c]-M03[c])*(M21[c]+M03[c])*(3.0*(M30[c]+M12[c])*(M30[c]+M12[c])-
1702 (M21[c]+M03[c])*(M21[c]+M03[c]));
1703 channel_moments[c].invariant[5]=(M20[c]-M02[c])*((M30[c]+M12[c])*
1704 (M30[c]+M12[c])-(M21[c]+M03[c])*(M21[c]+M03[c]))+4.0*M11[c]*
1705 (M30[c]+M12[c])*(M21[c]+M03[c]);
1706 channel_moments[c].invariant[6]=(3.0*M21[c]-M03[c])*(M30[c]+M12[c])*
1707 ((M30[c]+M12[c])*(M30[c]+M12[c])-3.0*(M21[c]+M03[c])*(M21[c]+M03[c]))-
1708 (M30[c]-3*M12[c])*(M21[c]+M03[c])*(3.0*(M30[c]+M12[c])*(M30[c]+M12[c])-
1709 (M21[c]+M03[c])*(M21[c]+M03[c]));
1710 channel_moments[c].invariant[7]=M11[c]*((M30[c]+M12[c])*(M30[c]+M12[c])-
1711 (M03[c]+M21[c])*(M03[c]+M21[c]))-(M20[c]-M02[c])*(M30[c]+M12[c])*
1712 (M03[c]+M21[c]);
1713 }
1714 if (y < (ssize_t) image->rows)
1715 channel_moments=(ChannelMoments *) RelinquishMagickMemory(channel_moments);
1716 return(channel_moments);
1717}
1718␌
1719/*
1720%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1721% %
1722% %
1723% %
1724% 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 %
1725% %
1726% %
1727% %
1728%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1729%
1730% GetImagePerceptualHash() returns the perceptual hash of one or more
1731% image channels.
1732%
1733% The format of the GetImagePerceptualHash method is:
1734%
1735% ChannelPerceptualHash *GetImagePerceptualHash(const Image *image,
1736% ExceptionInfo *exception)
1737%
1738% A description of each parameter follows:
1739%
1740% o image: the image.
1741%
1742% o exception: return any errors or warnings in this structure.
1743%
1744*/
1745MagickExport ChannelPerceptualHash *GetImagePerceptualHash(const Image *image,
1746 ExceptionInfo *exception)
1747{
1748 ChannelPerceptualHash
1749 *perceptual_hash;
1750
1751 char
1752 *colorspaces,
1753 *p,
1754 *q;
1755
1756 const char
1757 *artifact;
1758
1759 MagickBooleanType
1760 status;
1761
1762 ssize_t
1763 i;
1764
1765 perceptual_hash=(ChannelPerceptualHash *) AcquireQuantumMemory(
1766 MaxPixelChannels+1UL,sizeof(*perceptual_hash));
1767 if (perceptual_hash == (ChannelPerceptualHash *) NULL)
1768 return((ChannelPerceptualHash *) NULL);
1769 artifact=GetImageArtifact(image,"phash:colorspaces");
1770 if (artifact != (const char *) NULL)
1771 colorspaces=AcquireString(artifact);
1772 else
1773 colorspaces=AcquireString("xyY,HSB");
1774 perceptual_hash[0].number_colorspaces=0;
1775 perceptual_hash[0].number_channels=0;
1776 q=colorspaces;
1777 for (i=0; (p=StringToken(",",&q)) != (char *) NULL; i++)
1778 {
1779 ChannelMoments
1780 *moments;
1781
1782 Image
1783 *hash_image;
1784
1785 size_t
1786 j;
1787
1788 ssize_t
1789 channel,
1790 colorspace;
1791
1792 if (i >= MaximumNumberOfPerceptualColorspaces)
1793 break;
1794 colorspace=ParseCommandOption(MagickColorspaceOptions,MagickFalse,p);
1795 if (colorspace < 0)
1796 break;
1797 perceptual_hash[0].colorspace[i]=(ColorspaceType) colorspace;
1798 hash_image=BlurImage(image,0.0,1.0,exception);
1799 if (hash_image == (Image *) NULL)
1800 break;
1801 hash_image->depth=8;
1802 status=TransformImageColorspace(hash_image,(ColorspaceType) colorspace,
1803 exception);
1804 if (status == MagickFalse)
1805 break;
1806 moments=GetImageMoments(hash_image,exception);
1807 perceptual_hash[0].number_colorspaces++;
1808 perceptual_hash[0].number_channels+=GetImageChannels(hash_image);
1809 hash_image=DestroyImage(hash_image);
1810 if (moments == (ChannelMoments *) NULL)
1811 break;
1812 for (channel=0; channel <= MaxPixelChannels; channel++)
1813 for (j=0; j < MaximumNumberOfImageMoments; j++)
1814 perceptual_hash[channel].phash[i][j]=
1815 (-MagickLog10(moments[channel].invariant[j]));
1816 moments=(ChannelMoments *) RelinquishMagickMemory(moments);
1817 }
1818 colorspaces=DestroyString(colorspaces);
1819 return(perceptual_hash);
1820}
1821␌
1822/*
1823%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1824% %
1825% %
1826% %
1827% G e t I m a g e R a n g e %
1828% %
1829% %
1830% %
1831%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1832%
1833% GetImageRange() returns the range of one or more image channels.
1834%
1835% The format of the GetImageRange method is:
1836%
1837% MagickBooleanType GetImageRange(const Image *image,double *minima,
1838% double *maxima,ExceptionInfo *exception)
1839%
1840% A description of each parameter follows:
1841%
1842% o image: the image.
1843%
1844% o minima: the minimum value in the channel.
1845%
1846% o maxima: the maximum value in the channel.
1847%
1848% o exception: return any errors or warnings in this structure.
1849%
1850*/
1851MagickExport MagickBooleanType GetImageRange(const Image *image,double *minima,
1852 double *maxima,ExceptionInfo *exception)
1853{
1854 CacheView
1855 *image_view;
1856
1857 MagickBooleanType
1858 status;
1859
1860 ssize_t
1861 y;
1862
1863 assert(image != (Image *) NULL);
1864 assert(image->signature == MagickCoreSignature);
1865 if (IsEventLogging() != MagickFalse)
1866 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1867 status=MagickTrue;
1868 *maxima=MagickMinimumValue;
1869 *minima=MagickMaximumValue;
1870 image_view=AcquireVirtualCacheView(image,exception);
1871#if defined(MAGICKCORE_OPENMP_SUPPORT)
1872 #pragma omp parallel for schedule(static) shared(status) \
1873 magick_number_threads(image,image,image->rows,1)
1874#endif
1875 for (y=0; y < (ssize_t) image->rows; y++)
1876 {
1877 double
1878 row_maxima = MagickMinimumValue,
1879 row_minima = MagickMaximumValue;
1880
1881 const Quantum
1882 *magick_restrict p;
1883
1884 ssize_t
1885 x;
1886
1887 if (status == MagickFalse)
1888 continue;
1889 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1890 if (p == (const Quantum *) NULL)
1891 {
1892 status=MagickFalse;
1893 continue;
1894 }
1895 for (x=0; x < (ssize_t) image->columns; x++)
1896 {
1897 ssize_t
1898 i;
1899
1900 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1901 {
1902 PixelChannel channel = GetPixelChannelChannel(image,i);
1903 PixelTrait traits = GetPixelChannelTraits(image,channel);
1904 if (traits == UndefinedPixelTrait)
1905 continue;
1906 if ((traits & UpdatePixelTrait) == 0)
1907 continue;
1908 if ((double) p[i] < row_minima)
1909 row_minima=(double) p[i];
1910 if ((double) p[i] > row_maxima)
1911 row_maxima=(double) p[i];
1912 }
1913 p+=(ptrdiff_t) GetPixelChannels(image);
1914 }
1915#if defined(MAGICKCORE_OPENMP_SUPPORT)
1916 #pragma omp critical (MagickCore_GetImageRange)
1917#endif
1918 {
1919 if (row_minima < *minima)
1920 *minima=row_minima;
1921 if (row_maxima > *maxima)
1922 *maxima=row_maxima;
1923 }
1924 }
1925 image_view=DestroyCacheView(image_view);
1926 return(status);
1927}
1928␌
1929/*
1930%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1931% %
1932% %
1933% %
1934% G e t I m a g e S t a t i s t i c s %
1935% %
1936% %
1937% %
1938%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1939%
1940% GetImageStatistics() returns statistics for each channel in the image. The
1941% statistics include the channel depth, its minima, maxima, mean, standard
1942% deviation, kurtosis and skewness. You can access the red channel mean, for
1943% example, like this:
1944%
1945% channel_statistics=GetImageStatistics(image,exception);
1946% red_mean=channel_statistics[RedPixelChannel].mean;
1947%
1948% Use MagickRelinquishMemory() to free the statistics buffer.
1949%
1950% The format of the GetImageStatistics method is:
1951%
1952% ChannelStatistics *GetImageStatistics(const Image *image,
1953% ExceptionInfo *exception)
1954%
1955% A description of each parameter follows:
1956%
1957% o image: the image.
1958%
1959% o exception: return any errors or warnings in this structure.
1960%
1961*/
1962
1963static ssize_t GetMedianPixel(Quantum *pixels,const size_t n)
1964{
1965#define SwapPixels(alpha,beta) \
1966{ \
1967 Quantum gamma=(alpha); \
1968 (alpha)=(beta);(beta)=gamma; \
1969}
1970
1971 ssize_t
1972 low = 0,
1973 high = (ssize_t) n-1,
1974 median = (low+high)/2;
1975
1976 for ( ; ; )
1977 {
1978 ssize_t
1979 l = low+1,
1980 h = high,
1981 mid = (low+high)/2;
1982
1983 if (high <= low)
1984 return(median);
1985 if (high == (low+1))
1986 {
1987 if (pixels[low] > pixels[high])
1988 SwapPixels(pixels[low],pixels[high]);
1989 return(median);
1990 }
1991 if (pixels[mid] > pixels[high])
1992 SwapPixels(pixels[mid],pixels[high]);
1993 if (pixels[low] > pixels[high])
1994 SwapPixels(pixels[low], pixels[high]);
1995 if (pixels[mid] > pixels[low])
1996 SwapPixels(pixels[mid],pixels[low]);
1997 SwapPixels(pixels[mid],pixels[low+1]);
1998 for ( ; ; )
1999 {
2000 do l++; while (pixels[low] > pixels[l]);
2001 do h--; while (pixels[h] > pixels[low]);
2002 if (h < l)
2003 break;
2004 SwapPixels(pixels[l],pixels[h]);
2005 }
2006 SwapPixels(pixels[low],pixels[h]);
2007 if (h <= median)
2008 low=l;
2009 if (h >= median)
2010 high=h-1;
2011 }
2012}
2013
2014static inline long double PerceptibleReciprocalLD(const long double x)
2015{
2016 long double
2017 sign;
2018
2019 /*
2020 Return 1/x where x is perceptible (not unlimited or infinitesimal).
2021 */
2022 sign=x < 0.0 ? -1.0 : 1.0;
2023 if ((sign*x) >= MagickEpsilon)
2024 return(1.0/x);
2025 return(sign/MagickEpsilon);
2026}
2027
2028MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
2029 ExceptionInfo *exception)
2030{
2031 ChannelStatistics
2032 *channel_statistics;
2033
2034 double
2035 *histogram;
2036
2037 long double
2038 area,
2039 channels;
2040
2041 MagickStatusType
2042 status;
2043
2044 MemoryInfo
2045 *median_info;
2046
2047 Quantum
2048 *median;
2049
2050 QuantumAny
2051 range;
2052
2053 size_t
2054 depth;
2055
2056 ssize_t
2057 i,
2058 y;
2059
2060 assert(image != (Image *) NULL);
2061 assert(image->signature == MagickCoreSignature);
2062 if (IsEventLogging() != MagickFalse)
2063 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2064 histogram=(double *) AcquireQuantumMemory(MaxMap+1UL,
2065 MagickMax(GetPixelChannels(image),1)*sizeof(*histogram));
2066 channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(
2067 MaxPixelChannels+1,sizeof(*channel_statistics));
2068 if ((channel_statistics == (ChannelStatistics *) NULL) ||
2069 (histogram == (double *) NULL))
2070 {
2071 if (histogram != (double *) NULL)
2072 histogram=(double *) RelinquishMagickMemory(histogram);
2073 if (channel_statistics != (ChannelStatistics *) NULL)
2074 channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
2075 channel_statistics);
2076 (void) ThrowMagickException(exception,GetMagickModule(),
2077 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
2078 return(channel_statistics);
2079 }
2080 (void) memset(channel_statistics,0,(MaxPixelChannels+1)*
2081 sizeof(*channel_statistics));
2082 for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
2083 {
2084 ChannelStatistics *cs = channel_statistics+i;
2085 cs->area=0.0;
2086 cs->depth=1;
2087 cs->maxima=(-MagickMaximumValue);
2088 cs->minima=MagickMaximumValue;
2089 cs->sum=0.0;
2090 cs->sumLD=0.0;
2091 cs->mean=0.0;
2092 cs->standard_deviation=0.0;
2093 cs->variance=0.0;
2094 cs->skewness=0.0;
2095 cs->kurtosis=0.0;
2096 cs->entropy=0.0;
2097 }
2098 (void) memset(histogram,0,(MaxMap+1)*GetPixelChannels(image)*
2099 sizeof(*histogram));
2100 for (y=0; y < (ssize_t) image->rows; y++)
2101 {
2102 const Quantum
2103 *magick_restrict p;
2104
2105 ssize_t
2106 x;
2107
2108 /*
2109 Compute pixel statistics.
2110 */
2111 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
2112 if (p == (const Quantum *) NULL)
2113 break;
2114 for (x=0; x < (ssize_t) image->columns; x++)
2115 {
2116 if (GetPixelReadMask(image,p) <= (QuantumRange/2))
2117 {
2118 p+=(ptrdiff_t) GetPixelChannels(image);
2119 continue;
2120 }
2121 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2122 {
2123 ChannelStatistics
2124 *cs;
2125
2126 PixelChannel channel = GetPixelChannelChannel(image,i);
2127 PixelTrait traits = GetPixelChannelTraits(image,channel);
2128 if (traits == UndefinedPixelTrait)
2129 continue;
2130 cs=channel_statistics+channel;
2131 if (cs->depth != MAGICKCORE_QUANTUM_DEPTH)
2132 {
2133 depth=cs->depth;
2134 range=GetQuantumRange(depth);
2135 status=p[i] != ScaleAnyToQuantum(ScaleQuantumToAny(p[i],range),range) ?
2136 MagickTrue : MagickFalse;
2137 if (status != MagickFalse)
2138 {
2139 cs->depth++;
2140 if (cs->depth > channel_statistics[CompositePixelChannel].depth)
2141 channel_statistics[CompositePixelChannel].depth=cs->depth;
2142 i--;
2143 continue;
2144 }
2145 }
2146 cs->area++;
2147 if ((double) p[i] < cs->minima)
2148 cs->minima=(double) p[i];
2149 if ((double) p[i] > cs->maxima)
2150 cs->maxima=(double) p[i];
2151 histogram[(ssize_t) GetPixelChannels(image)*ScaleQuantumToMap(
2152 ClampToQuantum((double) p[i]))+i]++;
2153 cs->sumLD+=(long double) p[i];
2154 /*
2155 sum_squared, sum_cubed and sum_fourth_power are not used in
2156 MagickCore or MagickWand, but are made available in
2157 Magick++/lib/Statistic.cpp, so we need to calculate these.
2158 */
2159 cs->sum_squared+=(double) p[i]*(double) p[i];
2160 cs->sum_cubed+=(double) p[i]*(double) p[i]*(double) p[i];
2161 cs->sum_fourth_power+=(double) p[i]*(double) p[i]*(double) p[i]*
2162 (double) p[i];
2163 {
2164 /* Calculate running totals for Welford's method.
2165 */
2166 double
2167 n = cs->area,
2168 n1 = cs->area-1;
2169
2170 long double
2171 delta,
2172 delta_n,
2173 delta_n2,
2174 term1;
2175
2176 delta=(double) p[i]-cs->M1;
2177 delta_n=delta/n;
2178 delta_n2=delta_n*delta_n;
2179 term1=delta*delta_n*n1;
2180 cs->M4+=term1*delta_n2*(n*n-3.0*n+3.0)+6.0*delta_n2*cs->M2-4.0*
2181 delta_n*cs->M3;
2182 cs->M3+=term1*delta_n*(n-2.0)-3.0*delta_n*cs->M2;
2183 cs->M2+=term1;
2184 cs->M1+=delta_n;
2185 }
2186 }
2187 p+=(ptrdiff_t) GetPixelChannels(image);
2188 }
2189 }
2190 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2191 {
2192 ChannelStatistics
2193 *cs;
2194
2195 PixelChannel channel = GetPixelChannelChannel(image,i);
2196 PixelTrait traits = GetPixelChannelTraits(image,channel);
2197 double AdjArea = 1.0;
2198 if (traits == UndefinedPixelTrait)
2199 continue;
2200 cs=channel_statistics+channel;
2201 cs->mean=0.0;
2202 if (cs->area > 0)
2203 {
2204 cs->mean=(double) (cs->sumLD/(long double) cs->area);
2205 if (cs->area > 1.0)
2206 AdjArea=cs->area/(cs->area-1.0);
2207 }
2208 cs->sum=(double) cs->sum;
2209 if (cs->M2 == 0.0)
2210 {
2211 cs->standard_deviation=0.0;
2212 cs->variance=0.0;
2213 cs->skewness=0.0;
2214 cs->kurtosis=0.0;
2215 }
2216 else
2217 {
2218 if (cs->area > 1.0)
2219 cs->standard_deviation=(double) sqrtl(cs->M2/((long double) cs->area-1.0));
2220 else
2221 cs->standard_deviation=(double) sqrtl(cs->M2/((long double) cs->area));
2222 cs->variance=cs->standard_deviation*cs->standard_deviation;
2223 cs->skewness=(double) (sqrtl(cs->area)*cs->M3/powl(cs->M2*AdjArea,1.5));
2224 cs->kurtosis=(double) (cs->area*cs->M4/(cs->M2*cs->M2*AdjArea*AdjArea)-3.0);
2225 }
2226 }
2227
2228 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2229 {
2230 double
2231 number_bins;
2232
2233 ssize_t
2234 j;
2235
2236 PixelChannel channel = GetPixelChannelChannel(image,i);
2237 ChannelStatistics *cs = channel_statistics+channel;
2238 if (cs->area > 0.0)
2239 {
2240 cs->sum/=cs->area;
2241 cs->sum_squared/=cs->area;
2242 cs->sum_cubed/=cs->area;
2243 cs->sum_fourth_power/=cs->area;
2244 }
2245 /*
2246 Compute pixel entropy.
2247 */
2248 number_bins=0.0;
2249 for (j=0; j <= (ssize_t) MaxMap; j++)
2250 if (histogram[(ssize_t) GetPixelChannels(image)*j+i] > 0.0)
2251 number_bins++;
2252 area=PerceptibleReciprocalLD(channel_statistics[channel].area);
2253 for (j=0; j <= (ssize_t) MaxMap; j++)
2254 {
2255 double
2256 count;
2257
2258 count=(double) (area*histogram[(ssize_t) GetPixelChannels(image)*j+i]);
2259 channel_statistics[channel].entropy+=((long double) -count*
2260 MagickLog10(count)*PerceptibleReciprocalLD((long double)
2261 MagickLog10(number_bins)));
2262 channel_statistics[CompositePixelChannel].entropy+=((long double) -count*
2263 MagickLog10(count)*PerceptibleReciprocalLD((long double)
2264 MagickLog10(number_bins))/GetPixelChannels(image));
2265 }
2266 }
2267 histogram=(double *) RelinquishMagickMemory(histogram);
2268 median_info=AcquireVirtualMemory(image->columns,image->rows*sizeof(*median));
2269 if (median_info == (MemoryInfo *) NULL)
2270 (void) ThrowMagickException(exception,GetMagickModule(),
2271 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
2272 else
2273 {
2274 median=(Quantum *) GetVirtualMemoryBlob(median_info);
2275 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2276 {
2277 size_t
2278 n = 0;
2279
2280 /*
2281 Compute median statistics for each channel.
2282 */
2283 PixelChannel channel = GetPixelChannelChannel(image,i);
2284 PixelTrait traits = GetPixelChannelTraits(image,channel);
2285 if (traits == UndefinedPixelTrait)
2286 continue;
2287 if ((traits & UpdatePixelTrait) == 0)
2288 continue;
2289 for (y=0; y < (ssize_t) image->rows; y++)
2290 {
2291 const Quantum
2292 *magick_restrict p;
2293
2294 ssize_t
2295 x;
2296
2297 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
2298 if (p == (const Quantum *) NULL)
2299 break;
2300 for (x=0; x < (ssize_t) image->columns; x++)
2301 {
2302 if (GetPixelReadMask(image,p) <= (QuantumRange/2))
2303 {
2304 p+=(ptrdiff_t) GetPixelChannels(image);
2305 continue;
2306 }
2307 median[n++]=p[i];
2308 p+=(ptrdiff_t) GetPixelChannels(image);
2309 }
2310 }
2311 channel_statistics[channel].median=(double) median[
2312 GetMedianPixel(median,n)];
2313 }
2314 median_info=RelinquishVirtualMemory(median_info);
2315 }
2316 {
2317 ChannelStatistics *csComp = channel_statistics+CompositePixelChannel;
2318 csComp->sum=0.0;
2319 csComp->sum_squared=0.0;
2320 csComp->sum_cubed=0.0;
2321 csComp->sum_fourth_power=0.0;
2322 csComp->maxima=(-MagickMaximumValue);
2323 csComp->minima=MagickMaximumValue;
2324 csComp->area=0.0;
2325 csComp->mean=0.0;
2326 csComp->median=0.0;
2327 csComp->variance=0.0;
2328 csComp->standard_deviation=0.0;
2329 csComp->entropy=0.0;
2330 csComp->skewness=0.0;
2331 csComp->kurtosis=0.0;
2332 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2333 {
2334 ChannelStatistics
2335 *cs;
2336
2337 PixelChannel channel = GetPixelChannelChannel(image,i);
2338 PixelTrait traits = GetPixelChannelTraits(image,channel);
2339 if (traits == UndefinedPixelTrait)
2340 continue;
2341 if ((traits & UpdatePixelTrait) == 0)
2342 continue;
2343 cs=channel_statistics+channel;
2344 if (csComp->maxima < cs->maxima)
2345 csComp->maxima=cs->maxima;
2346 if (csComp->minima > cs->minima)
2347 csComp->minima=cs->minima;
2348 csComp->sum+=cs->sum;
2349 csComp->sum_squared+=cs->sum_squared;
2350 csComp->sum_cubed+=cs->sum_cubed;
2351 csComp->sum_fourth_power+=cs->sum_fourth_power;
2352 csComp->median+=cs->median;
2353 csComp->area+=cs->area;
2354 csComp->mean+=cs->mean;
2355 csComp->variance+=cs->variance;
2356 csComp->standard_deviation+=cs->standard_deviation;
2357 csComp->skewness+=cs->skewness;
2358 csComp->kurtosis+=cs->kurtosis;
2359 csComp->entropy+=cs->entropy;
2360 }
2361 channels=(double) GetImageChannels(image);
2362 csComp->sum/=channels;
2363 csComp->sum_squared/=channels;
2364 csComp->sum_cubed/=channels;
2365 csComp->sum_fourth_power/=channels;
2366 csComp->median/=channels;
2367 csComp->area/=channels;
2368 csComp->mean/=channels;
2369 csComp->variance/=channels;
2370 csComp->standard_deviation/=channels;
2371 csComp->skewness/=channels;
2372 csComp->kurtosis/=channels;
2373 csComp->entropy/=channels;
2374 }
2375 if (y < (ssize_t) image->rows)
2376 channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
2377 channel_statistics);
2378 return(channel_statistics);
2379}
2380␌
2381/*
2382%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2383% %
2384% %
2385% %
2386% P o l y n o m i a l I m a g e %
2387% %
2388% %
2389% %
2390%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2391%
2392% PolynomialImage() returns a new image where each pixel is the sum of the
2393% pixels in the image sequence after applying its corresponding terms
2394% (coefficient and degree pairs).
2395%
2396% The format of the PolynomialImage method is:
2397%
2398% Image *PolynomialImage(const Image *images,const size_t number_terms,
2399% const double *terms,ExceptionInfo *exception)
2400%
2401% A description of each parameter follows:
2402%
2403% o images: the image sequence.
2404%
2405% o number_terms: the number of terms in the list. The actual list length
2406% is 2 x number_terms + 1 (the constant).
2407%
2408% o terms: the list of polynomial coefficients and degree pairs and a
2409% constant.
2410%
2411% o exception: return any errors or warnings in this structure.
2412%
2413*/
2414MagickExport Image *PolynomialImage(const Image *images,
2415 const size_t number_terms,const double *terms,ExceptionInfo *exception)
2416{
2417#define PolynomialImageTag "Polynomial/Image"
2418
2419 CacheView
2420 *polynomial_view;
2421
2422 Image
2423 *image;
2424
2425 MagickBooleanType
2426 status;
2427
2428 MagickOffsetType
2429 progress;
2430
2431 PixelChannels
2432 **magick_restrict polynomial_pixels;
2433
2434 size_t
2435 number_images;
2436
2437 ssize_t
2438 y;
2439
2440 assert(images != (Image *) NULL);
2441 assert(images->signature == MagickCoreSignature);
2442 if (IsEventLogging() != MagickFalse)
2443 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
2444 assert(exception != (ExceptionInfo *) NULL);
2445 assert(exception->signature == MagickCoreSignature);
2446 image=AcquireImageCanvas(images,exception);
2447 if (image == (Image *) NULL)
2448 return((Image *) NULL);
2449 if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
2450 {
2451 image=DestroyImage(image);
2452 return((Image *) NULL);
2453 }
2454 number_images=GetImageListLength(images);
2455 polynomial_pixels=AcquirePixelTLS(images);
2456 if (polynomial_pixels == (PixelChannels **) NULL)
2457 {
2458 image=DestroyImage(image);
2459 (void) ThrowMagickException(exception,GetMagickModule(),
2460 ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
2461 return((Image *) NULL);
2462 }
2463 /*
2464 Polynomial image pixels.
2465 */
2466 status=MagickTrue;
2467 progress=0;
2468 polynomial_view=AcquireAuthenticCacheView(image,exception);
2469#if defined(MAGICKCORE_OPENMP_SUPPORT)
2470 #pragma omp parallel for schedule(static) shared(progress,status) \
2471 magick_number_threads(image,image,image->rows,1)
2472#endif
2473 for (y=0; y < (ssize_t) image->rows; y++)
2474 {
2475 CacheView
2476 *image_view;
2477
2478 const Image
2479 *next;
2480
2481 const int
2482 id = GetOpenMPThreadId();
2483
2484 PixelChannels
2485 *polynomial_pixel;
2486
2487 Quantum
2488 *magick_restrict q;
2489
2490 ssize_t
2491 i,
2492 j,
2493 x;
2494
2495 if (status == MagickFalse)
2496 continue;
2497 q=QueueCacheViewAuthenticPixels(polynomial_view,0,y,image->columns,1,
2498 exception);
2499 if (q == (Quantum *) NULL)
2500 {
2501 status=MagickFalse;
2502 continue;
2503 }
2504 polynomial_pixel=polynomial_pixels[id];
2505 for (j=0; j < (ssize_t) image->columns; j++)
2506 for (i=0; i < MaxPixelChannels; i++)
2507 polynomial_pixel[j].channel[i]=0.0;
2508 next=images;
2509 for (j=0; j < (ssize_t) number_images; j++)
2510 {
2511 const Quantum
2512 *p;
2513
2514 if (j >= (ssize_t) number_terms)
2515 continue;
2516 image_view=AcquireVirtualCacheView(next,exception);
2517 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2518 if (p == (const Quantum *) NULL)
2519 {
2520 image_view=DestroyCacheView(image_view);
2521 break;
2522 }
2523 for (x=0; x < (ssize_t) image->columns; x++)
2524 {
2525 for (i=0; i < (ssize_t) GetPixelChannels(next); i++)
2526 {
2527 MagickRealType
2528 coefficient,
2529 degree;
2530
2531 PixelChannel channel = GetPixelChannelChannel(image,i);
2532 PixelTrait traits = GetPixelChannelTraits(next,channel);
2533 PixelTrait polynomial_traits=GetPixelChannelTraits(image,channel);
2534 if ((traits == UndefinedPixelTrait) ||
2535 (polynomial_traits == UndefinedPixelTrait))
2536 continue;
2537 if ((traits & UpdatePixelTrait) == 0)
2538 continue;
2539 coefficient=(MagickRealType) terms[2*j];
2540 degree=(MagickRealType) terms[(j << 1)+1];
2541 polynomial_pixel[x].channel[i]+=coefficient*
2542 pow(QuantumScale*(double) GetPixelChannel(image,channel,p),degree);
2543 }
2544 p+=(ptrdiff_t) GetPixelChannels(next);
2545 }
2546 image_view=DestroyCacheView(image_view);
2547 next=GetNextImageInList(next);
2548 }
2549 for (x=0; x < (ssize_t) image->columns; x++)
2550 {
2551 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2552 {
2553 PixelChannel channel = GetPixelChannelChannel(image,i);
2554 PixelTrait traits = GetPixelChannelTraits(image,channel);
2555 if (traits == UndefinedPixelTrait)
2556 continue;
2557 if ((traits & UpdatePixelTrait) == 0)
2558 continue;
2559 q[i]=ClampToQuantum((double) QuantumRange*
2560 polynomial_pixel[x].channel[i]);
2561 }
2562 q+=(ptrdiff_t) GetPixelChannels(image);
2563 }
2564 if (SyncCacheViewAuthenticPixels(polynomial_view,exception) == MagickFalse)
2565 status=MagickFalse;
2566 if (images->progress_monitor != (MagickProgressMonitor) NULL)
2567 {
2568 MagickBooleanType
2569 proceed;
2570
2571#if defined(MAGICKCORE_OPENMP_SUPPORT)
2572 #pragma omp atomic
2573#endif
2574 progress++;
2575 proceed=SetImageProgress(images,PolynomialImageTag,progress,
2576 image->rows);
2577 if (proceed == MagickFalse)
2578 status=MagickFalse;
2579 }
2580 }
2581 polynomial_view=DestroyCacheView(polynomial_view);
2582 polynomial_pixels=DestroyPixelTLS(images,polynomial_pixels);
2583 if (status == MagickFalse)
2584 image=DestroyImage(image);
2585 return(image);
2586}
2587␌
2588/*
2589%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2590% %
2591% %
2592% %
2593% S t a t i s t i c I m a g e %
2594% %
2595% %
2596% %
2597%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2598%
2599% StatisticImage() makes each pixel the min / max / median / mode / etc. of
2600% the neighborhood of the specified width and height.
2601%
2602% The format of the StatisticImage method is:
2603%
2604% Image *StatisticImage(const Image *image,const StatisticType type,
2605% const size_t width,const size_t height,ExceptionInfo *exception)
2606%
2607% A description of each parameter follows:
2608%
2609% o image: the image.
2610%
2611% o type: the statistic type (median, mode, etc.).
2612%
2613% o width: the width of the pixel neighborhood.
2614%
2615% o height: the height of the pixel neighborhood.
2616%
2617% o exception: return any errors or warnings in this structure.
2618%
2619*/
2620
2621typedef struct _SkipNode
2622{
2623 size_t
2624 next[9],
2625 count,
2626 signature;
2627} SkipNode;
2628
2629typedef struct _SkipList
2630{
2631 ssize_t
2632 level;
2633
2634 SkipNode
2635 *nodes;
2636} SkipList;
2637
2638typedef struct _PixelList
2639{
2640 size_t
2641 length,
2642 seed;
2643
2644 SkipList
2645 skip_list;
2646
2647 size_t
2648 signature;
2649} PixelList;
2650
2651static PixelList *DestroyPixelList(PixelList *pixel_list)
2652{
2653 if (pixel_list == (PixelList *) NULL)
2654 return((PixelList *) NULL);
2655 if (pixel_list->skip_list.nodes != (SkipNode *) NULL)
2656 pixel_list->skip_list.nodes=(SkipNode *) RelinquishAlignedMemory(
2657 pixel_list->skip_list.nodes);
2658 pixel_list=(PixelList *) RelinquishMagickMemory(pixel_list);
2659 return(pixel_list);
2660}
2661
2662static PixelList **DestroyPixelListTLS(PixelList **pixel_list)
2663{
2664 ssize_t
2665 i;
2666
2667 assert(pixel_list != (PixelList **) NULL);
2668 for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
2669 if (pixel_list[i] != (PixelList *) NULL)
2670 pixel_list[i]=DestroyPixelList(pixel_list[i]);
2671 pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list);
2672 return(pixel_list);
2673}
2674
2675static PixelList *AcquirePixelList(const size_t width,const size_t height)
2676{
2677 PixelList
2678 *pixel_list;
2679
2680 pixel_list=(PixelList *) AcquireMagickMemory(sizeof(*pixel_list));
2681 if (pixel_list == (PixelList *) NULL)
2682 return(pixel_list);
2683 (void) memset((void *) pixel_list,0,sizeof(*pixel_list));
2684 pixel_list->length=width*height;
2685 pixel_list->skip_list.nodes=(SkipNode *) AcquireAlignedMemory(65537UL,
2686 sizeof(*pixel_list->skip_list.nodes));
2687 if (pixel_list->skip_list.nodes == (SkipNode *) NULL)
2688 return(DestroyPixelList(pixel_list));
2689 (void) memset(pixel_list->skip_list.nodes,0,65537UL*
2690 sizeof(*pixel_list->skip_list.nodes));
2691 pixel_list->signature=MagickCoreSignature;
2692 return(pixel_list);
2693}
2694
2695static PixelList **AcquirePixelListTLS(const size_t width,
2696 const size_t height)
2697{
2698 PixelList
2699 **pixel_list;
2700
2701 ssize_t
2702 i;
2703
2704 size_t
2705 number_threads;
2706
2707 number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
2708 pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
2709 sizeof(*pixel_list));
2710 if (pixel_list == (PixelList **) NULL)
2711 return((PixelList **) NULL);
2712 (void) memset(pixel_list,0,number_threads*sizeof(*pixel_list));
2713 for (i=0; i < (ssize_t) number_threads; i++)
2714 {
2715 pixel_list[i]=AcquirePixelList(width,height);
2716 if (pixel_list[i] == (PixelList *) NULL)
2717 return(DestroyPixelListTLS(pixel_list));
2718 }
2719 return(pixel_list);
2720}
2721
2722static void AddNodePixelList(PixelList *pixel_list,const size_t color)
2723{
2724 SkipList
2725 *p;
2726
2727 ssize_t
2728 level;
2729
2730 size_t
2731 search,
2732 update[9];
2733
2734 /*
2735 Initialize the node.
2736 */
2737 p=(&pixel_list->skip_list);
2738 p->nodes[color].signature=pixel_list->signature;
2739 p->nodes[color].count=1;
2740 /*
2741 Determine where it belongs in the list.
2742 */
2743 search=65536UL;
2744 (void) memset(update,0,sizeof(update));
2745 for (level=p->level; level >= 0; level--)
2746 {
2747 while (p->nodes[search].next[level] < color)
2748 search=p->nodes[search].next[level];
2749 update[level]=search;
2750 }
2751 /*
2752 Generate a pseudo-random level for this node.
2753 */
2754 for (level=0; ; level++)
2755 {
2756 pixel_list->seed=(pixel_list->seed*42893621L)+1L;
2757 if ((pixel_list->seed & 0x300) != 0x300)
2758 break;
2759 }
2760 if (level > 8)
2761 level=8;
2762 if (level > (p->level+2))
2763 level=p->level+2;
2764 /*
2765 If we're raising the list's level, link back to the root node.
2766 */
2767 while (level > p->level)
2768 {
2769 p->level++;
2770 update[p->level]=65536UL;
2771 }
2772 /*
2773 Link the node into the skip-list.
2774 */
2775 do
2776 {
2777 p->nodes[color].next[level]=p->nodes[update[level]].next[level];
2778 p->nodes[update[level]].next[level]=color;
2779 } while (level-- > 0);
2780}
2781
2782static inline void GetMedianPixelList(PixelList *pixel_list,Quantum *pixel)
2783{
2784 SkipList
2785 *p;
2786
2787 size_t
2788 color;
2789
2790 ssize_t
2791 count;
2792
2793 /*
2794 Find the median value for each of the color.
2795 */
2796 p=(&pixel_list->skip_list);
2797 color=65536L;
2798 count=0;
2799 do
2800 {
2801 color=p->nodes[color].next[0];
2802 count+=(ssize_t) p->nodes[color].count;
2803 } while (count <= (ssize_t) (pixel_list->length >> 1));
2804 *pixel=ScaleShortToQuantum((unsigned short) color);
2805}
2806
2807static inline void GetModePixelList(PixelList *pixel_list,Quantum *pixel)
2808{
2809 SkipList
2810 *p;
2811
2812 size_t
2813 color,
2814 max_count,
2815 mode;
2816
2817 ssize_t
2818 count;
2819
2820 /*
2821 Make each pixel the 'predominant color' of the specified neighborhood.
2822 */
2823 p=(&pixel_list->skip_list);
2824 color=65536L;
2825 mode=color;
2826 max_count=p->nodes[mode].count;
2827 count=0;
2828 do
2829 {
2830 color=p->nodes[color].next[0];
2831 if (p->nodes[color].count > max_count)
2832 {
2833 mode=color;
2834 max_count=p->nodes[mode].count;
2835 }
2836 count+=(ssize_t) p->nodes[color].count;
2837 } while (count < (ssize_t) pixel_list->length);
2838 *pixel=ScaleShortToQuantum((unsigned short) mode);
2839}
2840
2841static inline void GetNonpeakPixelList(PixelList *pixel_list,Quantum *pixel)
2842{
2843 SkipList
2844 *p;
2845
2846 size_t
2847 color,
2848 next,
2849 previous;
2850
2851 ssize_t
2852 count;
2853
2854 /*
2855 Finds the non peak value for each of the colors.
2856 */
2857 p=(&pixel_list->skip_list);
2858 color=65536L;
2859 next=p->nodes[color].next[0];
2860 count=0;
2861 do
2862 {
2863 previous=color;
2864 color=next;
2865 next=p->nodes[color].next[0];
2866 count+=(ssize_t) p->nodes[color].count;
2867 } while (count <= (ssize_t) (pixel_list->length >> 1));
2868 if ((previous == 65536UL) && (next != 65536UL))
2869 color=next;
2870 else
2871 if ((previous != 65536UL) && (next == 65536UL))
2872 color=previous;
2873 *pixel=ScaleShortToQuantum((unsigned short) color);
2874}
2875
2876static inline void InsertPixelList(const Quantum pixel,PixelList *pixel_list)
2877{
2878 size_t
2879 signature;
2880
2881 unsigned short
2882 index;
2883
2884 index=ScaleQuantumToShort(pixel);
2885 signature=pixel_list->skip_list.nodes[index].signature;
2886 if (signature == pixel_list->signature)
2887 {
2888 pixel_list->skip_list.nodes[index].count++;
2889 return;
2890 }
2891 AddNodePixelList(pixel_list,index);
2892}
2893
2894static void ResetPixelList(PixelList *pixel_list)
2895{
2896 int
2897 level;
2898
2899 SkipNode
2900 *root;
2901
2902 SkipList
2903 *p;
2904
2905 /*
2906 Reset the skip-list.
2907 */
2908 p=(&pixel_list->skip_list);
2909 root=p->nodes+65536UL;
2910 p->level=0;
2911 for (level=0; level < 9; level++)
2912 root->next[level]=65536UL;
2913 pixel_list->seed=pixel_list->signature++;
2914}
2915
2916MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
2917 const size_t width,const size_t height,ExceptionInfo *exception)
2918{
2919#define StatisticImageTag "Statistic/Image"
2920
2921 CacheView
2922 *image_view,
2923 *statistic_view;
2924
2925 Image
2926 *statistic_image;
2927
2928 MagickBooleanType
2929 status;
2930
2931 MagickOffsetType
2932 progress;
2933
2934 PixelList
2935 **magick_restrict pixel_list;
2936
2937 ssize_t
2938 center,
2939 y;
2940
2941 /*
2942 Initialize statistics image attributes.
2943 */
2944 assert(image != (Image *) NULL);
2945 assert(image->signature == MagickCoreSignature);
2946 if (IsEventLogging() != MagickFalse)
2947 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2948 assert(exception != (ExceptionInfo *) NULL);
2949 assert(exception->signature == MagickCoreSignature);
2950 statistic_image=CloneImage(image,0,0,MagickTrue,
2951 exception);
2952 if (statistic_image == (Image *) NULL)
2953 return((Image *) NULL);
2954 status=SetImageStorageClass(statistic_image,DirectClass,exception);
2955 if (status == MagickFalse)
2956 {
2957 statistic_image=DestroyImage(statistic_image);
2958 return((Image *) NULL);
2959 }
2960 pixel_list=AcquirePixelListTLS(MagickMax(width,1),MagickMax(height,1));
2961 if (pixel_list == (PixelList **) NULL)
2962 {
2963 statistic_image=DestroyImage(statistic_image);
2964 ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2965 }
2966 /*
2967 Make each pixel the min / max / median / mode / etc. of the neighborhood.
2968 */
2969 center=(ssize_t) GetPixelChannels(image)*((ssize_t) image->columns+
2970 MagickMax((ssize_t) width,1L))*(MagickMax((ssize_t) height,1)/2L)+(ssize_t)
2971 GetPixelChannels(image)*(MagickMax((ssize_t) width,1L)/2L);
2972 status=MagickTrue;
2973 progress=0;
2974 image_view=AcquireVirtualCacheView(image,exception);
2975 statistic_view=AcquireAuthenticCacheView(statistic_image,exception);
2976#if defined(MAGICKCORE_OPENMP_SUPPORT)
2977 #pragma omp parallel for schedule(static) shared(progress,status) \
2978 magick_number_threads(image,statistic_image,statistic_image->rows,1)
2979#endif
2980 for (y=0; y < (ssize_t) statistic_image->rows; y++)
2981 {
2982 const int
2983 id = GetOpenMPThreadId();
2984
2985 const Quantum
2986 *magick_restrict p;
2987
2988 Quantum
2989 *magick_restrict q;
2990
2991 ssize_t
2992 x;
2993
2994 if (status == MagickFalse)
2995 continue;
2996 p=GetCacheViewVirtualPixels(image_view,-((ssize_t) MagickMax(width,1)/2L),y-
2997 (ssize_t) (MagickMax(height,1)/2L),image->columns+MagickMax(width,1),
2998 MagickMax(height,1),exception);
2999 q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns, 1,exception);
3000 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3001 {
3002 status=MagickFalse;
3003 continue;
3004 }
3005 for (x=0; x < (ssize_t) statistic_image->columns; x++)
3006 {
3007 ssize_t
3008 i;
3009
3010 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3011 {
3012 double
3013 area,
3014 maximum,
3015 minimum,
3016 sum,
3017 sum_squared;
3018
3019 Quantum
3020 pixel;
3021
3022 const Quantum
3023 *magick_restrict pixels;
3024
3025 ssize_t
3026 u;
3027
3028 ssize_t
3029 v;
3030
3031 PixelChannel channel = GetPixelChannelChannel(image,i);
3032 PixelTrait traits = GetPixelChannelTraits(image,channel);
3033 PixelTrait statistic_traits=GetPixelChannelTraits(statistic_image,
3034 channel);
3035 if ((traits == UndefinedPixelTrait) ||
3036 (statistic_traits == UndefinedPixelTrait))
3037 continue;
3038 if (((statistic_traits & CopyPixelTrait) != 0) ||
3039 (GetPixelWriteMask(image,p) <= (QuantumRange/2)))
3040 {
3041 SetPixelChannel(statistic_image,channel,p[center+i],q);
3042 continue;
3043 }
3044 if ((statistic_traits & UpdatePixelTrait) == 0)
3045 continue;
3046 pixels=p;
3047 area=0.0;
3048 minimum=pixels[i];
3049 maximum=pixels[i];
3050 sum=0.0;
3051 sum_squared=0.0;
3052 ResetPixelList(pixel_list[id]);
3053 for (v=0; v < (ssize_t) MagickMax(height,1); v++)
3054 {
3055 for (u=0; u < (ssize_t) MagickMax(width,1); u++)
3056 {
3057 if ((type == MedianStatistic) || (type == ModeStatistic) ||
3058 (type == NonpeakStatistic))
3059 {
3060 InsertPixelList(pixels[i],pixel_list[id]);
3061 pixels+=(ptrdiff_t) GetPixelChannels(image);
3062 continue;
3063 }
3064 area++;
3065 if ((double) pixels[i] < minimum)
3066 minimum=(double) pixels[i];
3067 if ((double) pixels[i] > maximum)
3068 maximum=(double) pixels[i];
3069 sum+=(double) pixels[i];
3070 sum_squared+=(double) pixels[i]*(double) pixels[i];
3071 pixels+=(ptrdiff_t) GetPixelChannels(image);
3072 }
3073 pixels+=(ptrdiff_t) GetPixelChannels(image)*image->columns;
3074 }
3075 switch (type)
3076 {
3077 case ContrastStatistic:
3078 {
3079 pixel=ClampToQuantum(MagickAbsoluteValue((maximum-minimum)*
3080 PerceptibleReciprocal(maximum+minimum)));
3081 break;
3082 }
3083 case GradientStatistic:
3084 {
3085 pixel=ClampToQuantum(MagickAbsoluteValue(maximum-minimum));
3086 break;
3087 }
3088 case MaximumStatistic:
3089 {
3090 pixel=ClampToQuantum(maximum);
3091 break;
3092 }
3093 case MeanStatistic:
3094 default:
3095 {
3096 pixel=ClampToQuantum(sum/area);
3097 break;
3098 }
3099 case MedianStatistic:
3100 {
3101 GetMedianPixelList(pixel_list[id],&pixel);
3102 break;
3103 }
3104 case MinimumStatistic:
3105 {
3106 pixel=ClampToQuantum(minimum);
3107 break;
3108 }
3109 case ModeStatistic:
3110 {
3111 GetModePixelList(pixel_list[id],&pixel);
3112 break;
3113 }
3114 case NonpeakStatistic:
3115 {
3116 GetNonpeakPixelList(pixel_list[id],&pixel);
3117 break;
3118 }
3119 case RootMeanSquareStatistic:
3120 {
3121 pixel=ClampToQuantum(sqrt(sum_squared/area));
3122 break;
3123 }
3124 case StandardDeviationStatistic:
3125 {
3126 pixel=ClampToQuantum(sqrt(sum_squared/area-(sum/area*sum/area)));
3127 break;
3128 }
3129 }
3130 SetPixelChannel(statistic_image,channel,pixel,q);
3131 }
3132 p+=(ptrdiff_t) GetPixelChannels(image);
3133 q+=(ptrdiff_t) GetPixelChannels(statistic_image);
3134 }
3135 if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
3136 status=MagickFalse;
3137 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3138 {
3139 MagickBooleanType
3140 proceed;
3141
3142#if defined(MAGICKCORE_OPENMP_SUPPORT)
3143 #pragma omp atomic
3144#endif
3145 progress++;
3146 proceed=SetImageProgress(image,StatisticImageTag,progress,image->rows);
3147 if (proceed == MagickFalse)
3148 status=MagickFalse;
3149 }
3150 }
3151 statistic_view=DestroyCacheView(statistic_view);
3152 image_view=DestroyCacheView(image_view);
3153 pixel_list=DestroyPixelListTLS(pixel_list);
3154 if (status == MagickFalse)
3155 statistic_image=DestroyImage(statistic_image);
3156 return(statistic_image);
3157}
_CacheView
Definition cache-view.c:66
_ChannelMoments
Definition statistic.h:59
_ChannelPerceptualHash
Definition statistic.h:74
_ChannelStatistics
Definition statistic.h:30
_ExceptionInfo
Definition exception.h:102
_Image
Definition image.h:132
_MemoryInfo
Definition memory.c:164
_PixelChannels
Definition statistic.c:136
_PixelList
Definition statistic.c:2639
_PointInfo
Definition geometry.h:123
_RandomInfo
Definition random.c:84
_SkipList
Definition statistic.c:2630
_SkipNode
Definition statistic.c:2622