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