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