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