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