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