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