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