1 /*
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % %
4 % %
5 % %
6 % SSSSS TTTTT AAA TTTTT IIIII SSSSS TTTTT IIIII CCCC %
7 % SS T A A T I SS T I C %
8 % SSS T AAAAA T I SSS T I C %
9 % SS T A A T I SS T I C %
10 % SSSSS T A A T IIIII SSSSS T IIIII CCCC %
11 % %
12 % %
13 % MagickCore Image Statistical Methods %
14 % %
15 % Software Design %
16 % Cristy %
17 % July 1992 %
18 % %
19 % %
20 % Copyright 1999-2021 ImageMagick Studio LLC, a non-profit organization %
21 % dedicated to making software imaging solutions freely available. %
22 % %
23 % You may not use this file except in compliance with the License. You may %
24 % obtain a copy of the License at %
25 % %
26 % https://imagemagick.org/script/license.php %
27 % %
28 % Unless required by applicable law or agreed to in writing, software %
29 % distributed under the License is distributed on an "AS IS" BASIS, %
30 % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. %
31 % See the License for the specific language governing permissions and %
32 % limitations under the License. %
33 % %
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35 %
36 %
37 %
38 */
39
40
41 /*
42 Include declarations.
43 */
44 #include "magick/studio.h"
45 #include "magick/accelerate-private.h"
46 #include "magick/animate.h"
47 #include "magick/animate.h"
48 #include "magick/blob.h"
49 #include "magick/blob-private.h"
50 #include "magick/cache.h"
51 #include "magick/cache-private.h"
52 #include "magick/cache-view.h"
53 #include "magick/client.h"
54 #include "magick/color.h"
55 #include "magick/color-private.h"
56 #include "magick/colorspace.h"
57 #include "magick/colorspace-private.h"
58 #include "magick/composite.h"
59 #include "magick/composite-private.h"
60 #include "magick/compress.h"
61 #include "magick/constitute.h"
62 #include "magick/deprecate.h"
63 #include "magick/display.h"
64 #include "magick/draw.h"
65 #include "magick/enhance.h"
66 #include "magick/exception.h"
67 #include "magick/exception-private.h"
68 #include "magick/gem.h"
69 #include "magick/geometry.h"
70 #include "magick/list.h"
71 #include "magick/image-private.h"
72 #include "magick/magic.h"
73 #include "magick/magick.h"
74 #include "magick/memory_.h"
75 #include "magick/module.h"
76 #include "magick/monitor.h"
77 #include "magick/monitor-private.h"
78 #include "magick/option.h"
79 #include "magick/paint.h"
80 #include "magick/pixel-private.h"
81 #include "magick/profile.h"
82 #include "magick/property.h"
83 #include "magick/quantize.h"
84 #include "magick/random_.h"
85 #include "magick/random-private.h"
86 #include "magick/resource_.h"
87 #include "magick/segment.h"
88 #include "magick/semaphore.h"
89 #include "magick/signature-private.h"
90 #include "magick/statistic.h"
91 #include "magick/string_.h"
92 #include "magick/thread-private.h"
93 #include "magick/timer.h"
94 #include "magick/utility.h"
95 #include "magick/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 EvaluateImageChannel 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 % MagickBooleanType EvaluateImageChannel(Image *image,
122 % const ChannelType channel,const MagickEvaluateOperator op,
123 % const double value,ExceptionInfo *exception)
124 %
125 % A description of each parameter follows:
126 %
127 % o image: the image.
128 %
129 % o channel: the channel.
130 %
131 % o op: A channel op.
132 %
133 % o value: A value value.
134 %
135 % o exception: return any errors or warnings in this structure.
136 %
137 */
138
DestroyPixelThreadSet(const Image * images,MagickPixelPacket ** pixels)139 static MagickPixelPacket **DestroyPixelThreadSet(const Image *images,
140 MagickPixelPacket **pixels)
141 {
142 ssize_t
143 i;
144
145 size_t
146 rows;
147
148 assert(pixels != (MagickPixelPacket **) NULL);
149 rows=MagickMax(GetImageListLength(images),
150 (size_t) GetMagickResourceLimit(ThreadResource));
151 for (i=0; i < (ssize_t) rows; i++)
152 if (pixels[i] != (MagickPixelPacket *) NULL)
153 pixels[i]=(MagickPixelPacket *) RelinquishMagickMemory(pixels[i]);
154 pixels=(MagickPixelPacket **) RelinquishMagickMemory(pixels);
155 return(pixels);
156 }
157
AcquirePixelThreadSet(const Image * images)158 static MagickPixelPacket **AcquirePixelThreadSet(const Image *images)
159 {
160 const Image
161 *next;
162
163 MagickPixelPacket
164 **pixels;
165
166 ssize_t
167 i,
168 j;
169
170 size_t
171 columns,
172 rows;
173
174 rows=MagickMax(GetImageListLength(images),
175 (size_t) GetMagickResourceLimit(ThreadResource));
176 pixels=(MagickPixelPacket **) AcquireQuantumMemory(rows,sizeof(*pixels));
177 if (pixels == (MagickPixelPacket **) NULL)
178 return((MagickPixelPacket **) NULL);
179 (void) memset(pixels,0,rows*sizeof(*pixels));
180 columns=GetImageListLength(images);
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 pixels[i]=(MagickPixelPacket *) AcquireQuantumMemory(columns,
186 sizeof(**pixels));
187 if (pixels[i] == (MagickPixelPacket *) NULL)
188 return(DestroyPixelThreadSet(images,pixels));
189 for (j=0; j < (ssize_t) columns; j++)
190 GetMagickPixelPacket(images,&pixels[i][j]);
191 }
192 return(pixels);
193 }
194
EvaluateMax(const double x,const double y)195 static inline double EvaluateMax(const double x,const double y)
196 {
197 if (x > y)
198 return(x);
199 return(y);
200 }
201
202 #if defined(__cplusplus) || defined(c_plusplus)
203 extern "C" {
204 #endif
205
IntensityCompare(const void * x,const void * y)206 static int IntensityCompare(const void *x,const void *y)
207 {
208 const MagickPixelPacket
209 *color_1,
210 *color_2;
211
212 int
213 intensity;
214
215 color_1=(const MagickPixelPacket *) x;
216 color_2=(const MagickPixelPacket *) y;
217 intensity=(int) MagickPixelIntensity(color_2)-(int)
218 MagickPixelIntensity(color_1);
219 return(intensity);
220 }
221
222 #if defined(__cplusplus) || defined(c_plusplus)
223 }
224 #endif
225
ApplyEvaluateOperator(RandomInfo * random_info,const Quantum pixel,const MagickEvaluateOperator op,const MagickRealType value)226 static MagickRealType ApplyEvaluateOperator(RandomInfo *random_info,
227 const Quantum pixel,const MagickEvaluateOperator op,
228 const MagickRealType value)
229 {
230 MagickRealType
231 result;
232
233 ssize_t
234 i;
235
236 result=0.0;
237 switch (op)
238 {
239 case UndefinedEvaluateOperator:
240 break;
241 case AbsEvaluateOperator:
242 {
243 result=(MagickRealType) fabs((double) (pixel+value));
244 break;
245 }
246 case AddEvaluateOperator:
247 {
248 result=(MagickRealType) (pixel+value);
249 break;
250 }
251 case AddModulusEvaluateOperator:
252 {
253 /*
254 This returns a 'floored modulus' of the addition which is a
255 positive result. It differs from % or fmod() which returns a
256 'truncated modulus' result, where floor() is replaced by trunc()
257 and could return a negative result (which is clipped).
258 */
259 result=pixel+value;
260 result-=(QuantumRange+1.0)*floor((double) result/(QuantumRange+1.0));
261 break;
262 }
263 case AndEvaluateOperator:
264 {
265 result=(MagickRealType) ((ssize_t) pixel & (ssize_t) (value+0.5));
266 break;
267 }
268 case CosineEvaluateOperator:
269 {
270 result=(MagickRealType) (QuantumRange*(0.5*cos((double) (2.0*MagickPI*
271 QuantumScale*pixel*value))+0.5));
272 break;
273 }
274 case DivideEvaluateOperator:
275 {
276 result=pixel/(value == 0.0 ? 1.0 : value);
277 break;
278 }
279 case ExponentialEvaluateOperator:
280 {
281 result=(MagickRealType) (QuantumRange*exp((double) (value*QuantumScale*
282 pixel)));
283 break;
284 }
285 case GaussianNoiseEvaluateOperator:
286 {
287 result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
288 GaussianNoise,value);
289 break;
290 }
291 case ImpulseNoiseEvaluateOperator:
292 {
293 result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
294 ImpulseNoise,value);
295 break;
296 }
297 case InverseLogEvaluateOperator:
298 {
299 result=(QuantumRange*pow((value+1.0),QuantumScale*pixel)-1.0)*
300 PerceptibleReciprocal(value);
301 break;
302 }
303 case LaplacianNoiseEvaluateOperator:
304 {
305 result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
306 LaplacianNoise,value);
307 break;
308 }
309 case LeftShiftEvaluateOperator:
310 {
311 result=(double) pixel;
312 for (i=0; i < (ssize_t) value; i++)
313 result*=2.0;
314 break;
315 }
316 case LogEvaluateOperator:
317 {
318 if ((QuantumScale*pixel) >= MagickEpsilon)
319 result=(MagickRealType) (QuantumRange*log((double) (QuantumScale*value*
320 pixel+1.0))/log((double) (value+1.0)));
321 break;
322 }
323 case MaxEvaluateOperator:
324 {
325 result=(MagickRealType) EvaluateMax((double) pixel,value);
326 break;
327 }
328 case MeanEvaluateOperator:
329 {
330 result=(MagickRealType) (pixel+value);
331 break;
332 }
333 case MedianEvaluateOperator:
334 {
335 result=(MagickRealType) (pixel+value);
336 break;
337 }
338 case MinEvaluateOperator:
339 {
340 result=(MagickRealType) MagickMin((double) pixel,value);
341 break;
342 }
343 case MultiplicativeNoiseEvaluateOperator:
344 {
345 result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
346 MultiplicativeGaussianNoise,value);
347 break;
348 }
349 case MultiplyEvaluateOperator:
350 {
351 result=(MagickRealType) (value*pixel);
352 break;
353 }
354 case OrEvaluateOperator:
355 {
356 result=(MagickRealType) ((ssize_t) pixel | (ssize_t) (value+0.5));
357 break;
358 }
359 case PoissonNoiseEvaluateOperator:
360 {
361 result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
362 PoissonNoise,value);
363 break;
364 }
365 case PowEvaluateOperator:
366 {
367 if (pixel < 0)
368 result=(MagickRealType) -(QuantumRange*pow((double) -(QuantumScale*
369 pixel),(double) value));
370 else
371 result=(MagickRealType) (QuantumRange*pow((double) (QuantumScale*pixel),
372 (double) value));
373 break;
374 }
375 case RightShiftEvaluateOperator:
376 {
377 result=(MagickRealType) pixel;
378 for (i=0; i < (ssize_t) value; i++)
379 result/=2.0;
380 break;
381 }
382 case RootMeanSquareEvaluateOperator:
383 {
384 result=((MagickRealType) pixel*pixel+value);
385 break;
386 }
387 case SetEvaluateOperator:
388 {
389 result=value;
390 break;
391 }
392 case SineEvaluateOperator:
393 {
394 result=(MagickRealType) (QuantumRange*(0.5*sin((double) (2.0*MagickPI*
395 QuantumScale*pixel*value))+0.5));
396 break;
397 }
398 case SubtractEvaluateOperator:
399 {
400 result=(MagickRealType) (pixel-value);
401 break;
402 }
403 case SumEvaluateOperator:
404 {
405 result=(MagickRealType) (pixel+value);
406 break;
407 }
408 case ThresholdEvaluateOperator:
409 {
410 result=(MagickRealType) (((MagickRealType) pixel <= value) ? 0 :
411 QuantumRange);
412 break;
413 }
414 case ThresholdBlackEvaluateOperator:
415 {
416 result=(MagickRealType) (((MagickRealType) pixel <= value) ? 0 : pixel);
417 break;
418 }
419 case ThresholdWhiteEvaluateOperator:
420 {
421 result=(MagickRealType) (((MagickRealType) pixel > value) ? QuantumRange :
422 pixel);
423 break;
424 }
425 case UniformNoiseEvaluateOperator:
426 {
427 result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
428 UniformNoise,value);
429 break;
430 }
431 case XorEvaluateOperator:
432 {
433 result=(MagickRealType) ((ssize_t) pixel ^ (ssize_t) (value+0.5));
434 break;
435 }
436 }
437 return(result);
438 }
439
AcquireImageCanvas(const Image * images,ExceptionInfo * exception)440 static Image *AcquireImageCanvas(const Image *images,ExceptionInfo *exception)
441 {
442 const Image
443 *p,
444 *q;
445
446 size_t
447 columns,
448 number_channels,
449 rows;
450
451 q=images;
452 columns=images->columns;
453 rows=images->rows;
454 number_channels=0;
455 for (p=images; p != (Image *) NULL; p=p->next)
456 {
457 size_t
458 channels;
459
460 channels=3;
461 if (p->matte != MagickFalse)
462 channels+=1;
463 if (p->colorspace == CMYKColorspace)
464 channels+=1;
465 if (channels > number_channels)
466 {
467 number_channels=channels;
468 q=p;
469 }
470 if (p->columns > columns)
471 columns=p->columns;
472 if (p->rows > rows)
473 rows=p->rows;
474 }
475 return(CloneImage(q,columns,rows,MagickTrue,exception));
476 }
477
EvaluateImage(Image * image,const MagickEvaluateOperator op,const double value,ExceptionInfo * exception)478 MagickExport MagickBooleanType EvaluateImage(Image *image,
479 const MagickEvaluateOperator op,const double value,ExceptionInfo *exception)
480 {
481 MagickBooleanType
482 status;
483
484 status=EvaluateImageChannel(image,CompositeChannels,op,value,exception);
485 return(status);
486 }
487
EvaluateImages(const Image * images,const MagickEvaluateOperator op,ExceptionInfo * exception)488 MagickExport Image *EvaluateImages(const Image *images,
489 const MagickEvaluateOperator op,ExceptionInfo *exception)
490 {
491 #define EvaluateImageTag "Evaluate/Image"
492
493 CacheView
494 *evaluate_view;
495
496 Image
497 *image;
498
499 MagickBooleanType
500 status;
501
502 MagickOffsetType
503 progress;
504
505 MagickPixelPacket
506 **magick_restrict evaluate_pixels,
507 zero;
508
509 RandomInfo
510 **magick_restrict random_info;
511
512 size_t
513 number_images;
514
515 ssize_t
516 y;
517
518 #if defined(MAGICKCORE_OPENMP_SUPPORT)
519 unsigned long
520 key;
521 #endif
522
523 assert(images != (Image *) NULL);
524 assert(images->signature == MagickCoreSignature);
525 if (images->debug != MagickFalse)
526 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
527 assert(exception != (ExceptionInfo *) NULL);
528 assert(exception->signature == MagickCoreSignature);
529 image=AcquireImageCanvas(images,exception);
530 if (image == (Image *) NULL)
531 return((Image *) NULL);
532 if (SetImageStorageClass(image,DirectClass) == MagickFalse)
533 {
534 InheritException(exception,&image->exception);
535 image=DestroyImage(image);
536 return((Image *) NULL);
537 }
538 evaluate_pixels=AcquirePixelThreadSet(images);
539 if (evaluate_pixels == (MagickPixelPacket **) NULL)
540 {
541 image=DestroyImage(image);
542 (void) ThrowMagickException(exception,GetMagickModule(),
543 ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
544 return((Image *) NULL);
545 }
546 /*
547 Evaluate image pixels.
548 */
549 status=MagickTrue;
550 progress=0;
551 number_images=GetImageListLength(images);
552 GetMagickPixelPacket(images,&zero);
553 random_info=AcquireRandomInfoThreadSet();
554 evaluate_view=AcquireAuthenticCacheView(image,exception);
555 if (op == MedianEvaluateOperator)
556 {
557 #if defined(MAGICKCORE_OPENMP_SUPPORT)
558 key=GetRandomSecretKey(random_info[0]);
559 #pragma omp parallel for schedule(static) shared(progress,status) \
560 magick_number_threads(image,images,image->rows,key == ~0UL)
561 #endif
562 for (y=0; y < (ssize_t) image->rows; y++)
563 {
564 CacheView
565 *image_view;
566
567 const Image
568 *next;
569
570 const int
571 id = GetOpenMPThreadId();
572
573 IndexPacket
574 *magick_restrict evaluate_indexes;
575
576 MagickPixelPacket
577 *evaluate_pixel;
578
579 PixelPacket
580 *magick_restrict q;
581
582 ssize_t
583 x;
584
585 if (status == MagickFalse)
586 continue;
587 q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,image->columns,1,
588 exception);
589 if (q == (PixelPacket *) NULL)
590 {
591 status=MagickFalse;
592 continue;
593 }
594 evaluate_indexes=GetCacheViewAuthenticIndexQueue(evaluate_view);
595 evaluate_pixel=evaluate_pixels[id];
596 for (x=0; x < (ssize_t) image->columns; x++)
597 {
598 ssize_t
599 i;
600
601 for (i=0; i < (ssize_t) number_images; i++)
602 evaluate_pixel[i]=zero;
603 next=images;
604 for (i=0; i < (ssize_t) number_images; i++)
605 {
606 const IndexPacket
607 *indexes;
608
609 const PixelPacket
610 *p;
611
612 image_view=AcquireVirtualCacheView(next,exception);
613 p=GetCacheViewVirtualPixels(image_view,x,y,1,1,exception);
614 if (p == (const PixelPacket *) NULL)
615 {
616 image_view=DestroyCacheView(image_view);
617 break;
618 }
619 indexes=GetCacheViewVirtualIndexQueue(image_view);
620 evaluate_pixel[i].red=ApplyEvaluateOperator(random_info[id],
621 GetPixelRed(p),op,evaluate_pixel[i].red);
622 evaluate_pixel[i].green=ApplyEvaluateOperator(random_info[id],
623 GetPixelGreen(p),op,evaluate_pixel[i].green);
624 evaluate_pixel[i].blue=ApplyEvaluateOperator(random_info[id],
625 GetPixelBlue(p),op,evaluate_pixel[i].blue);
626 evaluate_pixel[i].opacity=ApplyEvaluateOperator(random_info[id],
627 GetPixelAlpha(p),op,evaluate_pixel[i].opacity);
628 if (image->colorspace == CMYKColorspace)
629 evaluate_pixel[i].index=ApplyEvaluateOperator(random_info[id],
630 *indexes,op,evaluate_pixel[i].index);
631 image_view=DestroyCacheView(image_view);
632 next=GetNextImageInList(next);
633 }
634 qsort((void *) evaluate_pixel,number_images,sizeof(*evaluate_pixel),
635 IntensityCompare);
636 SetPixelRed(q,ClampToQuantum(evaluate_pixel[i/2].red));
637 SetPixelGreen(q,ClampToQuantum(evaluate_pixel[i/2].green));
638 SetPixelBlue(q,ClampToQuantum(evaluate_pixel[i/2].blue));
639 SetPixelAlpha(q,ClampToQuantum(evaluate_pixel[i/2].opacity));
640 if (image->colorspace == CMYKColorspace)
641 SetPixelIndex(evaluate_indexes+i,ClampToQuantum(
642 evaluate_pixel[i/2].index));
643 q++;
644 }
645 if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
646 status=MagickFalse;
647 if (images->progress_monitor != (MagickProgressMonitor) NULL)
648 {
649 MagickBooleanType
650 proceed;
651
652 #if defined(MAGICKCORE_OPENMP_SUPPORT)
653 #pragma omp atomic
654 #endif
655 progress++;
656 proceed=SetImageProgress(images,EvaluateImageTag,progress,
657 image->rows);
658 if (proceed == MagickFalse)
659 status=MagickFalse;
660 }
661 }
662 }
663 else
664 {
665 #if defined(MAGICKCORE_OPENMP_SUPPORT)
666 key=GetRandomSecretKey(random_info[0]);
667 #pragma omp parallel for schedule(static) shared(progress,status) \
668 magick_number_threads(image,images,image->rows,key == ~0UL)
669 #endif
670 for (y=0; y < (ssize_t) image->rows; y++)
671 {
672 CacheView
673 *image_view;
674
675 const Image
676 *next;
677
678 const int
679 id = GetOpenMPThreadId();
680
681 IndexPacket
682 *magick_restrict evaluate_indexes;
683
684 ssize_t
685 i,
686 x;
687
688 MagickPixelPacket
689 *evaluate_pixel;
690
691 PixelPacket
692 *magick_restrict q;
693
694 if (status == MagickFalse)
695 continue;
696 q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,image->columns,1,
697 exception);
698 if (q == (PixelPacket *) NULL)
699 {
700 status=MagickFalse;
701 continue;
702 }
703 evaluate_indexes=GetCacheViewAuthenticIndexQueue(evaluate_view);
704 evaluate_pixel=evaluate_pixels[id];
705 for (x=0; x < (ssize_t) image->columns; x++)
706 evaluate_pixel[x]=zero;
707 next=images;
708 for (i=0; i < (ssize_t) number_images; i++)
709 {
710 const IndexPacket
711 *indexes;
712
713 const PixelPacket
714 *p;
715
716 image_view=AcquireVirtualCacheView(next,exception);
717 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,
718 exception);
719 if (p == (const PixelPacket *) NULL)
720 {
721 image_view=DestroyCacheView(image_view);
722 break;
723 }
724 indexes=GetCacheViewVirtualIndexQueue(image_view);
725 for (x=0; x < (ssize_t) image->columns; x++)
726 {
727 evaluate_pixel[x].red=ApplyEvaluateOperator(random_info[id],
728 GetPixelRed(p),i == 0 ? AddEvaluateOperator : op,
729 evaluate_pixel[x].red);
730 evaluate_pixel[x].green=ApplyEvaluateOperator(random_info[id],
731 GetPixelGreen(p),i == 0 ? AddEvaluateOperator : op,
732 evaluate_pixel[x].green);
733 evaluate_pixel[x].blue=ApplyEvaluateOperator(random_info[id],
734 GetPixelBlue(p),i == 0 ? AddEvaluateOperator : op,
735 evaluate_pixel[x].blue);
736 evaluate_pixel[x].opacity=ApplyEvaluateOperator(random_info[id],
737 GetPixelAlpha(p),i == 0 ? AddEvaluateOperator : op,
738 evaluate_pixel[x].opacity);
739 if (image->colorspace == CMYKColorspace)
740 evaluate_pixel[x].index=ApplyEvaluateOperator(random_info[id],
741 GetPixelIndex(indexes+x),i == 0 ? AddEvaluateOperator : op,
742 evaluate_pixel[x].index);
743 p++;
744 }
745 image_view=DestroyCacheView(image_view);
746 next=GetNextImageInList(next);
747 }
748 if (op == MeanEvaluateOperator)
749 for (x=0; x < (ssize_t) image->columns; x++)
750 {
751 evaluate_pixel[x].red/=number_images;
752 evaluate_pixel[x].green/=number_images;
753 evaluate_pixel[x].blue/=number_images;
754 evaluate_pixel[x].opacity/=number_images;
755 evaluate_pixel[x].index/=number_images;
756 }
757 if (op == RootMeanSquareEvaluateOperator)
758 for (x=0; x < (ssize_t) image->columns; x++)
759 {
760 evaluate_pixel[x].red=sqrt((double) evaluate_pixel[x].red/
761 number_images);
762 evaluate_pixel[x].green=sqrt((double) evaluate_pixel[x].green/
763 number_images);
764 evaluate_pixel[x].blue=sqrt((double) evaluate_pixel[x].blue/
765 number_images);
766 evaluate_pixel[x].opacity=sqrt((double) evaluate_pixel[x].opacity/
767 number_images);
768 evaluate_pixel[x].index=sqrt((double) evaluate_pixel[x].index/
769 number_images);
770 }
771 if (op == MultiplyEvaluateOperator)
772 for (x=0; x < (ssize_t) image->columns; x++)
773 {
774 ssize_t
775 j;
776
777 for (j=0; j < (ssize_t) (number_images-1); j++)
778 {
779 evaluate_pixel[x].red*=(MagickRealType) QuantumScale;
780 evaluate_pixel[x].green*=(MagickRealType) QuantumScale;
781 evaluate_pixel[x].blue*=(MagickRealType) QuantumScale;
782 evaluate_pixel[x].opacity*=(MagickRealType) QuantumScale;
783 evaluate_pixel[x].index*=(MagickRealType) QuantumScale;
784 }
785 }
786 for (x=0; x < (ssize_t) image->columns; x++)
787 {
788 SetPixelRed(q,ClampToQuantum(evaluate_pixel[x].red));
789 SetPixelGreen(q,ClampToQuantum(evaluate_pixel[x].green));
790 SetPixelBlue(q,ClampToQuantum(evaluate_pixel[x].blue));
791 SetPixelAlpha(q,ClampToQuantum(evaluate_pixel[x].opacity));
792 if (image->colorspace == CMYKColorspace)
793 SetPixelIndex(evaluate_indexes+x,ClampToQuantum(
794 evaluate_pixel[x].index));
795 q++;
796 }
797 if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
798 status=MagickFalse;
799 if (images->progress_monitor != (MagickProgressMonitor) NULL)
800 {
801 MagickBooleanType
802 proceed;
803
804 proceed=SetImageProgress(images,EvaluateImageTag,progress++,
805 image->rows);
806 if (proceed == MagickFalse)
807 status=MagickFalse;
808 }
809 }
810 }
811 evaluate_view=DestroyCacheView(evaluate_view);
812 evaluate_pixels=DestroyPixelThreadSet(images,evaluate_pixels);
813 random_info=DestroyRandomInfoThreadSet(random_info);
814 if (status == MagickFalse)
815 image=DestroyImage(image);
816 return(image);
817 }
818
EvaluateImageChannel(Image * image,const ChannelType channel,const MagickEvaluateOperator op,const double value,ExceptionInfo * exception)819 MagickExport MagickBooleanType EvaluateImageChannel(Image *image,
820 const ChannelType channel,const MagickEvaluateOperator op,const double value,
821 ExceptionInfo *exception)
822 {
823 CacheView
824 *image_view;
825
826 MagickBooleanType
827 status;
828
829 MagickOffsetType
830 progress;
831
832 RandomInfo
833 **magick_restrict random_info;
834
835 ssize_t
836 y;
837
838 #if defined(MAGICKCORE_OPENMP_SUPPORT)
839 unsigned long
840 key;
841 #endif
842
843 assert(image != (Image *) NULL);
844 assert(image->signature == MagickCoreSignature);
845 if (image->debug != MagickFalse)
846 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
847 assert(exception != (ExceptionInfo *) NULL);
848 assert(exception->signature == MagickCoreSignature);
849 if (SetImageStorageClass(image,DirectClass) == MagickFalse)
850 {
851 InheritException(exception,&image->exception);
852 return(MagickFalse);
853 }
854 status=MagickTrue;
855 progress=0;
856 random_info=AcquireRandomInfoThreadSet();
857 image_view=AcquireAuthenticCacheView(image,exception);
858 #if defined(MAGICKCORE_OPENMP_SUPPORT)
859 key=GetRandomSecretKey(random_info[0]);
860 #pragma omp parallel for schedule(static) shared(progress,status) \
861 magick_number_threads(image,image,image->rows,key == ~0UL)
862 #endif
863 for (y=0; y < (ssize_t) image->rows; y++)
864 {
865 const int
866 id = GetOpenMPThreadId();
867
868 IndexPacket
869 *magick_restrict indexes;
870
871 PixelPacket
872 *magick_restrict q;
873
874 ssize_t
875 x;
876
877 if (status == MagickFalse)
878 continue;
879 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
880 if (q == (PixelPacket *) NULL)
881 {
882 status=MagickFalse;
883 continue;
884 }
885 indexes=GetCacheViewAuthenticIndexQueue(image_view);
886 for (x=0; x < (ssize_t) image->columns; x++)
887 {
888 MagickRealType
889 result;
890
891 if ((channel & RedChannel) != 0)
892 {
893 result=ApplyEvaluateOperator(random_info[id],GetPixelRed(q),op,value);
894 if (op == MeanEvaluateOperator)
895 result/=2.0;
896 SetPixelRed(q,ClampToQuantum(result));
897 }
898 if ((channel & GreenChannel) != 0)
899 {
900 result=ApplyEvaluateOperator(random_info[id],GetPixelGreen(q),op,
901 value);
902 if (op == MeanEvaluateOperator)
903 result/=2.0;
904 SetPixelGreen(q,ClampToQuantum(result));
905 }
906 if ((channel & BlueChannel) != 0)
907 {
908 result=ApplyEvaluateOperator(random_info[id],GetPixelBlue(q),op,
909 value);
910 if (op == MeanEvaluateOperator)
911 result/=2.0;
912 SetPixelBlue(q,ClampToQuantum(result));
913 }
914 if ((channel & OpacityChannel) != 0)
915 {
916 if (image->matte == MagickFalse)
917 {
918 result=ApplyEvaluateOperator(random_info[id],GetPixelOpacity(q),
919 op,value);
920 if (op == MeanEvaluateOperator)
921 result/=2.0;
922 SetPixelOpacity(q,ClampToQuantum(result));
923 }
924 else
925 {
926 result=ApplyEvaluateOperator(random_info[id],GetPixelAlpha(q),
927 op,value);
928 if (op == MeanEvaluateOperator)
929 result/=2.0;
930 SetPixelAlpha(q,ClampToQuantum(result));
931 }
932 }
933 if (((channel & IndexChannel) != 0) && (indexes != (IndexPacket *) NULL))
934 {
935 result=ApplyEvaluateOperator(random_info[id],GetPixelIndex(indexes+x),
936 op,value);
937 if (op == MeanEvaluateOperator)
938 result/=2.0;
939 SetPixelIndex(indexes+x,ClampToQuantum(result));
940 }
941 q++;
942 }
943 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
944 status=MagickFalse;
945 if (image->progress_monitor != (MagickProgressMonitor) NULL)
946 {
947 MagickBooleanType
948 proceed;
949
950 proceed=SetImageProgress(image,EvaluateImageTag,progress++,image->rows);
951 if (proceed == MagickFalse)
952 status=MagickFalse;
953 }
954 }
955 image_view=DestroyCacheView(image_view);
956 random_info=DestroyRandomInfoThreadSet(random_info);
957 return(status);
958 }
959
960 /*
961 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
962 % %
963 % %
964 % %
965 % F u n c t i o n I m a g e %
966 % %
967 % %
968 % %
969 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
970 %
971 % FunctionImage() applies a value to the image with an arithmetic, relational,
972 % or logical operator to an image. Use these operations to lighten or darken
973 % an image, to increase or decrease contrast in an image, or to produce the
974 % "negative" of an image.
975 %
976 % The format of the FunctionImageChannel method is:
977 %
978 % MagickBooleanType FunctionImage(Image *image,
979 % const MagickFunction function,const ssize_t number_parameters,
980 % const double *parameters,ExceptionInfo *exception)
981 % MagickBooleanType FunctionImageChannel(Image *image,
982 % const ChannelType channel,const MagickFunction function,
983 % const ssize_t number_parameters,const double *argument,
984 % ExceptionInfo *exception)
985 %
986 % A description of each parameter follows:
987 %
988 % o image: the image.
989 %
990 % o channel: the channel.
991 %
992 % o function: A channel function.
993 %
994 % o parameters: one or more parameters.
995 %
996 % o exception: return any errors or warnings in this structure.
997 %
998 */
999
ApplyFunction(Quantum pixel,const MagickFunction function,const size_t number_parameters,const double * parameters,ExceptionInfo * exception)1000 static Quantum ApplyFunction(Quantum pixel,const MagickFunction function,
1001 const size_t number_parameters,const double *parameters,
1002 ExceptionInfo *exception)
1003 {
1004 MagickRealType
1005 result;
1006
1007 ssize_t
1008 i;
1009
1010 (void) exception;
1011 result=0.0;
1012 switch (function)
1013 {
1014 case PolynomialFunction:
1015 {
1016 /*
1017 * Polynomial
1018 * Parameters: polynomial constants, highest to lowest order
1019 * For example: c0*x^3 + c1*x^2 + c2*x + c3
1020 */
1021 result=0.0;
1022 for (i=0; i < (ssize_t) number_parameters; i++)
1023 result=result*QuantumScale*pixel + parameters[i];
1024 result*=QuantumRange;
1025 break;
1026 }
1027 case SinusoidFunction:
1028 {
1029 /* Sinusoid Function
1030 * Parameters: Freq, Phase, Ampl, bias
1031 */
1032 double freq,phase,ampl,bias;
1033 freq = ( number_parameters >= 1 ) ? parameters[0] : 1.0;
1034 phase = ( number_parameters >= 2 ) ? parameters[1] : 0.0;
1035 ampl = ( number_parameters >= 3 ) ? parameters[2] : 0.5;
1036 bias = ( number_parameters >= 4 ) ? parameters[3] : 0.5;
1037 result=(MagickRealType) (QuantumRange*(ampl*sin((double) (2.0*MagickPI*
1038 (freq*QuantumScale*pixel + phase/360.0) )) + bias ) );
1039 break;
1040 }
1041 case ArcsinFunction:
1042 {
1043 double
1044 bias,
1045 center,
1046 range,
1047 width;
1048
1049 /* Arcsin Function (peged at range limits for invalid results)
1050 * Parameters: Width, Center, Range, Bias
1051 */
1052 width=(number_parameters >= 1) ? parameters[0] : 1.0;
1053 center=(number_parameters >= 2) ? parameters[1] : 0.5;
1054 range=(number_parameters >= 3) ? parameters[2] : 1.0;
1055 bias=(number_parameters >= 4) ? parameters[3] : 0.5;
1056 result=2.0*PerceptibleReciprocal(width)*(QuantumScale*pixel-center);
1057 if (result <= -1.0)
1058 result=bias-range/2.0;
1059 else
1060 if (result >= 1.0)
1061 result=bias+range/2.0;
1062 else
1063 result=(MagickRealType) (range/MagickPI*asin((double) result)+bias);
1064 result*=QuantumRange;
1065 break;
1066 }
1067 case ArctanFunction:
1068 {
1069 /* Arctan Function
1070 * Parameters: Slope, Center, Range, Bias
1071 */
1072 double slope,range,center,bias;
1073 slope = ( number_parameters >= 1 ) ? parameters[0] : 1.0;
1074 center = ( number_parameters >= 2 ) ? parameters[1] : 0.5;
1075 range = ( number_parameters >= 3 ) ? parameters[2] : 1.0;
1076 bias = ( number_parameters >= 4 ) ? parameters[3] : 0.5;
1077 result=(MagickRealType) (MagickPI*slope*(QuantumScale*pixel-center));
1078 result=(MagickRealType) (QuantumRange*(range/MagickPI*atan((double)
1079 result) + bias ) );
1080 break;
1081 }
1082 case UndefinedFunction:
1083 break;
1084 }
1085 return(ClampToQuantum(result));
1086 }
1087
FunctionImage(Image * image,const MagickFunction function,const size_t number_parameters,const double * parameters,ExceptionInfo * exception)1088 MagickExport MagickBooleanType FunctionImage(Image *image,
1089 const MagickFunction function,const size_t number_parameters,
1090 const double *parameters,ExceptionInfo *exception)
1091 {
1092 MagickBooleanType
1093 status;
1094
1095 status=FunctionImageChannel(image,CompositeChannels,function,
1096 number_parameters,parameters,exception);
1097 return(status);
1098 }
1099
FunctionImageChannel(Image * image,const ChannelType channel,const MagickFunction function,const size_t number_parameters,const double * parameters,ExceptionInfo * exception)1100 MagickExport MagickBooleanType FunctionImageChannel(Image *image,
1101 const ChannelType channel,const MagickFunction function,
1102 const size_t number_parameters,const double *parameters,
1103 ExceptionInfo *exception)
1104 {
1105 #define FunctionImageTag "Function/Image "
1106
1107 CacheView
1108 *image_view;
1109
1110 MagickBooleanType
1111 status;
1112
1113 MagickOffsetType
1114 progress;
1115
1116 ssize_t
1117 y;
1118
1119 assert(image != (Image *) NULL);
1120 assert(image->signature == MagickCoreSignature);
1121 if (image->debug != MagickFalse)
1122 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1123 assert(exception != (ExceptionInfo *) NULL);
1124 assert(exception->signature == MagickCoreSignature);
1125 if (SetImageStorageClass(image,DirectClass) == MagickFalse)
1126 {
1127 InheritException(exception,&image->exception);
1128 return(MagickFalse);
1129 }
1130 #if defined(MAGICKCORE_OPENCL_SUPPORT)
1131 status=AccelerateFunctionImage(image,channel,function,number_parameters,
1132 parameters,exception);
1133 if (status != MagickFalse)
1134 return(status);
1135 #endif
1136 status=MagickTrue;
1137 progress=0;
1138 image_view=AcquireAuthenticCacheView(image,exception);
1139 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1140 #pragma omp parallel for schedule(static) shared(progress,status) \
1141 magick_number_threads(image,image,image->rows,1)
1142 #endif
1143 for (y=0; y < (ssize_t) image->rows; y++)
1144 {
1145 IndexPacket
1146 *magick_restrict indexes;
1147
1148 ssize_t
1149 x;
1150
1151 PixelPacket
1152 *magick_restrict q;
1153
1154 if (status == MagickFalse)
1155 continue;
1156 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
1157 if (q == (PixelPacket *) NULL)
1158 {
1159 status=MagickFalse;
1160 continue;
1161 }
1162 indexes=GetCacheViewAuthenticIndexQueue(image_view);
1163 for (x=0; x < (ssize_t) image->columns; x++)
1164 {
1165 if ((channel & RedChannel) != 0)
1166 SetPixelRed(q,ApplyFunction(GetPixelRed(q),function,
1167 number_parameters,parameters,exception));
1168 if ((channel & GreenChannel) != 0)
1169 SetPixelGreen(q,ApplyFunction(GetPixelGreen(q),function,
1170 number_parameters,parameters,exception));
1171 if ((channel & BlueChannel) != 0)
1172 SetPixelBlue(q,ApplyFunction(GetPixelBlue(q),function,
1173 number_parameters,parameters,exception));
1174 if ((channel & OpacityChannel) != 0)
1175 {
1176 if (image->matte == MagickFalse)
1177 SetPixelOpacity(q,ApplyFunction(GetPixelOpacity(q),function,
1178 number_parameters,parameters,exception));
1179 else
1180 SetPixelAlpha(q,ApplyFunction((Quantum) GetPixelAlpha(q),function,
1181 number_parameters,parameters,exception));
1182 }
1183 if (((channel & IndexChannel) != 0) && (indexes != (IndexPacket *) NULL))
1184 SetPixelIndex(indexes+x,ApplyFunction(GetPixelIndex(indexes+x),function,
1185 number_parameters,parameters,exception));
1186 q++;
1187 }
1188 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1189 status=MagickFalse;
1190 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1191 {
1192 MagickBooleanType
1193 proceed;
1194
1195 proceed=SetImageProgress(image,FunctionImageTag,progress++,image->rows);
1196 if (proceed == MagickFalse)
1197 status=MagickFalse;
1198 }
1199 }
1200 image_view=DestroyCacheView(image_view);
1201 return(status);
1202 }
1203
1204 /*
1205 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1206 % %
1207 % %
1208 % %
1209 % G e t I m a g e C h a n n e l E n t r o p y %
1210 % %
1211 % %
1212 % %
1213 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1214 %
1215 % GetImageChannelEntropy() returns the entropy of one or more image channels.
1216 %
1217 % The format of the GetImageChannelEntropy method is:
1218 %
1219 % MagickBooleanType GetImageChannelEntropy(const Image *image,
1220 % const ChannelType channel,double *entropy,ExceptionInfo *exception)
1221 %
1222 % A description of each parameter follows:
1223 %
1224 % o image: the image.
1225 %
1226 % o channel: the channel.
1227 %
1228 % o entropy: the average entropy of the selected channels.
1229 %
1230 % o exception: return any errors or warnings in this structure.
1231 %
1232 */
1233
GetImageEntropy(const Image * image,double * entropy,ExceptionInfo * exception)1234 MagickExport MagickBooleanType GetImageEntropy(const Image *image,
1235 double *entropy,ExceptionInfo *exception)
1236 {
1237 MagickBooleanType
1238 status;
1239
1240 status=GetImageChannelEntropy(image,CompositeChannels,entropy,exception);
1241 return(status);
1242 }
1243
GetImageChannelEntropy(const Image * image,const ChannelType channel,double * entropy,ExceptionInfo * exception)1244 MagickExport MagickBooleanType GetImageChannelEntropy(const Image *image,
1245 const ChannelType channel,double *entropy,ExceptionInfo *exception)
1246 {
1247 ChannelStatistics
1248 *channel_statistics;
1249
1250 size_t
1251 channels;
1252
1253 assert(image != (Image *) NULL);
1254 assert(image->signature == MagickCoreSignature);
1255 if (image->debug != MagickFalse)
1256 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1257 channel_statistics=GetImageChannelStatistics(image,exception);
1258 if (channel_statistics == (ChannelStatistics *) NULL)
1259 return(MagickFalse);
1260 channels=0;
1261 channel_statistics[CompositeChannels].entropy=0.0;
1262 if ((channel & RedChannel) != 0)
1263 {
1264 channel_statistics[CompositeChannels].entropy+=
1265 channel_statistics[RedChannel].entropy;
1266 channels++;
1267 }
1268 if ((channel & GreenChannel) != 0)
1269 {
1270 channel_statistics[CompositeChannels].entropy+=
1271 channel_statistics[GreenChannel].entropy;
1272 channels++;
1273 }
1274 if ((channel & BlueChannel) != 0)
1275 {
1276 channel_statistics[CompositeChannels].entropy+=
1277 channel_statistics[BlueChannel].entropy;
1278 channels++;
1279 }
1280 if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
1281 {
1282 channel_statistics[CompositeChannels].entropy+=
1283 channel_statistics[OpacityChannel].entropy;
1284 channels++;
1285 }
1286 if (((channel & IndexChannel) != 0) &&
1287 (image->colorspace == CMYKColorspace))
1288 {
1289 channel_statistics[CompositeChannels].entropy+=
1290 channel_statistics[BlackChannel].entropy;
1291 channels++;
1292 }
1293 channel_statistics[CompositeChannels].entropy/=channels;
1294 *entropy=channel_statistics[CompositeChannels].entropy;
1295 channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1296 channel_statistics);
1297 return(MagickTrue);
1298 }
1299
1300 /*
1301 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1302 % %
1303 % %
1304 % %
1305 + G e t I m a g e C h a n n e l E x t r e m a %
1306 % %
1307 % %
1308 % %
1309 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1310 %
1311 % GetImageChannelExtrema() returns the extrema of one or more image channels.
1312 %
1313 % The format of the GetImageChannelExtrema method is:
1314 %
1315 % MagickBooleanType GetImageChannelExtrema(const Image *image,
1316 % const ChannelType channel,size_t *minima,size_t *maxima,
1317 % ExceptionInfo *exception)
1318 %
1319 % A description of each parameter follows:
1320 %
1321 % o image: the image.
1322 %
1323 % o channel: the channel.
1324 %
1325 % o minima: the minimum value in the channel.
1326 %
1327 % o maxima: the maximum value in the channel.
1328 %
1329 % o exception: return any errors or warnings in this structure.
1330 %
1331 */
1332
GetImageExtrema(const Image * image,size_t * minima,size_t * maxima,ExceptionInfo * exception)1333 MagickExport MagickBooleanType GetImageExtrema(const Image *image,
1334 size_t *minima,size_t *maxima,ExceptionInfo *exception)
1335 {
1336 MagickBooleanType
1337 status;
1338
1339 status=GetImageChannelExtrema(image,CompositeChannels,minima,maxima,
1340 exception);
1341 return(status);
1342 }
1343
GetImageChannelExtrema(const Image * image,const ChannelType channel,size_t * minima,size_t * maxima,ExceptionInfo * exception)1344 MagickExport MagickBooleanType GetImageChannelExtrema(const Image *image,
1345 const ChannelType channel,size_t *minima,size_t *maxima,
1346 ExceptionInfo *exception)
1347 {
1348 double
1349 max,
1350 min;
1351
1352 MagickBooleanType
1353 status;
1354
1355 assert(image != (Image *) NULL);
1356 assert(image->signature == MagickCoreSignature);
1357 if (image->debug != MagickFalse)
1358 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1359 status=GetImageChannelRange(image,channel,&min,&max,exception);
1360 *minima=(size_t) ceil(min-0.5);
1361 *maxima=(size_t) floor(max+0.5);
1362 return(status);
1363 }
1364
1365 /*
1366 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1367 % %
1368 % %
1369 % %
1370 % G e t I m a g e C h a n n e l K u r t o s i s %
1371 % %
1372 % %
1373 % %
1374 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1375 %
1376 % GetImageChannelKurtosis() returns the kurtosis and skewness of one or more
1377 % image channels.
1378 %
1379 % The format of the GetImageChannelKurtosis method is:
1380 %
1381 % MagickBooleanType GetImageChannelKurtosis(const Image *image,
1382 % const ChannelType channel,double *kurtosis,double *skewness,
1383 % ExceptionInfo *exception)
1384 %
1385 % A description of each parameter follows:
1386 %
1387 % o image: the image.
1388 %
1389 % o channel: the channel.
1390 %
1391 % o kurtosis: the kurtosis of the channel.
1392 %
1393 % o skewness: the skewness of the channel.
1394 %
1395 % o exception: return any errors or warnings in this structure.
1396 %
1397 */
1398
GetImageKurtosis(const Image * image,double * kurtosis,double * skewness,ExceptionInfo * exception)1399 MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
1400 double *kurtosis,double *skewness,ExceptionInfo *exception)
1401 {
1402 MagickBooleanType
1403 status;
1404
1405 status=GetImageChannelKurtosis(image,CompositeChannels,kurtosis,skewness,
1406 exception);
1407 return(status);
1408 }
1409
GetImageChannelKurtosis(const Image * image,const ChannelType channel,double * kurtosis,double * skewness,ExceptionInfo * exception)1410 MagickExport MagickBooleanType GetImageChannelKurtosis(const Image *image,
1411 const ChannelType channel,double *kurtosis,double *skewness,
1412 ExceptionInfo *exception)
1413 {
1414 double
1415 area,
1416 mean,
1417 standard_deviation,
1418 sum_squares,
1419 sum_cubes,
1420 sum_fourth_power;
1421
1422 ssize_t
1423 y;
1424
1425 assert(image != (Image *) NULL);
1426 assert(image->signature == MagickCoreSignature);
1427 if (image->debug != MagickFalse)
1428 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1429 *kurtosis=0.0;
1430 *skewness=0.0;
1431 area=0.0;
1432 mean=0.0;
1433 standard_deviation=0.0;
1434 sum_squares=0.0;
1435 sum_cubes=0.0;
1436 sum_fourth_power=0.0;
1437 for (y=0; y < (ssize_t) image->rows; y++)
1438 {
1439 const IndexPacket
1440 *magick_restrict indexes;
1441
1442 const PixelPacket
1443 *magick_restrict p;
1444
1445 ssize_t
1446 x;
1447
1448 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1449 if (p == (const PixelPacket *) NULL)
1450 break;
1451 indexes=GetVirtualIndexQueue(image);
1452 for (x=0; x < (ssize_t) image->columns; x++)
1453 {
1454 if ((channel & RedChannel) != 0)
1455 {
1456 mean+=GetPixelRed(p);
1457 sum_squares+=(double) GetPixelRed(p)*GetPixelRed(p);
1458 sum_cubes+=(double) GetPixelRed(p)*GetPixelRed(p)*GetPixelRed(p);
1459 sum_fourth_power+=(double) GetPixelRed(p)*GetPixelRed(p)*
1460 GetPixelRed(p)*GetPixelRed(p);
1461 area++;
1462 }
1463 if ((channel & GreenChannel) != 0)
1464 {
1465 mean+=GetPixelGreen(p);
1466 sum_squares+=(double) GetPixelGreen(p)*GetPixelGreen(p);
1467 sum_cubes+=(double) GetPixelGreen(p)*GetPixelGreen(p)*
1468 GetPixelGreen(p);
1469 sum_fourth_power+=(double) GetPixelGreen(p)*GetPixelGreen(p)*
1470 GetPixelGreen(p)*GetPixelGreen(p);
1471 area++;
1472 }
1473 if ((channel & BlueChannel) != 0)
1474 {
1475 mean+=GetPixelBlue(p);
1476 sum_squares+=(double) GetPixelBlue(p)*GetPixelBlue(p);
1477 sum_cubes+=(double) GetPixelBlue(p)*GetPixelBlue(p)*GetPixelBlue(p);
1478 sum_fourth_power+=(double) GetPixelBlue(p)*GetPixelBlue(p)*
1479 GetPixelBlue(p)*GetPixelBlue(p);
1480 area++;
1481 }
1482 if ((channel & OpacityChannel) != 0)
1483 {
1484 mean+=GetPixelAlpha(p);
1485 sum_squares+=(double) GetPixelOpacity(p)*GetPixelAlpha(p);
1486 sum_cubes+=(double) GetPixelOpacity(p)*GetPixelAlpha(p)*
1487 GetPixelAlpha(p);
1488 sum_fourth_power+=(double) GetPixelAlpha(p)*GetPixelAlpha(p)*
1489 GetPixelAlpha(p)*GetPixelAlpha(p);
1490 area++;
1491 }
1492 if (((channel & IndexChannel) != 0) &&
1493 (image->colorspace == CMYKColorspace))
1494 {
1495 double
1496 index;
1497
1498 index=(double) GetPixelIndex(indexes+x);
1499 mean+=index;
1500 sum_squares+=index*index;
1501 sum_cubes+=index*index*index;
1502 sum_fourth_power+=index*index*index*index;
1503 area++;
1504 }
1505 p++;
1506 }
1507 }
1508 if (y < (ssize_t) image->rows)
1509 return(MagickFalse);
1510 if (area != 0.0)
1511 {
1512 mean/=area;
1513 sum_squares/=area;
1514 sum_cubes/=area;
1515 sum_fourth_power/=area;
1516 }
1517 standard_deviation=sqrt(sum_squares-(mean*mean));
1518 if (standard_deviation != 0.0)
1519 {
1520 *kurtosis=sum_fourth_power-4.0*mean*sum_cubes+6.0*mean*mean*sum_squares-
1521 3.0*mean*mean*mean*mean;
1522 *kurtosis/=standard_deviation*standard_deviation*standard_deviation*
1523 standard_deviation;
1524 *kurtosis-=3.0;
1525 *skewness=sum_cubes-3.0*mean*sum_squares+2.0*mean*mean*mean;
1526 *skewness/=standard_deviation*standard_deviation*standard_deviation;
1527 }
1528 return(y == (ssize_t) image->rows ? MagickTrue : MagickFalse);
1529 }
1530
1531 /*
1532 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1533 % %
1534 % %
1535 % %
1536 % G e t I m a g e C h a n n e l M e a n %
1537 % %
1538 % %
1539 % %
1540 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1541 %
1542 % GetImageChannelMean() returns the mean and standard deviation of one or more
1543 % image channels.
1544 %
1545 % The format of the GetImageChannelMean method is:
1546 %
1547 % MagickBooleanType GetImageChannelMean(const Image *image,
1548 % const ChannelType channel,double *mean,double *standard_deviation,
1549 % ExceptionInfo *exception)
1550 %
1551 % A description of each parameter follows:
1552 %
1553 % o image: the image.
1554 %
1555 % o channel: the channel.
1556 %
1557 % o mean: the average value in the channel.
1558 %
1559 % o standard_deviation: the standard deviation of the channel.
1560 %
1561 % o exception: return any errors or warnings in this structure.
1562 %
1563 */
1564
GetImageMean(const Image * image,double * mean,double * standard_deviation,ExceptionInfo * exception)1565 MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean,
1566 double *standard_deviation,ExceptionInfo *exception)
1567 {
1568 MagickBooleanType
1569 status;
1570
1571 status=GetImageChannelMean(image,CompositeChannels,mean,standard_deviation,
1572 exception);
1573 return(status);
1574 }
1575
GetImageChannelMean(const Image * image,const ChannelType channel,double * mean,double * standard_deviation,ExceptionInfo * exception)1576 MagickExport MagickBooleanType GetImageChannelMean(const Image *image,
1577 const ChannelType channel,double *mean,double *standard_deviation,
1578 ExceptionInfo *exception)
1579 {
1580 ChannelStatistics
1581 *channel_statistics;
1582
1583 size_t
1584 channels;
1585
1586 assert(image != (Image *) NULL);
1587 assert(image->signature == MagickCoreSignature);
1588 if (image->debug != MagickFalse)
1589 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1590 channel_statistics=GetImageChannelStatistics(image,exception);
1591 if (channel_statistics == (ChannelStatistics *) NULL)
1592 return(MagickFalse);
1593 channels=0;
1594 channel_statistics[CompositeChannels].mean=0.0;
1595 channel_statistics[CompositeChannels].standard_deviation=0.0;
1596 if ((channel & RedChannel) != 0)
1597 {
1598 channel_statistics[CompositeChannels].mean+=
1599 channel_statistics[RedChannel].mean;
1600 channel_statistics[CompositeChannels].standard_deviation+=
1601 channel_statistics[RedChannel].standard_deviation;
1602 channels++;
1603 }
1604 if ((channel & GreenChannel) != 0)
1605 {
1606 channel_statistics[CompositeChannels].mean+=
1607 channel_statistics[GreenChannel].mean;
1608 channel_statistics[CompositeChannels].standard_deviation+=
1609 channel_statistics[GreenChannel].standard_deviation;
1610 channels++;
1611 }
1612 if ((channel & BlueChannel) != 0)
1613 {
1614 channel_statistics[CompositeChannels].mean+=
1615 channel_statistics[BlueChannel].mean;
1616 channel_statistics[CompositeChannels].standard_deviation+=
1617 channel_statistics[BlueChannel].standard_deviation;
1618 channels++;
1619 }
1620 if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
1621 {
1622 channel_statistics[CompositeChannels].mean+=
1623 channel_statistics[OpacityChannel].mean;
1624 channel_statistics[CompositeChannels].standard_deviation+=
1625 channel_statistics[OpacityChannel].standard_deviation;
1626 channels++;
1627 }
1628 if (((channel & IndexChannel) != 0) && (image->colorspace == CMYKColorspace))
1629 {
1630 channel_statistics[CompositeChannels].mean+=
1631 channel_statistics[BlackChannel].mean;
1632 channel_statistics[CompositeChannels].standard_deviation+=
1633 channel_statistics[CompositeChannels].standard_deviation;
1634 channels++;
1635 }
1636 channel_statistics[CompositeChannels].mean/=channels;
1637 channel_statistics[CompositeChannels].standard_deviation/=channels;
1638 *mean=channel_statistics[CompositeChannels].mean;
1639 *standard_deviation=channel_statistics[CompositeChannels].standard_deviation;
1640 channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1641 channel_statistics);
1642 return(MagickTrue);
1643 }
1644
1645 /*
1646 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1647 % %
1648 % %
1649 % %
1650 % G e t I m a g e C h a n n e l M o m e n t s %
1651 % %
1652 % %
1653 % %
1654 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1655 %
1656 % GetImageChannelMoments() returns the normalized moments of one or more image
1657 % channels.
1658 %
1659 % The format of the GetImageChannelMoments method is:
1660 %
1661 % ChannelMoments *GetImageChannelMoments(const Image *image,
1662 % ExceptionInfo *exception)
1663 %
1664 % A description of each parameter follows:
1665 %
1666 % o image: the image.
1667 %
1668 % o exception: return any errors or warnings in this structure.
1669 %
1670 */
GetImageChannelMoments(const Image * image,ExceptionInfo * exception)1671 MagickExport ChannelMoments *GetImageChannelMoments(const Image *image,
1672 ExceptionInfo *exception)
1673 {
1674 #define MaxNumberImageMoments 8
1675
1676 ChannelMoments
1677 *channel_moments;
1678
1679 double
1680 M00[CompositeChannels+1],
1681 M01[CompositeChannels+1],
1682 M02[CompositeChannels+1],
1683 M03[CompositeChannels+1],
1684 M10[CompositeChannels+1],
1685 M11[CompositeChannels+1],
1686 M12[CompositeChannels+1],
1687 M20[CompositeChannels+1],
1688 M21[CompositeChannels+1],
1689 M22[CompositeChannels+1],
1690 M30[CompositeChannels+1];
1691
1692 MagickPixelPacket
1693 pixel;
1694
1695 PointInfo
1696 centroid[CompositeChannels+1];
1697
1698 ssize_t
1699 channel,
1700 channels,
1701 y;
1702
1703 size_t
1704 length;
1705
1706 assert(image != (Image *) NULL);
1707 assert(image->signature == MagickCoreSignature);
1708 if (image->debug != MagickFalse)
1709 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1710 length=CompositeChannels+1UL;
1711 channel_moments=(ChannelMoments *) AcquireQuantumMemory(length,
1712 sizeof(*channel_moments));
1713 if (channel_moments == (ChannelMoments *) NULL)
1714 return(channel_moments);
1715 (void) memset(channel_moments,0,length*sizeof(*channel_moments));
1716 (void) memset(centroid,0,sizeof(centroid));
1717 (void) memset(M00,0,sizeof(M00));
1718 (void) memset(M01,0,sizeof(M01));
1719 (void) memset(M02,0,sizeof(M02));
1720 (void) memset(M03,0,sizeof(M03));
1721 (void) memset(M10,0,sizeof(M10));
1722 (void) memset(M11,0,sizeof(M11));
1723 (void) memset(M12,0,sizeof(M12));
1724 (void) memset(M20,0,sizeof(M20));
1725 (void) memset(M21,0,sizeof(M21));
1726 (void) memset(M22,0,sizeof(M22));
1727 (void) memset(M30,0,sizeof(M30));
1728 GetMagickPixelPacket(image,&pixel);
1729 for (y=0; y < (ssize_t) image->rows; y++)
1730 {
1731 const IndexPacket
1732 *magick_restrict indexes;
1733
1734 const PixelPacket
1735 *magick_restrict p;
1736
1737 ssize_t
1738 x;
1739
1740 /*
1741 Compute center of mass (centroid).
1742 */
1743 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1744 if (p == (const PixelPacket *) NULL)
1745 break;
1746 indexes=GetVirtualIndexQueue(image);
1747 for (x=0; x < (ssize_t) image->columns; x++)
1748 {
1749 SetMagickPixelPacket(image,p,indexes+x,&pixel);
1750 M00[RedChannel]+=QuantumScale*pixel.red;
1751 M10[RedChannel]+=x*QuantumScale*pixel.red;
1752 M01[RedChannel]+=y*QuantumScale*pixel.red;
1753 M00[GreenChannel]+=QuantumScale*pixel.green;
1754 M10[GreenChannel]+=x*QuantumScale*pixel.green;
1755 M01[GreenChannel]+=y*QuantumScale*pixel.green;
1756 M00[BlueChannel]+=QuantumScale*pixel.blue;
1757 M10[BlueChannel]+=x*QuantumScale*pixel.blue;
1758 M01[BlueChannel]+=y*QuantumScale*pixel.blue;
1759 if (image->matte != MagickFalse)
1760 {
1761 M00[OpacityChannel]+=QuantumScale*pixel.opacity;
1762 M10[OpacityChannel]+=x*QuantumScale*pixel.opacity;
1763 M01[OpacityChannel]+=y*QuantumScale*pixel.opacity;
1764 }
1765 if (image->colorspace == CMYKColorspace)
1766 {
1767 M00[IndexChannel]+=QuantumScale*pixel.index;
1768 M10[IndexChannel]+=x*QuantumScale*pixel.index;
1769 M01[IndexChannel]+=y*QuantumScale*pixel.index;
1770 }
1771 p++;
1772 }
1773 }
1774 for (channel=0; channel <= CompositeChannels; channel++)
1775 {
1776 /*
1777 Compute center of mass (centroid).
1778 */
1779 if (M00[channel] < MagickEpsilon)
1780 {
1781 M00[channel]+=MagickEpsilon;
1782 centroid[channel].x=(double) image->columns/2.0;
1783 centroid[channel].y=(double) image->rows/2.0;
1784 continue;
1785 }
1786 M00[channel]+=MagickEpsilon;
1787 centroid[channel].x=M10[channel]/M00[channel];
1788 centroid[channel].y=M01[channel]/M00[channel];
1789 }
1790 for (y=0; y < (ssize_t) image->rows; y++)
1791 {
1792 const IndexPacket
1793 *magick_restrict indexes;
1794
1795 const PixelPacket
1796 *magick_restrict p;
1797
1798 ssize_t
1799 x;
1800
1801 /*
1802 Compute the image moments.
1803 */
1804 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1805 if (p == (const PixelPacket *) NULL)
1806 break;
1807 indexes=GetVirtualIndexQueue(image);
1808 for (x=0; x < (ssize_t) image->columns; x++)
1809 {
1810 SetMagickPixelPacket(image,p,indexes+x,&pixel);
1811 M11[RedChannel]+=(x-centroid[RedChannel].x)*(y-
1812 centroid[RedChannel].y)*QuantumScale*pixel.red;
1813 M20[RedChannel]+=(x-centroid[RedChannel].x)*(x-
1814 centroid[RedChannel].x)*QuantumScale*pixel.red;
1815 M02[RedChannel]+=(y-centroid[RedChannel].y)*(y-
1816 centroid[RedChannel].y)*QuantumScale*pixel.red;
1817 M21[RedChannel]+=(x-centroid[RedChannel].x)*(x-
1818 centroid[RedChannel].x)*(y-centroid[RedChannel].y)*QuantumScale*
1819 pixel.red;
1820 M12[RedChannel]+=(x-centroid[RedChannel].x)*(y-
1821 centroid[RedChannel].y)*(y-centroid[RedChannel].y)*QuantumScale*
1822 pixel.red;
1823 M22[RedChannel]+=(x-centroid[RedChannel].x)*(x-
1824 centroid[RedChannel].x)*(y-centroid[RedChannel].y)*(y-
1825 centroid[RedChannel].y)*QuantumScale*pixel.red;
1826 M30[RedChannel]+=(x-centroid[RedChannel].x)*(x-
1827 centroid[RedChannel].x)*(x-centroid[RedChannel].x)*QuantumScale*
1828 pixel.red;
1829 M03[RedChannel]+=(y-centroid[RedChannel].y)*(y-
1830 centroid[RedChannel].y)*(y-centroid[RedChannel].y)*QuantumScale*
1831 pixel.red;
1832 M11[GreenChannel]+=(x-centroid[GreenChannel].x)*(y-
1833 centroid[GreenChannel].y)*QuantumScale*pixel.green;
1834 M20[GreenChannel]+=(x-centroid[GreenChannel].x)*(x-
1835 centroid[GreenChannel].x)*QuantumScale*pixel.green;
1836 M02[GreenChannel]+=(y-centroid[GreenChannel].y)*(y-
1837 centroid[GreenChannel].y)*QuantumScale*pixel.green;
1838 M21[GreenChannel]+=(x-centroid[GreenChannel].x)*(x-
1839 centroid[GreenChannel].x)*(y-centroid[GreenChannel].y)*QuantumScale*
1840 pixel.green;
1841 M12[GreenChannel]+=(x-centroid[GreenChannel].x)*(y-
1842 centroid[GreenChannel].y)*(y-centroid[GreenChannel].y)*QuantumScale*
1843 pixel.green;
1844 M22[GreenChannel]+=(x-centroid[GreenChannel].x)*(x-
1845 centroid[GreenChannel].x)*(y-centroid[GreenChannel].y)*(y-
1846 centroid[GreenChannel].y)*QuantumScale*pixel.green;
1847 M30[GreenChannel]+=(x-centroid[GreenChannel].x)*(x-
1848 centroid[GreenChannel].x)*(x-centroid[GreenChannel].x)*QuantumScale*
1849 pixel.green;
1850 M03[GreenChannel]+=(y-centroid[GreenChannel].y)*(y-
1851 centroid[GreenChannel].y)*(y-centroid[GreenChannel].y)*QuantumScale*
1852 pixel.green;
1853 M11[BlueChannel]+=(x-centroid[BlueChannel].x)*(y-
1854 centroid[BlueChannel].y)*QuantumScale*pixel.blue;
1855 M20[BlueChannel]+=(x-centroid[BlueChannel].x)*(x-
1856 centroid[BlueChannel].x)*QuantumScale*pixel.blue;
1857 M02[BlueChannel]+=(y-centroid[BlueChannel].y)*(y-
1858 centroid[BlueChannel].y)*QuantumScale*pixel.blue;
1859 M21[BlueChannel]+=(x-centroid[BlueChannel].x)*(x-
1860 centroid[BlueChannel].x)*(y-centroid[BlueChannel].y)*QuantumScale*
1861 pixel.blue;
1862 M12[BlueChannel]+=(x-centroid[BlueChannel].x)*(y-
1863 centroid[BlueChannel].y)*(y-centroid[BlueChannel].y)*QuantumScale*
1864 pixel.blue;
1865 M22[BlueChannel]+=(x-centroid[BlueChannel].x)*(x-
1866 centroid[BlueChannel].x)*(y-centroid[BlueChannel].y)*(y-
1867 centroid[BlueChannel].y)*QuantumScale*pixel.blue;
1868 M30[BlueChannel]+=(x-centroid[BlueChannel].x)*(x-
1869 centroid[BlueChannel].x)*(x-centroid[BlueChannel].x)*QuantumScale*
1870 pixel.blue;
1871 M03[BlueChannel]+=(y-centroid[BlueChannel].y)*(y-
1872 centroid[BlueChannel].y)*(y-centroid[BlueChannel].y)*QuantumScale*
1873 pixel.blue;
1874 if (image->matte != MagickFalse)
1875 {
1876 M11[OpacityChannel]+=(x-centroid[OpacityChannel].x)*(y-
1877 centroid[OpacityChannel].y)*QuantumScale*pixel.opacity;
1878 M20[OpacityChannel]+=(x-centroid[OpacityChannel].x)*(x-
1879 centroid[OpacityChannel].x)*QuantumScale*pixel.opacity;
1880 M02[OpacityChannel]+=(y-centroid[OpacityChannel].y)*(y-
1881 centroid[OpacityChannel].y)*QuantumScale*pixel.opacity;
1882 M21[OpacityChannel]+=(x-centroid[OpacityChannel].x)*(x-
1883 centroid[OpacityChannel].x)*(y-centroid[OpacityChannel].y)*
1884 QuantumScale*pixel.opacity;
1885 M12[OpacityChannel]+=(x-centroid[OpacityChannel].x)*(y-
1886 centroid[OpacityChannel].y)*(y-centroid[OpacityChannel].y)*
1887 QuantumScale*pixel.opacity;
1888 M22[OpacityChannel]+=(x-centroid[OpacityChannel].x)*(x-
1889 centroid[OpacityChannel].x)*(y-centroid[OpacityChannel].y)*(y-
1890 centroid[OpacityChannel].y)*QuantumScale*pixel.opacity;
1891 M30[OpacityChannel]+=(x-centroid[OpacityChannel].x)*(x-
1892 centroid[OpacityChannel].x)*(x-centroid[OpacityChannel].x)*
1893 QuantumScale*pixel.opacity;
1894 M03[OpacityChannel]+=(y-centroid[OpacityChannel].y)*(y-
1895 centroid[OpacityChannel].y)*(y-centroid[OpacityChannel].y)*
1896 QuantumScale*pixel.opacity;
1897 }
1898 if (image->colorspace == CMYKColorspace)
1899 {
1900 M11[IndexChannel]+=(x-centroid[IndexChannel].x)*(y-
1901 centroid[IndexChannel].y)*QuantumScale*pixel.index;
1902 M20[IndexChannel]+=(x-centroid[IndexChannel].x)*(x-
1903 centroid[IndexChannel].x)*QuantumScale*pixel.index;
1904 M02[IndexChannel]+=(y-centroid[IndexChannel].y)*(y-
1905 centroid[IndexChannel].y)*QuantumScale*pixel.index;
1906 M21[IndexChannel]+=(x-centroid[IndexChannel].x)*(x-
1907 centroid[IndexChannel].x)*(y-centroid[IndexChannel].y)*
1908 QuantumScale*pixel.index;
1909 M12[IndexChannel]+=(x-centroid[IndexChannel].x)*(y-
1910 centroid[IndexChannel].y)*(y-centroid[IndexChannel].y)*
1911 QuantumScale*pixel.index;
1912 M22[IndexChannel]+=(x-centroid[IndexChannel].x)*(x-
1913 centroid[IndexChannel].x)*(y-centroid[IndexChannel].y)*(y-
1914 centroid[IndexChannel].y)*QuantumScale*pixel.index;
1915 M30[IndexChannel]+=(x-centroid[IndexChannel].x)*(x-
1916 centroid[IndexChannel].x)*(x-centroid[IndexChannel].x)*
1917 QuantumScale*pixel.index;
1918 M03[IndexChannel]+=(y-centroid[IndexChannel].y)*(y-
1919 centroid[IndexChannel].y)*(y-centroid[IndexChannel].y)*
1920 QuantumScale*pixel.index;
1921 }
1922 p++;
1923 }
1924 }
1925 channels=3;
1926 M00[CompositeChannels]+=(M00[RedChannel]+M00[GreenChannel]+M00[BlueChannel]);
1927 M01[CompositeChannels]+=(M01[RedChannel]+M01[GreenChannel]+M01[BlueChannel]);
1928 M02[CompositeChannels]+=(M02[RedChannel]+M02[GreenChannel]+M02[BlueChannel]);
1929 M03[CompositeChannels]+=(M03[RedChannel]+M03[GreenChannel]+M03[BlueChannel]);
1930 M10[CompositeChannels]+=(M10[RedChannel]+M10[GreenChannel]+M10[BlueChannel]);
1931 M11[CompositeChannels]+=(M11[RedChannel]+M11[GreenChannel]+M11[BlueChannel]);
1932 M12[CompositeChannels]+=(M12[RedChannel]+M12[GreenChannel]+M12[BlueChannel]);
1933 M20[CompositeChannels]+=(M20[RedChannel]+M20[GreenChannel]+M20[BlueChannel]);
1934 M21[CompositeChannels]+=(M21[RedChannel]+M21[GreenChannel]+M21[BlueChannel]);
1935 M22[CompositeChannels]+=(M22[RedChannel]+M22[GreenChannel]+M22[BlueChannel]);
1936 M30[CompositeChannels]+=(M30[RedChannel]+M30[GreenChannel]+M30[BlueChannel]);
1937 if (image->matte != MagickFalse)
1938 {
1939 channels+=1;
1940 M00[CompositeChannels]+=M00[OpacityChannel];
1941 M01[CompositeChannels]+=M01[OpacityChannel];
1942 M02[CompositeChannels]+=M02[OpacityChannel];
1943 M03[CompositeChannels]+=M03[OpacityChannel];
1944 M10[CompositeChannels]+=M10[OpacityChannel];
1945 M11[CompositeChannels]+=M11[OpacityChannel];
1946 M12[CompositeChannels]+=M12[OpacityChannel];
1947 M20[CompositeChannels]+=M20[OpacityChannel];
1948 M21[CompositeChannels]+=M21[OpacityChannel];
1949 M22[CompositeChannels]+=M22[OpacityChannel];
1950 M30[CompositeChannels]+=M30[OpacityChannel];
1951 }
1952 if (image->colorspace == CMYKColorspace)
1953 {
1954 channels+=1;
1955 M00[CompositeChannels]+=M00[IndexChannel];
1956 M01[CompositeChannels]+=M01[IndexChannel];
1957 M02[CompositeChannels]+=M02[IndexChannel];
1958 M03[CompositeChannels]+=M03[IndexChannel];
1959 M10[CompositeChannels]+=M10[IndexChannel];
1960 M11[CompositeChannels]+=M11[IndexChannel];
1961 M12[CompositeChannels]+=M12[IndexChannel];
1962 M20[CompositeChannels]+=M20[IndexChannel];
1963 M21[CompositeChannels]+=M21[IndexChannel];
1964 M22[CompositeChannels]+=M22[IndexChannel];
1965 M30[CompositeChannels]+=M30[IndexChannel];
1966 }
1967 M00[CompositeChannels]/=(double) channels;
1968 M01[CompositeChannels]/=(double) channels;
1969 M02[CompositeChannels]/=(double) channels;
1970 M03[CompositeChannels]/=(double) channels;
1971 M10[CompositeChannels]/=(double) channels;
1972 M11[CompositeChannels]/=(double) channels;
1973 M12[CompositeChannels]/=(double) channels;
1974 M20[CompositeChannels]/=(double) channels;
1975 M21[CompositeChannels]/=(double) channels;
1976 M22[CompositeChannels]/=(double) channels;
1977 M30[CompositeChannels]/=(double) channels;
1978 for (channel=0; channel <= CompositeChannels; channel++)
1979 {
1980 /*
1981 Compute elliptical angle, major and minor axes, eccentricity, & intensity.
1982 */
1983 channel_moments[channel].centroid=centroid[channel];
1984 channel_moments[channel].ellipse_axis.x=sqrt((2.0*
1985 PerceptibleReciprocal(M00[channel]))*((M20[channel]+M02[channel])+
1986 sqrt(4.0*M11[channel]*M11[channel]+(M20[channel]-M02[channel])*
1987 (M20[channel]-M02[channel]))));
1988 channel_moments[channel].ellipse_axis.y=sqrt((2.0*
1989 PerceptibleReciprocal(M00[channel]))*((M20[channel]+M02[channel])-
1990 sqrt(4.0*M11[channel]*M11[channel]+(M20[channel]-M02[channel])*
1991 (M20[channel]-M02[channel]))));
1992 channel_moments[channel].ellipse_angle=RadiansToDegrees(1.0/2.0*atan(2.0*
1993 M11[channel]*PerceptibleReciprocal(M20[channel]-M02[channel])));
1994 if (fabs(M11[channel]) < 0.0)
1995 {
1996 if ((fabs(M20[channel]-M02[channel]) >= 0.0) &&
1997 ((M20[channel]-M02[channel]) < 0.0))
1998 channel_moments[channel].ellipse_angle+=90.0;
1999 }
2000 else
2001 if (M11[channel] < 0.0)
2002 {
2003 if (fabs(M20[channel]-M02[channel]) >= 0.0)
2004 {
2005 if ((M20[channel]-M02[channel]) < 0.0)
2006 channel_moments[channel].ellipse_angle+=90.0;
2007 else
2008 channel_moments[channel].ellipse_angle+=180.0;
2009 }
2010 }
2011 else
2012 if ((fabs(M20[channel]-M02[channel]) >= 0.0) &&
2013 ((M20[channel]-M02[channel]) < 0.0))
2014 channel_moments[channel].ellipse_angle+=90.0;
2015 channel_moments[channel].ellipse_eccentricity=sqrt(1.0-(
2016 channel_moments[channel].ellipse_axis.y*
2017 channel_moments[channel].ellipse_axis.y*PerceptibleReciprocal(
2018 channel_moments[channel].ellipse_axis.x*
2019 channel_moments[channel].ellipse_axis.x)));
2020 channel_moments[channel].ellipse_intensity=M00[channel]/
2021 (MagickPI*channel_moments[channel].ellipse_axis.x*
2022 channel_moments[channel].ellipse_axis.y+MagickEpsilon);
2023 }
2024 for (channel=0; channel <= CompositeChannels; channel++)
2025 {
2026 /*
2027 Normalize image moments.
2028 */
2029 M10[channel]=0.0;
2030 M01[channel]=0.0;
2031 M11[channel]/=pow(M00[channel],1.0+(1.0+1.0)/2.0);
2032 M20[channel]/=pow(M00[channel],1.0+(2.0+0.0)/2.0);
2033 M02[channel]/=pow(M00[channel],1.0+(0.0+2.0)/2.0);
2034 M21[channel]/=pow(M00[channel],1.0+(2.0+1.0)/2.0);
2035 M12[channel]/=pow(M00[channel],1.0+(1.0+2.0)/2.0);
2036 M22[channel]/=pow(M00[channel],1.0+(2.0+2.0)/2.0);
2037 M30[channel]/=pow(M00[channel],1.0+(3.0+0.0)/2.0);
2038 M03[channel]/=pow(M00[channel],1.0+(0.0+3.0)/2.0);
2039 M00[channel]=1.0;
2040 }
2041 for (channel=0; channel <= CompositeChannels; channel++)
2042 {
2043 /*
2044 Compute Hu invariant moments.
2045 */
2046 channel_moments[channel].I[0]=M20[channel]+M02[channel];
2047 channel_moments[channel].I[1]=(M20[channel]-M02[channel])*
2048 (M20[channel]-M02[channel])+4.0*M11[channel]*M11[channel];
2049 channel_moments[channel].I[2]=(M30[channel]-3.0*M12[channel])*
2050 (M30[channel]-3.0*M12[channel])+(3.0*M21[channel]-M03[channel])*
2051 (3.0*M21[channel]-M03[channel]);
2052 channel_moments[channel].I[3]=(M30[channel]+M12[channel])*
2053 (M30[channel]+M12[channel])+(M21[channel]+M03[channel])*
2054 (M21[channel]+M03[channel]);
2055 channel_moments[channel].I[4]=(M30[channel]-3.0*M12[channel])*
2056 (M30[channel]+M12[channel])*((M30[channel]+M12[channel])*
2057 (M30[channel]+M12[channel])-3.0*(M21[channel]+M03[channel])*
2058 (M21[channel]+M03[channel]))+(3.0*M21[channel]-M03[channel])*
2059 (M21[channel]+M03[channel])*(3.0*(M30[channel]+M12[channel])*
2060 (M30[channel]+M12[channel])-(M21[channel]+M03[channel])*
2061 (M21[channel]+M03[channel]));
2062 channel_moments[channel].I[5]=(M20[channel]-M02[channel])*
2063 ((M30[channel]+M12[channel])*(M30[channel]+M12[channel])-
2064 (M21[channel]+M03[channel])*(M21[channel]+M03[channel]))+
2065 4.0*M11[channel]*(M30[channel]+M12[channel])*(M21[channel]+M03[channel]);
2066 channel_moments[channel].I[6]=(3.0*M21[channel]-M03[channel])*
2067 (M30[channel]+M12[channel])*((M30[channel]+M12[channel])*
2068 (M30[channel]+M12[channel])-3.0*(M21[channel]+M03[channel])*
2069 (M21[channel]+M03[channel]))-(M30[channel]-3*M12[channel])*
2070 (M21[channel]+M03[channel])*(3.0*(M30[channel]+M12[channel])*
2071 (M30[channel]+M12[channel])-(M21[channel]+M03[channel])*
2072 (M21[channel]+M03[channel]));
2073 channel_moments[channel].I[7]=M11[channel]*((M30[channel]+M12[channel])*
2074 (M30[channel]+M12[channel])-(M03[channel]+M21[channel])*
2075 (M03[channel]+M21[channel]))-(M20[channel]-M02[channel])*
2076 (M30[channel]+M12[channel])*(M03[channel]+M21[channel]);
2077 }
2078 if (y < (ssize_t) image->rows)
2079 channel_moments=(ChannelMoments *) RelinquishMagickMemory(channel_moments);
2080 return(channel_moments);
2081 }
2082
2083 /*
2084 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2085 % %
2086 % %
2087 % %
2088 % 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 %
2089 % %
2090 % %
2091 % %
2092 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2093 %
2094 % GetImageChannelPerceptualHash() returns the perceptual hash of one or more
2095 % image channels.
2096 %
2097 % The format of the GetImageChannelPerceptualHash method is:
2098 %
2099 % ChannelPerceptualHash *GetImageChannelPerceptualHash(const Image *image,
2100 % ExceptionInfo *exception)
2101 %
2102 % A description of each parameter follows:
2103 %
2104 % o image: the image.
2105 %
2106 % o exception: return any errors or warnings in this structure.
2107 %
2108 */
2109
MagickLog10(const double x)2110 static inline double MagickLog10(const double x)
2111 {
2112 #define Log10Epsilon (1.0e-11)
2113
2114 if (fabs(x) < Log10Epsilon)
2115 return(log10(Log10Epsilon));
2116 return(log10(fabs(x)));
2117 }
2118
GetImageChannelPerceptualHash(const Image * image,ExceptionInfo * exception)2119 MagickExport ChannelPerceptualHash *GetImageChannelPerceptualHash(
2120 const Image *image,ExceptionInfo *exception)
2121 {
2122 ChannelMoments
2123 *moments;
2124
2125 ChannelPerceptualHash
2126 *perceptual_hash;
2127
2128 Image
2129 *hash_image;
2130
2131 MagickBooleanType
2132 status;
2133
2134 ssize_t
2135 i;
2136
2137 ssize_t
2138 channel;
2139
2140 /*
2141 Blur then transform to sRGB colorspace.
2142 */
2143 hash_image=BlurImage(image,0.0,1.0,exception);
2144 if (hash_image == (Image *) NULL)
2145 return((ChannelPerceptualHash *) NULL);
2146 hash_image->depth=8;
2147 status=TransformImageColorspace(hash_image,sRGBColorspace);
2148 if (status == MagickFalse)
2149 return((ChannelPerceptualHash *) NULL);
2150 moments=GetImageChannelMoments(hash_image,exception);
2151 hash_image=DestroyImage(hash_image);
2152 if (moments == (ChannelMoments *) NULL)
2153 return((ChannelPerceptualHash *) NULL);
2154 perceptual_hash=(ChannelPerceptualHash *) AcquireQuantumMemory(
2155 CompositeChannels+1UL,sizeof(*perceptual_hash));
2156 if (perceptual_hash == (ChannelPerceptualHash *) NULL)
2157 return((ChannelPerceptualHash *) NULL);
2158 for (channel=0; channel <= CompositeChannels; channel++)
2159 for (i=0; i < MaximumNumberOfImageMoments; i++)
2160 perceptual_hash[channel].P[i]=(-MagickLog10(moments[channel].I[i]));
2161 moments=(ChannelMoments *) RelinquishMagickMemory(moments);
2162 /*
2163 Blur then transform to HCLp colorspace.
2164 */
2165 hash_image=BlurImage(image,0.0,1.0,exception);
2166 if (hash_image == (Image *) NULL)
2167 {
2168 perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
2169 perceptual_hash);
2170 return((ChannelPerceptualHash *) NULL);
2171 }
2172 hash_image->depth=8;
2173 status=TransformImageColorspace(hash_image,HCLpColorspace);
2174 if (status == MagickFalse)
2175 {
2176 perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
2177 perceptual_hash);
2178 return((ChannelPerceptualHash *) NULL);
2179 }
2180 moments=GetImageChannelMoments(hash_image,exception);
2181 hash_image=DestroyImage(hash_image);
2182 if (moments == (ChannelMoments *) NULL)
2183 {
2184 perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
2185 perceptual_hash);
2186 return((ChannelPerceptualHash *) NULL);
2187 }
2188 for (channel=0; channel <= CompositeChannels; channel++)
2189 for (i=0; i < MaximumNumberOfImageMoments; i++)
2190 perceptual_hash[channel].Q[i]=(-MagickLog10(moments[channel].I[i]));
2191 moments=(ChannelMoments *) RelinquishMagickMemory(moments);
2192 return(perceptual_hash);
2193 }
2194
2195 /*
2196 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2197 % %
2198 % %
2199 % %
2200 % G e t I m a g e C h a n n e l R a n g e %
2201 % %
2202 % %
2203 % %
2204 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2205 %
2206 % GetImageChannelRange() returns the range of one or more image channels.
2207 %
2208 % The format of the GetImageChannelRange method is:
2209 %
2210 % MagickBooleanType GetImageChannelRange(const Image *image,
2211 % const ChannelType channel,double *minima,double *maxima,
2212 % ExceptionInfo *exception)
2213 %
2214 % A description of each parameter follows:
2215 %
2216 % o image: the image.
2217 %
2218 % o channel: the channel.
2219 %
2220 % o minima: the minimum value in the channel.
2221 %
2222 % o maxima: the maximum value in the channel.
2223 %
2224 % o exception: return any errors or warnings in this structure.
2225 %
2226 */
2227
GetImageRange(const Image * image,double * minima,double * maxima,ExceptionInfo * exception)2228 MagickExport MagickBooleanType GetImageRange(const Image *image,
2229 double *minima,double *maxima,ExceptionInfo *exception)
2230 {
2231 return(GetImageChannelRange(image,CompositeChannels,minima,maxima,exception));
2232 }
2233
GetImageChannelRange(const Image * image,const ChannelType channel,double * minima,double * maxima,ExceptionInfo * exception)2234 MagickExport MagickBooleanType GetImageChannelRange(const Image *image,
2235 const ChannelType channel,double *minima,double *maxima,
2236 ExceptionInfo *exception)
2237 {
2238 MagickPixelPacket
2239 pixel;
2240
2241 ssize_t
2242 y;
2243
2244 assert(image != (Image *) NULL);
2245 assert(image->signature == MagickCoreSignature);
2246 if (image->debug != MagickFalse)
2247 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2248 *maxima=(-MagickMaximumValue);
2249 *minima=MagickMaximumValue;
2250 GetMagickPixelPacket(image,&pixel);
2251 for (y=0; y < (ssize_t) image->rows; y++)
2252 {
2253 const IndexPacket
2254 *magick_restrict indexes;
2255
2256 const PixelPacket
2257 *magick_restrict p;
2258
2259 ssize_t
2260 x;
2261
2262 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
2263 if (p == (const PixelPacket *) NULL)
2264 break;
2265 indexes=GetVirtualIndexQueue(image);
2266 for (x=0; x < (ssize_t) image->columns; x++)
2267 {
2268 SetMagickPixelPacket(image,p,indexes+x,&pixel);
2269 if ((channel & RedChannel) != 0)
2270 {
2271 if (pixel.red < *minima)
2272 *minima=(double) pixel.red;
2273 if (pixel.red > *maxima)
2274 *maxima=(double) pixel.red;
2275 }
2276 if ((channel & GreenChannel) != 0)
2277 {
2278 if (pixel.green < *minima)
2279 *minima=(double) pixel.green;
2280 if (pixel.green > *maxima)
2281 *maxima=(double) pixel.green;
2282 }
2283 if ((channel & BlueChannel) != 0)
2284 {
2285 if (pixel.blue < *minima)
2286 *minima=(double) pixel.blue;
2287 if (pixel.blue > *maxima)
2288 *maxima=(double) pixel.blue;
2289 }
2290 if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
2291 {
2292 if ((QuantumRange-pixel.opacity) < *minima)
2293 *minima=(double) (QuantumRange-pixel.opacity);
2294 if ((QuantumRange-pixel.opacity) > *maxima)
2295 *maxima=(double) (QuantumRange-pixel.opacity);
2296 }
2297 if (((channel & IndexChannel) != 0) &&
2298 (image->colorspace == CMYKColorspace))
2299 {
2300 if ((double) pixel.index < *minima)
2301 *minima=(double) pixel.index;
2302 if ((double) pixel.index > *maxima)
2303 *maxima=(double) pixel.index;
2304 }
2305 p++;
2306 }
2307 }
2308 return(y == (ssize_t) image->rows ? MagickTrue : MagickFalse);
2309 }
2310
2311 /*
2312 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2313 % %
2314 % %
2315 % %
2316 % G e t I m a g e C h a n n e l S t a t i s t i c s %
2317 % %
2318 % %
2319 % %
2320 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2321 %
2322 % GetImageChannelStatistics() returns statistics for each channel in the
2323 % image. The statistics include the channel depth, its minima, maxima, mean,
2324 % standard deviation, kurtosis and skewness. You can access the red channel
2325 % mean, for example, like this:
2326 %
2327 % channel_statistics=GetImageChannelStatistics(image,exception);
2328 % red_mean=channel_statistics[RedChannel].mean;
2329 %
2330 % Use MagickRelinquishMemory() to free the statistics buffer.
2331 %
2332 % The format of the GetImageChannelStatistics method is:
2333 %
2334 % ChannelStatistics *GetImageChannelStatistics(const Image *image,
2335 % ExceptionInfo *exception)
2336 %
2337 % A description of each parameter follows:
2338 %
2339 % o image: the image.
2340 %
2341 % o exception: return any errors or warnings in this structure.
2342 %
2343 */
GetImageChannelStatistics(const Image * image,ExceptionInfo * exception)2344 MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
2345 ExceptionInfo *exception)
2346 {
2347 ChannelStatistics
2348 *channel_statistics;
2349
2350 double
2351 area,
2352 standard_deviation;
2353
2354 MagickPixelPacket
2355 number_bins,
2356 *histogram;
2357
2358 QuantumAny
2359 range;
2360
2361 ssize_t
2362 i;
2363
2364 size_t
2365 channels,
2366 depth,
2367 length;
2368
2369 ssize_t
2370 y;
2371
2372 assert(image != (Image *) NULL);
2373 assert(image->signature == MagickCoreSignature);
2374 if (image->debug != MagickFalse)
2375 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2376 length=CompositeChannels+1UL;
2377 channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(length,
2378 sizeof(*channel_statistics));
2379 histogram=(MagickPixelPacket *) AcquireQuantumMemory(MaxMap+1U,
2380 sizeof(*histogram));
2381 if ((channel_statistics == (ChannelStatistics *) NULL) ||
2382 (histogram == (MagickPixelPacket *) NULL))
2383 {
2384 if (histogram != (MagickPixelPacket *) NULL)
2385 histogram=(MagickPixelPacket *) RelinquishMagickMemory(histogram);
2386 if (channel_statistics != (ChannelStatistics *) NULL)
2387 channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
2388 channel_statistics);
2389 return(channel_statistics);
2390 }
2391 (void) memset(channel_statistics,0,length*
2392 sizeof(*channel_statistics));
2393 for (i=0; i <= (ssize_t) CompositeChannels; i++)
2394 {
2395 channel_statistics[i].depth=1;
2396 channel_statistics[i].maxima=(-MagickMaximumValue);
2397 channel_statistics[i].minima=MagickMaximumValue;
2398 }
2399 (void) memset(histogram,0,(MaxMap+1U)*sizeof(*histogram));
2400 (void) memset(&number_bins,0,sizeof(number_bins));
2401 for (y=0; y < (ssize_t) image->rows; y++)
2402 {
2403 const IndexPacket
2404 *magick_restrict indexes;
2405
2406 const PixelPacket
2407 *magick_restrict p;
2408
2409 ssize_t
2410 x;
2411
2412 /*
2413 Compute pixel statistics.
2414 */
2415 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
2416 if (p == (const PixelPacket *) NULL)
2417 break;
2418 indexes=GetVirtualIndexQueue(image);
2419 for (x=0; x < (ssize_t) image->columns; )
2420 {
2421 if (channel_statistics[RedChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2422 {
2423 depth=channel_statistics[RedChannel].depth;
2424 range=GetQuantumRange(depth);
2425 if (IsPixelAtDepth(GetPixelRed(p),range) == MagickFalse)
2426 {
2427 channel_statistics[RedChannel].depth++;
2428 continue;
2429 }
2430 }
2431 if (channel_statistics[GreenChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2432 {
2433 depth=channel_statistics[GreenChannel].depth;
2434 range=GetQuantumRange(depth);
2435 if (IsPixelAtDepth(GetPixelGreen(p),range) == MagickFalse)
2436 {
2437 channel_statistics[GreenChannel].depth++;
2438 continue;
2439 }
2440 }
2441 if (channel_statistics[BlueChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2442 {
2443 depth=channel_statistics[BlueChannel].depth;
2444 range=GetQuantumRange(depth);
2445 if (IsPixelAtDepth(GetPixelBlue(p),range) == MagickFalse)
2446 {
2447 channel_statistics[BlueChannel].depth++;
2448 continue;
2449 }
2450 }
2451 if (image->matte != MagickFalse)
2452 {
2453 if (channel_statistics[OpacityChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2454 {
2455 depth=channel_statistics[OpacityChannel].depth;
2456 range=GetQuantumRange(depth);
2457 if (IsPixelAtDepth(GetPixelAlpha(p),range) == MagickFalse)
2458 {
2459 channel_statistics[OpacityChannel].depth++;
2460 continue;
2461 }
2462 }
2463 }
2464 if (image->colorspace == CMYKColorspace)
2465 {
2466 if (channel_statistics[BlackChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2467 {
2468 depth=channel_statistics[BlackChannel].depth;
2469 range=GetQuantumRange(depth);
2470 if (IsPixelAtDepth(GetPixelIndex(indexes+x),range) == MagickFalse)
2471 {
2472 channel_statistics[BlackChannel].depth++;
2473 continue;
2474 }
2475 }
2476 }
2477 if ((double) GetPixelRed(p) < channel_statistics[RedChannel].minima)
2478 channel_statistics[RedChannel].minima=(double) GetPixelRed(p);
2479 if ((double) GetPixelRed(p) > channel_statistics[RedChannel].maxima)
2480 channel_statistics[RedChannel].maxima=(double) GetPixelRed(p);
2481 channel_statistics[RedChannel].sum+=GetPixelRed(p);
2482 channel_statistics[RedChannel].sum_squared+=(double) GetPixelRed(p)*
2483 GetPixelRed(p);
2484 channel_statistics[RedChannel].sum_cubed+=(double)
2485 GetPixelRed(p)*GetPixelRed(p)*GetPixelRed(p);
2486 channel_statistics[RedChannel].sum_fourth_power+=(double)
2487 GetPixelRed(p)*GetPixelRed(p)*GetPixelRed(p)*GetPixelRed(p);
2488 if ((double) GetPixelGreen(p) < channel_statistics[GreenChannel].minima)
2489 channel_statistics[GreenChannel].minima=(double) GetPixelGreen(p);
2490 if ((double) GetPixelGreen(p) > channel_statistics[GreenChannel].maxima)
2491 channel_statistics[GreenChannel].maxima=(double) GetPixelGreen(p);
2492 channel_statistics[GreenChannel].sum+=GetPixelGreen(p);
2493 channel_statistics[GreenChannel].sum_squared+=(double) GetPixelGreen(p)*
2494 GetPixelGreen(p);
2495 channel_statistics[GreenChannel].sum_cubed+=(double) GetPixelGreen(p)*
2496 GetPixelGreen(p)*GetPixelGreen(p);
2497 channel_statistics[GreenChannel].sum_fourth_power+=(double)
2498 GetPixelGreen(p)*GetPixelGreen(p)*GetPixelGreen(p)*GetPixelGreen(p);
2499 if ((double) GetPixelBlue(p) < channel_statistics[BlueChannel].minima)
2500 channel_statistics[BlueChannel].minima=(double) GetPixelBlue(p);
2501 if ((double) GetPixelBlue(p) > channel_statistics[BlueChannel].maxima)
2502 channel_statistics[BlueChannel].maxima=(double) GetPixelBlue(p);
2503 channel_statistics[BlueChannel].sum+=GetPixelBlue(p);
2504 channel_statistics[BlueChannel].sum_squared+=(double) GetPixelBlue(p)*
2505 GetPixelBlue(p);
2506 channel_statistics[BlueChannel].sum_cubed+=(double) GetPixelBlue(p)*
2507 GetPixelBlue(p)*GetPixelBlue(p);
2508 channel_statistics[BlueChannel].sum_fourth_power+=(double)
2509 GetPixelBlue(p)*GetPixelBlue(p)*GetPixelBlue(p)*GetPixelBlue(p);
2510 histogram[ScaleQuantumToMap(GetPixelRed(p))].red++;
2511 histogram[ScaleQuantumToMap(GetPixelGreen(p))].green++;
2512 histogram[ScaleQuantumToMap(GetPixelBlue(p))].blue++;
2513 if (image->matte != MagickFalse)
2514 {
2515 if ((double) GetPixelAlpha(p) < channel_statistics[OpacityChannel].minima)
2516 channel_statistics[OpacityChannel].minima=(double) GetPixelAlpha(p);
2517 if ((double) GetPixelAlpha(p) > channel_statistics[OpacityChannel].maxima)
2518 channel_statistics[OpacityChannel].maxima=(double) GetPixelAlpha(p);
2519 channel_statistics[OpacityChannel].sum+=GetPixelAlpha(p);
2520 channel_statistics[OpacityChannel].sum_squared+=(double)
2521 GetPixelAlpha(p)*GetPixelAlpha(p);
2522 channel_statistics[OpacityChannel].sum_cubed+=(double)
2523 GetPixelAlpha(p)*GetPixelAlpha(p)*GetPixelAlpha(p);
2524 channel_statistics[OpacityChannel].sum_fourth_power+=(double)
2525 GetPixelAlpha(p)*GetPixelAlpha(p)*GetPixelAlpha(p)*GetPixelAlpha(p);
2526 histogram[ScaleQuantumToMap(GetPixelAlpha(p))].opacity++;
2527 }
2528 if (image->colorspace == CMYKColorspace)
2529 {
2530 if ((double) GetPixelIndex(indexes+x) < channel_statistics[BlackChannel].minima)
2531 channel_statistics[BlackChannel].minima=(double)
2532 GetPixelIndex(indexes+x);
2533 if ((double) GetPixelIndex(indexes+x) > channel_statistics[BlackChannel].maxima)
2534 channel_statistics[BlackChannel].maxima=(double)
2535 GetPixelIndex(indexes+x);
2536 channel_statistics[BlackChannel].sum+=GetPixelIndex(indexes+x);
2537 channel_statistics[BlackChannel].sum_squared+=(double)
2538 GetPixelIndex(indexes+x)*GetPixelIndex(indexes+x);
2539 channel_statistics[BlackChannel].sum_cubed+=(double)
2540 GetPixelIndex(indexes+x)*GetPixelIndex(indexes+x)*
2541 GetPixelIndex(indexes+x);
2542 channel_statistics[BlackChannel].sum_fourth_power+=(double)
2543 GetPixelIndex(indexes+x)*GetPixelIndex(indexes+x)*
2544 GetPixelIndex(indexes+x)*GetPixelIndex(indexes+x);
2545 histogram[ScaleQuantumToMap(GetPixelIndex(indexes+x))].index++;
2546 }
2547 x++;
2548 p++;
2549 }
2550 }
2551 for (i=0; i < (ssize_t) CompositeChannels; i++)
2552 {
2553 double
2554 area,
2555 mean,
2556 standard_deviation;
2557
2558 /*
2559 Normalize pixel statistics.
2560 */
2561 area=PerceptibleReciprocal((double) image->columns*image->rows);
2562 mean=channel_statistics[i].sum*area;
2563 channel_statistics[i].sum=mean;
2564 channel_statistics[i].sum_squared*=area;
2565 channel_statistics[i].sum_cubed*=area;
2566 channel_statistics[i].sum_fourth_power*=area;
2567 channel_statistics[i].mean=mean;
2568 channel_statistics[i].variance=channel_statistics[i].sum_squared;
2569 standard_deviation=sqrt(channel_statistics[i].variance-(mean*mean));
2570 area=PerceptibleReciprocal((double) image->columns*image->rows-1.0)*
2571 ((double) image->columns*image->rows);
2572 standard_deviation=sqrt(area*standard_deviation*standard_deviation);
2573 channel_statistics[i].standard_deviation=standard_deviation;
2574 }
2575 for (i=0; i < (ssize_t) (MaxMap+1U); i++)
2576 {
2577 if (histogram[i].red > 0.0)
2578 number_bins.red++;
2579 if (histogram[i].green > 0.0)
2580 number_bins.green++;
2581 if (histogram[i].blue > 0.0)
2582 number_bins.blue++;
2583 if ((image->matte != MagickFalse) && (histogram[i].opacity > 0.0))
2584 number_bins.opacity++;
2585 if ((image->colorspace == CMYKColorspace) && (histogram[i].index > 0.0))
2586 number_bins.index++;
2587 }
2588 area=PerceptibleReciprocal((double) image->columns*image->rows);
2589 for (i=0; i < (ssize_t) (MaxMap+1U); i++)
2590 {
2591 /*
2592 Compute pixel entropy.
2593 */
2594 histogram[i].red*=area;
2595 channel_statistics[RedChannel].entropy+=-histogram[i].red*
2596 MagickLog10(histogram[i].red)*
2597 PerceptibleReciprocal(MagickLog10((double) number_bins.red));
2598 histogram[i].green*=area;
2599 channel_statistics[GreenChannel].entropy+=-histogram[i].green*
2600 MagickLog10(histogram[i].green)*
2601 PerceptibleReciprocal(MagickLog10((double) number_bins.green));
2602 histogram[i].blue*=area;
2603 channel_statistics[BlueChannel].entropy+=-histogram[i].blue*
2604 MagickLog10(histogram[i].blue)*
2605 PerceptibleReciprocal(MagickLog10((double) number_bins.blue));
2606 if (image->matte != MagickFalse)
2607 {
2608 histogram[i].opacity*=area;
2609 channel_statistics[OpacityChannel].entropy+=-histogram[i].opacity*
2610 MagickLog10(histogram[i].opacity)*
2611 PerceptibleReciprocal(MagickLog10((double) number_bins.opacity));
2612 }
2613 if (image->colorspace == CMYKColorspace)
2614 {
2615 histogram[i].index*=area;
2616 channel_statistics[IndexChannel].entropy+=-histogram[i].index*
2617 MagickLog10(histogram[i].index)*
2618 PerceptibleReciprocal(MagickLog10((double) number_bins.index));
2619 }
2620 }
2621 /*
2622 Compute overall statistics.
2623 */
2624 for (i=0; i < (ssize_t) CompositeChannels; i++)
2625 {
2626 channel_statistics[CompositeChannels].depth=(size_t) EvaluateMax((double)
2627 channel_statistics[CompositeChannels].depth,(double)
2628 channel_statistics[i].depth);
2629 channel_statistics[CompositeChannels].minima=MagickMin(
2630 channel_statistics[CompositeChannels].minima,
2631 channel_statistics[i].minima);
2632 channel_statistics[CompositeChannels].maxima=EvaluateMax(
2633 channel_statistics[CompositeChannels].maxima,
2634 channel_statistics[i].maxima);
2635 channel_statistics[CompositeChannels].sum+=channel_statistics[i].sum;
2636 channel_statistics[CompositeChannels].sum_squared+=
2637 channel_statistics[i].sum_squared;
2638 channel_statistics[CompositeChannels].sum_cubed+=
2639 channel_statistics[i].sum_cubed;
2640 channel_statistics[CompositeChannels].sum_fourth_power+=
2641 channel_statistics[i].sum_fourth_power;
2642 channel_statistics[CompositeChannels].mean+=channel_statistics[i].mean;
2643 channel_statistics[CompositeChannels].variance+=
2644 channel_statistics[i].variance-channel_statistics[i].mean*
2645 channel_statistics[i].mean;
2646 standard_deviation=sqrt(channel_statistics[i].variance-
2647 (channel_statistics[i].mean*channel_statistics[i].mean));
2648 area=PerceptibleReciprocal((double) image->columns*image->rows-1.0)*
2649 ((double) image->columns*image->rows);
2650 standard_deviation=sqrt(area*standard_deviation*standard_deviation);
2651 channel_statistics[CompositeChannels].standard_deviation=standard_deviation;
2652 channel_statistics[CompositeChannels].entropy+=
2653 channel_statistics[i].entropy;
2654 }
2655 channels=3;
2656 if (image->matte != MagickFalse)
2657 channels++;
2658 if (image->colorspace == CMYKColorspace)
2659 channels++;
2660 channel_statistics[CompositeChannels].sum/=channels;
2661 channel_statistics[CompositeChannels].sum_squared/=channels;
2662 channel_statistics[CompositeChannels].sum_cubed/=channels;
2663 channel_statistics[CompositeChannels].sum_fourth_power/=channels;
2664 channel_statistics[CompositeChannels].mean/=channels;
2665 channel_statistics[CompositeChannels].kurtosis/=channels;
2666 channel_statistics[CompositeChannels].skewness/=channels;
2667 channel_statistics[CompositeChannels].entropy/=channels;
2668 i=CompositeChannels;
2669 area=PerceptibleReciprocal((double) channels*image->columns*image->rows);
2670 channel_statistics[i].variance=channel_statistics[i].sum_squared;
2671 channel_statistics[i].mean=channel_statistics[i].sum;
2672 standard_deviation=sqrt(channel_statistics[i].variance-
2673 (channel_statistics[i].mean*channel_statistics[i].mean));
2674 standard_deviation=sqrt(PerceptibleReciprocal((double) channels*
2675 image->columns*image->rows-1.0)*channels*image->columns*image->rows*
2676 standard_deviation*standard_deviation);
2677 channel_statistics[i].standard_deviation=standard_deviation;
2678 for (i=0; i <= (ssize_t) CompositeChannels; i++)
2679 {
2680 /*
2681 Compute kurtosis & skewness statistics.
2682 */
2683 standard_deviation=PerceptibleReciprocal(
2684 channel_statistics[i].standard_deviation);
2685 channel_statistics[i].skewness=(channel_statistics[i].sum_cubed-3.0*
2686 channel_statistics[i].mean*channel_statistics[i].sum_squared+2.0*
2687 channel_statistics[i].mean*channel_statistics[i].mean*
2688 channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2689 standard_deviation);
2690 channel_statistics[i].kurtosis=(channel_statistics[i].sum_fourth_power-4.0*
2691 channel_statistics[i].mean*channel_statistics[i].sum_cubed+6.0*
2692 channel_statistics[i].mean*channel_statistics[i].mean*
2693 channel_statistics[i].sum_squared-3.0*channel_statistics[i].mean*
2694 channel_statistics[i].mean*1.0*channel_statistics[i].mean*
2695 channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2696 standard_deviation*standard_deviation)-3.0;
2697 }
2698 channel_statistics[CompositeChannels].mean=0.0;
2699 channel_statistics[CompositeChannels].standard_deviation=0.0;
2700 for (i=0; i < (ssize_t) CompositeChannels; i++)
2701 {
2702 channel_statistics[CompositeChannels].mean+=
2703 channel_statistics[i].mean;
2704 channel_statistics[CompositeChannels].standard_deviation+=
2705 channel_statistics[i].standard_deviation;
2706 }
2707 channel_statistics[CompositeChannels].mean/=(double) channels;
2708 channel_statistics[CompositeChannels].standard_deviation/=(double) channels;
2709 histogram=(MagickPixelPacket *) RelinquishMagickMemory(histogram);
2710 if (y < (ssize_t) image->rows)
2711 channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
2712 channel_statistics);
2713 return(channel_statistics);
2714 }
2715
2716 /*
2717 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2718 % %
2719 % %
2720 % %
2721 % P o l y n o m i a l I m a g e %
2722 % %
2723 % %
2724 % %
2725 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2726 %
2727 % PolynomialImage() returns a new image where each pixel is the sum of the
2728 % pixels in the image sequence after applying its corresponding terms
2729 % (coefficient and degree pairs).
2730 %
2731 % The format of the PolynomialImage method is:
2732 %
2733 % Image *PolynomialImage(const Image *images,const size_t number_terms,
2734 % const double *terms,ExceptionInfo *exception)
2735 % Image *PolynomialImageChannel(const Image *images,
2736 % const size_t number_terms,const ChannelType channel,
2737 % const double *terms,ExceptionInfo *exception)
2738 %
2739 % A description of each parameter follows:
2740 %
2741 % o images: the image sequence.
2742 %
2743 % o channel: the channel.
2744 %
2745 % o number_terms: the number of terms in the list. The actual list length
2746 % is 2 x number_terms + 1 (the constant).
2747 %
2748 % o terms: the list of polynomial coefficients and degree pairs and a
2749 % constant.
2750 %
2751 % o exception: return any errors or warnings in this structure.
2752 %
2753 */
PolynomialImage(const Image * images,const size_t number_terms,const double * terms,ExceptionInfo * exception)2754 MagickExport Image *PolynomialImage(const Image *images,
2755 const size_t number_terms,const double *terms,ExceptionInfo *exception)
2756 {
2757 Image
2758 *polynomial_image;
2759
2760 polynomial_image=PolynomialImageChannel(images,DefaultChannels,number_terms,
2761 terms,exception);
2762 return(polynomial_image);
2763 }
2764
PolynomialImageChannel(const Image * images,const ChannelType channel,const size_t number_terms,const double * terms,ExceptionInfo * exception)2765 MagickExport Image *PolynomialImageChannel(const Image *images,
2766 const ChannelType channel,const size_t number_terms,const double *terms,
2767 ExceptionInfo *exception)
2768 {
2769 #define PolynomialImageTag "Polynomial/Image"
2770
2771 CacheView
2772 *polynomial_view;
2773
2774 Image
2775 *image;
2776
2777 MagickBooleanType
2778 status;
2779
2780 MagickOffsetType
2781 progress;
2782
2783 MagickPixelPacket
2784 **magick_restrict polynomial_pixels,
2785 zero;
2786
2787 ssize_t
2788 y;
2789
2790 assert(images != (Image *) NULL);
2791 assert(images->signature == MagickCoreSignature);
2792 if (images->debug != MagickFalse)
2793 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
2794 assert(exception != (ExceptionInfo *) NULL);
2795 assert(exception->signature == MagickCoreSignature);
2796 image=AcquireImageCanvas(images,exception);
2797 if (image == (Image *) NULL)
2798 return((Image *) NULL);
2799 if (SetImageStorageClass(image,DirectClass) == MagickFalse)
2800 {
2801 InheritException(exception,&image->exception);
2802 image=DestroyImage(image);
2803 return((Image *) NULL);
2804 }
2805 polynomial_pixels=AcquirePixelThreadSet(images);
2806 if (polynomial_pixels == (MagickPixelPacket **) NULL)
2807 {
2808 image=DestroyImage(image);
2809 (void) ThrowMagickException(exception,GetMagickModule(),
2810 ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
2811 return((Image *) NULL);
2812 }
2813 /*
2814 Polynomial image pixels.
2815 */
2816 status=MagickTrue;
2817 progress=0;
2818 GetMagickPixelPacket(images,&zero);
2819 polynomial_view=AcquireAuthenticCacheView(image,exception);
2820 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2821 #pragma omp parallel for schedule(static) shared(progress,status) \
2822 magick_number_threads(image,image,image->rows,1)
2823 #endif
2824 for (y=0; y < (ssize_t) image->rows; y++)
2825 {
2826 CacheView
2827 *image_view;
2828
2829 const Image
2830 *next;
2831
2832 const int
2833 id = GetOpenMPThreadId();
2834
2835 IndexPacket
2836 *magick_restrict polynomial_indexes;
2837
2838 MagickPixelPacket
2839 *polynomial_pixel;
2840
2841 PixelPacket
2842 *magick_restrict q;
2843
2844 ssize_t
2845 i,
2846 x;
2847
2848 size_t
2849 number_images;
2850
2851 if (status == MagickFalse)
2852 continue;
2853 q=QueueCacheViewAuthenticPixels(polynomial_view,0,y,image->columns,1,
2854 exception);
2855 if (q == (PixelPacket *) NULL)
2856 {
2857 status=MagickFalse;
2858 continue;
2859 }
2860 polynomial_indexes=GetCacheViewAuthenticIndexQueue(polynomial_view);
2861 polynomial_pixel=polynomial_pixels[id];
2862 for (x=0; x < (ssize_t) image->columns; x++)
2863 polynomial_pixel[x]=zero;
2864 next=images;
2865 number_images=GetImageListLength(images);
2866 for (i=0; i < (ssize_t) number_images; i++)
2867 {
2868 const IndexPacket
2869 *indexes;
2870
2871 const PixelPacket
2872 *p;
2873
2874 if (i >= (ssize_t) number_terms)
2875 break;
2876 image_view=AcquireVirtualCacheView(next,exception);
2877 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2878 if (p == (const PixelPacket *) NULL)
2879 {
2880 image_view=DestroyCacheView(image_view);
2881 break;
2882 }
2883 indexes=GetCacheViewVirtualIndexQueue(image_view);
2884 for (x=0; x < (ssize_t) image->columns; x++)
2885 {
2886 double
2887 coefficient,
2888 degree;
2889
2890 coefficient=terms[i << 1];
2891 degree=terms[(i << 1)+1];
2892 if ((channel & RedChannel) != 0)
2893 polynomial_pixel[x].red+=coefficient*pow(QuantumScale*p->red,degree);
2894 if ((channel & GreenChannel) != 0)
2895 polynomial_pixel[x].green+=coefficient*pow(QuantumScale*p->green,
2896 degree);
2897 if ((channel & BlueChannel) != 0)
2898 polynomial_pixel[x].blue+=coefficient*pow(QuantumScale*p->blue,
2899 degree);
2900 if ((channel & OpacityChannel) != 0)
2901 polynomial_pixel[x].opacity+=coefficient*pow(QuantumScale*
2902 (QuantumRange-p->opacity),degree);
2903 if (((channel & IndexChannel) != 0) &&
2904 (image->colorspace == CMYKColorspace))
2905 polynomial_pixel[x].index+=coefficient*pow(QuantumScale*indexes[x],
2906 degree);
2907 p++;
2908 }
2909 image_view=DestroyCacheView(image_view);
2910 next=GetNextImageInList(next);
2911 }
2912 for (x=0; x < (ssize_t) image->columns; x++)
2913 {
2914 SetPixelRed(q,ClampToQuantum(QuantumRange*polynomial_pixel[x].red));
2915 SetPixelGreen(q,ClampToQuantum(QuantumRange*polynomial_pixel[x].green));
2916 SetPixelBlue(q,ClampToQuantum(QuantumRange*polynomial_pixel[x].blue));
2917 if (image->matte == MagickFalse)
2918 SetPixelOpacity(q,ClampToQuantum(QuantumRange-QuantumRange*
2919 polynomial_pixel[x].opacity));
2920 else
2921 SetPixelAlpha(q,ClampToQuantum(QuantumRange-QuantumRange*
2922 polynomial_pixel[x].opacity));
2923 if (image->colorspace == CMYKColorspace)
2924 SetPixelIndex(polynomial_indexes+x,ClampToQuantum(QuantumRange*
2925 polynomial_pixel[x].index));
2926 q++;
2927 }
2928 if (SyncCacheViewAuthenticPixels(polynomial_view,exception) == MagickFalse)
2929 status=MagickFalse;
2930 if (images->progress_monitor != (MagickProgressMonitor) NULL)
2931 {
2932 MagickBooleanType
2933 proceed;
2934
2935 proceed=SetImageProgress(images,PolynomialImageTag,progress++,
2936 image->rows);
2937 if (proceed == MagickFalse)
2938 status=MagickFalse;
2939 }
2940 }
2941 polynomial_view=DestroyCacheView(polynomial_view);
2942 polynomial_pixels=DestroyPixelThreadSet(images,polynomial_pixels);
2943 if (status == MagickFalse)
2944 image=DestroyImage(image);
2945 return(image);
2946 }
2947
2948 /*
2949 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2950 % %
2951 % %
2952 % %
2953 % S t a t i s t i c I m a g e %
2954 % %
2955 % %
2956 % %
2957 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2958 %
2959 % StatisticImage() makes each pixel the min / max / median / mode / etc. of
2960 % the neighborhood of the specified width and height.
2961 %
2962 % The format of the StatisticImage method is:
2963 %
2964 % Image *StatisticImage(const Image *image,const StatisticType type,
2965 % const size_t width,const size_t height,ExceptionInfo *exception)
2966 % Image *StatisticImageChannel(const Image *image,
2967 % const ChannelType channel,const StatisticType type,
2968 % const size_t width,const size_t height,ExceptionInfo *exception)
2969 %
2970 % A description of each parameter follows:
2971 %
2972 % o image: the image.
2973 %
2974 % o channel: the image channel.
2975 %
2976 % o type: the statistic type (median, mode, etc.).
2977 %
2978 % o width: the width of the pixel neighborhood.
2979 %
2980 % o height: the height of the pixel neighborhood.
2981 %
2982 % o exception: return any errors or warnings in this structure.
2983 %
2984 */
2985
2986 #define ListChannels 5
2987
2988 typedef struct _ListNode
2989 {
2990 size_t
2991 next[9],
2992 count,
2993 signature;
2994 } ListNode;
2995
2996 typedef struct _SkipList
2997 {
2998 ssize_t
2999 level;
3000
3001 ListNode
3002 *nodes;
3003 } SkipList;
3004
3005 typedef struct _PixelList
3006 {
3007 size_t
3008 length,
3009 seed,
3010 signature;
3011
3012 SkipList
3013 lists[ListChannels];
3014 } PixelList;
3015
DestroyPixelList(PixelList * pixel_list)3016 static PixelList *DestroyPixelList(PixelList *pixel_list)
3017 {
3018 ssize_t
3019 i;
3020
3021 if (pixel_list == (PixelList *) NULL)
3022 return((PixelList *) NULL);
3023 for (i=0; i < ListChannels; i++)
3024 if (pixel_list->lists[i].nodes != (ListNode *) NULL)
3025 pixel_list->lists[i].nodes=(ListNode *) RelinquishAlignedMemory(
3026 pixel_list->lists[i].nodes);
3027 pixel_list=(PixelList *) RelinquishMagickMemory(pixel_list);
3028 return(pixel_list);
3029 }
3030
DestroyPixelListThreadSet(PixelList ** pixel_list)3031 static PixelList **DestroyPixelListThreadSet(PixelList **pixel_list)
3032 {
3033 ssize_t
3034 i;
3035
3036 assert(pixel_list != (PixelList **) NULL);
3037 for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
3038 if (pixel_list[i] != (PixelList *) NULL)
3039 pixel_list[i]=DestroyPixelList(pixel_list[i]);
3040 pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list);
3041 return(pixel_list);
3042 }
3043
AcquirePixelList(const size_t width,const size_t height)3044 static PixelList *AcquirePixelList(const size_t width,const size_t height)
3045 {
3046 PixelList
3047 *pixel_list;
3048
3049 ssize_t
3050 i;
3051
3052 pixel_list=(PixelList *) AcquireMagickMemory(sizeof(*pixel_list));
3053 if (pixel_list == (PixelList *) NULL)
3054 return(pixel_list);
3055 (void) memset((void *) pixel_list,0,sizeof(*pixel_list));
3056 pixel_list->length=width*height;
3057 for (i=0; i < ListChannels; i++)
3058 {
3059 pixel_list->lists[i].nodes=(ListNode *) AcquireAlignedMemory(65537UL,
3060 sizeof(*pixel_list->lists[i].nodes));
3061 if (pixel_list->lists[i].nodes == (ListNode *) NULL)
3062 return(DestroyPixelList(pixel_list));
3063 (void) memset(pixel_list->lists[i].nodes,0,65537UL*
3064 sizeof(*pixel_list->lists[i].nodes));
3065 }
3066 pixel_list->signature=MagickCoreSignature;
3067 return(pixel_list);
3068 }
3069
AcquirePixelListThreadSet(const size_t width,const size_t height)3070 static PixelList **AcquirePixelListThreadSet(const size_t width,
3071 const size_t height)
3072 {
3073 PixelList
3074 **pixel_list;
3075
3076 ssize_t
3077 i;
3078
3079 size_t
3080 number_threads;
3081
3082 number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
3083 pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
3084 sizeof(*pixel_list));
3085 if (pixel_list == (PixelList **) NULL)
3086 return((PixelList **) NULL);
3087 (void) memset(pixel_list,0,number_threads*sizeof(*pixel_list));
3088 for (i=0; i < (ssize_t) number_threads; i++)
3089 {
3090 pixel_list[i]=AcquirePixelList(width,height);
3091 if (pixel_list[i] == (PixelList *) NULL)
3092 return(DestroyPixelListThreadSet(pixel_list));
3093 }
3094 return(pixel_list);
3095 }
3096
AddNodePixelList(PixelList * pixel_list,const ssize_t channel,const size_t color)3097 static void AddNodePixelList(PixelList *pixel_list,const ssize_t channel,
3098 const size_t color)
3099 {
3100 SkipList
3101 *list;
3102
3103 ssize_t
3104 level;
3105
3106 size_t
3107 search,
3108 update[9];
3109
3110 /*
3111 Initialize the node.
3112 */
3113 list=pixel_list->lists+channel;
3114 list->nodes[color].signature=pixel_list->signature;
3115 list->nodes[color].count=1;
3116 /*
3117 Determine where it belongs in the list.
3118 */
3119 search=65536UL;
3120 for (level=list->level; level >= 0; level--)
3121 {
3122 while (list->nodes[search].next[level] < color)
3123 search=list->nodes[search].next[level];
3124 update[level]=search;
3125 }
3126 /*
3127 Generate a pseudo-random level for this node.
3128 */
3129 for (level=0; ; level++)
3130 {
3131 pixel_list->seed=(pixel_list->seed*42893621L)+1L;
3132 if ((pixel_list->seed & 0x300) != 0x300)
3133 break;
3134 }
3135 if (level > 8)
3136 level=8;
3137 if (level > (list->level+2))
3138 level=list->level+2;
3139 /*
3140 If we're raising the list's level, link back to the root node.
3141 */
3142 while (level > list->level)
3143 {
3144 list->level++;
3145 update[list->level]=65536UL;
3146 }
3147 /*
3148 Link the node into the skip-list.
3149 */
3150 do
3151 {
3152 list->nodes[color].next[level]=list->nodes[update[level]].next[level];
3153 list->nodes[update[level]].next[level]=color;
3154 } while (level-- > 0);
3155 }
3156
GetMaximumPixelList(PixelList * pixel_list,MagickPixelPacket * pixel)3157 static void GetMaximumPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3158 {
3159 SkipList
3160 *list;
3161
3162 ssize_t
3163 channel;
3164
3165 size_t
3166 color,
3167 maximum;
3168
3169 ssize_t
3170 count;
3171
3172 unsigned short
3173 channels[ListChannels];
3174
3175 /*
3176 Find the maximum value for each of the color.
3177 */
3178 for (channel=0; channel < 5; channel++)
3179 {
3180 list=pixel_list->lists+channel;
3181 color=65536L;
3182 count=0;
3183 maximum=list->nodes[color].next[0];
3184 do
3185 {
3186 color=list->nodes[color].next[0];
3187 if (color > maximum)
3188 maximum=color;
3189 count+=list->nodes[color].count;
3190 } while (count < (ssize_t) pixel_list->length);
3191 channels[channel]=(unsigned short) maximum;
3192 }
3193 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3194 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3195 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3196 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3197 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3198 }
3199
GetMeanPixelList(PixelList * pixel_list,MagickPixelPacket * pixel)3200 static void GetMeanPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3201 {
3202 MagickRealType
3203 sum;
3204
3205 SkipList
3206 *list;
3207
3208 ssize_t
3209 channel;
3210
3211 size_t
3212 color;
3213
3214 ssize_t
3215 count;
3216
3217 unsigned short
3218 channels[ListChannels];
3219
3220 /*
3221 Find the mean value for each of the color.
3222 */
3223 for (channel=0; channel < 5; channel++)
3224 {
3225 list=pixel_list->lists+channel;
3226 color=65536L;
3227 count=0;
3228 sum=0.0;
3229 do
3230 {
3231 color=list->nodes[color].next[0];
3232 sum+=(MagickRealType) list->nodes[color].count*color;
3233 count+=list->nodes[color].count;
3234 } while (count < (ssize_t) pixel_list->length);
3235 sum/=pixel_list->length;
3236 channels[channel]=(unsigned short) sum;
3237 }
3238 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3239 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3240 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3241 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3242 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3243 }
3244
GetMedianPixelList(PixelList * pixel_list,MagickPixelPacket * pixel)3245 static void GetMedianPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3246 {
3247 SkipList
3248 *list;
3249
3250 ssize_t
3251 channel;
3252
3253 size_t
3254 color;
3255
3256 ssize_t
3257 count;
3258
3259 unsigned short
3260 channels[ListChannels];
3261
3262 /*
3263 Find the median value for each of the color.
3264 */
3265 for (channel=0; channel < 5; channel++)
3266 {
3267 list=pixel_list->lists+channel;
3268 color=65536L;
3269 count=0;
3270 do
3271 {
3272 color=list->nodes[color].next[0];
3273 count+=list->nodes[color].count;
3274 } while (count <= (ssize_t) (pixel_list->length >> 1));
3275 channels[channel]=(unsigned short) color;
3276 }
3277 GetMagickPixelPacket((const Image *) NULL,pixel);
3278 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3279 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3280 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3281 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3282 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3283 }
3284
GetMinimumPixelList(PixelList * pixel_list,MagickPixelPacket * pixel)3285 static void GetMinimumPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3286 {
3287 SkipList
3288 *list;
3289
3290 ssize_t
3291 channel;
3292
3293 size_t
3294 color,
3295 minimum;
3296
3297 ssize_t
3298 count;
3299
3300 unsigned short
3301 channels[ListChannels];
3302
3303 /*
3304 Find the minimum value for each of the color.
3305 */
3306 for (channel=0; channel < 5; channel++)
3307 {
3308 list=pixel_list->lists+channel;
3309 count=0;
3310 color=65536UL;
3311 minimum=list->nodes[color].next[0];
3312 do
3313 {
3314 color=list->nodes[color].next[0];
3315 if (color < minimum)
3316 minimum=color;
3317 count+=list->nodes[color].count;
3318 } while (count < (ssize_t) pixel_list->length);
3319 channels[channel]=(unsigned short) minimum;
3320 }
3321 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3322 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3323 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3324 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3325 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3326 }
3327
GetModePixelList(PixelList * pixel_list,MagickPixelPacket * pixel)3328 static void GetModePixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3329 {
3330 SkipList
3331 *list;
3332
3333 ssize_t
3334 channel;
3335
3336 size_t
3337 color,
3338 max_count,
3339 mode;
3340
3341 ssize_t
3342 count;
3343
3344 unsigned short
3345 channels[5];
3346
3347 /*
3348 Make each pixel the 'predominant color' of the specified neighborhood.
3349 */
3350 for (channel=0; channel < 5; channel++)
3351 {
3352 list=pixel_list->lists+channel;
3353 color=65536L;
3354 mode=color;
3355 max_count=list->nodes[mode].count;
3356 count=0;
3357 do
3358 {
3359 color=list->nodes[color].next[0];
3360 if (list->nodes[color].count > max_count)
3361 {
3362 mode=color;
3363 max_count=list->nodes[mode].count;
3364 }
3365 count+=list->nodes[color].count;
3366 } while (count < (ssize_t) pixel_list->length);
3367 channels[channel]=(unsigned short) mode;
3368 }
3369 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3370 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3371 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3372 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3373 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3374 }
3375
GetNonpeakPixelList(PixelList * pixel_list,MagickPixelPacket * pixel)3376 static void GetNonpeakPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3377 {
3378 SkipList
3379 *list;
3380
3381 ssize_t
3382 channel;
3383
3384 size_t
3385 color,
3386 next,
3387 previous;
3388
3389 ssize_t
3390 count;
3391
3392 unsigned short
3393 channels[5];
3394
3395 /*
3396 Finds the non peak value for each of the colors.
3397 */
3398 for (channel=0; channel < 5; channel++)
3399 {
3400 list=pixel_list->lists+channel;
3401 color=65536L;
3402 next=list->nodes[color].next[0];
3403 count=0;
3404 do
3405 {
3406 previous=color;
3407 color=next;
3408 next=list->nodes[color].next[0];
3409 count+=list->nodes[color].count;
3410 } while (count <= (ssize_t) (pixel_list->length >> 1));
3411 if ((previous == 65536UL) && (next != 65536UL))
3412 color=next;
3413 else
3414 if ((previous != 65536UL) && (next == 65536UL))
3415 color=previous;
3416 channels[channel]=(unsigned short) color;
3417 }
3418 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3419 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3420 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3421 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3422 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3423 }
3424
GetRootMeanSquarePixelList(PixelList * pixel_list,MagickPixelPacket * pixel)3425 static void GetRootMeanSquarePixelList(PixelList *pixel_list,
3426 MagickPixelPacket *pixel)
3427 {
3428 MagickRealType
3429 sum;
3430
3431 SkipList
3432 *list;
3433
3434 ssize_t
3435 channel;
3436
3437 size_t
3438 color;
3439
3440 ssize_t
3441 count;
3442
3443 unsigned short
3444 channels[ListChannels];
3445
3446 /*
3447 Find the root mean square value for each of the color.
3448 */
3449 for (channel=0; channel < 5; channel++)
3450 {
3451 list=pixel_list->lists+channel;
3452 color=65536L;
3453 count=0;
3454 sum=0.0;
3455 do
3456 {
3457 color=list->nodes[color].next[0];
3458 sum+=(MagickRealType) (list->nodes[color].count*color*color);
3459 count+=list->nodes[color].count;
3460 } while (count < (ssize_t) pixel_list->length);
3461 sum/=pixel_list->length;
3462 channels[channel]=(unsigned short) sqrt(sum);
3463 }
3464 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3465 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3466 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3467 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3468 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3469 }
3470
GetStandardDeviationPixelList(PixelList * pixel_list,MagickPixelPacket * pixel)3471 static void GetStandardDeviationPixelList(PixelList *pixel_list,
3472 MagickPixelPacket *pixel)
3473 {
3474 MagickRealType
3475 sum,
3476 sum_squared;
3477
3478 SkipList
3479 *list;
3480
3481 ssize_t
3482 channel;
3483
3484 size_t
3485 color;
3486
3487 ssize_t
3488 count;
3489
3490 unsigned short
3491 channels[ListChannels];
3492
3493 /*
3494 Find the standard-deviation value for each of the color.
3495 */
3496 for (channel=0; channel < 5; channel++)
3497 {
3498 list=pixel_list->lists+channel;
3499 color=65536L;
3500 count=0;
3501 sum=0.0;
3502 sum_squared=0.0;
3503 do
3504 {
3505 ssize_t
3506 i;
3507
3508 color=list->nodes[color].next[0];
3509 sum+=(MagickRealType) list->nodes[color].count*color;
3510 for (i=0; i < (ssize_t) list->nodes[color].count; i++)
3511 sum_squared+=((MagickRealType) color)*((MagickRealType) color);
3512 count+=list->nodes[color].count;
3513 } while (count < (ssize_t) pixel_list->length);
3514 sum/=pixel_list->length;
3515 sum_squared/=pixel_list->length;
3516 channels[channel]=(unsigned short) sqrt(sum_squared-(sum*sum));
3517 }
3518 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3519 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3520 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3521 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3522 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3523 }
3524
InsertPixelList(const Image * image,const PixelPacket * pixel,const IndexPacket * indexes,PixelList * pixel_list)3525 static inline void InsertPixelList(const Image *image,const PixelPacket *pixel,
3526 const IndexPacket *indexes,PixelList *pixel_list)
3527 {
3528 size_t
3529 signature;
3530
3531 unsigned short
3532 index;
3533
3534 index=ScaleQuantumToShort(GetPixelRed(pixel));
3535 signature=pixel_list->lists[0].nodes[index].signature;
3536 if (signature == pixel_list->signature)
3537 pixel_list->lists[0].nodes[index].count++;
3538 else
3539 AddNodePixelList(pixel_list,0,index);
3540 index=ScaleQuantumToShort(GetPixelGreen(pixel));
3541 signature=pixel_list->lists[1].nodes[index].signature;
3542 if (signature == pixel_list->signature)
3543 pixel_list->lists[1].nodes[index].count++;
3544 else
3545 AddNodePixelList(pixel_list,1,index);
3546 index=ScaleQuantumToShort(GetPixelBlue(pixel));
3547 signature=pixel_list->lists[2].nodes[index].signature;
3548 if (signature == pixel_list->signature)
3549 pixel_list->lists[2].nodes[index].count++;
3550 else
3551 AddNodePixelList(pixel_list,2,index);
3552 index=ScaleQuantumToShort(GetPixelOpacity(pixel));
3553 signature=pixel_list->lists[3].nodes[index].signature;
3554 if (signature == pixel_list->signature)
3555 pixel_list->lists[3].nodes[index].count++;
3556 else
3557 AddNodePixelList(pixel_list,3,index);
3558 if (image->colorspace == CMYKColorspace)
3559 index=ScaleQuantumToShort(GetPixelIndex(indexes));
3560 signature=pixel_list->lists[4].nodes[index].signature;
3561 if (signature == pixel_list->signature)
3562 pixel_list->lists[4].nodes[index].count++;
3563 else
3564 AddNodePixelList(pixel_list,4,index);
3565 }
3566
ResetPixelList(PixelList * pixel_list)3567 static void ResetPixelList(PixelList *pixel_list)
3568 {
3569 int
3570 level;
3571
3572 ListNode
3573 *root;
3574
3575 SkipList
3576 *list;
3577
3578 ssize_t
3579 channel;
3580
3581 /*
3582 Reset the skip-list.
3583 */
3584 for (channel=0; channel < 5; channel++)
3585 {
3586 list=pixel_list->lists+channel;
3587 root=list->nodes+65536UL;
3588 list->level=0;
3589 for (level=0; level < 9; level++)
3590 root->next[level]=65536UL;
3591 }
3592 pixel_list->seed=pixel_list->signature++;
3593 }
3594
StatisticImage(const Image * image,const StatisticType type,const size_t width,const size_t height,ExceptionInfo * exception)3595 MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
3596 const size_t width,const size_t height,ExceptionInfo *exception)
3597 {
3598 Image
3599 *statistic_image;
3600
3601 statistic_image=StatisticImageChannel(image,DefaultChannels,type,width,
3602 height,exception);
3603 return(statistic_image);
3604 }
3605
StatisticImageChannel(const Image * image,const ChannelType channel,const StatisticType type,const size_t width,const size_t height,ExceptionInfo * exception)3606 MagickExport Image *StatisticImageChannel(const Image *image,
3607 const ChannelType channel,const StatisticType type,const size_t width,
3608 const size_t height,ExceptionInfo *exception)
3609 {
3610 #define StatisticImageTag "Statistic/Image"
3611
3612 CacheView
3613 *image_view,
3614 *statistic_view;
3615
3616 Image
3617 *statistic_image;
3618
3619 MagickBooleanType
3620 status;
3621
3622 MagickOffsetType
3623 progress;
3624
3625 PixelList
3626 **magick_restrict pixel_list;
3627
3628 size_t
3629 neighbor_height,
3630 neighbor_width;
3631
3632 ssize_t
3633 y;
3634
3635 /*
3636 Initialize statistics image attributes.
3637 */
3638 assert(image != (Image *) NULL);
3639 assert(image->signature == MagickCoreSignature);
3640 if (image->debug != MagickFalse)
3641 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3642 assert(exception != (ExceptionInfo *) NULL);
3643 assert(exception->signature == MagickCoreSignature);
3644 statistic_image=CloneImage(image,0,0,MagickTrue,exception);
3645 if (statistic_image == (Image *) NULL)
3646 return((Image *) NULL);
3647 if (SetImageStorageClass(statistic_image,DirectClass) == MagickFalse)
3648 {
3649 InheritException(exception,&statistic_image->exception);
3650 statistic_image=DestroyImage(statistic_image);
3651 return((Image *) NULL);
3652 }
3653 neighbor_width=width == 0 ? GetOptimalKernelWidth2D((double) width,0.5) :
3654 width;
3655 neighbor_height=height == 0 ? GetOptimalKernelWidth2D((double) height,0.5) :
3656 height;
3657 pixel_list=AcquirePixelListThreadSet(neighbor_width,neighbor_height);
3658 if (pixel_list == (PixelList **) NULL)
3659 {
3660 statistic_image=DestroyImage(statistic_image);
3661 ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3662 }
3663 /*
3664 Make each pixel the min / max / median / mode / etc. of the neighborhood.
3665 */
3666 status=MagickTrue;
3667 progress=0;
3668 image_view=AcquireVirtualCacheView(image,exception);
3669 statistic_view=AcquireAuthenticCacheView(statistic_image,exception);
3670 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3671 #pragma omp parallel for schedule(static) shared(progress,status) \
3672 magick_number_threads(image,statistic_image,statistic_image->rows,1)
3673 #endif
3674 for (y=0; y < (ssize_t) statistic_image->rows; y++)
3675 {
3676 const int
3677 id = GetOpenMPThreadId();
3678
3679 const IndexPacket
3680 *magick_restrict indexes;
3681
3682 const PixelPacket
3683 *magick_restrict p;
3684
3685 IndexPacket
3686 *magick_restrict statistic_indexes;
3687
3688 PixelPacket
3689 *magick_restrict q;
3690
3691 ssize_t
3692 x;
3693
3694 if (status == MagickFalse)
3695 continue;
3696 p=GetCacheViewVirtualPixels(image_view,-((ssize_t) neighbor_width/2L),y-
3697 (ssize_t) (neighbor_height/2L),image->columns+neighbor_width,
3698 neighbor_height,exception);
3699 q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns, 1,exception);
3700 if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
3701 {
3702 status=MagickFalse;
3703 continue;
3704 }
3705 indexes=GetCacheViewVirtualIndexQueue(image_view);
3706 statistic_indexes=GetCacheViewAuthenticIndexQueue(statistic_view);
3707 for (x=0; x < (ssize_t) statistic_image->columns; x++)
3708 {
3709 MagickPixelPacket
3710 pixel;
3711
3712 const IndexPacket
3713 *magick_restrict s;
3714
3715 const PixelPacket
3716 *magick_restrict r;
3717
3718 ssize_t
3719 u,
3720 v;
3721
3722 r=p;
3723 s=indexes+x;
3724 ResetPixelList(pixel_list[id]);
3725 for (v=0; v < (ssize_t) neighbor_height; v++)
3726 {
3727 for (u=0; u < (ssize_t) neighbor_width; u++)
3728 InsertPixelList(image,r+u,s+u,pixel_list[id]);
3729 r+=image->columns+neighbor_width;
3730 s+=image->columns+neighbor_width;
3731 }
3732 GetMagickPixelPacket(image,&pixel);
3733 SetMagickPixelPacket(image,p+neighbor_width*neighbor_height/2,indexes+x+
3734 neighbor_width*neighbor_height/2,&pixel);
3735 switch (type)
3736 {
3737 case GradientStatistic:
3738 {
3739 MagickPixelPacket
3740 maximum,
3741 minimum;
3742
3743 GetMinimumPixelList(pixel_list[id],&pixel);
3744 minimum=pixel;
3745 GetMaximumPixelList(pixel_list[id],&pixel);
3746 maximum=pixel;
3747 pixel.red=MagickAbsoluteValue(maximum.red-minimum.red);
3748 pixel.green=MagickAbsoluteValue(maximum.green-minimum.green);
3749 pixel.blue=MagickAbsoluteValue(maximum.blue-minimum.blue);
3750 pixel.opacity=MagickAbsoluteValue(maximum.opacity-minimum.opacity);
3751 if (image->colorspace == CMYKColorspace)
3752 pixel.index=MagickAbsoluteValue(maximum.index-minimum.index);
3753 break;
3754 }
3755 case MaximumStatistic:
3756 {
3757 GetMaximumPixelList(pixel_list[id],&pixel);
3758 break;
3759 }
3760 case MeanStatistic:
3761 {
3762 GetMeanPixelList(pixel_list[id],&pixel);
3763 break;
3764 }
3765 case MedianStatistic:
3766 default:
3767 {
3768 GetMedianPixelList(pixel_list[id],&pixel);
3769 break;
3770 }
3771 case MinimumStatistic:
3772 {
3773 GetMinimumPixelList(pixel_list[id],&pixel);
3774 break;
3775 }
3776 case ModeStatistic:
3777 {
3778 GetModePixelList(pixel_list[id],&pixel);
3779 break;
3780 }
3781 case NonpeakStatistic:
3782 {
3783 GetNonpeakPixelList(pixel_list[id],&pixel);
3784 break;
3785 }
3786 case RootMeanSquareStatistic:
3787 {
3788 GetRootMeanSquarePixelList(pixel_list[id],&pixel);
3789 break;
3790 }
3791 case StandardDeviationStatistic:
3792 {
3793 GetStandardDeviationPixelList(pixel_list[id],&pixel);
3794 break;
3795 }
3796 }
3797 if ((channel & RedChannel) != 0)
3798 SetPixelRed(q,ClampToQuantum(pixel.red));
3799 if ((channel & GreenChannel) != 0)
3800 SetPixelGreen(q,ClampToQuantum(pixel.green));
3801 if ((channel & BlueChannel) != 0)
3802 SetPixelBlue(q,ClampToQuantum(pixel.blue));
3803 if ((channel & OpacityChannel) != 0)
3804 SetPixelOpacity(q,ClampToQuantum(pixel.opacity));
3805 if (((channel & IndexChannel) != 0) &&
3806 (image->colorspace == CMYKColorspace))
3807 SetPixelIndex(statistic_indexes+x,ClampToQuantum(pixel.index));
3808 p++;
3809 q++;
3810 }
3811 if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
3812 status=MagickFalse;
3813 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3814 {
3815 MagickBooleanType
3816 proceed;
3817
3818 proceed=SetImageProgress(image,StatisticImageTag,progress++,
3819 image->rows);
3820 if (proceed == MagickFalse)
3821 status=MagickFalse;
3822 }
3823 }
3824 statistic_view=DestroyCacheView(statistic_view);
3825 image_view=DestroyCacheView(image_view);
3826 pixel_list=DestroyPixelListThreadSet(pixel_list);
3827 if (status == MagickFalse)
3828 statistic_image=DestroyImage(statistic_image);
3829 return(statistic_image);
3830 }
3831