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