MagickCore 6.9.13
Loading...
Searching...
No Matches
statistic.c
1/*
2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3% %
4% %
5% %
6% SSSSS TTTTT AAA TTTTT IIIII SSSSS TTTTT IIIII CCCC %
7% SS T A A T I SS T I C %
8% SSS T AAAAA T I SSS T I C %
9% SS T A A T I SS T I C %
10% SSSSS T A A T IIIII SSSSS T IIIII CCCC %
11% %
12% %
13% MagickCore Image Statistical Methods %
14% %
15% Software Design %
16% Cristy %
17% July 1992 %
18% %
19% %
20% Copyright 1999 ImageMagick Studio LLC, a non-profit organization %
21% dedicated to making software imaging solutions freely available. %
22% %
23% You may not use this file except in compliance with the License. You may %
24% obtain a copy of the License at %
25% %
26% https://imagemagick.org/script/license.php %
27% %
28% Unless required by applicable law or agreed to in writing, software %
29% distributed under the License is distributed on an "AS IS" BASIS, %
30% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. %
31% See the License for the specific language governing permissions and %
32% limitations under the License. %
33% %
34%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35%
36%
37%
38*/
39␌
40/*
41 Include declarations.
42*/
43#include "magick/studio.h"
44#include "magick/accelerate-private.h"
45#include "magick/animate.h"
46#include "magick/animate.h"
47#include "magick/blob.h"
48#include "magick/blob-private.h"
49#include "magick/cache.h"
50#include "magick/cache-private.h"
51#include "magick/cache-view.h"
52#include "magick/client.h"
53#include "magick/color.h"
54#include "magick/color-private.h"
55#include "magick/colorspace.h"
56#include "magick/colorspace-private.h"
57#include "magick/composite.h"
58#include "magick/composite-private.h"
59#include "magick/compress.h"
60#include "magick/constitute.h"
61#include "magick/deprecate.h"
62#include "magick/display.h"
63#include "magick/draw.h"
64#include "magick/enhance.h"
65#include "magick/exception.h"
66#include "magick/exception-private.h"
67#include "magick/gem.h"
68#include "magick/geometry.h"
69#include "magick/list.h"
70#include "magick/image-private.h"
71#include "magick/magic.h"
72#include "magick/magick.h"
73#include "magick/memory_.h"
74#include "magick/module.h"
75#include "magick/monitor.h"
76#include "magick/monitor-private.h"
77#include "magick/option.h"
78#include "magick/paint.h"
79#include "magick/pixel-private.h"
80#include "magick/profile.h"
81#include "magick/property.h"
82#include "magick/quantize.h"
83#include "magick/random_.h"
84#include "magick/random-private.h"
85#include "magick/resource_.h"
86#include "magick/segment.h"
87#include "magick/semaphore.h"
88#include "magick/signature-private.h"
89#include "magick/statistic.h"
90#include "magick/statistic-private.h"
91#include "magick/string_.h"
92#include "magick/thread-private.h"
93#include "magick/timer.h"
94#include "magick/utility.h"
95#include "magick/version.h"
96␌
97/*
98%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
99% %
100% %
101% %
102% E v a l u a t e I m a g e %
103% %
104% %
105% %
106%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
107%
108% EvaluateImage() applies a value to the image with an arithmetic, relational,
109% or logical operator to an image. Use these operations to lighten or darken
110% an image, to increase or decrease contrast in an image, or to produce the
111% "negative" of an image.
112%
113% The format of the EvaluateImageChannel method is:
114%
115% MagickBooleanType EvaluateImage(Image *image,
116% const MagickEvaluateOperator op,const double value,
117% ExceptionInfo *exception)
118% MagickBooleanType EvaluateImages(Image *images,
119% const MagickEvaluateOperator op,const double value,
120% ExceptionInfo *exception)
121% MagickBooleanType EvaluateImageChannel(Image *image,
122% const ChannelType channel,const MagickEvaluateOperator op,
123% const double value,ExceptionInfo *exception)
124%
125% A description of each parameter follows:
126%
127% o image: the image.
128%
129% o channel: the channel.
130%
131% o op: A channel op.
132%
133% o value: A value value.
134%
135% o exception: return any errors or warnings in this structure.
136%
137*/
138
139static MagickPixelPacket **DestroyPixelTLS(const Image *images,
140 MagickPixelPacket **pixels)
141{
142 ssize_t
143 i;
144
145 size_t
146 rows;
147
148 assert(pixels != (MagickPixelPacket **) NULL);
149 rows=MagickMax(GetImageListLength(images),
150 (size_t) GetMagickResourceLimit(ThreadResource));
151 for (i=0; i < (ssize_t) rows; i++)
152 if (pixels[i] != (MagickPixelPacket *) NULL)
153 pixels[i]=(MagickPixelPacket *) RelinquishMagickMemory(pixels[i]);
154 pixels=(MagickPixelPacket **) RelinquishMagickMemory(pixels);
155 return(pixels);
156}
157
158static MagickPixelPacket **AcquirePixelTLS(const Image *images)
159{
160 const Image
161 *next;
162
163 MagickPixelPacket
164 **pixels;
165
166 ssize_t
167 i,
168 j;
169
170 size_t
171 columns,
172 rows;
173
174 rows=MagickMax(GetImageListLength(images),
175 (size_t) GetMagickResourceLimit(ThreadResource));
176 pixels=(MagickPixelPacket **) AcquireQuantumMemory(rows,sizeof(*pixels));
177 if (pixels == (MagickPixelPacket **) NULL)
178 return((MagickPixelPacket **) NULL);
179 (void) memset(pixels,0,rows*sizeof(*pixels));
180 columns=GetImageListLength(images);
181 for (next=images; next != (Image *) NULL; next=next->next)
182 columns=MagickMax(next->columns,columns);
183 for (i=0; i < (ssize_t) rows; i++)
184 {
185 pixels[i]=(MagickPixelPacket *) AcquireQuantumMemory(columns,
186 sizeof(**pixels));
187 if (pixels[i] == (MagickPixelPacket *) NULL)
188 return(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)*PerceptibleReciprocal(value);
302 break;
303 }
304 case LaplacianNoiseEvaluateOperator:
305 {
306 result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
307 LaplacianNoise,value);
308 break;
309 }
310 case LeftShiftEvaluateOperator:
311 {
312 result=(double) pixel;
313 for (i=0; i < (ssize_t) value; i++)
314 result*=2.0;
315 break;
316 }
317 case LogEvaluateOperator:
318 {
319 if ((QuantumScale*(MagickRealType) pixel) >= MagickEpsilon)
320 result=(MagickRealType) QuantumRange*log((double) (QuantumScale*value*
321 (MagickRealType) pixel+1.0))/log((double) (value+1.0));
322 break;
323 }
324 case MaxEvaluateOperator:
325 {
326 result=(MagickRealType) EvaluateMax((double) pixel,value);
327 break;
328 }
329 case MeanEvaluateOperator:
330 {
331 result=(MagickRealType) pixel+value;
332 break;
333 }
334 case MedianEvaluateOperator:
335 {
336 result=(MagickRealType) pixel+value;
337 break;
338 }
339 case MinEvaluateOperator:
340 {
341 result=(MagickRealType) MagickMin((double) pixel,value);
342 break;
343 }
344 case MultiplicativeNoiseEvaluateOperator:
345 {
346 result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
347 MultiplicativeGaussianNoise,value);
348 break;
349 }
350 case MultiplyEvaluateOperator:
351 {
352 result=(MagickRealType) pixel*value;
353 break;
354 }
355 case OrEvaluateOperator:
356 {
357 result=(MagickRealType) ((ssize_t) pixel | (ssize_t) (value+0.5));
358 break;
359 }
360 case PoissonNoiseEvaluateOperator:
361 {
362 result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
363 PoissonNoise,value);
364 break;
365 }
366 case PowEvaluateOperator:
367 {
368 if (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*PerceptibleReciprocal(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 PerceptibleReciprocal(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 PerceptibleReciprocal(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]*PerceptibleReciprocal(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*PerceptibleReciprocal(
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]=(-MagickLog10(moments[channel].I[i]));
2168 moments=(ChannelMoments *) RelinquishMagickMemory(moments);
2169 /*
2170 Blur then transform to HCLp colorspace.
2171 */
2172 hash_image=BlurImage(image,0.0,1.0,exception);
2173 if (hash_image == (Image *) NULL)
2174 {
2175 perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
2176 perceptual_hash);
2177 return((ChannelPerceptualHash *) NULL);
2178 }
2179 hash_image->depth=8;
2180 status=TransformImageColorspace(hash_image,HSBColorspace);
2181 if (status == MagickFalse)
2182 {
2183 perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
2184 perceptual_hash);
2185 return((ChannelPerceptualHash *) NULL);
2186 }
2187 moments=GetImageChannelMoments(hash_image,exception);
2188 hash_image=DestroyImage(hash_image);
2189 if (moments == (ChannelMoments *) NULL)
2190 {
2191 perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
2192 perceptual_hash);
2193 return((ChannelPerceptualHash *) NULL);
2194 }
2195 for (channel=0; channel <= CompositeChannels; channel++)
2196 for (i=0; i < MaximumNumberOfImageMoments; i++)
2197 perceptual_hash[channel].Q[i]=(-MagickLog10(moments[channel].I[i]));
2198 moments=(ChannelMoments *) RelinquishMagickMemory(moments);
2199 return(perceptual_hash);
2200}
2201␌
2202/*
2203%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2204% %
2205% %
2206% %
2207% G e t I m a g e C h a n n e l R a n g e %
2208% %
2209% %
2210% %
2211%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2212%
2213% GetImageChannelRange() returns the range of one or more image channels.
2214%
2215% The format of the GetImageChannelRange method is:
2216%
2217% MagickBooleanType GetImageChannelRange(const Image *image,
2218% const ChannelType channel,double *minima,double *maxima,
2219% ExceptionInfo *exception)
2220%
2221% A description of each parameter follows:
2222%
2223% o image: the image.
2224%
2225% o channel: the channel.
2226%
2227% o minima: the minimum value in the channel.
2228%
2229% o maxima: the maximum value in the channel.
2230%
2231% o exception: return any errors or warnings in this structure.
2232%
2233*/
2234
2235MagickExport MagickBooleanType GetImageRange(const Image *image,
2236 double *minima,double *maxima,ExceptionInfo *exception)
2237{
2238 return(GetImageChannelRange(image,CompositeChannels,minima,maxima,exception));
2239}
2240
2241MagickExport MagickBooleanType GetImageChannelRange(const Image *image,
2242 const ChannelType channel,double *minima,double *maxima,
2243 ExceptionInfo *exception)
2244{
2245 MagickPixelPacket
2246 pixel;
2247
2248 ssize_t
2249 y;
2250
2251 assert(image != (Image *) NULL);
2252 assert(image->signature == MagickCoreSignature);
2253 if (IsEventLogging() != MagickFalse)
2254 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2255 *maxima=(-MagickMaximumValue);
2256 *minima=MagickMaximumValue;
2257 GetMagickPixelPacket(image,&pixel);
2258 for (y=0; y < (ssize_t) image->rows; y++)
2259 {
2260 const IndexPacket
2261 *magick_restrict indexes;
2262
2263 const PixelPacket
2264 *magick_restrict p;
2265
2266 ssize_t
2267 x;
2268
2269 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
2270 if (p == (const PixelPacket *) NULL)
2271 break;
2272 indexes=GetVirtualIndexQueue(image);
2273 for (x=0; x < (ssize_t) image->columns; x++)
2274 {
2275 SetMagickPixelPacket(image,p,indexes+x,&pixel);
2276 if ((channel & RedChannel) != 0)
2277 {
2278 if (pixel.red < *minima)
2279 *minima=(double) pixel.red;
2280 if (pixel.red > *maxima)
2281 *maxima=(double) pixel.red;
2282 }
2283 if ((channel & GreenChannel) != 0)
2284 {
2285 if (pixel.green < *minima)
2286 *minima=(double) pixel.green;
2287 if (pixel.green > *maxima)
2288 *maxima=(double) pixel.green;
2289 }
2290 if ((channel & BlueChannel) != 0)
2291 {
2292 if (pixel.blue < *minima)
2293 *minima=(double) pixel.blue;
2294 if (pixel.blue > *maxima)
2295 *maxima=(double) pixel.blue;
2296 }
2297 if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
2298 {
2299 if (((MagickRealType) QuantumRange-(MagickRealType) pixel.opacity) < *minima)
2300 *minima=(double) ((MagickRealType) QuantumRange-(MagickRealType)
2301 pixel.opacity);
2302 if (((MagickRealType) QuantumRange-(MagickRealType) pixel.opacity) > *maxima)
2303 *maxima=(double) ((MagickRealType) QuantumRange-(MagickRealType)
2304 pixel.opacity);
2305 }
2306 if (((channel & IndexChannel) != 0) &&
2307 (image->colorspace == CMYKColorspace))
2308 {
2309 if ((double) pixel.index < *minima)
2310 *minima=(double) pixel.index;
2311 if ((double) pixel.index > *maxima)
2312 *maxima=(double) pixel.index;
2313 }
2314 p++;
2315 }
2316 }
2317 return(y == (ssize_t) image->rows ? MagickTrue : MagickFalse);
2318}
2319␌
2320/*
2321%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2322% %
2323% %
2324% %
2325% 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 %
2326% %
2327% %
2328% %
2329%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2330%
2331% GetImageChannelStatistics() returns statistics for each channel in the
2332% image. The statistics include the channel depth, its minima, maxima, mean,
2333% standard deviation, kurtosis and skewness. You can access the red channel
2334% mean, for example, like this:
2335%
2336% channel_statistics=GetImageChannelStatistics(image,exception);
2337% red_mean=channel_statistics[RedChannel].mean;
2338%
2339% Use MagickRelinquishMemory() to free the statistics buffer.
2340%
2341% The format of the GetImageChannelStatistics method is:
2342%
2343% ChannelStatistics *GetImageChannelStatistics(const Image *image,
2344% ExceptionInfo *exception)
2345%
2346% A description of each parameter follows:
2347%
2348% o image: the image.
2349%
2350% o exception: return any errors or warnings in this structure.
2351%
2352*/
2353MagickExport ChannelStatistics *GetImageChannelStatistics(const Image *image,
2354 ExceptionInfo *exception)
2355{
2356 ChannelStatistics
2357 *channel_statistics;
2358
2359 double
2360 area,
2361 standard_deviation;
2362
2363 MagickPixelPacket
2364 number_bins,
2365 *histogram;
2366
2367 QuantumAny
2368 range;
2369
2370 size_t
2371 channels,
2372 depth,
2373 length;
2374
2375 ssize_t
2376 i,
2377 y;
2378
2379 assert(image != (Image *) NULL);
2380 assert(image->signature == MagickCoreSignature);
2381 if (IsEventLogging() != MagickFalse)
2382 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2383 length=CompositeChannels+1UL;
2384 channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(length,
2385 sizeof(*channel_statistics));
2386 histogram=(MagickPixelPacket *) AcquireQuantumMemory(MaxMap+1U,
2387 sizeof(*histogram));
2388 if ((channel_statistics == (ChannelStatistics *) NULL) ||
2389 (histogram == (MagickPixelPacket *) NULL))
2390 {
2391 if (histogram != (MagickPixelPacket *) NULL)
2392 histogram=(MagickPixelPacket *) RelinquishMagickMemory(histogram);
2393 if (channel_statistics != (ChannelStatistics *) NULL)
2394 channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
2395 channel_statistics);
2396 return(channel_statistics);
2397 }
2398 (void) memset(channel_statistics,0,length*
2399 sizeof(*channel_statistics));
2400 for (i=0; i <= (ssize_t) CompositeChannels; i++)
2401 {
2402 channel_statistics[i].depth=1;
2403 channel_statistics[i].maxima=(-MagickMaximumValue);
2404 channel_statistics[i].minima=MagickMaximumValue;
2405 }
2406 (void) memset(histogram,0,(MaxMap+1U)*sizeof(*histogram));
2407 (void) memset(&number_bins,0,sizeof(number_bins));
2408 for (y=0; y < (ssize_t) image->rows; y++)
2409 {
2410 const IndexPacket
2411 *magick_restrict indexes;
2412
2413 const PixelPacket
2414 *magick_restrict p;
2415
2416 ssize_t
2417 x;
2418
2419 /*
2420 Compute pixel statistics.
2421 */
2422 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
2423 if (p == (const PixelPacket *) NULL)
2424 break;
2425 indexes=GetVirtualIndexQueue(image);
2426 for (x=0; x < (ssize_t) image->columns; )
2427 {
2428 if (channel_statistics[RedChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2429 {
2430 depth=channel_statistics[RedChannel].depth;
2431 range=GetQuantumRange(depth);
2432 if (IsPixelAtDepth(GetPixelRed(p),range) == MagickFalse)
2433 {
2434 channel_statistics[RedChannel].depth++;
2435 continue;
2436 }
2437 }
2438 if (channel_statistics[GreenChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2439 {
2440 depth=channel_statistics[GreenChannel].depth;
2441 range=GetQuantumRange(depth);
2442 if (IsPixelAtDepth(GetPixelGreen(p),range) == MagickFalse)
2443 {
2444 channel_statistics[GreenChannel].depth++;
2445 continue;
2446 }
2447 }
2448 if (channel_statistics[BlueChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2449 {
2450 depth=channel_statistics[BlueChannel].depth;
2451 range=GetQuantumRange(depth);
2452 if (IsPixelAtDepth(GetPixelBlue(p),range) == MagickFalse)
2453 {
2454 channel_statistics[BlueChannel].depth++;
2455 continue;
2456 }
2457 }
2458 if (image->matte != MagickFalse)
2459 {
2460 if (channel_statistics[OpacityChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2461 {
2462 depth=channel_statistics[OpacityChannel].depth;
2463 range=GetQuantumRange(depth);
2464 if (IsPixelAtDepth(GetPixelAlpha(p),range) == MagickFalse)
2465 {
2466 channel_statistics[OpacityChannel].depth++;
2467 continue;
2468 }
2469 }
2470 }
2471 if (image->colorspace == CMYKColorspace)
2472 {
2473 if (channel_statistics[BlackChannel].depth != MAGICKCORE_QUANTUM_DEPTH)
2474 {
2475 depth=channel_statistics[BlackChannel].depth;
2476 range=GetQuantumRange(depth);
2477 if (IsPixelAtDepth(GetPixelIndex(indexes+x),range) == MagickFalse)
2478 {
2479 channel_statistics[BlackChannel].depth++;
2480 continue;
2481 }
2482 }
2483 }
2484 if ((double) GetPixelRed(p) < channel_statistics[RedChannel].minima)
2485 channel_statistics[RedChannel].minima=(double) GetPixelRed(p);
2486 if ((double) GetPixelRed(p) > channel_statistics[RedChannel].maxima)
2487 channel_statistics[RedChannel].maxima=(double) GetPixelRed(p);
2488 channel_statistics[RedChannel].sum+=QuantumScale*GetPixelRed(p);
2489 channel_statistics[RedChannel].sum_squared+=QuantumScale*GetPixelRed(p)*
2490 QuantumScale*GetPixelRed(p);
2491 channel_statistics[RedChannel].sum_cubed+=QuantumScale*GetPixelRed(p)*
2492 QuantumScale*GetPixelRed(p)*QuantumScale*GetPixelRed(p);
2493 channel_statistics[RedChannel].sum_fourth_power+=QuantumScale*
2494 GetPixelRed(p)*QuantumScale*GetPixelRed(p)*QuantumScale*GetPixelRed(p)*
2495 QuantumScale*GetPixelRed(p);
2496 if ((double) GetPixelGreen(p) < channel_statistics[GreenChannel].minima)
2497 channel_statistics[GreenChannel].minima=(double) GetPixelGreen(p);
2498 if ((double) GetPixelGreen(p) > channel_statistics[GreenChannel].maxima)
2499 channel_statistics[GreenChannel].maxima=(double) GetPixelGreen(p);
2500 channel_statistics[GreenChannel].sum+=QuantumScale*GetPixelGreen(p);
2501 channel_statistics[GreenChannel].sum_squared+=QuantumScale*GetPixelGreen(p)*
2502 QuantumScale*GetPixelGreen(p);
2503 channel_statistics[GreenChannel].sum_cubed+=QuantumScale*GetPixelGreen(p)*
2504 QuantumScale*GetPixelGreen(p)*QuantumScale*GetPixelGreen(p);
2505 channel_statistics[GreenChannel].sum_fourth_power+=QuantumScale*
2506 GetPixelGreen(p)*QuantumScale*GetPixelGreen(p)*QuantumScale*
2507 GetPixelGreen(p)*QuantumScale*GetPixelGreen(p);
2508 if ((double) GetPixelBlue(p) < channel_statistics[BlueChannel].minima)
2509 channel_statistics[BlueChannel].minima=(double) GetPixelBlue(p);
2510 if ((double) GetPixelBlue(p) > channel_statistics[BlueChannel].maxima)
2511 channel_statistics[BlueChannel].maxima=(double) GetPixelBlue(p);
2512 channel_statistics[BlueChannel].sum+=QuantumScale*GetPixelBlue(p);
2513 channel_statistics[BlueChannel].sum_squared+=QuantumScale*GetPixelBlue(p)*
2514 QuantumScale*GetPixelBlue(p);
2515 channel_statistics[BlueChannel].sum_cubed+=QuantumScale*GetPixelBlue(p)*
2516 QuantumScale*GetPixelBlue(p)*QuantumScale*GetPixelBlue(p);
2517 channel_statistics[BlueChannel].sum_fourth_power+=QuantumScale*
2518 GetPixelBlue(p)*QuantumScale*GetPixelBlue(p)*QuantumScale*
2519 GetPixelBlue(p)*QuantumScale*GetPixelBlue(p);
2520 histogram[ScaleQuantumToMap(GetPixelRed(p))].red++;
2521 histogram[ScaleQuantumToMap(GetPixelGreen(p))].green++;
2522 histogram[ScaleQuantumToMap(GetPixelBlue(p))].blue++;
2523 if (image->matte != MagickFalse)
2524 {
2525 if ((double) GetPixelAlpha(p) < channel_statistics[OpacityChannel].minima)
2526 channel_statistics[OpacityChannel].minima=(double) GetPixelAlpha(p);
2527 if ((double) GetPixelAlpha(p) > channel_statistics[OpacityChannel].maxima)
2528 channel_statistics[OpacityChannel].maxima=(double) GetPixelAlpha(p);
2529 channel_statistics[OpacityChannel].sum+=QuantumScale*GetPixelAlpha(p);
2530 channel_statistics[OpacityChannel].sum_squared+=QuantumScale*
2531 GetPixelAlpha(p)*QuantumScale*GetPixelAlpha(p);
2532 channel_statistics[OpacityChannel].sum_cubed+=QuantumScale*
2533 GetPixelAlpha(p)*QuantumScale*GetPixelAlpha(p)*QuantumScale*
2534 GetPixelAlpha(p);
2535 channel_statistics[OpacityChannel].sum_fourth_power+=QuantumScale*
2536 GetPixelAlpha(p)*QuantumScale*GetPixelAlpha(p)*QuantumScale*
2537 GetPixelAlpha(p)*QuantumScale*GetPixelAlpha(p);
2538 histogram[ScaleQuantumToMap(GetPixelAlpha(p))].opacity++;
2539 }
2540 if (image->colorspace == CMYKColorspace)
2541 {
2542 if ((double) GetPixelIndex(indexes+x) < channel_statistics[BlackChannel].minima)
2543 channel_statistics[BlackChannel].minima=(double)
2544 GetPixelIndex(indexes+x);
2545 if ((double) GetPixelIndex(indexes+x) > channel_statistics[BlackChannel].maxima)
2546 channel_statistics[BlackChannel].maxima=(double)
2547 GetPixelIndex(indexes+x);
2548 channel_statistics[BlackChannel].sum+=QuantumScale*
2549 GetPixelIndex(indexes+x);
2550 channel_statistics[BlackChannel].sum_squared+=QuantumScale*
2551 GetPixelIndex(indexes+x)*QuantumScale*GetPixelIndex(indexes+x);
2552 channel_statistics[BlackChannel].sum_cubed+=QuantumScale*
2553 GetPixelIndex(indexes+x)*QuantumScale*GetPixelIndex(indexes+x)*
2554 QuantumScale*GetPixelIndex(indexes+x);
2555 channel_statistics[BlackChannel].sum_fourth_power+=QuantumScale*
2556 GetPixelIndex(indexes+x)*QuantumScale*GetPixelIndex(indexes+x)*
2557 QuantumScale*GetPixelIndex(indexes+x)*QuantumScale*
2558 GetPixelIndex(indexes+x);
2559 histogram[ScaleQuantumToMap(GetPixelIndex(indexes+x))].index++;
2560 }
2561 x++;
2562 p++;
2563 }
2564 }
2565 for (i=0; i < (ssize_t) CompositeChannels; i++)
2566 {
2567 double
2568 area,
2569 mean,
2570 standard_deviation;
2571
2572 /*
2573 Normalize pixel statistics.
2574 */
2575 area=PerceptibleReciprocal((double) image->columns*image->rows);
2576 mean=channel_statistics[i].sum*area;
2577 channel_statistics[i].sum=mean;
2578 channel_statistics[i].sum_squared*=area;
2579 channel_statistics[i].sum_cubed*=area;
2580 channel_statistics[i].sum_fourth_power*=area;
2581 channel_statistics[i].mean=mean;
2582 channel_statistics[i].variance=channel_statistics[i].sum_squared;
2583 standard_deviation=sqrt(channel_statistics[i].variance-(mean*mean));
2584 area=PerceptibleReciprocal((double) image->columns*image->rows-1.0)*
2585 ((double) image->columns*image->rows);
2586 standard_deviation=sqrt(area*standard_deviation*standard_deviation);
2587 channel_statistics[i].standard_deviation=standard_deviation;
2588 }
2589 for (i=0; i < (ssize_t) (MaxMap+1U); i++)
2590 {
2591 if (histogram[i].red > 0.0)
2592 number_bins.red++;
2593 if (histogram[i].green > 0.0)
2594 number_bins.green++;
2595 if (histogram[i].blue > 0.0)
2596 number_bins.blue++;
2597 if ((image->matte != MagickFalse) && (histogram[i].opacity > 0.0))
2598 number_bins.opacity++;
2599 if ((image->colorspace == CMYKColorspace) && (histogram[i].index > 0.0))
2600 number_bins.index++;
2601 }
2602 area=PerceptibleReciprocal((double) image->columns*image->rows);
2603 for (i=0; i < (ssize_t) (MaxMap+1U); i++)
2604 {
2605 /*
2606 Compute pixel entropy.
2607 */
2608 histogram[i].red*=area;
2609 channel_statistics[RedChannel].entropy+=-histogram[i].red*
2610 MagickLog10(histogram[i].red)*
2611 PerceptibleReciprocal(MagickLog10((double) number_bins.red));
2612 histogram[i].green*=area;
2613 channel_statistics[GreenChannel].entropy+=-histogram[i].green*
2614 MagickLog10(histogram[i].green)*
2615 PerceptibleReciprocal(MagickLog10((double) number_bins.green));
2616 histogram[i].blue*=area;
2617 channel_statistics[BlueChannel].entropy+=-histogram[i].blue*
2618 MagickLog10(histogram[i].blue)*
2619 PerceptibleReciprocal(MagickLog10((double) number_bins.blue));
2620 if (image->matte != MagickFalse)
2621 {
2622 histogram[i].opacity*=area;
2623 channel_statistics[OpacityChannel].entropy+=-histogram[i].opacity*
2624 MagickLog10(histogram[i].opacity)*
2625 PerceptibleReciprocal(MagickLog10((double) number_bins.opacity));
2626 }
2627 if (image->colorspace == CMYKColorspace)
2628 {
2629 histogram[i].index*=area;
2630 channel_statistics[IndexChannel].entropy+=-histogram[i].index*
2631 MagickLog10(histogram[i].index)*
2632 PerceptibleReciprocal(MagickLog10((double) number_bins.index));
2633 }
2634 }
2635 /*
2636 Compute overall statistics.
2637 */
2638 for (i=0; i < (ssize_t) CompositeChannels; i++)
2639 {
2640 channel_statistics[CompositeChannels].depth=(size_t) EvaluateMax((double)
2641 channel_statistics[CompositeChannels].depth,(double)
2642 channel_statistics[i].depth);
2643 channel_statistics[CompositeChannels].minima=MagickMin(
2644 channel_statistics[CompositeChannels].minima,
2645 channel_statistics[i].minima);
2646 channel_statistics[CompositeChannels].maxima=EvaluateMax(
2647 channel_statistics[CompositeChannels].maxima,
2648 channel_statistics[i].maxima);
2649 channel_statistics[CompositeChannels].sum+=channel_statistics[i].sum;
2650 channel_statistics[CompositeChannels].sum_squared+=
2651 channel_statistics[i].sum_squared;
2652 channel_statistics[CompositeChannels].sum_cubed+=
2653 channel_statistics[i].sum_cubed;
2654 channel_statistics[CompositeChannels].sum_fourth_power+=
2655 channel_statistics[i].sum_fourth_power;
2656 channel_statistics[CompositeChannels].mean+=channel_statistics[i].mean;
2657 channel_statistics[CompositeChannels].variance+=
2658 channel_statistics[i].variance-channel_statistics[i].mean*
2659 channel_statistics[i].mean;
2660 standard_deviation=sqrt(channel_statistics[i].variance-
2661 (channel_statistics[i].mean*channel_statistics[i].mean));
2662 area=PerceptibleReciprocal((double) image->columns*image->rows-1.0)*
2663 ((double) image->columns*image->rows);
2664 standard_deviation=sqrt(area*standard_deviation*standard_deviation);
2665 channel_statistics[CompositeChannels].standard_deviation=standard_deviation;
2666 channel_statistics[CompositeChannels].entropy+=
2667 channel_statistics[i].entropy;
2668 }
2669 channels=3;
2670 if (image->matte != MagickFalse)
2671 channels++;
2672 if (image->colorspace == CMYKColorspace)
2673 channels++;
2674 channel_statistics[CompositeChannels].sum/=channels;
2675 channel_statistics[CompositeChannels].sum_squared/=channels;
2676 channel_statistics[CompositeChannels].sum_cubed/=channels;
2677 channel_statistics[CompositeChannels].sum_fourth_power/=channels;
2678 channel_statistics[CompositeChannels].mean/=channels;
2679 channel_statistics[CompositeChannels].kurtosis/=channels;
2680 channel_statistics[CompositeChannels].skewness/=channels;
2681 channel_statistics[CompositeChannels].entropy/=channels;
2682 i=CompositeChannels;
2683 area=PerceptibleReciprocal((double) channels*image->columns*image->rows);
2684 channel_statistics[i].variance=channel_statistics[i].sum_squared;
2685 channel_statistics[i].mean=channel_statistics[i].sum;
2686 standard_deviation=sqrt(channel_statistics[i].variance-
2687 (channel_statistics[i].mean*channel_statistics[i].mean));
2688 standard_deviation=sqrt(PerceptibleReciprocal((double) channels*
2689 image->columns*image->rows-1.0)*channels*image->columns*image->rows*
2690 standard_deviation*standard_deviation);
2691 channel_statistics[i].standard_deviation=standard_deviation;
2692 for (i=0; i <= (ssize_t) CompositeChannels; i++)
2693 {
2694 /*
2695 Compute kurtosis & skewness statistics.
2696 */
2697 standard_deviation=PerceptibleReciprocal(
2698 channel_statistics[i].standard_deviation);
2699 channel_statistics[i].skewness=(channel_statistics[i].sum_cubed-3.0*
2700 channel_statistics[i].mean*channel_statistics[i].sum_squared+2.0*
2701 channel_statistics[i].mean*channel_statistics[i].mean*
2702 channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2703 standard_deviation);
2704 channel_statistics[i].kurtosis=(channel_statistics[i].sum_fourth_power-4.0*
2705 channel_statistics[i].mean*channel_statistics[i].sum_cubed+6.0*
2706 channel_statistics[i].mean*channel_statistics[i].mean*
2707 channel_statistics[i].sum_squared-3.0*channel_statistics[i].mean*
2708 channel_statistics[i].mean*1.0*channel_statistics[i].mean*
2709 channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2710 standard_deviation*standard_deviation)-3.0;
2711 }
2712 for (i=0; i <= (ssize_t) CompositeChannels; i++)
2713 {
2714 channel_statistics[i].mean*=QuantumRange;
2715 channel_statistics[i].variance*=QuantumRange;
2716 channel_statistics[i].standard_deviation*=QuantumRange;
2717 channel_statistics[i].sum*=QuantumRange;
2718 channel_statistics[i].sum_squared*=QuantumRange;
2719 channel_statistics[i].sum_cubed*=QuantumRange;
2720 channel_statistics[i].sum_fourth_power*=QuantumRange;
2721 }
2722 channel_statistics[CompositeChannels].mean=0.0;
2723 channel_statistics[CompositeChannels].standard_deviation=0.0;
2724 for (i=0; i < (ssize_t) CompositeChannels; i++)
2725 {
2726 channel_statistics[CompositeChannels].mean+=
2727 channel_statistics[i].mean;
2728 channel_statistics[CompositeChannels].standard_deviation+=
2729 channel_statistics[i].standard_deviation;
2730 }
2731 channel_statistics[CompositeChannels].mean/=(double) channels;
2732 channel_statistics[CompositeChannels].standard_deviation/=(double) channels;
2733 histogram=(MagickPixelPacket *) RelinquishMagickMemory(histogram);
2734 if (y < (ssize_t) image->rows)
2735 channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
2736 channel_statistics);
2737 return(channel_statistics);
2738}
2739␌
2740/*
2741%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2742% %
2743% %
2744% %
2745% P o l y n o m i a l I m a g e %
2746% %
2747% %
2748% %
2749%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2750%
2751% PolynomialImage() returns a new image where each pixel is the sum of the
2752% pixels in the image sequence after applying its corresponding terms
2753% (coefficient and degree pairs).
2754%
2755% The format of the PolynomialImage method is:
2756%
2757% Image *PolynomialImage(const Image *images,const size_t number_terms,
2758% const double *terms,ExceptionInfo *exception)
2759% Image *PolynomialImageChannel(const Image *images,
2760% const size_t number_terms,const ChannelType channel,
2761% const double *terms,ExceptionInfo *exception)
2762%
2763% A description of each parameter follows:
2764%
2765% o images: the image sequence.
2766%
2767% o channel: the channel.
2768%
2769% o number_terms: the number of terms in the list. The actual list length
2770% is 2 x number_terms + 1 (the constant).
2771%
2772% o terms: the list of polynomial coefficients and degree pairs and a
2773% constant.
2774%
2775% o exception: return any errors or warnings in this structure.
2776%
2777*/
2778MagickExport Image *PolynomialImage(const Image *images,
2779 const size_t number_terms,const double *terms,ExceptionInfo *exception)
2780{
2781 Image
2782 *polynomial_image;
2783
2784 polynomial_image=PolynomialImageChannel(images,DefaultChannels,number_terms,
2785 terms,exception);
2786 return(polynomial_image);
2787}
2788
2789MagickExport Image *PolynomialImageChannel(const Image *images,
2790 const ChannelType channel,const size_t number_terms,const double *terms,
2791 ExceptionInfo *exception)
2792{
2793#define PolynomialImageTag "Polynomial/Image"
2794
2795 CacheView
2796 *polynomial_view;
2797
2798 Image
2799 *image;
2800
2801 MagickBooleanType
2802 status;
2803
2804 MagickOffsetType
2805 progress;
2806
2807 MagickPixelPacket
2808 **magick_restrict polynomial_pixels,
2809 zero;
2810
2811 ssize_t
2812 y;
2813
2814 assert(images != (Image *) NULL);
2815 assert(images->signature == MagickCoreSignature);
2816 if (IsEventLogging() != MagickFalse)
2817 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
2818 assert(exception != (ExceptionInfo *) NULL);
2819 assert(exception->signature == MagickCoreSignature);
2820 image=AcquireImageCanvas(images,exception);
2821 if (image == (Image *) NULL)
2822 return((Image *) NULL);
2823 if (SetImageStorageClass(image,DirectClass) == MagickFalse)
2824 {
2825 InheritException(exception,&image->exception);
2826 image=DestroyImage(image);
2827 return((Image *) NULL);
2828 }
2829 polynomial_pixels=AcquirePixelTLS(images);
2830 if (polynomial_pixels == (MagickPixelPacket **) NULL)
2831 {
2832 image=DestroyImage(image);
2833 (void) ThrowMagickException(exception,GetMagickModule(),
2834 ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
2835 return((Image *) NULL);
2836 }
2837 /*
2838 Polynomial image pixels.
2839 */
2840 status=MagickTrue;
2841 progress=0;
2842 GetMagickPixelPacket(images,&zero);
2843 polynomial_view=AcquireAuthenticCacheView(image,exception);
2844#if defined(MAGICKCORE_OPENMP_SUPPORT)
2845 #pragma omp parallel for schedule(static) shared(progress,status) \
2846 magick_number_threads(image,image,image->rows,1)
2847#endif
2848 for (y=0; y < (ssize_t) image->rows; y++)
2849 {
2850 CacheView
2851 *image_view;
2852
2853 const Image
2854 *next;
2855
2856 const int
2857 id = GetOpenMPThreadId();
2858
2859 IndexPacket
2860 *magick_restrict polynomial_indexes;
2861
2862 MagickPixelPacket
2863 *polynomial_pixel;
2864
2865 PixelPacket
2866 *magick_restrict q;
2867
2868 ssize_t
2869 i,
2870 x;
2871
2872 size_t
2873 number_images;
2874
2875 if (status == MagickFalse)
2876 continue;
2877 q=QueueCacheViewAuthenticPixels(polynomial_view,0,y,image->columns,1,
2878 exception);
2879 if (q == (PixelPacket *) NULL)
2880 {
2881 status=MagickFalse;
2882 continue;
2883 }
2884 polynomial_indexes=GetCacheViewAuthenticIndexQueue(polynomial_view);
2885 polynomial_pixel=polynomial_pixels[id];
2886 for (x=0; x < (ssize_t) image->columns; x++)
2887 polynomial_pixel[x]=zero;
2888 next=images;
2889 number_images=GetImageListLength(images);
2890 for (i=0; i < (ssize_t) number_images; i++)
2891 {
2892 const IndexPacket
2893 *indexes;
2894
2895 const PixelPacket
2896 *p;
2897
2898 if (i >= (ssize_t) number_terms)
2899 break;
2900 image_view=AcquireVirtualCacheView(next,exception);
2901 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2902 if (p == (const PixelPacket *) NULL)
2903 {
2904 image_view=DestroyCacheView(image_view);
2905 break;
2906 }
2907 indexes=GetCacheViewVirtualIndexQueue(image_view);
2908 for (x=0; x < (ssize_t) image->columns; x++)
2909 {
2910 double
2911 coefficient,
2912 degree;
2913
2914 coefficient=terms[i << 1];
2915 degree=terms[(i << 1)+1];
2916 if ((channel & RedChannel) != 0)
2917 polynomial_pixel[x].red+=coefficient*pow(QuantumScale*(double)
2918 p->red,degree);
2919 if ((channel & GreenChannel) != 0)
2920 polynomial_pixel[x].green+=coefficient*pow(QuantumScale*(double)
2921 p->green,
2922 degree);
2923 if ((channel & BlueChannel) != 0)
2924 polynomial_pixel[x].blue+=coefficient*pow(QuantumScale*(double)
2925 p->blue,degree);
2926 if ((channel & OpacityChannel) != 0)
2927 polynomial_pixel[x].opacity+=coefficient*pow(QuantumScale*
2928 ((double) QuantumRange-(double) p->opacity),degree);
2929 if (((channel & IndexChannel) != 0) &&
2930 (image->colorspace == CMYKColorspace))
2931 polynomial_pixel[x].index+=coefficient*pow(QuantumScale*(double)
2932 indexes[x],degree);
2933 p++;
2934 }
2935 image_view=DestroyCacheView(image_view);
2936 next=GetNextImageInList(next);
2937 }
2938 for (x=0; x < (ssize_t) image->columns; x++)
2939 {
2940 SetPixelRed(q,ClampToQuantum((MagickRealType) QuantumRange*
2941 polynomial_pixel[x].red));
2942 SetPixelGreen(q,ClampToQuantum((MagickRealType) QuantumRange*
2943 polynomial_pixel[x].green));
2944 SetPixelBlue(q,ClampToQuantum((MagickRealType) QuantumRange*
2945 polynomial_pixel[x].blue));
2946 if (image->matte == MagickFalse)
2947 SetPixelOpacity(q,ClampToQuantum((MagickRealType) QuantumRange-
2948 (MagickRealType) QuantumRange*polynomial_pixel[x].opacity));
2949 else
2950 SetPixelAlpha(q,ClampToQuantum((MagickRealType) QuantumRange-
2951 (MagickRealType) QuantumRange*polynomial_pixel[x].opacity));
2952 if (image->colorspace == CMYKColorspace)
2953 SetPixelIndex(polynomial_indexes+x,ClampToQuantum((MagickRealType)
2954 QuantumRange*polynomial_pixel[x].index));
2955 q++;
2956 }
2957 if (SyncCacheViewAuthenticPixels(polynomial_view,exception) == MagickFalse)
2958 status=MagickFalse;
2959 if (images->progress_monitor != (MagickProgressMonitor) NULL)
2960 {
2961 MagickBooleanType
2962 proceed;
2963
2964 proceed=SetImageProgress(images,PolynomialImageTag,progress++,
2965 image->rows);
2966 if (proceed == MagickFalse)
2967 status=MagickFalse;
2968 }
2969 }
2970 polynomial_view=DestroyCacheView(polynomial_view);
2971 polynomial_pixels=DestroyPixelTLS(images,polynomial_pixels);
2972 if (status == MagickFalse)
2973 image=DestroyImage(image);
2974 return(image);
2975}
2976␌
2977/*
2978%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2979% %
2980% %
2981% %
2982% S t a t i s t i c I m a g e %
2983% %
2984% %
2985% %
2986%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2987%
2988% StatisticImage() makes each pixel the min / max / median / mode / etc. of
2989% the neighborhood of the specified width and height.
2990%
2991% The format of the StatisticImage method is:
2992%
2993% Image *StatisticImage(const Image *image,const StatisticType type,
2994% const size_t width,const size_t height,ExceptionInfo *exception)
2995% Image *StatisticImageChannel(const Image *image,
2996% const ChannelType channel,const StatisticType type,
2997% const size_t width,const size_t height,ExceptionInfo *exception)
2998%
2999% A description of each parameter follows:
3000%
3001% o image: the image.
3002%
3003% o channel: the image channel.
3004%
3005% o type: the statistic type (median, mode, etc.).
3006%
3007% o width: the width of the pixel neighborhood.
3008%
3009% o height: the height of the pixel neighborhood.
3010%
3011% o exception: return any errors or warnings in this structure.
3012%
3013*/
3014
3015#define ListChannels 5
3016
3017typedef struct _ListNode
3018{
3019 size_t
3020 next[9],
3021 count,
3022 signature;
3023} ListNode;
3024
3025typedef struct _SkipList
3026{
3027 ssize_t
3028 level;
3029
3030 ListNode
3031 *nodes;
3032} SkipList;
3033
3034typedef struct _PixelList
3035{
3036 size_t
3037 length,
3038 seed,
3039 signature;
3040
3041 SkipList
3042 lists[ListChannels];
3043} PixelList;
3044
3045static PixelList *DestroyPixelList(PixelList *pixel_list)
3046{
3047 ssize_t
3048 i;
3049
3050 if (pixel_list == (PixelList *) NULL)
3051 return((PixelList *) NULL);
3052 for (i=0; i < ListChannels; i++)
3053 if (pixel_list->lists[i].nodes != (ListNode *) NULL)
3054 pixel_list->lists[i].nodes=(ListNode *) RelinquishAlignedMemory(
3055 pixel_list->lists[i].nodes);
3056 pixel_list=(PixelList *) RelinquishMagickMemory(pixel_list);
3057 return(pixel_list);
3058}
3059
3060static PixelList **DestroyPixelListTLS(PixelList **pixel_list)
3061{
3062 ssize_t
3063 i;
3064
3065 assert(pixel_list != (PixelList **) NULL);
3066 for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
3067 if (pixel_list[i] != (PixelList *) NULL)
3068 pixel_list[i]=DestroyPixelList(pixel_list[i]);
3069 pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list);
3070 return(pixel_list);
3071}
3072
3073static PixelList *AcquirePixelList(const size_t width,const size_t height)
3074{
3075 PixelList
3076 *pixel_list;
3077
3078 ssize_t
3079 i;
3080
3081 pixel_list=(PixelList *) AcquireMagickMemory(sizeof(*pixel_list));
3082 if (pixel_list == (PixelList *) NULL)
3083 return(pixel_list);
3084 (void) memset((void *) pixel_list,0,sizeof(*pixel_list));
3085 pixel_list->length=width*height;
3086 for (i=0; i < ListChannels; i++)
3087 {
3088 pixel_list->lists[i].nodes=(ListNode *) AcquireAlignedMemory(65537UL,
3089 sizeof(*pixel_list->lists[i].nodes));
3090 if (pixel_list->lists[i].nodes == (ListNode *) NULL)
3091 return(DestroyPixelList(pixel_list));
3092 (void) memset(pixel_list->lists[i].nodes,0,65537UL*
3093 sizeof(*pixel_list->lists[i].nodes));
3094 }
3095 pixel_list->signature=MagickCoreSignature;
3096 return(pixel_list);
3097}
3098
3099static PixelList **AcquirePixelListTLS(const size_t width,
3100 const size_t height)
3101{
3102 PixelList
3103 **pixel_list;
3104
3105 ssize_t
3106 i;
3107
3108 size_t
3109 number_threads;
3110
3111 number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
3112 pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
3113 sizeof(*pixel_list));
3114 if (pixel_list == (PixelList **) NULL)
3115 return((PixelList **) NULL);
3116 (void) memset(pixel_list,0,number_threads*sizeof(*pixel_list));
3117 for (i=0; i < (ssize_t) number_threads; i++)
3118 {
3119 pixel_list[i]=AcquirePixelList(width,height);
3120 if (pixel_list[i] == (PixelList *) NULL)
3121 return(DestroyPixelListTLS(pixel_list));
3122 }
3123 return(pixel_list);
3124}
3125
3126static void AddNodePixelList(PixelList *pixel_list,const ssize_t channel,
3127 const size_t color)
3128{
3129 SkipList
3130 *list;
3131
3132 ssize_t
3133 level;
3134
3135 size_t
3136 search,
3137 update[9];
3138
3139 /*
3140 Initialize the node.
3141 */
3142 list=pixel_list->lists+channel;
3143 list->nodes[color].signature=pixel_list->signature;
3144 list->nodes[color].count=1;
3145 /*
3146 Determine where it belongs in the list.
3147 */
3148 search=65536UL;
3149 (void) memset(update,0,sizeof(update));
3150 for (level=list->level; level >= 0; level--)
3151 {
3152 while (list->nodes[search].next[level] < color)
3153 search=list->nodes[search].next[level];
3154 update[level]=search;
3155 }
3156 /*
3157 Generate a pseudo-random level for this node.
3158 */
3159 for (level=0; ; level++)
3160 {
3161 pixel_list->seed=(pixel_list->seed*42893621L)+1L;
3162 if ((pixel_list->seed & 0x300) != 0x300)
3163 break;
3164 }
3165 if (level > 8)
3166 level=8;
3167 if (level > (list->level+2))
3168 level=list->level+2;
3169 /*
3170 If we're raising the list's level, link back to the root node.
3171 */
3172 while (level > list->level)
3173 {
3174 list->level++;
3175 update[list->level]=65536UL;
3176 }
3177 /*
3178 Link the node into the skip-list.
3179 */
3180 do
3181 {
3182 list->nodes[color].next[level]=list->nodes[update[level]].next[level];
3183 list->nodes[update[level]].next[level]=color;
3184 } while (level-- > 0);
3185}
3186
3187static void GetMaximumPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3188{
3189 SkipList
3190 *list;
3191
3192 ssize_t
3193 channel;
3194
3195 size_t
3196 color,
3197 maximum;
3198
3199 ssize_t
3200 count;
3201
3202 unsigned short
3203 channels[ListChannels];
3204
3205 /*
3206 Find the maximum value for each of the color.
3207 */
3208 for (channel=0; channel < 5; channel++)
3209 {
3210 list=pixel_list->lists+channel;
3211 color=65536L;
3212 count=0;
3213 maximum=list->nodes[color].next[0];
3214 do
3215 {
3216 color=list->nodes[color].next[0];
3217 if (color > maximum)
3218 maximum=color;
3219 count+=list->nodes[color].count;
3220 } while (count < (ssize_t) pixel_list->length);
3221 channels[channel]=(unsigned short) maximum;
3222 }
3223 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3224 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3225 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3226 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3227 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3228}
3229
3230static void GetMeanPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3231{
3232 MagickRealType
3233 sum;
3234
3235 SkipList
3236 *list;
3237
3238 ssize_t
3239 channel;
3240
3241 size_t
3242 color;
3243
3244 ssize_t
3245 count;
3246
3247 unsigned short
3248 channels[ListChannels];
3249
3250 /*
3251 Find the mean value for each of the color.
3252 */
3253 for (channel=0; channel < 5; channel++)
3254 {
3255 list=pixel_list->lists+channel;
3256 color=65536L;
3257 count=0;
3258 sum=0.0;
3259 do
3260 {
3261 color=list->nodes[color].next[0];
3262 sum+=(MagickRealType) list->nodes[color].count*color;
3263 count+=list->nodes[color].count;
3264 } while (count < (ssize_t) pixel_list->length);
3265 sum/=pixel_list->length;
3266 channels[channel]=(unsigned short) sum;
3267 }
3268 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3269 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3270 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3271 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3272 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3273}
3274
3275static void GetMedianPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3276{
3277 SkipList
3278 *list;
3279
3280 ssize_t
3281 channel;
3282
3283 size_t
3284 color;
3285
3286 ssize_t
3287 count;
3288
3289 unsigned short
3290 channels[ListChannels];
3291
3292 /*
3293 Find the median value for each of the color.
3294 */
3295 for (channel=0; channel < 5; channel++)
3296 {
3297 list=pixel_list->lists+channel;
3298 color=65536L;
3299 count=0;
3300 do
3301 {
3302 color=list->nodes[color].next[0];
3303 count+=list->nodes[color].count;
3304 } while (count <= (ssize_t) (pixel_list->length >> 1));
3305 channels[channel]=(unsigned short) color;
3306 }
3307 GetMagickPixelPacket((const Image *) NULL,pixel);
3308 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3309 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3310 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3311 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3312 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3313}
3314
3315static void GetMinimumPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3316{
3317 SkipList
3318 *list;
3319
3320 ssize_t
3321 channel;
3322
3323 size_t
3324 color,
3325 minimum;
3326
3327 ssize_t
3328 count;
3329
3330 unsigned short
3331 channels[ListChannels];
3332
3333 /*
3334 Find the minimum value for each of the color.
3335 */
3336 for (channel=0; channel < 5; channel++)
3337 {
3338 list=pixel_list->lists+channel;
3339 count=0;
3340 color=65536UL;
3341 minimum=list->nodes[color].next[0];
3342 do
3343 {
3344 color=list->nodes[color].next[0];
3345 if (color < minimum)
3346 minimum=color;
3347 count+=list->nodes[color].count;
3348 } while (count < (ssize_t) pixel_list->length);
3349 channels[channel]=(unsigned short) minimum;
3350 }
3351 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3352 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3353 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3354 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3355 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3356}
3357
3358static void GetModePixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3359{
3360 SkipList
3361 *list;
3362
3363 ssize_t
3364 channel;
3365
3366 size_t
3367 color,
3368 max_count,
3369 mode;
3370
3371 ssize_t
3372 count;
3373
3374 unsigned short
3375 channels[5];
3376
3377 /*
3378 Make each pixel the 'predominant color' of the specified neighborhood.
3379 */
3380 for (channel=0; channel < 5; channel++)
3381 {
3382 list=pixel_list->lists+channel;
3383 color=65536L;
3384 mode=color;
3385 max_count=list->nodes[mode].count;
3386 count=0;
3387 do
3388 {
3389 color=list->nodes[color].next[0];
3390 if (list->nodes[color].count > max_count)
3391 {
3392 mode=color;
3393 max_count=list->nodes[mode].count;
3394 }
3395 count+=list->nodes[color].count;
3396 } while (count < (ssize_t) pixel_list->length);
3397 channels[channel]=(unsigned short) mode;
3398 }
3399 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3400 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3401 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3402 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3403 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3404}
3405
3406static void GetNonpeakPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3407{
3408 SkipList
3409 *list;
3410
3411 ssize_t
3412 channel;
3413
3414 size_t
3415 color,
3416 next,
3417 previous;
3418
3419 ssize_t
3420 count;
3421
3422 unsigned short
3423 channels[5];
3424
3425 /*
3426 Finds the non peak value for each of the colors.
3427 */
3428 for (channel=0; channel < 5; channel++)
3429 {
3430 list=pixel_list->lists+channel;
3431 color=65536L;
3432 next=list->nodes[color].next[0];
3433 count=0;
3434 do
3435 {
3436 previous=color;
3437 color=next;
3438 next=list->nodes[color].next[0];
3439 count+=list->nodes[color].count;
3440 } while (count <= (ssize_t) (pixel_list->length >> 1));
3441 if ((previous == 65536UL) && (next != 65536UL))
3442 color=next;
3443 else
3444 if ((previous != 65536UL) && (next == 65536UL))
3445 color=previous;
3446 channels[channel]=(unsigned short) color;
3447 }
3448 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3449 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3450 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3451 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3452 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3453}
3454
3455static void GetRootMeanSquarePixelList(PixelList *pixel_list,
3456 MagickPixelPacket *pixel)
3457{
3458 MagickRealType
3459 sum;
3460
3461 SkipList
3462 *list;
3463
3464 ssize_t
3465 channel;
3466
3467 size_t
3468 color;
3469
3470 ssize_t
3471 count;
3472
3473 unsigned short
3474 channels[ListChannels];
3475
3476 /*
3477 Find the root mean square value for each of the color.
3478 */
3479 for (channel=0; channel < 5; channel++)
3480 {
3481 list=pixel_list->lists+channel;
3482 color=65536L;
3483 count=0;
3484 sum=0.0;
3485 do
3486 {
3487 color=list->nodes[color].next[0];
3488 sum+=(MagickRealType) (list->nodes[color].count*color*color);
3489 count+=list->nodes[color].count;
3490 } while (count < (ssize_t) pixel_list->length);
3491 sum/=pixel_list->length;
3492 channels[channel]=(unsigned short) sqrt(sum);
3493 }
3494 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3495 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3496 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3497 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3498 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3499}
3500
3501static void GetStandardDeviationPixelList(PixelList *pixel_list,
3502 MagickPixelPacket *pixel)
3503{
3504 MagickRealType
3505 sum,
3506 sum_squared;
3507
3508 SkipList
3509 *list;
3510
3511 size_t
3512 color;
3513
3514 ssize_t
3515 channel,
3516 count;
3517
3518 unsigned short
3519 channels[ListChannels];
3520
3521 /*
3522 Find the standard-deviation value for each of the color.
3523 */
3524 for (channel=0; channel < 5; channel++)
3525 {
3526 list=pixel_list->lists+channel;
3527 color=65536L;
3528 count=0;
3529 sum=0.0;
3530 sum_squared=0.0;
3531 do
3532 {
3533 ssize_t
3534 i;
3535
3536 color=list->nodes[color].next[0];
3537 sum+=(MagickRealType) list->nodes[color].count*color;
3538 for (i=0; i < (ssize_t) list->nodes[color].count; i++)
3539 sum_squared+=((MagickRealType) color)*((MagickRealType) color);
3540 count+=list->nodes[color].count;
3541 } while (count < (ssize_t) pixel_list->length);
3542 sum/=pixel_list->length;
3543 sum_squared/=pixel_list->length;
3544 channels[channel]=(unsigned short) sqrt(sum_squared-(sum*sum));
3545 }
3546 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3547 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3548 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3549 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3550 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3551}
3552
3553static inline void InsertPixelList(const Image *image,const PixelPacket *pixel,
3554 const IndexPacket *indexes,PixelList *pixel_list)
3555{
3556 size_t
3557 signature;
3558
3559 unsigned short
3560 index;
3561
3562 index=ScaleQuantumToShort(GetPixelRed(pixel));
3563 signature=pixel_list->lists[0].nodes[index].signature;
3564 if (signature == pixel_list->signature)
3565 pixel_list->lists[0].nodes[index].count++;
3566 else
3567 AddNodePixelList(pixel_list,0,index);
3568 index=ScaleQuantumToShort(GetPixelGreen(pixel));
3569 signature=pixel_list->lists[1].nodes[index].signature;
3570 if (signature == pixel_list->signature)
3571 pixel_list->lists[1].nodes[index].count++;
3572 else
3573 AddNodePixelList(pixel_list,1,index);
3574 index=ScaleQuantumToShort(GetPixelBlue(pixel));
3575 signature=pixel_list->lists[2].nodes[index].signature;
3576 if (signature == pixel_list->signature)
3577 pixel_list->lists[2].nodes[index].count++;
3578 else
3579 AddNodePixelList(pixel_list,2,index);
3580 index=ScaleQuantumToShort(GetPixelOpacity(pixel));
3581 signature=pixel_list->lists[3].nodes[index].signature;
3582 if (signature == pixel_list->signature)
3583 pixel_list->lists[3].nodes[index].count++;
3584 else
3585 AddNodePixelList(pixel_list,3,index);
3586 if (image->colorspace == CMYKColorspace)
3587 index=ScaleQuantumToShort(GetPixelIndex(indexes));
3588 signature=pixel_list->lists[4].nodes[index].signature;
3589 if (signature == pixel_list->signature)
3590 pixel_list->lists[4].nodes[index].count++;
3591 else
3592 AddNodePixelList(pixel_list,4,index);
3593}
3594
3595static void ResetPixelList(PixelList *pixel_list)
3596{
3597 int
3598 level;
3599
3600 ListNode
3601 *root;
3602
3603 SkipList
3604 *list;
3605
3606 ssize_t
3607 channel;
3608
3609 /*
3610 Reset the skip-list.
3611 */
3612 for (channel=0; channel < 5; channel++)
3613 {
3614 list=pixel_list->lists+channel;
3615 root=list->nodes+65536UL;
3616 list->level=0;
3617 for (level=0; level < 9; level++)
3618 root->next[level]=65536UL;
3619 }
3620 pixel_list->seed=pixel_list->signature++;
3621}
3622
3623MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
3624 const size_t width,const size_t height,ExceptionInfo *exception)
3625{
3626 Image
3627 *statistic_image;
3628
3629 statistic_image=StatisticImageChannel(image,DefaultChannels,type,width,
3630 height,exception);
3631 return(statistic_image);
3632}
3633
3634MagickExport Image *StatisticImageChannel(const Image *image,
3635 const ChannelType channel,const StatisticType type,const size_t width,
3636 const size_t height,ExceptionInfo *exception)
3637{
3638#define StatisticImageTag "Statistic/Image"
3639
3640 CacheView
3641 *image_view,
3642 *statistic_view;
3643
3644 Image
3645 *statistic_image;
3646
3647 MagickBooleanType
3648 status;
3649
3650 MagickOffsetType
3651 progress;
3652
3653 PixelList
3654 **magick_restrict pixel_list;
3655
3656 size_t
3657 neighbor_height,
3658 neighbor_width;
3659
3660 ssize_t
3661 y;
3662
3663 /*
3664 Initialize statistics image attributes.
3665 */
3666 assert(image != (Image *) NULL);
3667 assert(image->signature == MagickCoreSignature);
3668 if (IsEventLogging() != MagickFalse)
3669 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3670 assert(exception != (ExceptionInfo *) NULL);
3671 assert(exception->signature == MagickCoreSignature);
3672 statistic_image=CloneImage(image,0,0,MagickTrue,exception);
3673 if (statistic_image == (Image *) NULL)
3674 return((Image *) NULL);
3675 if (SetImageStorageClass(statistic_image,DirectClass) == MagickFalse)
3676 {
3677 InheritException(exception,&statistic_image->exception);
3678 statistic_image=DestroyImage(statistic_image);
3679 return((Image *) NULL);
3680 }
3681 neighbor_width=width == 0 ? GetOptimalKernelWidth2D((double) width,0.5) :
3682 width;
3683 neighbor_height=height == 0 ? GetOptimalKernelWidth2D((double) height,0.5) :
3684 height;
3685 pixel_list=AcquirePixelListTLS(neighbor_width,neighbor_height);
3686 if (pixel_list == (PixelList **) NULL)
3687 {
3688 statistic_image=DestroyImage(statistic_image);
3689 ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3690 }
3691 /*
3692 Make each pixel the min / max / median / mode / etc. of the neighborhood.
3693 */
3694 status=MagickTrue;
3695 progress=0;
3696 image_view=AcquireVirtualCacheView(image,exception);
3697 statistic_view=AcquireAuthenticCacheView(statistic_image,exception);
3698#if defined(MAGICKCORE_OPENMP_SUPPORT)
3699 #pragma omp parallel for schedule(static) shared(progress,status) \
3700 magick_number_threads(image,statistic_image,statistic_image->rows,1)
3701#endif
3702 for (y=0; y < (ssize_t) statistic_image->rows; y++)
3703 {
3704 const int
3705 id = GetOpenMPThreadId();
3706
3707 const IndexPacket
3708 *magick_restrict indexes;
3709
3710 const PixelPacket
3711 *magick_restrict p;
3712
3713 IndexPacket
3714 *magick_restrict statistic_indexes;
3715
3716 PixelPacket
3717 *magick_restrict q;
3718
3719 ssize_t
3720 x;
3721
3722 if (status == MagickFalse)
3723 continue;
3724 p=GetCacheViewVirtualPixels(image_view,-((ssize_t) neighbor_width/2L),y-
3725 (ssize_t) (neighbor_height/2L),image->columns+neighbor_width,
3726 neighbor_height,exception);
3727 q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns, 1,exception);
3728 if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
3729 {
3730 status=MagickFalse;
3731 continue;
3732 }
3733 indexes=GetCacheViewVirtualIndexQueue(image_view);
3734 statistic_indexes=GetCacheViewAuthenticIndexQueue(statistic_view);
3735 for (x=0; x < (ssize_t) statistic_image->columns; x++)
3736 {
3737 MagickPixelPacket
3738 pixel;
3739
3740 const IndexPacket
3741 *magick_restrict s;
3742
3743 const PixelPacket
3744 *magick_restrict r;
3745
3746 ssize_t
3747 u,
3748 v;
3749
3750 r=p;
3751 s=indexes+x;
3752 ResetPixelList(pixel_list[id]);
3753 for (v=0; v < (ssize_t) neighbor_height; v++)
3754 {
3755 for (u=0; u < (ssize_t) neighbor_width; u++)
3756 InsertPixelList(image,r+u,s+u,pixel_list[id]);
3757 r+=(ptrdiff_t) image->columns+neighbor_width;
3758 s+=(ptrdiff_t) image->columns+neighbor_width;
3759 }
3760 GetMagickPixelPacket(image,&pixel);
3761 SetMagickPixelPacket(image,p+neighbor_width*neighbor_height/2,indexes+x+
3762 neighbor_width*neighbor_height/2,&pixel);
3763 switch (type)
3764 {
3765 case GradientStatistic:
3766 {
3767 MagickPixelPacket
3768 maximum,
3769 minimum;
3770
3771 GetMinimumPixelList(pixel_list[id],&pixel);
3772 minimum=pixel;
3773 GetMaximumPixelList(pixel_list[id],&pixel);
3774 maximum=pixel;
3775 pixel.red=MagickAbsoluteValue(maximum.red-minimum.red);
3776 pixel.green=MagickAbsoluteValue(maximum.green-minimum.green);
3777 pixel.blue=MagickAbsoluteValue(maximum.blue-minimum.blue);
3778 pixel.opacity=MagickAbsoluteValue(maximum.opacity-minimum.opacity);
3779 if (image->colorspace == CMYKColorspace)
3780 pixel.index=MagickAbsoluteValue(maximum.index-minimum.index);
3781 break;
3782 }
3783 case MaximumStatistic:
3784 {
3785 GetMaximumPixelList(pixel_list[id],&pixel);
3786 break;
3787 }
3788 case MeanStatistic:
3789 {
3790 GetMeanPixelList(pixel_list[id],&pixel);
3791 break;
3792 }
3793 case MedianStatistic:
3794 default:
3795 {
3796 GetMedianPixelList(pixel_list[id],&pixel);
3797 break;
3798 }
3799 case MinimumStatistic:
3800 {
3801 GetMinimumPixelList(pixel_list[id],&pixel);
3802 break;
3803 }
3804 case ModeStatistic:
3805 {
3806 GetModePixelList(pixel_list[id],&pixel);
3807 break;
3808 }
3809 case NonpeakStatistic:
3810 {
3811 GetNonpeakPixelList(pixel_list[id],&pixel);
3812 break;
3813 }
3814 case RootMeanSquareStatistic:
3815 {
3816 GetRootMeanSquarePixelList(pixel_list[id],&pixel);
3817 break;
3818 }
3819 case StandardDeviationStatistic:
3820 {
3821 GetStandardDeviationPixelList(pixel_list[id],&pixel);
3822 break;
3823 }
3824 }
3825 if ((channel & RedChannel) != 0)
3826 SetPixelRed(q,ClampToQuantum(pixel.red));
3827 if ((channel & GreenChannel) != 0)
3828 SetPixelGreen(q,ClampToQuantum(pixel.green));
3829 if ((channel & BlueChannel) != 0)
3830 SetPixelBlue(q,ClampToQuantum(pixel.blue));
3831 if ((channel & OpacityChannel) != 0)
3832 SetPixelOpacity(q,ClampToQuantum(pixel.opacity));
3833 if (((channel & IndexChannel) != 0) &&
3834 (image->colorspace == CMYKColorspace))
3835 SetPixelIndex(statistic_indexes+x,ClampToQuantum(pixel.index));
3836 p++;
3837 q++;
3838 }
3839 if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
3840 status=MagickFalse;
3841 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3842 {
3843 MagickBooleanType
3844 proceed;
3845
3846 proceed=SetImageProgress(image,StatisticImageTag,progress++,
3847 image->rows);
3848 if (proceed == MagickFalse)
3849 status=MagickFalse;
3850 }
3851 }
3852 statistic_view=DestroyCacheView(statistic_view);
3853 image_view=DestroyCacheView(image_view);
3854 pixel_list=DestroyPixelListTLS(pixel_list);
3855 if (status == MagickFalse)
3856 statistic_image=DestroyImage(statistic_image);
3857 return(statistic_image);
3858}
_CacheView
Definition cache-view.c:66
_ChannelMoments
Definition statistic.h:53
_ChannelPerceptualHash
Definition statistic.h:68
_ChannelStatistics
Definition statistic.h:31
_ExceptionInfo
Definition exception.h:103
_Image
Definition image.h:134
_ListNode
Definition statistic.c:3018
_MagickPixelPacket
Definition pixel.h:105
_PixelList
Definition statistic.c:3035
_PixelPacket
Definition pixel.h:132
_PointInfo
Definition draw.h:140
_RandomInfo
Definition random.c:82
_SkipList
Definition statistic.c:3026