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