MagickCore 6.9.13-10
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 ? 0 : 1)
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 ? 0 : 1)
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 ? 0 : 1)
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,2)
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+=QuantumScale*GetPixelRed(p);
1460 sum_squares+=QuantumScale*GetPixelRed(p)*QuantumScale*GetPixelRed(p);
1461 sum_cubes+=QuantumScale*GetPixelRed(p)*QuantumScale*GetPixelRed(p)*
1462 QuantumScale*GetPixelRed(p);
1463 sum_fourth_power+=QuantumScale*GetPixelRed(p)*QuantumScale*
1464 GetPixelRed(p)*QuantumScale*GetPixelRed(p)*QuantumScale*
1465 GetPixelRed(p);
1466 area++;
1467 }
1468 if ((channel & GreenChannel) != 0)
1469 {
1470 mean+=QuantumScale*GetPixelGreen(p);
1471 sum_squares+=QuantumScale*GetPixelGreen(p)*QuantumScale*
1472 GetPixelGreen(p);
1473 sum_cubes+=QuantumScale*GetPixelGreen(p)*QuantumScale*
1474 GetPixelGreen(p)*QuantumScale*GetPixelGreen(p);
1475 sum_fourth_power+=QuantumScale*GetPixelGreen(p)*QuantumScale*
1476 GetPixelGreen(p)*QuantumScale*GetPixelGreen(p)*QuantumScale*
1477 GetPixelGreen(p);
1478 area++;
1479 }
1480 if ((channel & BlueChannel) != 0)
1481 {
1482 mean+=QuantumScale*GetPixelBlue(p);
1483 sum_squares+=QuantumScale*GetPixelBlue(p)*QuantumScale*
1484 GetPixelBlue(p);
1485 sum_cubes+=QuantumScale*GetPixelBlue(p)*QuantumScale*GetPixelBlue(p)*
1486 QuantumScale*GetPixelBlue(p);
1487 sum_fourth_power+=QuantumScale*GetPixelBlue(p)*QuantumScale*
1488 GetPixelBlue(p)*QuantumScale*GetPixelBlue(p)*QuantumScale*
1489 GetPixelBlue(p);
1490 area++;
1491 }
1492 if ((channel & OpacityChannel) != 0)
1493 {
1494 mean+=QuantumScale*GetPixelAlpha(p);
1495 sum_squares+=QuantumScale*GetPixelOpacity(p)*QuantumScale*
1496 GetPixelAlpha(p);
1497 sum_cubes+=QuantumScale*GetPixelOpacity(p)*QuantumScale*
1498 GetPixelAlpha(p)*QuantumScale*GetPixelAlpha(p);
1499 sum_fourth_power+=QuantumScale*GetPixelAlpha(p)*QuantumScale*
1500 GetPixelAlpha(p)*QuantumScale*GetPixelAlpha(p)*GetPixelAlpha(p);
1501 area++;
1502 }
1503 if (((channel & IndexChannel) != 0) &&
1504 (image->colorspace == CMYKColorspace))
1505 {
1506 double
1507 index;
1508
1509 index=QuantumScale*GetPixelIndex(indexes+x);
1510 mean+=index;
1511 sum_squares+=index*index;
1512 sum_cubes+=index*index*index;
1513 sum_fourth_power+=index*index*index*index;
1514 area++;
1515 }
1516 p++;
1517 }
1518 }
1519 if (y < (ssize_t) image->rows)
1520 return(MagickFalse);
1521 if (area != 0.0)
1522 {
1523 mean/=area;
1524 sum_squares/=area;
1525 sum_cubes/=area;
1526 sum_fourth_power/=area;
1527 }
1528 standard_deviation=sqrt(sum_squares-(mean*mean));
1529 if (standard_deviation != 0.0)
1530 {
1531 *kurtosis=sum_fourth_power-4.0*mean*sum_cubes+6.0*mean*mean*sum_squares-
1532 3.0*mean*mean*mean*mean;
1533 *kurtosis/=standard_deviation*standard_deviation*standard_deviation*
1534 standard_deviation;
1535 *kurtosis-=3.0;
1536 *skewness=sum_cubes-3.0*mean*sum_squares+2.0*mean*mean*mean;
1537 *skewness/=standard_deviation*standard_deviation*standard_deviation;
1538 }
1539 return(y == (ssize_t) image->rows ? MagickTrue : MagickFalse);
1540}
1541
1542/*
1543%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1544% %
1545% %
1546% %
1547% G e t I m a g e C h a n n e l M e a n %
1548% %
1549% %
1550% %
1551%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1552%
1553% GetImageChannelMean() returns the mean and standard deviation of one or more
1554% image channels.
1555%
1556% The format of the GetImageChannelMean method is:
1557%
1558% MagickBooleanType GetImageChannelMean(const Image *image,
1559% const ChannelType channel,double *mean,double *standard_deviation,
1560% ExceptionInfo *exception)
1561%
1562% A description of each parameter follows:
1563%
1564% o image: the image.
1565%
1566% o channel: the channel.
1567%
1568% o mean: the average value in the channel.
1569%
1570% o standard_deviation: the standard deviation of the channel.
1571%
1572% o exception: return any errors or warnings in this structure.
1573%
1574*/
1575
1576MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean,
1577 double *standard_deviation,ExceptionInfo *exception)
1578{
1579 MagickBooleanType
1580 status;
1581
1582 status=GetImageChannelMean(image,CompositeChannels,mean,standard_deviation,
1583 exception);
1584 return(status);
1585}
1586
1587MagickExport MagickBooleanType GetImageChannelMean(const Image *image,
1588 const ChannelType channel,double *mean,double *standard_deviation,
1589 ExceptionInfo *exception)
1590{
1592 *channel_statistics;
1593
1594 size_t
1595 channels;
1596
1597 assert(image != (Image *) NULL);
1598 assert(image->signature == MagickCoreSignature);
1599 if (IsEventLogging() != MagickFalse)
1600 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1601 channel_statistics=GetImageChannelStatistics(image,exception);
1602 if (channel_statistics == (ChannelStatistics *) NULL)
1603 return(MagickFalse);
1604 channels=0;
1605 channel_statistics[CompositeChannels].mean=0.0;
1606 channel_statistics[CompositeChannels].standard_deviation=0.0;
1607 if ((channel & RedChannel) != 0)
1608 {
1609 channel_statistics[CompositeChannels].mean+=
1610 channel_statistics[RedChannel].mean;
1611 channel_statistics[CompositeChannels].standard_deviation+=
1612 channel_statistics[RedChannel].standard_deviation;
1613 channels++;
1614 }
1615 if ((channel & GreenChannel) != 0)
1616 {
1617 channel_statistics[CompositeChannels].mean+=
1618 channel_statistics[GreenChannel].mean;
1619 channel_statistics[CompositeChannels].standard_deviation+=
1620 channel_statistics[GreenChannel].standard_deviation;
1621 channels++;
1622 }
1623 if ((channel & BlueChannel) != 0)
1624 {
1625 channel_statistics[CompositeChannels].mean+=
1626 channel_statistics[BlueChannel].mean;
1627 channel_statistics[CompositeChannels].standard_deviation+=
1628 channel_statistics[BlueChannel].standard_deviation;
1629 channels++;
1630 }
1631 if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
1632 {
1633 channel_statistics[CompositeChannels].mean+=
1634 channel_statistics[OpacityChannel].mean;
1635 channel_statistics[CompositeChannels].standard_deviation+=
1636 channel_statistics[OpacityChannel].standard_deviation;
1637 channels++;
1638 }
1639 if (((channel & IndexChannel) != 0) && (image->colorspace == CMYKColorspace))
1640 {
1641 channel_statistics[CompositeChannels].mean+=
1642 channel_statistics[BlackChannel].mean;
1643 channel_statistics[CompositeChannels].standard_deviation+=
1644 channel_statistics[CompositeChannels].standard_deviation;
1645 channels++;
1646 }
1647 channel_statistics[CompositeChannels].mean/=channels;
1648 channel_statistics[CompositeChannels].standard_deviation/=channels;
1649 *mean=channel_statistics[CompositeChannels].mean;
1650 *standard_deviation=channel_statistics[CompositeChannels].standard_deviation;
1651 channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1652 channel_statistics);
1653 return(MagickTrue);
1654}
1655
1656/*
1657%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1658% %
1659% %
1660% %
1661% G e t I m a g e C h a n n e l M o m e n t s %
1662% %
1663% %
1664% %
1665%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1666%
1667% GetImageChannelMoments() returns the normalized moments of one or more image
1668% channels.
1669%
1670% The format of the GetImageChannelMoments method is:
1671%
1672% ChannelMoments *GetImageChannelMoments(const Image *image,
1673% ExceptionInfo *exception)
1674%
1675% A description of each parameter follows:
1676%
1677% o image: the image.
1678%
1679% o exception: return any errors or warnings in this structure.
1680%
1681*/
1682MagickExport ChannelMoments *GetImageChannelMoments(const Image *image,
1683 ExceptionInfo *exception)
1684{
1685#define MaxNumberImageMoments 8
1686
1688 *channel_moments;
1689
1690 double
1691 M00[CompositeChannels+1],
1692 M01[CompositeChannels+1],
1693 M02[CompositeChannels+1],
1694 M03[CompositeChannels+1],
1695 M10[CompositeChannels+1],
1696 M11[CompositeChannels+1],
1697 M12[CompositeChannels+1],
1698 M20[CompositeChannels+1],
1699 M21[CompositeChannels+1],
1700 M22[CompositeChannels+1],
1701 M30[CompositeChannels+1];
1702
1704 pixel;
1705
1706 PointInfo
1707 centroid[CompositeChannels+1];
1708
1709 ssize_t
1710 channel,
1711 channels,
1712 y;
1713
1714 size_t
1715 length;
1716
1717 assert(image != (Image *) NULL);
1718 assert(image->signature == MagickCoreSignature);
1719 if (IsEventLogging() != MagickFalse)
1720 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1721 length=CompositeChannels+1UL;
1722 channel_moments=(ChannelMoments *) AcquireQuantumMemory(length,
1723 sizeof(*channel_moments));
1724 if (channel_moments == (ChannelMoments *) NULL)
1725 return(channel_moments);
1726 (void) memset(channel_moments,0,length*sizeof(*channel_moments));
1727 (void) memset(centroid,0,sizeof(centroid));
1728 (void) memset(M00,0,sizeof(M00));
1729 (void) memset(M01,0,sizeof(M01));
1730 (void) memset(M02,0,sizeof(M02));
1731 (void) memset(M03,0,sizeof(M03));
1732 (void) memset(M10,0,sizeof(M10));
1733 (void) memset(M11,0,sizeof(M11));
1734 (void) memset(M12,0,sizeof(M12));
1735 (void) memset(M20,0,sizeof(M20));
1736 (void) memset(M21,0,sizeof(M21));
1737 (void) memset(M22,0,sizeof(M22));
1738 (void) memset(M30,0,sizeof(M30));
1739 GetMagickPixelPacket(image,&pixel);
1740 for (y=0; y < (ssize_t) image->rows; y++)
1741 {
1742 const IndexPacket
1743 *magick_restrict indexes;
1744
1745 const PixelPacket
1746 *magick_restrict p;
1747
1748 ssize_t
1749 x;
1750
1751 /*
1752 Compute center of mass (centroid).
1753 */
1754 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1755 if (p == (const PixelPacket *) NULL)
1756 break;
1757 indexes=GetVirtualIndexQueue(image);
1758 for (x=0; x < (ssize_t) image->columns; x++)
1759 {
1760 SetMagickPixelPacket(image,p,indexes+x,&pixel);
1761 M00[RedChannel]+=QuantumScale*pixel.red;
1762 M10[RedChannel]+=x*QuantumScale*pixel.red;
1763 M01[RedChannel]+=y*QuantumScale*pixel.red;
1764 M00[GreenChannel]+=QuantumScale*pixel.green;
1765 M10[GreenChannel]+=x*QuantumScale*pixel.green;
1766 M01[GreenChannel]+=y*QuantumScale*pixel.green;
1767 M00[BlueChannel]+=QuantumScale*pixel.blue;
1768 M10[BlueChannel]+=x*QuantumScale*pixel.blue;
1769 M01[BlueChannel]+=y*QuantumScale*pixel.blue;
1770 if (image->matte != MagickFalse)
1771 {
1772 M00[OpacityChannel]+=QuantumScale*pixel.opacity;
1773 M10[OpacityChannel]+=x*QuantumScale*pixel.opacity;
1774 M01[OpacityChannel]+=y*QuantumScale*pixel.opacity;
1775 }
1776 if (image->colorspace == CMYKColorspace)
1777 {
1778 M00[IndexChannel]+=QuantumScale*pixel.index;
1779 M10[IndexChannel]+=x*QuantumScale*pixel.index;
1780 M01[IndexChannel]+=y*QuantumScale*pixel.index;
1781 }
1782 p++;
1783 }
1784 }
1785 for (channel=0; channel <= CompositeChannels; channel++)
1786 {
1787 /*
1788 Compute center of mass (centroid).
1789 */
1790 if (M00[channel] < MagickEpsilon)
1791 {
1792 M00[channel]+=MagickEpsilon;
1793 centroid[channel].x=(double) image->columns/2.0;
1794 centroid[channel].y=(double) image->rows/2.0;
1795 continue;
1796 }
1797 M00[channel]+=MagickEpsilon;
1798 centroid[channel].x=M10[channel]/M00[channel];
1799 centroid[channel].y=M01[channel]/M00[channel];
1800 }
1801 for (y=0; y < (ssize_t) image->rows; y++)
1802 {
1803 const IndexPacket
1804 *magick_restrict indexes;
1805
1806 const PixelPacket
1807 *magick_restrict p;
1808
1809 ssize_t
1810 x;
1811
1812 /*
1813 Compute the image moments.
1814 */
1815 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1816 if (p == (const PixelPacket *) NULL)
1817 break;
1818 indexes=GetVirtualIndexQueue(image);
1819 for (x=0; x < (ssize_t) image->columns; x++)
1820 {
1821 SetMagickPixelPacket(image,p,indexes+x,&pixel);
1822 M11[RedChannel]+=(x-centroid[RedChannel].x)*(y-
1823 centroid[RedChannel].y)*QuantumScale*pixel.red;
1824 M20[RedChannel]+=(x-centroid[RedChannel].x)*(x-
1825 centroid[RedChannel].x)*QuantumScale*pixel.red;
1826 M02[RedChannel]+=(y-centroid[RedChannel].y)*(y-
1827 centroid[RedChannel].y)*QuantumScale*pixel.red;
1828 M21[RedChannel]+=(x-centroid[RedChannel].x)*(x-
1829 centroid[RedChannel].x)*(y-centroid[RedChannel].y)*QuantumScale*
1830 pixel.red;
1831 M12[RedChannel]+=(x-centroid[RedChannel].x)*(y-
1832 centroid[RedChannel].y)*(y-centroid[RedChannel].y)*QuantumScale*
1833 pixel.red;
1834 M22[RedChannel]+=(x-centroid[RedChannel].x)*(x-
1835 centroid[RedChannel].x)*(y-centroid[RedChannel].y)*(y-
1836 centroid[RedChannel].y)*QuantumScale*pixel.red;
1837 M30[RedChannel]+=(x-centroid[RedChannel].x)*(x-
1838 centroid[RedChannel].x)*(x-centroid[RedChannel].x)*QuantumScale*
1839 pixel.red;
1840 M03[RedChannel]+=(y-centroid[RedChannel].y)*(y-
1841 centroid[RedChannel].y)*(y-centroid[RedChannel].y)*QuantumScale*
1842 pixel.red;
1843 M11[GreenChannel]+=(x-centroid[GreenChannel].x)*(y-
1844 centroid[GreenChannel].y)*QuantumScale*pixel.green;
1845 M20[GreenChannel]+=(x-centroid[GreenChannel].x)*(x-
1846 centroid[GreenChannel].x)*QuantumScale*pixel.green;
1847 M02[GreenChannel]+=(y-centroid[GreenChannel].y)*(y-
1848 centroid[GreenChannel].y)*QuantumScale*pixel.green;
1849 M21[GreenChannel]+=(x-centroid[GreenChannel].x)*(x-
1850 centroid[GreenChannel].x)*(y-centroid[GreenChannel].y)*QuantumScale*
1851 pixel.green;
1852 M12[GreenChannel]+=(x-centroid[GreenChannel].x)*(y-
1853 centroid[GreenChannel].y)*(y-centroid[GreenChannel].y)*QuantumScale*
1854 pixel.green;
1855 M22[GreenChannel]+=(x-centroid[GreenChannel].x)*(x-
1856 centroid[GreenChannel].x)*(y-centroid[GreenChannel].y)*(y-
1857 centroid[GreenChannel].y)*QuantumScale*pixel.green;
1858 M30[GreenChannel]+=(x-centroid[GreenChannel].x)*(x-
1859 centroid[GreenChannel].x)*(x-centroid[GreenChannel].x)*QuantumScale*
1860 pixel.green;
1861 M03[GreenChannel]+=(y-centroid[GreenChannel].y)*(y-
1862 centroid[GreenChannel].y)*(y-centroid[GreenChannel].y)*QuantumScale*
1863 pixel.green;
1864 M11[BlueChannel]+=(x-centroid[BlueChannel].x)*(y-
1865 centroid[BlueChannel].y)*QuantumScale*pixel.blue;
1866 M20[BlueChannel]+=(x-centroid[BlueChannel].x)*(x-
1867 centroid[BlueChannel].x)*QuantumScale*pixel.blue;
1868 M02[BlueChannel]+=(y-centroid[BlueChannel].y)*(y-
1869 centroid[BlueChannel].y)*QuantumScale*pixel.blue;
1870 M21[BlueChannel]+=(x-centroid[BlueChannel].x)*(x-
1871 centroid[BlueChannel].x)*(y-centroid[BlueChannel].y)*QuantumScale*
1872 pixel.blue;
1873 M12[BlueChannel]+=(x-centroid[BlueChannel].x)*(y-
1874 centroid[BlueChannel].y)*(y-centroid[BlueChannel].y)*QuantumScale*
1875 pixel.blue;
1876 M22[BlueChannel]+=(x-centroid[BlueChannel].x)*(x-
1877 centroid[BlueChannel].x)*(y-centroid[BlueChannel].y)*(y-
1878 centroid[BlueChannel].y)*QuantumScale*pixel.blue;
1879 M30[BlueChannel]+=(x-centroid[BlueChannel].x)*(x-
1880 centroid[BlueChannel].x)*(x-centroid[BlueChannel].x)*QuantumScale*
1881 pixel.blue;
1882 M03[BlueChannel]+=(y-centroid[BlueChannel].y)*(y-
1883 centroid[BlueChannel].y)*(y-centroid[BlueChannel].y)*QuantumScale*
1884 pixel.blue;
1885 if (image->matte != MagickFalse)
1886 {
1887 M11[OpacityChannel]+=(x-centroid[OpacityChannel].x)*(y-
1888 centroid[OpacityChannel].y)*QuantumScale*pixel.opacity;
1889 M20[OpacityChannel]+=(x-centroid[OpacityChannel].x)*(x-
1890 centroid[OpacityChannel].x)*QuantumScale*pixel.opacity;
1891 M02[OpacityChannel]+=(y-centroid[OpacityChannel].y)*(y-
1892 centroid[OpacityChannel].y)*QuantumScale*pixel.opacity;
1893 M21[OpacityChannel]+=(x-centroid[OpacityChannel].x)*(x-
1894 centroid[OpacityChannel].x)*(y-centroid[OpacityChannel].y)*
1895 QuantumScale*pixel.opacity;
1896 M12[OpacityChannel]+=(x-centroid[OpacityChannel].x)*(y-
1897 centroid[OpacityChannel].y)*(y-centroid[OpacityChannel].y)*
1898 QuantumScale*pixel.opacity;
1899 M22[OpacityChannel]+=(x-centroid[OpacityChannel].x)*(x-
1900 centroid[OpacityChannel].x)*(y-centroid[OpacityChannel].y)*(y-
1901 centroid[OpacityChannel].y)*QuantumScale*pixel.opacity;
1902 M30[OpacityChannel]+=(x-centroid[OpacityChannel].x)*(x-
1903 centroid[OpacityChannel].x)*(x-centroid[OpacityChannel].x)*
1904 QuantumScale*pixel.opacity;
1905 M03[OpacityChannel]+=(y-centroid[OpacityChannel].y)*(y-
1906 centroid[OpacityChannel].y)*(y-centroid[OpacityChannel].y)*
1907 QuantumScale*pixel.opacity;
1908 }
1909 if (image->colorspace == CMYKColorspace)
1910 {
1911 M11[IndexChannel]+=(x-centroid[IndexChannel].x)*(y-
1912 centroid[IndexChannel].y)*QuantumScale*pixel.index;
1913 M20[IndexChannel]+=(x-centroid[IndexChannel].x)*(x-
1914 centroid[IndexChannel].x)*QuantumScale*pixel.index;
1915 M02[IndexChannel]+=(y-centroid[IndexChannel].y)*(y-
1916 centroid[IndexChannel].y)*QuantumScale*pixel.index;
1917 M21[IndexChannel]+=(x-centroid[IndexChannel].x)*(x-
1918 centroid[IndexChannel].x)*(y-centroid[IndexChannel].y)*
1919 QuantumScale*pixel.index;
1920 M12[IndexChannel]+=(x-centroid[IndexChannel].x)*(y-
1921 centroid[IndexChannel].y)*(y-centroid[IndexChannel].y)*
1922 QuantumScale*pixel.index;
1923 M22[IndexChannel]+=(x-centroid[IndexChannel].x)*(x-
1924 centroid[IndexChannel].x)*(y-centroid[IndexChannel].y)*(y-
1925 centroid[IndexChannel].y)*QuantumScale*pixel.index;
1926 M30[IndexChannel]+=(x-centroid[IndexChannel].x)*(x-
1927 centroid[IndexChannel].x)*(x-centroid[IndexChannel].x)*
1928 QuantumScale*pixel.index;
1929 M03[IndexChannel]+=(y-centroid[IndexChannel].y)*(y-
1930 centroid[IndexChannel].y)*(y-centroid[IndexChannel].y)*
1931 QuantumScale*pixel.index;
1932 }
1933 p++;
1934 }
1935 }
1936 channels=3;
1937 M00[CompositeChannels]+=(M00[RedChannel]+M00[GreenChannel]+M00[BlueChannel]);
1938 M01[CompositeChannels]+=(M01[RedChannel]+M01[GreenChannel]+M01[BlueChannel]);
1939 M02[CompositeChannels]+=(M02[RedChannel]+M02[GreenChannel]+M02[BlueChannel]);
1940 M03[CompositeChannels]+=(M03[RedChannel]+M03[GreenChannel]+M03[BlueChannel]);
1941 M10[CompositeChannels]+=(M10[RedChannel]+M10[GreenChannel]+M10[BlueChannel]);
1942 M11[CompositeChannels]+=(M11[RedChannel]+M11[GreenChannel]+M11[BlueChannel]);
1943 M12[CompositeChannels]+=(M12[RedChannel]+M12[GreenChannel]+M12[BlueChannel]);
1944 M20[CompositeChannels]+=(M20[RedChannel]+M20[GreenChannel]+M20[BlueChannel]);
1945 M21[CompositeChannels]+=(M21[RedChannel]+M21[GreenChannel]+M21[BlueChannel]);
1946 M22[CompositeChannels]+=(M22[RedChannel]+M22[GreenChannel]+M22[BlueChannel]);
1947 M30[CompositeChannels]+=(M30[RedChannel]+M30[GreenChannel]+M30[BlueChannel]);
1948 if (image->matte != MagickFalse)
1949 {
1950 channels+=1;
1951 M00[CompositeChannels]+=M00[OpacityChannel];
1952 M01[CompositeChannels]+=M01[OpacityChannel];
1953 M02[CompositeChannels]+=M02[OpacityChannel];
1954 M03[CompositeChannels]+=M03[OpacityChannel];
1955 M10[CompositeChannels]+=M10[OpacityChannel];
1956 M11[CompositeChannels]+=M11[OpacityChannel];
1957 M12[CompositeChannels]+=M12[OpacityChannel];
1958 M20[CompositeChannels]+=M20[OpacityChannel];
1959 M21[CompositeChannels]+=M21[OpacityChannel];
1960 M22[CompositeChannels]+=M22[OpacityChannel];
1961 M30[CompositeChannels]+=M30[OpacityChannel];
1962 }
1963 if (image->colorspace == CMYKColorspace)
1964 {
1965 channels+=1;
1966 M00[CompositeChannels]+=M00[IndexChannel];
1967 M01[CompositeChannels]+=M01[IndexChannel];
1968 M02[CompositeChannels]+=M02[IndexChannel];
1969 M03[CompositeChannels]+=M03[IndexChannel];
1970 M10[CompositeChannels]+=M10[IndexChannel];
1971 M11[CompositeChannels]+=M11[IndexChannel];
1972 M12[CompositeChannels]+=M12[IndexChannel];
1973 M20[CompositeChannels]+=M20[IndexChannel];
1974 M21[CompositeChannels]+=M21[IndexChannel];
1975 M22[CompositeChannels]+=M22[IndexChannel];
1976 M30[CompositeChannels]+=M30[IndexChannel];
1977 }
1978 M00[CompositeChannels]/=(double) channels;
1979 M01[CompositeChannels]/=(double) channels;
1980 M02[CompositeChannels]/=(double) channels;
1981 M03[CompositeChannels]/=(double) channels;
1982 M10[CompositeChannels]/=(double) channels;
1983 M11[CompositeChannels]/=(double) channels;
1984 M12[CompositeChannels]/=(double) channels;
1985 M20[CompositeChannels]/=(double) channels;
1986 M21[CompositeChannels]/=(double) channels;
1987 M22[CompositeChannels]/=(double) channels;
1988 M30[CompositeChannels]/=(double) channels;
1989 for (channel=0; channel <= CompositeChannels; channel++)
1990 {
1991 /*
1992 Compute elliptical angle, major and minor axes, eccentricity, & intensity.
1993 */
1994 channel_moments[channel].centroid=centroid[channel];
1995 channel_moments[channel].ellipse_axis.x=sqrt((2.0*
1996 PerceptibleReciprocal(M00[channel]))*((M20[channel]+M02[channel])+
1997 sqrt(4.0*M11[channel]*M11[channel]+(M20[channel]-M02[channel])*
1998 (M20[channel]-M02[channel]))));
1999 channel_moments[channel].ellipse_axis.y=sqrt((2.0*
2000 PerceptibleReciprocal(M00[channel]))*((M20[channel]+M02[channel])-
2001 sqrt(4.0*M11[channel]*M11[channel]+(M20[channel]-M02[channel])*
2002 (M20[channel]-M02[channel]))));
2003 channel_moments[channel].ellipse_angle=RadiansToDegrees(1.0/2.0*atan(2.0*
2004 M11[channel]*PerceptibleReciprocal(M20[channel]-M02[channel])));
2005 if (fabs(M11[channel]) < 0.0)
2006 {
2007 if ((fabs(M20[channel]-M02[channel]) >= 0.0) &&
2008 ((M20[channel]-M02[channel]) < 0.0))
2009 channel_moments[channel].ellipse_angle+=90.0;
2010 }
2011 else
2012 if (M11[channel] < 0.0)
2013 {
2014 if (fabs(M20[channel]-M02[channel]) >= 0.0)
2015 {
2016 if ((M20[channel]-M02[channel]) < 0.0)
2017 channel_moments[channel].ellipse_angle+=90.0;
2018 else
2019 channel_moments[channel].ellipse_angle+=180.0;
2020 }
2021 }
2022 else
2023 if ((fabs(M20[channel]-M02[channel]) >= 0.0) &&
2024 ((M20[channel]-M02[channel]) < 0.0))
2025 channel_moments[channel].ellipse_angle+=90.0;
2026 channel_moments[channel].ellipse_eccentricity=sqrt(1.0-(
2027 channel_moments[channel].ellipse_axis.y*
2028 channel_moments[channel].ellipse_axis.y*PerceptibleReciprocal(
2029 channel_moments[channel].ellipse_axis.x*
2030 channel_moments[channel].ellipse_axis.x)));
2031 channel_moments[channel].ellipse_intensity=M00[channel]/
2032 (MagickPI*channel_moments[channel].ellipse_axis.x*
2033 channel_moments[channel].ellipse_axis.y+MagickEpsilon);
2034 }
2035 for (channel=0; channel <= CompositeChannels; channel++)
2036 {
2037 /*
2038 Normalize image moments.
2039 */
2040 M10[channel]=0.0;
2041 M01[channel]=0.0;
2042 M11[channel]/=pow(M00[channel],1.0+(1.0+1.0)/2.0);
2043 M20[channel]/=pow(M00[channel],1.0+(2.0+0.0)/2.0);
2044 M02[channel]/=pow(M00[channel],1.0+(0.0+2.0)/2.0);
2045 M21[channel]/=pow(M00[channel],1.0+(2.0+1.0)/2.0);
2046 M12[channel]/=pow(M00[channel],1.0+(1.0+2.0)/2.0);
2047 M22[channel]/=pow(M00[channel],1.0+(2.0+2.0)/2.0);
2048 M30[channel]/=pow(M00[channel],1.0+(3.0+0.0)/2.0);
2049 M03[channel]/=pow(M00[channel],1.0+(0.0+3.0)/2.0);
2050 M00[channel]=1.0;
2051 }
2052 for (channel=0; channel <= CompositeChannels; channel++)
2053 {
2054 /*
2055 Compute Hu invariant moments.
2056 */
2057 channel_moments[channel].I[0]=M20[channel]+M02[channel];
2058 channel_moments[channel].I[1]=(M20[channel]-M02[channel])*
2059 (M20[channel]-M02[channel])+4.0*M11[channel]*M11[channel];
2060 channel_moments[channel].I[2]=(M30[channel]-3.0*M12[channel])*
2061 (M30[channel]-3.0*M12[channel])+(3.0*M21[channel]-M03[channel])*
2062 (3.0*M21[channel]-M03[channel]);
2063 channel_moments[channel].I[3]=(M30[channel]+M12[channel])*
2064 (M30[channel]+M12[channel])+(M21[channel]+M03[channel])*
2065 (M21[channel]+M03[channel]);
2066 channel_moments[channel].I[4]=(M30[channel]-3.0*M12[channel])*
2067 (M30[channel]+M12[channel])*((M30[channel]+M12[channel])*
2068 (M30[channel]+M12[channel])-3.0*(M21[channel]+M03[channel])*
2069 (M21[channel]+M03[channel]))+(3.0*M21[channel]-M03[channel])*
2070 (M21[channel]+M03[channel])*(3.0*(M30[channel]+M12[channel])*
2071 (M30[channel]+M12[channel])-(M21[channel]+M03[channel])*
2072 (M21[channel]+M03[channel]));
2073 channel_moments[channel].I[5]=(M20[channel]-M02[channel])*
2074 ((M30[channel]+M12[channel])*(M30[channel]+M12[channel])-
2075 (M21[channel]+M03[channel])*(M21[channel]+M03[channel]))+
2076 4.0*M11[channel]*(M30[channel]+M12[channel])*(M21[channel]+M03[channel]);
2077 channel_moments[channel].I[6]=(3.0*M21[channel]-M03[channel])*
2078 (M30[channel]+M12[channel])*((M30[channel]+M12[channel])*
2079 (M30[channel]+M12[channel])-3.0*(M21[channel]+M03[channel])*
2080 (M21[channel]+M03[channel]))-(M30[channel]-3*M12[channel])*
2081 (M21[channel]+M03[channel])*(3.0*(M30[channel]+M12[channel])*
2082 (M30[channel]+M12[channel])-(M21[channel]+M03[channel])*
2083 (M21[channel]+M03[channel]));
2084 channel_moments[channel].I[7]=M11[channel]*((M30[channel]+M12[channel])*
2085 (M30[channel]+M12[channel])-(M03[channel]+M21[channel])*
2086 (M03[channel]+M21[channel]))-(M20[channel]-M02[channel])*
2087 (M30[channel]+M12[channel])*(M03[channel]+M21[channel]);
2088 }
2089 if (y < (ssize_t) image->rows)
2090 channel_moments=(ChannelMoments *) RelinquishMagickMemory(channel_moments);
2091 return(channel_moments);
2092}
2093
2094/*
2095%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2096% %
2097% %
2098% %
2099% 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 %
2100% %
2101% %
2102% %
2103%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2104%
2105% GetImageChannelPerceptualHash() returns the perceptual hash of one or more
2106% image channels.
2107%
2108% The format of the GetImageChannelPerceptualHash method is:
2109%
2110% ChannelPerceptualHash *GetImageChannelPerceptualHash(const Image *image,
2111% ExceptionInfo *exception)
2112%
2113% A description of each parameter follows:
2114%
2115% o image: the image.
2116%
2117% o exception: return any errors or warnings in this structure.
2118%
2119*/
2120MagickExport ChannelPerceptualHash *GetImageChannelPerceptualHash(
2121 const Image *image,ExceptionInfo *exception)
2122{
2124 *moments;
2125
2127 *perceptual_hash;
2128
2129 Image
2130 *hash_image;
2131
2132 MagickBooleanType
2133 status;
2134
2135 ssize_t
2136 i;
2137
2138 ssize_t
2139 channel;
2140
2141 /*
2142 Blur then transform to xyY colorspace.
2143 */
2144 hash_image=BlurImage(image,0.0,1.0,exception);
2145 if (hash_image == (Image *) NULL)
2146 return((ChannelPerceptualHash *) NULL);
2147 hash_image->depth=8;
2148 status=TransformImageColorspace(hash_image,xyYColorspace);
2149 if (status == MagickFalse)
2150 return((ChannelPerceptualHash *) NULL);
2151 moments=GetImageChannelMoments(hash_image,exception);
2152 hash_image=DestroyImage(hash_image);
2153 if (moments == (ChannelMoments *) NULL)
2154 return((ChannelPerceptualHash *) NULL);
2155 perceptual_hash=(ChannelPerceptualHash *) AcquireQuantumMemory(
2156 CompositeChannels+1UL,sizeof(*perceptual_hash));
2157 if (perceptual_hash == (ChannelPerceptualHash *) NULL)
2158 return((ChannelPerceptualHash *) NULL);
2159 for (channel=0; channel <= CompositeChannels; channel++)
2160 for (i=0; i < MaximumNumberOfImageMoments; i++)
2161 perceptual_hash[channel].P[i]=(-MagickLog10(moments[channel].I[i]));
2162 moments=(ChannelMoments *) RelinquishMagickMemory(moments);
2163 /*
2164 Blur then transform to HCLp colorspace.
2165 */
2166 hash_image=BlurImage(image,0.0,1.0,exception);
2167 if (hash_image == (Image *) NULL)
2168 {
2169 perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
2170 perceptual_hash);
2171 return((ChannelPerceptualHash *) NULL);
2172 }
2173 hash_image->depth=8;
2174 status=TransformImageColorspace(hash_image,HSBColorspace);
2175 if (status == MagickFalse)
2176 {
2177 perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
2178 perceptual_hash);
2179 return((ChannelPerceptualHash *) NULL);
2180 }
2181 moments=GetImageChannelMoments(hash_image,exception);
2182 hash_image=DestroyImage(hash_image);
2183 if (moments == (ChannelMoments *) NULL)
2184 {
2185 perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
2186 perceptual_hash);
2187 return((ChannelPerceptualHash *) NULL);
2188 }
2189 for (channel=0; channel <= CompositeChannels; channel++)
2190 for (i=0; i < MaximumNumberOfImageMoments; i++)
2191 perceptual_hash[channel].Q[i]=(-MagickLog10(moments[channel].I[i]));
2192 moments=(ChannelMoments *) RelinquishMagickMemory(moments);
2193 return(perceptual_hash);
2194}
2195
2196/*
2197%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2198% %
2199% %
2200% %
2201% G e t I m a g e C h a n n e l R a n g e %
2202% %
2203% %
2204% %
2205%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2206%
2207% GetImageChannelRange() returns the range of one or more image channels.
2208%
2209% The format of the GetImageChannelRange method is:
2210%
2211% MagickBooleanType GetImageChannelRange(const Image *image,
2212% const ChannelType channel,double *minima,double *maxima,
2213% ExceptionInfo *exception)
2214%
2215% A description of each parameter follows:
2216%
2217% o image: the image.
2218%
2219% o channel: the channel.
2220%
2221% o minima: the minimum value in the channel.
2222%
2223% o maxima: the maximum value in the channel.
2224%
2225% o exception: return any errors or warnings in this structure.
2226%
2227*/
2228
2229MagickExport MagickBooleanType GetImageRange(const Image *image,
2230 double *minima,double *maxima,ExceptionInfo *exception)
2231{
2232 return(GetImageChannelRange(image,CompositeChannels,minima,maxima,exception));
2233}
2234
2235MagickExport MagickBooleanType GetImageChannelRange(const Image *image,
2236 const ChannelType channel,double *minima,double *maxima,
2237 ExceptionInfo *exception)
2238{
2240 pixel;
2241
2242 ssize_t
2243 y;
2244
2245 assert(image != (Image *) NULL);
2246 assert(image->signature == MagickCoreSignature);
2247 if (IsEventLogging() != MagickFalse)
2248 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2249 *maxima=(-MagickMaximumValue);
2250 *minima=MagickMaximumValue;
2251 GetMagickPixelPacket(image,&pixel);
2252 for (y=0; y < (ssize_t) image->rows; y++)
2253 {
2254 const IndexPacket
2255 *magick_restrict indexes;
2256
2257 const PixelPacket
2258 *magick_restrict p;
2259
2260 ssize_t
2261 x;
2262
2263 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
2264 if (p == (const PixelPacket *) NULL)
2265 break;
2266 indexes=GetVirtualIndexQueue(image);
2267 for (x=0; x < (ssize_t) image->columns; x++)
2268 {
2269 SetMagickPixelPacket(image,p,indexes+x,&pixel);
2270 if ((channel & RedChannel) != 0)
2271 {
2272 if (pixel.red < *minima)
2273 *minima=(double) pixel.red;
2274 if (pixel.red > *maxima)
2275 *maxima=(double) pixel.red;
2276 }
2277 if ((channel & GreenChannel) != 0)
2278 {
2279 if (pixel.green < *minima)
2280 *minima=(double) pixel.green;
2281 if (pixel.green > *maxima)
2282 *maxima=(double) pixel.green;
2283 }
2284 if ((channel & BlueChannel) != 0)
2285 {
2286 if (pixel.blue < *minima)
2287 *minima=(double) pixel.blue;
2288 if (pixel.blue > *maxima)
2289 *maxima=(double) pixel.blue;
2290 }
2291 if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
2292 {
2293 if (((MagickRealType) QuantumRange-(MagickRealType) pixel.opacity) < *minima)
2294 *minima=(double) ((MagickRealType) QuantumRange-(MagickRealType)
2295 pixel.opacity);
2296 if (((MagickRealType) QuantumRange-(MagickRealType) pixel.opacity) > *maxima)
2297 *maxima=(double) ((MagickRealType) QuantumRange-(MagickRealType)
2298 pixel.opacity);
2299 }
2300 if (((channel & IndexChannel) != 0) &&
2301 (image->colorspace == CMYKColorspace))
2302 {
2303 if ((double) pixel.index < *minima)
2304 *minima=(double) pixel.index;
2305 if ((double) pixel.index > *maxima)
2306 *maxima=(double) pixel.index;
2307 }
2308 p++;
2309 }
2310 }
2311 return(y == (ssize_t) image->rows ? MagickTrue : MagickFalse);
2312}
2313
2314/*
2315%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2316% %
2317% %
2318% %
2319% 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 %
2320% %
2321% %
2322% %
2323%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2324%
2325% GetImageChannelStatistics() returns statistics for each channel in the
2326% image. The statistics include the channel depth, its minima, maxima, mean,
2327% standard deviation, kurtosis and skewness. You can access the red channel
2328% mean, for example, like this:
2329%
2330% channel_statistics=GetImageChannelStatistics(image,exception);
2331% red_mean=channel_statistics[RedChannel].mean;
2332%
2333% Use MagickRelinquishMemory() to free the statistics buffer.
2334%
2335% The format of the GetImageChannelStatistics method is:
2336%
2337% ChannelStatistics *GetImageChannelStatistics(const Image *image,
2338% ExceptionInfo *exception)
2339%
2340% A description of each parameter follows:
2341%
2342% o image: the image.
2343%
2344% o exception: return any errors or warnings in this structure.
2345%
2346*/
2347MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
2348 ExceptionInfo *exception)
2349{
2351 *channel_statistics;
2352
2353 double
2354 area,
2355 standard_deviation;
2356
2358 number_bins,
2359 *histogram;
2360
2361 QuantumAny
2362 range;
2363
2364 size_t
2365 channels,
2366 depth,
2367 length;
2368
2369 ssize_t
2370 i,
2371 y;
2372
2373 assert(image != (Image *) NULL);
2374 assert(image->signature == MagickCoreSignature);
2375 if (IsEventLogging() != MagickFalse)
2376 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2377 length=CompositeChannels+1UL;
2378 channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(length,
2379 sizeof(*channel_statistics));
2380 histogram=(MagickPixelPacket *) AcquireQuantumMemory(MaxMap+1U,
2381 sizeof(*histogram));
2382 if ((channel_statistics == (ChannelStatistics *) NULL) ||
2383 (histogram == (MagickPixelPacket *) NULL))
2384 {
2385 if (histogram != (MagickPixelPacket *) NULL)
2386 histogram=(MagickPixelPacket *) RelinquishMagickMemory(histogram);
2387 if (channel_statistics != (ChannelStatistics *) NULL)
2388 channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
2389 channel_statistics);
2390 return(channel_statistics);
2391 }
2392 (void) memset(channel_statistics,0,length*
2393 sizeof(*channel_statistics));
2394 for (i=0; i <= (ssize_t) CompositeChannels; i++)
2395 {
2396 channel_statistics[i].depth=1;
2397 channel_statistics[i].maxima=(-MagickMaximumValue);
2398 channel_statistics[i].minima=MagickMaximumValue;
2399 }
2400 (void) memset(histogram,0,(MaxMap+1U)*sizeof(*histogram));
2401 (void) memset(&number_bins,0,sizeof(number_bins));
2402 for (y=0; y < (ssize_t) image->rows; y++)
2403 {
2404 const IndexPacket
2405 *magick_restrict indexes;
2406
2407 const PixelPacket
2408 *magick_restrict p;
2409
2410 ssize_t
2411 x;
2412
2413 /*
2414 Compute pixel statistics.
2415 */
2416 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
2417 if (p == (const PixelPacket *) NULL)
2418 break;
2419 indexes=GetVirtualIndexQueue(image);
2420 for (x=0; x < (ssize_t) image->columns; )
2421 {
2422 if (channel_statistics[RedChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2423 {
2424 depth=channel_statistics[RedChannel].depth;
2425 range=GetQuantumRange(depth);
2426 if (IsPixelAtDepth(GetPixelRed(p),range) == MagickFalse)
2427 {
2428 channel_statistics[RedChannel].depth++;
2429 continue;
2430 }
2431 }
2432 if (channel_statistics[GreenChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2433 {
2434 depth=channel_statistics[GreenChannel].depth;
2435 range=GetQuantumRange(depth);
2436 if (IsPixelAtDepth(GetPixelGreen(p),range) == MagickFalse)
2437 {
2438 channel_statistics[GreenChannel].depth++;
2439 continue;
2440 }
2441 }
2442 if (channel_statistics[BlueChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2443 {
2444 depth=channel_statistics[BlueChannel].depth;
2445 range=GetQuantumRange(depth);
2446 if (IsPixelAtDepth(GetPixelBlue(p),range) == MagickFalse)
2447 {
2448 channel_statistics[BlueChannel].depth++;
2449 continue;
2450 }
2451 }
2452 if (image->matte != MagickFalse)
2453 {
2454 if (channel_statistics[OpacityChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2455 {
2456 depth=channel_statistics[OpacityChannel].depth;
2457 range=GetQuantumRange(depth);
2458 if (IsPixelAtDepth(GetPixelAlpha(p),range) == MagickFalse)
2459 {
2460 channel_statistics[OpacityChannel].depth++;
2461 continue;
2462 }
2463 }
2464 }
2465 if (image->colorspace == CMYKColorspace)
2466 {
2467 if (channel_statistics[BlackChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2468 {
2469 depth=channel_statistics[BlackChannel].depth;
2470 range=GetQuantumRange(depth);
2471 if (IsPixelAtDepth(GetPixelIndex(indexes+x),range) == MagickFalse)
2472 {
2473 channel_statistics[BlackChannel].depth++;
2474 continue;
2475 }
2476 }
2477 }
2478 if ((double) GetPixelRed(p) < channel_statistics[RedChannel].minima)
2479 channel_statistics[RedChannel].minima=(double) GetPixelRed(p);
2480 if ((double) GetPixelRed(p) > channel_statistics[RedChannel].maxima)
2481 channel_statistics[RedChannel].maxima=(double) GetPixelRed(p);
2482 channel_statistics[RedChannel].sum+=QuantumScale*GetPixelRed(p);
2483 channel_statistics[RedChannel].sum_squared+=QuantumScale*GetPixelRed(p)*
2484 QuantumScale*GetPixelRed(p);
2485 channel_statistics[RedChannel].sum_cubed+=QuantumScale*GetPixelRed(p)*
2486 QuantumScale*GetPixelRed(p)*QuantumScale*GetPixelRed(p);
2487 channel_statistics[RedChannel].sum_fourth_power+=QuantumScale*
2488 GetPixelRed(p)*QuantumScale*GetPixelRed(p)*QuantumScale*GetPixelRed(p)*
2489 QuantumScale*GetPixelRed(p);
2490 if ((double) GetPixelGreen(p) < channel_statistics[GreenChannel].minima)
2491 channel_statistics[GreenChannel].minima=(double) GetPixelGreen(p);
2492 if ((double) GetPixelGreen(p) > channel_statistics[GreenChannel].maxima)
2493 channel_statistics[GreenChannel].maxima=(double) GetPixelGreen(p);
2494 channel_statistics[GreenChannel].sum+=QuantumScale*GetPixelGreen(p);
2495 channel_statistics[GreenChannel].sum_squared+=QuantumScale*GetPixelGreen(p)*
2496 QuantumScale*GetPixelGreen(p);
2497 channel_statistics[GreenChannel].sum_cubed+=QuantumScale*GetPixelGreen(p)*
2498 QuantumScale*GetPixelGreen(p)*QuantumScale*GetPixelGreen(p);
2499 channel_statistics[GreenChannel].sum_fourth_power+=QuantumScale*
2500 GetPixelGreen(p)*QuantumScale*GetPixelGreen(p)*QuantumScale*
2501 GetPixelGreen(p)*QuantumScale*GetPixelGreen(p);
2502 if ((double) GetPixelBlue(p) < channel_statistics[BlueChannel].minima)
2503 channel_statistics[BlueChannel].minima=(double) GetPixelBlue(p);
2504 if ((double) GetPixelBlue(p) > channel_statistics[BlueChannel].maxima)
2505 channel_statistics[BlueChannel].maxima=(double) GetPixelBlue(p);
2506 channel_statistics[BlueChannel].sum+=QuantumScale*GetPixelBlue(p);
2507 channel_statistics[BlueChannel].sum_squared+=QuantumScale*GetPixelBlue(p)*
2508 QuantumScale*GetPixelBlue(p);
2509 channel_statistics[BlueChannel].sum_cubed+=QuantumScale*GetPixelBlue(p)*
2510 QuantumScale*GetPixelBlue(p)*QuantumScale*GetPixelBlue(p);
2511 channel_statistics[BlueChannel].sum_fourth_power+=QuantumScale*
2512 GetPixelBlue(p)*QuantumScale*GetPixelBlue(p)*QuantumScale*
2513 GetPixelBlue(p)*QuantumScale*GetPixelBlue(p);
2514 histogram[ScaleQuantumToMap(GetPixelRed(p))].red++;
2515 histogram[ScaleQuantumToMap(GetPixelGreen(p))].green++;
2516 histogram[ScaleQuantumToMap(GetPixelBlue(p))].blue++;
2517 if (image->matte != MagickFalse)
2518 {
2519 if ((double) GetPixelAlpha(p) < channel_statistics[OpacityChannel].minima)
2520 channel_statistics[OpacityChannel].minima=(double) GetPixelAlpha(p);
2521 if ((double) GetPixelAlpha(p) > channel_statistics[OpacityChannel].maxima)
2522 channel_statistics[OpacityChannel].maxima=(double) GetPixelAlpha(p);
2523 channel_statistics[OpacityChannel].sum+=QuantumScale*GetPixelAlpha(p);
2524 channel_statistics[OpacityChannel].sum_squared+=QuantumScale*
2525 GetPixelAlpha(p)*QuantumScale*GetPixelAlpha(p);
2526 channel_statistics[OpacityChannel].sum_cubed+=QuantumScale*
2527 GetPixelAlpha(p)*QuantumScale*GetPixelAlpha(p)*QuantumScale*
2528 GetPixelAlpha(p);
2529 channel_statistics[OpacityChannel].sum_fourth_power+=QuantumScale*
2530 GetPixelAlpha(p)*QuantumScale*GetPixelAlpha(p)*QuantumScale*
2531 GetPixelAlpha(p)*QuantumScale*GetPixelAlpha(p);
2532 histogram[ScaleQuantumToMap(GetPixelAlpha(p))].opacity++;
2533 }
2534 if (image->colorspace == CMYKColorspace)
2535 {
2536 if ((double) GetPixelIndex(indexes+x) < channel_statistics[BlackChannel].minima)
2537 channel_statistics[BlackChannel].minima=(double)
2538 GetPixelIndex(indexes+x);
2539 if ((double) GetPixelIndex(indexes+x) > channel_statistics[BlackChannel].maxima)
2540 channel_statistics[BlackChannel].maxima=(double)
2541 GetPixelIndex(indexes+x);
2542 channel_statistics[BlackChannel].sum+=QuantumScale*
2543 GetPixelIndex(indexes+x);
2544 channel_statistics[BlackChannel].sum_squared+=QuantumScale*
2545 GetPixelIndex(indexes+x)*QuantumScale*GetPixelIndex(indexes+x);
2546 channel_statistics[BlackChannel].sum_cubed+=QuantumScale*
2547 GetPixelIndex(indexes+x)*QuantumScale*GetPixelIndex(indexes+x)*
2548 QuantumScale*GetPixelIndex(indexes+x);
2549 channel_statistics[BlackChannel].sum_fourth_power+=QuantumScale*
2550 GetPixelIndex(indexes+x)*QuantumScale*GetPixelIndex(indexes+x)*
2551 QuantumScale*GetPixelIndex(indexes+x)*QuantumScale*
2552 GetPixelIndex(indexes+x);
2553 histogram[ScaleQuantumToMap(GetPixelIndex(indexes+x))].index++;
2554 }
2555 x++;
2556 p++;
2557 }
2558 }
2559 for (i=0; i < (ssize_t) CompositeChannels; i++)
2560 {
2561 double
2562 area,
2563 mean,
2564 standard_deviation;
2565
2566 /*
2567 Normalize pixel statistics.
2568 */
2569 area=PerceptibleReciprocal((double) image->columns*image->rows);
2570 mean=channel_statistics[i].sum*area;
2571 channel_statistics[i].sum=mean;
2572 channel_statistics[i].sum_squared*=area;
2573 channel_statistics[i].sum_cubed*=area;
2574 channel_statistics[i].sum_fourth_power*=area;
2575 channel_statistics[i].mean=mean;
2576 channel_statistics[i].variance=channel_statistics[i].sum_squared;
2577 standard_deviation=sqrt(channel_statistics[i].variance-(mean*mean));
2578 area=PerceptibleReciprocal((double) image->columns*image->rows-1.0)*
2579 ((double) image->columns*image->rows);
2580 standard_deviation=sqrt(area*standard_deviation*standard_deviation);
2581 channel_statistics[i].standard_deviation=standard_deviation;
2582 }
2583 for (i=0; i < (ssize_t) (MaxMap+1U); i++)
2584 {
2585 if (histogram[i].red > 0.0)
2586 number_bins.red++;
2587 if (histogram[i].green > 0.0)
2588 number_bins.green++;
2589 if (histogram[i].blue > 0.0)
2590 number_bins.blue++;
2591 if ((image->matte != MagickFalse) && (histogram[i].opacity > 0.0))
2592 number_bins.opacity++;
2593 if ((image->colorspace == CMYKColorspace) && (histogram[i].index > 0.0))
2594 number_bins.index++;
2595 }
2596 area=PerceptibleReciprocal((double) image->columns*image->rows);
2597 for (i=0; i < (ssize_t) (MaxMap+1U); i++)
2598 {
2599 /*
2600 Compute pixel entropy.
2601 */
2602 histogram[i].red*=area;
2603 channel_statistics[RedChannel].entropy+=-histogram[i].red*
2604 MagickLog10(histogram[i].red)*
2605 PerceptibleReciprocal(MagickLog10((double) number_bins.red));
2606 histogram[i].green*=area;
2607 channel_statistics[GreenChannel].entropy+=-histogram[i].green*
2608 MagickLog10(histogram[i].green)*
2609 PerceptibleReciprocal(MagickLog10((double) number_bins.green));
2610 histogram[i].blue*=area;
2611 channel_statistics[BlueChannel].entropy+=-histogram[i].blue*
2612 MagickLog10(histogram[i].blue)*
2613 PerceptibleReciprocal(MagickLog10((double) number_bins.blue));
2614 if (image->matte != MagickFalse)
2615 {
2616 histogram[i].opacity*=area;
2617 channel_statistics[OpacityChannel].entropy+=-histogram[i].opacity*
2618 MagickLog10(histogram[i].opacity)*
2619 PerceptibleReciprocal(MagickLog10((double) number_bins.opacity));
2620 }
2621 if (image->colorspace == CMYKColorspace)
2622 {
2623 histogram[i].index*=area;
2624 channel_statistics[IndexChannel].entropy+=-histogram[i].index*
2625 MagickLog10(histogram[i].index)*
2626 PerceptibleReciprocal(MagickLog10((double) number_bins.index));
2627 }
2628 }
2629 /*
2630 Compute overall statistics.
2631 */
2632 for (i=0; i < (ssize_t) CompositeChannels; i++)
2633 {
2634 channel_statistics[CompositeChannels].depth=(size_t) EvaluateMax((double)
2635 channel_statistics[CompositeChannels].depth,(double)
2636 channel_statistics[i].depth);
2637 channel_statistics[CompositeChannels].minima=MagickMin(
2638 channel_statistics[CompositeChannels].minima,
2639 channel_statistics[i].minima);
2640 channel_statistics[CompositeChannels].maxima=EvaluateMax(
2641 channel_statistics[CompositeChannels].maxima,
2642 channel_statistics[i].maxima);
2643 channel_statistics[CompositeChannels].sum+=channel_statistics[i].sum;
2644 channel_statistics[CompositeChannels].sum_squared+=
2645 channel_statistics[i].sum_squared;
2646 channel_statistics[CompositeChannels].sum_cubed+=
2647 channel_statistics[i].sum_cubed;
2648 channel_statistics[CompositeChannels].sum_fourth_power+=
2649 channel_statistics[i].sum_fourth_power;
2650 channel_statistics[CompositeChannels].mean+=channel_statistics[i].mean;
2651 channel_statistics[CompositeChannels].variance+=
2652 channel_statistics[i].variance-channel_statistics[i].mean*
2653 channel_statistics[i].mean;
2654 standard_deviation=sqrt(channel_statistics[i].variance-
2655 (channel_statistics[i].mean*channel_statistics[i].mean));
2656 area=PerceptibleReciprocal((double) image->columns*image->rows-1.0)*
2657 ((double) image->columns*image->rows);
2658 standard_deviation=sqrt(area*standard_deviation*standard_deviation);
2659 channel_statistics[CompositeChannels].standard_deviation=standard_deviation;
2660 channel_statistics[CompositeChannels].entropy+=
2661 channel_statistics[i].entropy;
2662 }
2663 channels=3;
2664 if (image->matte != MagickFalse)
2665 channels++;
2666 if (image->colorspace == CMYKColorspace)
2667 channels++;
2668 channel_statistics[CompositeChannels].sum/=channels;
2669 channel_statistics[CompositeChannels].sum_squared/=channels;
2670 channel_statistics[CompositeChannels].sum_cubed/=channels;
2671 channel_statistics[CompositeChannels].sum_fourth_power/=channels;
2672 channel_statistics[CompositeChannels].mean/=channels;
2673 channel_statistics[CompositeChannels].kurtosis/=channels;
2674 channel_statistics[CompositeChannels].skewness/=channels;
2675 channel_statistics[CompositeChannels].entropy/=channels;
2676 i=CompositeChannels;
2677 area=PerceptibleReciprocal((double) channels*image->columns*image->rows);
2678 channel_statistics[i].variance=channel_statistics[i].sum_squared;
2679 channel_statistics[i].mean=channel_statistics[i].sum;
2680 standard_deviation=sqrt(channel_statistics[i].variance-
2681 (channel_statistics[i].mean*channel_statistics[i].mean));
2682 standard_deviation=sqrt(PerceptibleReciprocal((double) channels*
2683 image->columns*image->rows-1.0)*channels*image->columns*image->rows*
2684 standard_deviation*standard_deviation);
2685 channel_statistics[i].standard_deviation=standard_deviation;
2686 for (i=0; i <= (ssize_t) CompositeChannels; i++)
2687 {
2688 /*
2689 Compute kurtosis & skewness statistics.
2690 */
2691 standard_deviation=PerceptibleReciprocal(
2692 channel_statistics[i].standard_deviation);
2693 channel_statistics[i].skewness=(channel_statistics[i].sum_cubed-3.0*
2694 channel_statistics[i].mean*channel_statistics[i].sum_squared+2.0*
2695 channel_statistics[i].mean*channel_statistics[i].mean*
2696 channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2697 standard_deviation);
2698 channel_statistics[i].kurtosis=(channel_statistics[i].sum_fourth_power-4.0*
2699 channel_statistics[i].mean*channel_statistics[i].sum_cubed+6.0*
2700 channel_statistics[i].mean*channel_statistics[i].mean*
2701 channel_statistics[i].sum_squared-3.0*channel_statistics[i].mean*
2702 channel_statistics[i].mean*1.0*channel_statistics[i].mean*
2703 channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2704 standard_deviation*standard_deviation)-3.0;
2705 }
2706 for (i=0; i <= (ssize_t) CompositeChannels; i++)
2707 {
2708 channel_statistics[i].mean*=QuantumRange;
2709 channel_statistics[i].variance*=QuantumRange;
2710 channel_statistics[i].standard_deviation*=QuantumRange;
2711 channel_statistics[i].sum*=QuantumRange;
2712 channel_statistics[i].sum_squared*=QuantumRange;
2713 channel_statistics[i].sum_cubed*=QuantumRange;
2714 channel_statistics[i].sum_fourth_power*=QuantumRange;
2715 }
2716 channel_statistics[CompositeChannels].mean=0.0;
2717 channel_statistics[CompositeChannels].standard_deviation=0.0;
2718 for (i=0; i < (ssize_t) CompositeChannels; i++)
2719 {
2720 channel_statistics[CompositeChannels].mean+=
2721 channel_statistics[i].mean;
2722 channel_statistics[CompositeChannels].standard_deviation+=
2723 channel_statistics[i].standard_deviation;
2724 }
2725 channel_statistics[CompositeChannels].mean/=(double) channels;
2726 channel_statistics[CompositeChannels].standard_deviation/=(double) channels;
2727 histogram=(MagickPixelPacket *) RelinquishMagickMemory(histogram);
2728 if (y < (ssize_t) image->rows)
2729 channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
2730 channel_statistics);
2731 return(channel_statistics);
2732}
2733
2734/*
2735%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2736% %
2737% %
2738% %
2739% P o l y n o m i a l I m a g e %
2740% %
2741% %
2742% %
2743%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2744%
2745% PolynomialImage() returns a new image where each pixel is the sum of the
2746% pixels in the image sequence after applying its corresponding terms
2747% (coefficient and degree pairs).
2748%
2749% The format of the PolynomialImage method is:
2750%
2751% Image *PolynomialImage(const Image *images,const size_t number_terms,
2752% const double *terms,ExceptionInfo *exception)
2753% Image *PolynomialImageChannel(const Image *images,
2754% const size_t number_terms,const ChannelType channel,
2755% const double *terms,ExceptionInfo *exception)
2756%
2757% A description of each parameter follows:
2758%
2759% o images: the image sequence.
2760%
2761% o channel: the channel.
2762%
2763% o number_terms: the number of terms in the list. The actual list length
2764% is 2 x number_terms + 1 (the constant).
2765%
2766% o terms: the list of polynomial coefficients and degree pairs and a
2767% constant.
2768%
2769% o exception: return any errors or warnings in this structure.
2770%
2771*/
2772MagickExport Image *PolynomialImage(const Image *images,
2773 const size_t number_terms,const double *terms,ExceptionInfo *exception)
2774{
2775 Image
2776 *polynomial_image;
2777
2778 polynomial_image=PolynomialImageChannel(images,DefaultChannels,number_terms,
2779 terms,exception);
2780 return(polynomial_image);
2781}
2782
2783MagickExport Image *PolynomialImageChannel(const Image *images,
2784 const ChannelType channel,const size_t number_terms,const double *terms,
2785 ExceptionInfo *exception)
2786{
2787#define PolynomialImageTag "Polynomial/Image"
2788
2789 CacheView
2790 *polynomial_view;
2791
2792 Image
2793 *image;
2794
2795 MagickBooleanType
2796 status;
2797
2798 MagickOffsetType
2799 progress;
2800
2802 **magick_restrict polynomial_pixels,
2803 zero;
2804
2805 ssize_t
2806 y;
2807
2808 assert(images != (Image *) NULL);
2809 assert(images->signature == MagickCoreSignature);
2810 if (IsEventLogging() != MagickFalse)
2811 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
2812 assert(exception != (ExceptionInfo *) NULL);
2813 assert(exception->signature == MagickCoreSignature);
2814 image=AcquireImageCanvas(images,exception);
2815 if (image == (Image *) NULL)
2816 return((Image *) NULL);
2817 if (SetImageStorageClass(image,DirectClass) == MagickFalse)
2818 {
2819 InheritException(exception,&image->exception);
2820 image=DestroyImage(image);
2821 return((Image *) NULL);
2822 }
2823 polynomial_pixels=AcquirePixelTLS(images);
2824 if (polynomial_pixels == (MagickPixelPacket **) NULL)
2825 {
2826 image=DestroyImage(image);
2827 (void) ThrowMagickException(exception,GetMagickModule(),
2828 ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
2829 return((Image *) NULL);
2830 }
2831 /*
2832 Polynomial image pixels.
2833 */
2834 status=MagickTrue;
2835 progress=0;
2836 GetMagickPixelPacket(images,&zero);
2837 polynomial_view=AcquireAuthenticCacheView(image,exception);
2838#if defined(MAGICKCORE_OPENMP_SUPPORT)
2839 #pragma omp parallel for schedule(static) shared(progress,status) \
2840 magick_number_threads(image,image,image->rows,1)
2841#endif
2842 for (y=0; y < (ssize_t) image->rows; y++)
2843 {
2844 CacheView
2845 *image_view;
2846
2847 const Image
2848 *next;
2849
2850 const int
2851 id = GetOpenMPThreadId();
2852
2853 IndexPacket
2854 *magick_restrict polynomial_indexes;
2855
2857 *polynomial_pixel;
2858
2860 *magick_restrict q;
2861
2862 ssize_t
2863 i,
2864 x;
2865
2866 size_t
2867 number_images;
2868
2869 if (status == MagickFalse)
2870 continue;
2871 q=QueueCacheViewAuthenticPixels(polynomial_view,0,y,image->columns,1,
2872 exception);
2873 if (q == (PixelPacket *) NULL)
2874 {
2875 status=MagickFalse;
2876 continue;
2877 }
2878 polynomial_indexes=GetCacheViewAuthenticIndexQueue(polynomial_view);
2879 polynomial_pixel=polynomial_pixels[id];
2880 for (x=0; x < (ssize_t) image->columns; x++)
2881 polynomial_pixel[x]=zero;
2882 next=images;
2883 number_images=GetImageListLength(images);
2884 for (i=0; i < (ssize_t) number_images; i++)
2885 {
2886 const IndexPacket
2887 *indexes;
2888
2889 const PixelPacket
2890 *p;
2891
2892 if (i >= (ssize_t) number_terms)
2893 break;
2894 image_view=AcquireVirtualCacheView(next,exception);
2895 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2896 if (p == (const PixelPacket *) NULL)
2897 {
2898 image_view=DestroyCacheView(image_view);
2899 break;
2900 }
2901 indexes=GetCacheViewVirtualIndexQueue(image_view);
2902 for (x=0; x < (ssize_t) image->columns; x++)
2903 {
2904 double
2905 coefficient,
2906 degree;
2907
2908 coefficient=terms[i << 1];
2909 degree=terms[(i << 1)+1];
2910 if ((channel & RedChannel) != 0)
2911 polynomial_pixel[x].red+=coefficient*pow(QuantumScale*(double)
2912 p->red,degree);
2913 if ((channel & GreenChannel) != 0)
2914 polynomial_pixel[x].green+=coefficient*pow(QuantumScale*(double)
2915 p->green,
2916 degree);
2917 if ((channel & BlueChannel) != 0)
2918 polynomial_pixel[x].blue+=coefficient*pow(QuantumScale*(double)
2919 p->blue,degree);
2920 if ((channel & OpacityChannel) != 0)
2921 polynomial_pixel[x].opacity+=coefficient*pow(QuantumScale*
2922 ((double) QuantumRange-(double) p->opacity),degree);
2923 if (((channel & IndexChannel) != 0) &&
2924 (image->colorspace == CMYKColorspace))
2925 polynomial_pixel[x].index+=coefficient*pow(QuantumScale*(double)
2926 indexes[x],degree);
2927 p++;
2928 }
2929 image_view=DestroyCacheView(image_view);
2930 next=GetNextImageInList(next);
2931 }
2932 for (x=0; x < (ssize_t) image->columns; x++)
2933 {
2934 SetPixelRed(q,ClampToQuantum((MagickRealType) QuantumRange*
2935 polynomial_pixel[x].red));
2936 SetPixelGreen(q,ClampToQuantum((MagickRealType) QuantumRange*
2937 polynomial_pixel[x].green));
2938 SetPixelBlue(q,ClampToQuantum((MagickRealType) QuantumRange*
2939 polynomial_pixel[x].blue));
2940 if (image->matte == MagickFalse)
2941 SetPixelOpacity(q,ClampToQuantum((MagickRealType) QuantumRange-
2942 (MagickRealType) QuantumRange*polynomial_pixel[x].opacity));
2943 else
2944 SetPixelAlpha(q,ClampToQuantum((MagickRealType) QuantumRange-
2945 (MagickRealType) QuantumRange*polynomial_pixel[x].opacity));
2946 if (image->colorspace == CMYKColorspace)
2947 SetPixelIndex(polynomial_indexes+x,ClampToQuantum((MagickRealType)
2948 QuantumRange*polynomial_pixel[x].index));
2949 q++;
2950 }
2951 if (SyncCacheViewAuthenticPixels(polynomial_view,exception) == MagickFalse)
2952 status=MagickFalse;
2953 if (images->progress_monitor != (MagickProgressMonitor) NULL)
2954 {
2955 MagickBooleanType
2956 proceed;
2957
2958 proceed=SetImageProgress(images,PolynomialImageTag,progress++,
2959 image->rows);
2960 if (proceed == MagickFalse)
2961 status=MagickFalse;
2962 }
2963 }
2964 polynomial_view=DestroyCacheView(polynomial_view);
2965 polynomial_pixels=DestroyPixelTLS(images,polynomial_pixels);
2966 if (status == MagickFalse)
2967 image=DestroyImage(image);
2968 return(image);
2969}
2970
2971/*
2972%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2973% %
2974% %
2975% %
2976% S t a t i s t i c I m a g e %
2977% %
2978% %
2979% %
2980%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2981%
2982% StatisticImage() makes each pixel the min / max / median / mode / etc. of
2983% the neighborhood of the specified width and height.
2984%
2985% The format of the StatisticImage method is:
2986%
2987% Image *StatisticImage(const Image *image,const StatisticType type,
2988% const size_t width,const size_t height,ExceptionInfo *exception)
2989% Image *StatisticImageChannel(const Image *image,
2990% const ChannelType channel,const StatisticType type,
2991% const size_t width,const size_t height,ExceptionInfo *exception)
2992%
2993% A description of each parameter follows:
2994%
2995% o image: the image.
2996%
2997% o channel: the image channel.
2998%
2999% o type: the statistic type (median, mode, etc.).
3000%
3001% o width: the width of the pixel neighborhood.
3002%
3003% o height: the height of the pixel neighborhood.
3004%
3005% o exception: return any errors or warnings in this structure.
3006%
3007*/
3008
3009#define ListChannels 5
3010
3011typedef struct _ListNode
3012{
3013 size_t
3014 next[9],
3015 count,
3016 signature;
3017} ListNode;
3018
3019typedef struct _SkipList
3020{
3021 ssize_t
3022 level;
3023
3024 ListNode
3025 *nodes;
3026} SkipList;
3027
3028typedef struct _PixelList
3029{
3030 size_t
3031 length,
3032 seed,
3033 signature;
3034
3035 SkipList
3036 lists[ListChannels];
3037} PixelList;
3038
3039static PixelList *DestroyPixelList(PixelList *pixel_list)
3040{
3041 ssize_t
3042 i;
3043
3044 if (pixel_list == (PixelList *) NULL)
3045 return((PixelList *) NULL);
3046 for (i=0; i < ListChannels; i++)
3047 if (pixel_list->lists[i].nodes != (ListNode *) NULL)
3048 pixel_list->lists[i].nodes=(ListNode *) RelinquishAlignedMemory(
3049 pixel_list->lists[i].nodes);
3050 pixel_list=(PixelList *) RelinquishMagickMemory(pixel_list);
3051 return(pixel_list);
3052}
3053
3054static PixelList **DestroyPixelListTLS(PixelList **pixel_list)
3055{
3056 ssize_t
3057 i;
3058
3059 assert(pixel_list != (PixelList **) NULL);
3060 for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
3061 if (pixel_list[i] != (PixelList *) NULL)
3062 pixel_list[i]=DestroyPixelList(pixel_list[i]);
3063 pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list);
3064 return(pixel_list);
3065}
3066
3067static PixelList *AcquirePixelList(const size_t width,const size_t height)
3068{
3069 PixelList
3070 *pixel_list;
3071
3072 ssize_t
3073 i;
3074
3075 pixel_list=(PixelList *) AcquireMagickMemory(sizeof(*pixel_list));
3076 if (pixel_list == (PixelList *) NULL)
3077 return(pixel_list);
3078 (void) memset((void *) pixel_list,0,sizeof(*pixel_list));
3079 pixel_list->length=width*height;
3080 for (i=0; i < ListChannels; i++)
3081 {
3082 pixel_list->lists[i].nodes=(ListNode *) AcquireAlignedMemory(65537UL,
3083 sizeof(*pixel_list->lists[i].nodes));
3084 if (pixel_list->lists[i].nodes == (ListNode *) NULL)
3085 return(DestroyPixelList(pixel_list));
3086 (void) memset(pixel_list->lists[i].nodes,0,65537UL*
3087 sizeof(*pixel_list->lists[i].nodes));
3088 }
3089 pixel_list->signature=MagickCoreSignature;
3090 return(pixel_list);
3091}
3092
3093static PixelList **AcquirePixelListTLS(const size_t width,
3094 const size_t height)
3095{
3096 PixelList
3097 **pixel_list;
3098
3099 ssize_t
3100 i;
3101
3102 size_t
3103 number_threads;
3104
3105 number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
3106 pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
3107 sizeof(*pixel_list));
3108 if (pixel_list == (PixelList **) NULL)
3109 return((PixelList **) NULL);
3110 (void) memset(pixel_list,0,number_threads*sizeof(*pixel_list));
3111 for (i=0; i < (ssize_t) number_threads; i++)
3112 {
3113 pixel_list[i]=AcquirePixelList(width,height);
3114 if (pixel_list[i] == (PixelList *) NULL)
3115 return(DestroyPixelListTLS(pixel_list));
3116 }
3117 return(pixel_list);
3118}
3119
3120static void AddNodePixelList(PixelList *pixel_list,const ssize_t channel,
3121 const size_t color)
3122{
3123 SkipList
3124 *list;
3125
3126 ssize_t
3127 level;
3128
3129 size_t
3130 search,
3131 update[9];
3132
3133 /*
3134 Initialize the node.
3135 */
3136 list=pixel_list->lists+channel;
3137 list->nodes[color].signature=pixel_list->signature;
3138 list->nodes[color].count=1;
3139 /*
3140 Determine where it belongs in the list.
3141 */
3142 search=65536UL;
3143 (void) memset(update,0,sizeof(update));
3144 for (level=list->level; level >= 0; level--)
3145 {
3146 while (list->nodes[search].next[level] < color)
3147 search=list->nodes[search].next[level];
3148 update[level]=search;
3149 }
3150 /*
3151 Generate a pseudo-random level for this node.
3152 */
3153 for (level=0; ; level++)
3154 {
3155 pixel_list->seed=(pixel_list->seed*42893621L)+1L;
3156 if ((pixel_list->seed & 0x300) != 0x300)
3157 break;
3158 }
3159 if (level > 8)
3160 level=8;
3161 if (level > (list->level+2))
3162 level=list->level+2;
3163 /*
3164 If we're raising the list's level, link back to the root node.
3165 */
3166 while (level > list->level)
3167 {
3168 list->level++;
3169 update[list->level]=65536UL;
3170 }
3171 /*
3172 Link the node into the skip-list.
3173 */
3174 do
3175 {
3176 list->nodes[color].next[level]=list->nodes[update[level]].next[level];
3177 list->nodes[update[level]].next[level]=color;
3178 } while (level-- > 0);
3179}
3180
3181static void GetMaximumPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3182{
3183 SkipList
3184 *list;
3185
3186 ssize_t
3187 channel;
3188
3189 size_t
3190 color,
3191 maximum;
3192
3193 ssize_t
3194 count;
3195
3196 unsigned short
3197 channels[ListChannels];
3198
3199 /*
3200 Find the maximum value for each of the color.
3201 */
3202 for (channel=0; channel < 5; channel++)
3203 {
3204 list=pixel_list->lists+channel;
3205 color=65536L;
3206 count=0;
3207 maximum=list->nodes[color].next[0];
3208 do
3209 {
3210 color=list->nodes[color].next[0];
3211 if (color > maximum)
3212 maximum=color;
3213 count+=list->nodes[color].count;
3214 } while (count < (ssize_t) pixel_list->length);
3215 channels[channel]=(unsigned short) maximum;
3216 }
3217 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3218 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3219 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3220 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3221 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3222}
3223
3224static void GetMeanPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3225{
3226 MagickRealType
3227 sum;
3228
3229 SkipList
3230 *list;
3231
3232 ssize_t
3233 channel;
3234
3235 size_t
3236 color;
3237
3238 ssize_t
3239 count;
3240
3241 unsigned short
3242 channels[ListChannels];
3243
3244 /*
3245 Find the mean value for each of the color.
3246 */
3247 for (channel=0; channel < 5; channel++)
3248 {
3249 list=pixel_list->lists+channel;
3250 color=65536L;
3251 count=0;
3252 sum=0.0;
3253 do
3254 {
3255 color=list->nodes[color].next[0];
3256 sum+=(MagickRealType) list->nodes[color].count*color;
3257 count+=list->nodes[color].count;
3258 } while (count < (ssize_t) pixel_list->length);
3259 sum/=pixel_list->length;
3260 channels[channel]=(unsigned short) sum;
3261 }
3262 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3263 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3264 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3265 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3266 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3267}
3268
3269static void GetMedianPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3270{
3271 SkipList
3272 *list;
3273
3274 ssize_t
3275 channel;
3276
3277 size_t
3278 color;
3279
3280 ssize_t
3281 count;
3282
3283 unsigned short
3284 channels[ListChannels];
3285
3286 /*
3287 Find the median value for each of the color.
3288 */
3289 for (channel=0; channel < 5; channel++)
3290 {
3291 list=pixel_list->lists+channel;
3292 color=65536L;
3293 count=0;
3294 do
3295 {
3296 color=list->nodes[color].next[0];
3297 count+=list->nodes[color].count;
3298 } while (count <= (ssize_t) (pixel_list->length >> 1));
3299 channels[channel]=(unsigned short) color;
3300 }
3301 GetMagickPixelPacket((const Image *) NULL,pixel);
3302 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3303 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3304 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3305 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3306 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3307}
3308
3309static void GetMinimumPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3310{
3311 SkipList
3312 *list;
3313
3314 ssize_t
3315 channel;
3316
3317 size_t
3318 color,
3319 minimum;
3320
3321 ssize_t
3322 count;
3323
3324 unsigned short
3325 channels[ListChannels];
3326
3327 /*
3328 Find the minimum value for each of the color.
3329 */
3330 for (channel=0; channel < 5; channel++)
3331 {
3332 list=pixel_list->lists+channel;
3333 count=0;
3334 color=65536UL;
3335 minimum=list->nodes[color].next[0];
3336 do
3337 {
3338 color=list->nodes[color].next[0];
3339 if (color < minimum)
3340 minimum=color;
3341 count+=list->nodes[color].count;
3342 } while (count < (ssize_t) pixel_list->length);
3343 channels[channel]=(unsigned short) minimum;
3344 }
3345 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3346 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3347 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3348 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3349 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3350}
3351
3352static void GetModePixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3353{
3354 SkipList
3355 *list;
3356
3357 ssize_t
3358 channel;
3359
3360 size_t
3361 color,
3362 max_count,
3363 mode;
3364
3365 ssize_t
3366 count;
3367
3368 unsigned short
3369 channels[5];
3370
3371 /*
3372 Make each pixel the 'predominant color' of the specified neighborhood.
3373 */
3374 for (channel=0; channel < 5; channel++)
3375 {
3376 list=pixel_list->lists+channel;
3377 color=65536L;
3378 mode=color;
3379 max_count=list->nodes[mode].count;
3380 count=0;
3381 do
3382 {
3383 color=list->nodes[color].next[0];
3384 if (list->nodes[color].count > max_count)
3385 {
3386 mode=color;
3387 max_count=list->nodes[mode].count;
3388 }
3389 count+=list->nodes[color].count;
3390 } while (count < (ssize_t) pixel_list->length);
3391 channels[channel]=(unsigned short) mode;
3392 }
3393 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3394 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3395 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3396 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3397 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3398}
3399
3400static void GetNonpeakPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3401{
3402 SkipList
3403 *list;
3404
3405 ssize_t
3406 channel;
3407
3408 size_t
3409 color,
3410 next,
3411 previous;
3412
3413 ssize_t
3414 count;
3415
3416 unsigned short
3417 channels[5];
3418
3419 /*
3420 Finds the non peak value for each of the colors.
3421 */
3422 for (channel=0; channel < 5; channel++)
3423 {
3424 list=pixel_list->lists+channel;
3425 color=65536L;
3426 next=list->nodes[color].next[0];
3427 count=0;
3428 do
3429 {
3430 previous=color;
3431 color=next;
3432 next=list->nodes[color].next[0];
3433 count+=list->nodes[color].count;
3434 } while (count <= (ssize_t) (pixel_list->length >> 1));
3435 if ((previous == 65536UL) && (next != 65536UL))
3436 color=next;
3437 else
3438 if ((previous != 65536UL) && (next == 65536UL))
3439 color=previous;
3440 channels[channel]=(unsigned short) color;
3441 }
3442 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3443 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3444 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3445 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3446 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3447}
3448
3449static void GetRootMeanSquarePixelList(PixelList *pixel_list,
3450 MagickPixelPacket *pixel)
3451{
3452 MagickRealType
3453 sum;
3454
3455 SkipList
3456 *list;
3457
3458 ssize_t
3459 channel;
3460
3461 size_t
3462 color;
3463
3464 ssize_t
3465 count;
3466
3467 unsigned short
3468 channels[ListChannels];
3469
3470 /*
3471 Find the root mean square value for each of the color.
3472 */
3473 for (channel=0; channel < 5; channel++)
3474 {
3475 list=pixel_list->lists+channel;
3476 color=65536L;
3477 count=0;
3478 sum=0.0;
3479 do
3480 {
3481 color=list->nodes[color].next[0];
3482 sum+=(MagickRealType) (list->nodes[color].count*color*color);
3483 count+=list->nodes[color].count;
3484 } while (count < (ssize_t) pixel_list->length);
3485 sum/=pixel_list->length;
3486 channels[channel]=(unsigned short) sqrt(sum);
3487 }
3488 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3489 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3490 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3491 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3492 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3493}
3494
3495static void GetStandardDeviationPixelList(PixelList *pixel_list,
3496 MagickPixelPacket *pixel)
3497{
3498 MagickRealType
3499 sum,
3500 sum_squared;
3501
3502 SkipList
3503 *list;
3504
3505 size_t
3506 color;
3507
3508 ssize_t
3509 channel,
3510 count;
3511
3512 unsigned short
3513 channels[ListChannels];
3514
3515 /*
3516 Find the standard-deviation value for each of the color.
3517 */
3518 for (channel=0; channel < 5; channel++)
3519 {
3520 list=pixel_list->lists+channel;
3521 color=65536L;
3522 count=0;
3523 sum=0.0;
3524 sum_squared=0.0;
3525 do
3526 {
3527 ssize_t
3528 i;
3529
3530 color=list->nodes[color].next[0];
3531 sum+=(MagickRealType) list->nodes[color].count*color;
3532 for (i=0; i < (ssize_t) list->nodes[color].count; i++)
3533 sum_squared+=((MagickRealType) color)*((MagickRealType) color);
3534 count+=list->nodes[color].count;
3535 } while (count < (ssize_t) pixel_list->length);
3536 sum/=pixel_list->length;
3537 sum_squared/=pixel_list->length;
3538 channels[channel]=(unsigned short) sqrt(sum_squared-(sum*sum));
3539 }
3540 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3541 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3542 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3543 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3544 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3545}
3546
3547static inline void InsertPixelList(const Image *image,const PixelPacket *pixel,
3548 const IndexPacket *indexes,PixelList *pixel_list)
3549{
3550 size_t
3551 signature;
3552
3553 unsigned short
3554 index;
3555
3556 index=ScaleQuantumToShort(GetPixelRed(pixel));
3557 signature=pixel_list->lists[0].nodes[index].signature;
3558 if (signature == pixel_list->signature)
3559 pixel_list->lists[0].nodes[index].count++;
3560 else
3561 AddNodePixelList(pixel_list,0,index);
3562 index=ScaleQuantumToShort(GetPixelGreen(pixel));
3563 signature=pixel_list->lists[1].nodes[index].signature;
3564 if (signature == pixel_list->signature)
3565 pixel_list->lists[1].nodes[index].count++;
3566 else
3567 AddNodePixelList(pixel_list,1,index);
3568 index=ScaleQuantumToShort(GetPixelBlue(pixel));
3569 signature=pixel_list->lists[2].nodes[index].signature;
3570 if (signature == pixel_list->signature)
3571 pixel_list->lists[2].nodes[index].count++;
3572 else
3573 AddNodePixelList(pixel_list,2,index);
3574 index=ScaleQuantumToShort(GetPixelOpacity(pixel));
3575 signature=pixel_list->lists[3].nodes[index].signature;
3576 if (signature == pixel_list->signature)
3577 pixel_list->lists[3].nodes[index].count++;
3578 else
3579 AddNodePixelList(pixel_list,3,index);
3580 if (image->colorspace == CMYKColorspace)
3581 index=ScaleQuantumToShort(GetPixelIndex(indexes));
3582 signature=pixel_list->lists[4].nodes[index].signature;
3583 if (signature == pixel_list->signature)
3584 pixel_list->lists[4].nodes[index].count++;
3585 else
3586 AddNodePixelList(pixel_list,4,index);
3587}
3588
3589static void ResetPixelList(PixelList *pixel_list)
3590{
3591 int
3592 level;
3593
3594 ListNode
3595 *root;
3596
3597 SkipList
3598 *list;
3599
3600 ssize_t
3601 channel;
3602
3603 /*
3604 Reset the skip-list.
3605 */
3606 for (channel=0; channel < 5; channel++)
3607 {
3608 list=pixel_list->lists+channel;
3609 root=list->nodes+65536UL;
3610 list->level=0;
3611 for (level=0; level < 9; level++)
3612 root->next[level]=65536UL;
3613 }
3614 pixel_list->seed=pixel_list->signature++;
3615}
3616
3617MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
3618 const size_t width,const size_t height,ExceptionInfo *exception)
3619{
3620 Image
3621 *statistic_image;
3622
3623 statistic_image=StatisticImageChannel(image,DefaultChannels,type,width,
3624 height,exception);
3625 return(statistic_image);
3626}
3627
3628MagickExport Image *StatisticImageChannel(const Image *image,
3629 const ChannelType channel,const StatisticType type,const size_t width,
3630 const size_t height,ExceptionInfo *exception)
3631{
3632#define StatisticImageTag "Statistic/Image"
3633
3634 CacheView
3635 *image_view,
3636 *statistic_view;
3637
3638 Image
3639 *statistic_image;
3640
3641 MagickBooleanType
3642 status;
3643
3644 MagickOffsetType
3645 progress;
3646
3647 PixelList
3648 **magick_restrict pixel_list;
3649
3650 size_t
3651 neighbor_height,
3652 neighbor_width;
3653
3654 ssize_t
3655 y;
3656
3657 /*
3658 Initialize statistics image attributes.
3659 */
3660 assert(image != (Image *) NULL);
3661 assert(image->signature == MagickCoreSignature);
3662 if (IsEventLogging() != MagickFalse)
3663 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3664 assert(exception != (ExceptionInfo *) NULL);
3665 assert(exception->signature == MagickCoreSignature);
3666 statistic_image=CloneImage(image,0,0,MagickTrue,exception);
3667 if (statistic_image == (Image *) NULL)
3668 return((Image *) NULL);
3669 if (SetImageStorageClass(statistic_image,DirectClass) == MagickFalse)
3670 {
3671 InheritException(exception,&statistic_image->exception);
3672 statistic_image=DestroyImage(statistic_image);
3673 return((Image *) NULL);
3674 }
3675 neighbor_width=width == 0 ? GetOptimalKernelWidth2D((double) width,0.5) :
3676 width;
3677 neighbor_height=height == 0 ? GetOptimalKernelWidth2D((double) height,0.5) :
3678 height;
3679 pixel_list=AcquirePixelListTLS(neighbor_width,neighbor_height);
3680 if (pixel_list == (PixelList **) NULL)
3681 {
3682 statistic_image=DestroyImage(statistic_image);
3683 ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3684 }
3685 /*
3686 Make each pixel the min / max / median / mode / etc. of the neighborhood.
3687 */
3688 status=MagickTrue;
3689 progress=0;
3690 image_view=AcquireVirtualCacheView(image,exception);
3691 statistic_view=AcquireAuthenticCacheView(statistic_image,exception);
3692#if defined(MAGICKCORE_OPENMP_SUPPORT)
3693 #pragma omp parallel for schedule(static) shared(progress,status) \
3694 magick_number_threads(image,statistic_image,statistic_image->rows,1)
3695#endif
3696 for (y=0; y < (ssize_t) statistic_image->rows; y++)
3697 {
3698 const int
3699 id = GetOpenMPThreadId();
3700
3701 const IndexPacket
3702 *magick_restrict indexes;
3703
3704 const PixelPacket
3705 *magick_restrict p;
3706
3707 IndexPacket
3708 *magick_restrict statistic_indexes;
3709
3711 *magick_restrict q;
3712
3713 ssize_t
3714 x;
3715
3716 if (status == MagickFalse)
3717 continue;
3718 p=GetCacheViewVirtualPixels(image_view,-((ssize_t) neighbor_width/2L),y-
3719 (ssize_t) (neighbor_height/2L),image->columns+neighbor_width,
3720 neighbor_height,exception);
3721 q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns, 1,exception);
3722 if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
3723 {
3724 status=MagickFalse;
3725 continue;
3726 }
3727 indexes=GetCacheViewVirtualIndexQueue(image_view);
3728 statistic_indexes=GetCacheViewAuthenticIndexQueue(statistic_view);
3729 for (x=0; x < (ssize_t) statistic_image->columns; x++)
3730 {
3732 pixel;
3733
3734 const IndexPacket
3735 *magick_restrict s;
3736
3737 const PixelPacket
3738 *magick_restrict r;
3739
3740 ssize_t
3741 u,
3742 v;
3743
3744 r=p;
3745 s=indexes+x;
3746 ResetPixelList(pixel_list[id]);
3747 for (v=0; v < (ssize_t) neighbor_height; v++)
3748 {
3749 for (u=0; u < (ssize_t) neighbor_width; u++)
3750 InsertPixelList(image,r+u,s+u,pixel_list[id]);
3751 r+=image->columns+neighbor_width;
3752 s+=image->columns+neighbor_width;
3753 }
3754 GetMagickPixelPacket(image,&pixel);
3755 SetMagickPixelPacket(image,p+neighbor_width*neighbor_height/2,indexes+x+
3756 neighbor_width*neighbor_height/2,&pixel);
3757 switch (type)
3758 {
3759 case GradientStatistic:
3760 {
3762 maximum,
3763 minimum;
3764
3765 GetMinimumPixelList(pixel_list[id],&pixel);
3766 minimum=pixel;
3767 GetMaximumPixelList(pixel_list[id],&pixel);
3768 maximum=pixel;
3769 pixel.red=MagickAbsoluteValue(maximum.red-minimum.red);
3770 pixel.green=MagickAbsoluteValue(maximum.green-minimum.green);
3771 pixel.blue=MagickAbsoluteValue(maximum.blue-minimum.blue);
3772 pixel.opacity=MagickAbsoluteValue(maximum.opacity-minimum.opacity);
3773 if (image->colorspace == CMYKColorspace)
3774 pixel.index=MagickAbsoluteValue(maximum.index-minimum.index);
3775 break;
3776 }
3777 case MaximumStatistic:
3778 {
3779 GetMaximumPixelList(pixel_list[id],&pixel);
3780 break;
3781 }
3782 case MeanStatistic:
3783 {
3784 GetMeanPixelList(pixel_list[id],&pixel);
3785 break;
3786 }
3787 case MedianStatistic:
3788 default:
3789 {
3790 GetMedianPixelList(pixel_list[id],&pixel);
3791 break;
3792 }
3793 case MinimumStatistic:
3794 {
3795 GetMinimumPixelList(pixel_list[id],&pixel);
3796 break;
3797 }
3798 case ModeStatistic:
3799 {
3800 GetModePixelList(pixel_list[id],&pixel);
3801 break;
3802 }
3803 case NonpeakStatistic:
3804 {
3805 GetNonpeakPixelList(pixel_list[id],&pixel);
3806 break;
3807 }
3808 case RootMeanSquareStatistic:
3809 {
3810 GetRootMeanSquarePixelList(pixel_list[id],&pixel);
3811 break;
3812 }
3813 case StandardDeviationStatistic:
3814 {
3815 GetStandardDeviationPixelList(pixel_list[id],&pixel);
3816 break;
3817 }
3818 }
3819 if ((channel & RedChannel) != 0)
3820 SetPixelRed(q,ClampToQuantum(pixel.red));
3821 if ((channel & GreenChannel) != 0)
3822 SetPixelGreen(q,ClampToQuantum(pixel.green));
3823 if ((channel & BlueChannel) != 0)
3824 SetPixelBlue(q,ClampToQuantum(pixel.blue));
3825 if ((channel & OpacityChannel) != 0)
3826 SetPixelOpacity(q,ClampToQuantum(pixel.opacity));
3827 if (((channel & IndexChannel) != 0) &&
3828 (image->colorspace == CMYKColorspace))
3829 SetPixelIndex(statistic_indexes+x,ClampToQuantum(pixel.index));
3830 p++;
3831 q++;
3832 }
3833 if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
3834 status=MagickFalse;
3835 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3836 {
3837 MagickBooleanType
3838 proceed;
3839
3840 proceed=SetImageProgress(image,StatisticImageTag,progress++,
3841 image->rows);
3842 if (proceed == MagickFalse)
3843 status=MagickFalse;
3844 }
3845 }
3846 statistic_view=DestroyCacheView(statistic_view);
3847 image_view=DestroyCacheView(image_view);
3848 pixel_list=DestroyPixelListTLS(pixel_list);
3849 if (status == MagickFalse)
3850 statistic_image=DestroyImage(statistic_image);
3851 return(statistic_image);
3852}
Definition: image.h:134