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