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