MagickCore 6.9.13
Loading...
Searching...
No Matches
statistic.c
1/*
2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3% %
4% %
5% %
6% SSSSS TTTTT AAA TTTTT IIIII SSSSS TTTTT IIIII CCCC %
7% SS T A A T I SS T I C %
8% SSS T AAAAA T I SSS T I C %
9% SS T A A T I SS T I C %
10% SSSSS T A A T IIIII SSSSS T IIIII CCCC %
11% %
12% %
13% MagickCore Image Statistical Methods %
14% %
15% Software Design %
16% Cristy %
17% July 1992 %
18% %
19% %
20% Copyright 1999 ImageMagick Studio LLC, a non-profit organization %
21% dedicated to making software imaging solutions freely available. %
22% %
23% You may not use this file except in compliance with the License. You may %
24% obtain a copy of the License at %
25% %
26% https://imagemagick.org/script/license.php %
27% %
28% Unless required by applicable law or agreed to in writing, software %
29% distributed under the License is distributed on an "AS IS" BASIS, %
30% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. %
31% See the License for the specific language governing permissions and %
32% limitations under the License. %
33% %
34%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35%
36%
37%
38*/
39␌
40/*
41 Include declarations.
42*/
43#include "magick/studio.h"
44#include "magick/accelerate-private.h"
45#include "magick/animate.h"
46#include "magick/animate.h"
47#include "magick/blob.h"
48#include "magick/blob-private.h"
49#include "magick/cache.h"
50#include "magick/cache-private.h"
51#include "magick/cache-view.h"
52#include "magick/client.h"
53#include "magick/color.h"
54#include "magick/color-private.h"
55#include "magick/colorspace.h"
56#include "magick/colorspace-private.h"
57#include "magick/composite.h"
58#include "magick/composite-private.h"
59#include "magick/compress.h"
60#include "magick/constitute.h"
61#include "magick/deprecate.h"
62#include "magick/display.h"
63#include "magick/draw.h"
64#include "magick/enhance.h"
65#include "magick/exception.h"
66#include "magick/exception-private.h"
67#include "magick/gem.h"
68#include "magick/geometry.h"
69#include "magick/list.h"
70#include "magick/image-private.h"
71#include "magick/magic.h"
72#include "magick/magick.h"
73#include "magick/memory_.h"
74#include "magick/module.h"
75#include "magick/monitor.h"
76#include "magick/monitor-private.h"
77#include "magick/option.h"
78#include "magick/paint.h"
79#include "magick/pixel-private.h"
80#include "magick/profile.h"
81#include "magick/property.h"
82#include "magick/quantize.h"
83#include "magick/random_.h"
84#include "magick/random-private.h"
85#include "magick/resource_.h"
86#include "magick/segment.h"
87#include "magick/semaphore.h"
88#include "magick/signature-private.h"
89#include "magick/statistic.h"
90#include "magick/statistic-private.h"
91#include "magick/string_.h"
92#include "magick/thread-private.h"
93#include "magick/timer.h"
94#include "magick/utility.h"
95#include "magick/version.h"
96␌
97/*
98%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
99% %
100% %
101% %
102% E v a l u a t e I m a g e %
103% %
104% %
105% %
106%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
107%
108% EvaluateImage() applies a value to the image with an arithmetic, relational,
109% or logical operator to an image. Use these operations to lighten or darken
110% an image, to increase or decrease contrast in an image, or to produce the
111% "negative" of an image.
112%
113% The format of the EvaluateImageChannel method is:
114%
115% MagickBooleanType EvaluateImage(Image *image,
116% const MagickEvaluateOperator op,const double value,
117% ExceptionInfo *exception)
118% MagickBooleanType EvaluateImages(Image *images,
119% const MagickEvaluateOperator op,const double value,
120% ExceptionInfo *exception)
121% MagickBooleanType EvaluateImageChannel(Image *image,
122% const ChannelType channel,const MagickEvaluateOperator op,
123% const double value,ExceptionInfo *exception)
124%
125% A description of each parameter follows:
126%
127% o image: the image.
128%
129% o channel: the channel.
130%
131% o op: A channel op.
132%
133% o value: A value value.
134%
135% o exception: return any errors or warnings in this structure.
136%
137*/
138
139static MagickPixelPacket **DestroyPixelTLS(const Image *images,
140 MagickPixelPacket **pixels)
141{
142 ssize_t
143 i;
144
145 size_t
146 rows;
147
148 assert(pixels != (MagickPixelPacket **) NULL);
149 rows=MagickMax(GetImageListLength(images),
150 (size_t) GetMagickResourceLimit(ThreadResource));
151 for (i=0; i < (ssize_t) rows; i++)
152 if (pixels[i] != (MagickPixelPacket *) NULL)
153 pixels[i]=(MagickPixelPacket *) RelinquishMagickMemory(pixels[i]);
154 pixels=(MagickPixelPacket **) RelinquishMagickMemory(pixels);
155 return(pixels);
156}
157
158static MagickPixelPacket **AcquirePixelTLS(const Image *images)
159{
160 const Image
161 *next;
162
163 MagickPixelPacket
164 **pixels;
165
166 ssize_t
167 i,
168 j;
169
170 size_t
171 columns,
172 rows;
173
174 rows=MagickMax(GetImageListLength(images),
175 (size_t) GetMagickResourceLimit(ThreadResource));
176 pixels=(MagickPixelPacket **) AcquireQuantumMemory(rows,sizeof(*pixels));
177 if (pixels == (MagickPixelPacket **) NULL)
178 return((MagickPixelPacket **) NULL);
179 (void) memset(pixels,0,rows*sizeof(*pixels));
180 columns=GetImageListLength(images);
181 for (next=images; next != (Image *) NULL; next=next->next)
182 columns=MagickMax(next->columns,columns);
183 for (i=0; i < (ssize_t) rows; i++)
184 {
185 pixels[i]=(MagickPixelPacket *) AcquireQuantumMemory(columns,
186 sizeof(**pixels));
187 if (pixels[i] == (MagickPixelPacket *) NULL)
188 return(DestroyPixelTLS(images,pixels));
189 for (j=0; j < (ssize_t) columns; j++)
190 GetMagickPixelPacket(images,&pixels[i][j]);
191 }
192 return(pixels);
193}
194
195static inline double EvaluateMax(const double x,const double y)
196{
197 if (x > y)
198 return(x);
199 return(y);
200}
201
202#if defined(__cplusplus) || defined(c_plusplus)
203extern "C" {
204#endif
205
206static int IntensityCompare(const void *x,const void *y)
207{
208 const MagickPixelPacket
209 *color_1,
210 *color_2;
211
212 int
213 intensity;
214
215 color_1=(const MagickPixelPacket *) x;
216 color_2=(const MagickPixelPacket *) y;
217 intensity=(int) MagickPixelIntensity(color_2)-(int)
218 MagickPixelIntensity(color_1);
219 return(intensity);
220}
221
222#if defined(__cplusplus) || defined(c_plusplus)
223}
224#endif
225
226static MagickRealType ApplyEvaluateOperator(RandomInfo *random_info,
227 const Quantum pixel,const MagickEvaluateOperator op,
228 const MagickRealType value)
229{
230 MagickRealType
231 result;
232
233 ssize_t
234 i;
235
236 result=0.0;
237 switch (op)
238 {
239 case UndefinedEvaluateOperator:
240 break;
241 case AbsEvaluateOperator:
242 {
243 result=(MagickRealType) fabs((double) pixel+value);
244 break;
245 }
246 case AddEvaluateOperator:
247 {
248 result=(MagickRealType) pixel+value;
249 break;
250 }
251 case AddModulusEvaluateOperator:
252 {
253 /*
254 This returns a 'floored modulus' of the addition which is a
255 positive result. It differs from % or fmod() which returns a
256 'truncated modulus' result, where floor() is replaced by trunc()
257 and could return a negative result (which is clipped).
258 */
259 result=(MagickRealType) pixel+value;
260 result-=((MagickRealType) QuantumRange+1.0)*floor((double) result/
261 ((MagickRealType) QuantumRange+1.0));
262 break;
263 }
264 case AndEvaluateOperator:
265 {
266 result=(MagickRealType) ((ssize_t) pixel & (ssize_t) (value+0.5));
267 break;
268 }
269 case CosineEvaluateOperator:
270 {
271 result=(MagickRealType) QuantumRange*(0.5*cos((double) (2.0*MagickPI*
272 QuantumScale*(MagickRealType) pixel*value))+0.5);
273 break;
274 }
275 case DivideEvaluateOperator:
276 {
277 result=(MagickRealType) pixel/(value == 0.0 ? 1.0 : value);
278 break;
279 }
280 case ExponentialEvaluateOperator:
281 {
282 result=(MagickRealType) QuantumRange*exp(value*QuantumScale*(double)
283 pixel);
284 break;
285 }
286 case GaussianNoiseEvaluateOperator:
287 {
288 result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
289 GaussianNoise,value);
290 break;
291 }
292 case ImpulseNoiseEvaluateOperator:
293 {
294 result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel,
295 ImpulseNoise,value);
296 break;
297 }
298 case InverseLogEvaluateOperator:
299 {
300 result=((MagickRealType) QuantumRange*pow((value+1.0),
301 QuantumScale*(MagickRealType) pixel)-1.0)*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(log10((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(log10((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(log10((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(log10((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(log10((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,
3116 const size_t height)
3117{
3118 PixelList
3119 **pixel_list;
3120
3121 ssize_t
3122 i;
3123
3124 size_t
3125 number_threads;
3126
3127 number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
3128 pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
3129 sizeof(*pixel_list));
3130 if (pixel_list == (PixelList **) NULL)
3131 return((PixelList **) NULL);
3132 (void) memset(pixel_list,0,number_threads*sizeof(*pixel_list));
3133 for (i=0; i < (ssize_t) number_threads; i++)
3134 {
3135 pixel_list[i]=AcquirePixelList(width,height);
3136 if (pixel_list[i] == (PixelList *) NULL)
3137 return(DestroyPixelListTLS(pixel_list));
3138 }
3139 return(pixel_list);
3140}
3141
3142static void AddNodePixelList(PixelList *pixel_list,const ssize_t channel,
3143 const size_t color)
3144{
3145 SkipList
3146 *list;
3147
3148 ssize_t
3149 level;
3150
3151 size_t
3152 search,
3153 update[9];
3154
3155 /*
3156 Initialize the node.
3157 */
3158 list=pixel_list->lists+channel;
3159 list->nodes[color].signature=pixel_list->signature;
3160 list->nodes[color].count=1;
3161 /*
3162 Determine where it belongs in the list.
3163 */
3164 search=65536UL;
3165 (void) memset(update,0,sizeof(update));
3166 for (level=list->level; level >= 0; level--)
3167 {
3168 while (list->nodes[search].next[level] < color)
3169 search=list->nodes[search].next[level];
3170 update[level]=search;
3171 }
3172 /*
3173 Generate a pseudo-random level for this node.
3174 */
3175 for (level=0; ; level++)
3176 {
3177 pixel_list->seed=(pixel_list->seed*42893621L)+1L;
3178 if ((pixel_list->seed & 0x300) != 0x300)
3179 break;
3180 }
3181 if (level > 8)
3182 level=8;
3183 if (level > (list->level+2))
3184 level=list->level+2;
3185 /*
3186 If we're raising the list's level, link back to the root node.
3187 */
3188 while (level > list->level)
3189 {
3190 list->level++;
3191 update[list->level]=65536UL;
3192 }
3193 /*
3194 Link the node into the skip-list.
3195 */
3196 do
3197 {
3198 list->nodes[color].next[level]=list->nodes[update[level]].next[level];
3199 list->nodes[update[level]].next[level]=color;
3200 } while (level-- > 0);
3201}
3202
3203static void GetMaximumPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3204{
3205 SkipList
3206 *list;
3207
3208 ssize_t
3209 channel;
3210
3211 size_t
3212 color,
3213 maximum;
3214
3215 ssize_t
3216 count;
3217
3218 unsigned short
3219 channels[ListChannels];
3220
3221 /*
3222 Find the maximum value for each of the color.
3223 */
3224 for (channel=0; channel < 5; channel++)
3225 {
3226 list=pixel_list->lists+channel;
3227 color=65536L;
3228 count=0;
3229 maximum=list->nodes[color].next[0];
3230 do
3231 {
3232 color=list->nodes[color].next[0];
3233 if (color > maximum)
3234 maximum=color;
3235 count+=list->nodes[color].count;
3236 } while (count < (ssize_t) pixel_list->length);
3237 channels[channel]=(unsigned short) maximum;
3238 }
3239 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3240 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3241 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3242 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3243 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3244}
3245
3246static void GetMeanPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3247{
3248 MagickRealType
3249 sum;
3250
3251 SkipList
3252 *list;
3253
3254 ssize_t
3255 channel;
3256
3257 size_t
3258 color;
3259
3260 ssize_t
3261 count;
3262
3263 unsigned short
3264 channels[ListChannels];
3265
3266 /*
3267 Find the mean value for each of the color.
3268 */
3269 for (channel=0; channel < 5; channel++)
3270 {
3271 list=pixel_list->lists+channel;
3272 color=65536L;
3273 count=0;
3274 sum=0.0;
3275 do
3276 {
3277 color=list->nodes[color].next[0];
3278 sum+=(MagickRealType) list->nodes[color].count*color;
3279 count+=list->nodes[color].count;
3280 } while (count < (ssize_t) pixel_list->length);
3281 sum/=pixel_list->length;
3282 channels[channel]=(unsigned short) sum;
3283 }
3284 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3285 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3286 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3287 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3288 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3289}
3290
3291static void GetMedianPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3292{
3293 SkipList
3294 *list;
3295
3296 ssize_t
3297 channel;
3298
3299 size_t
3300 color;
3301
3302 ssize_t
3303 count;
3304
3305 unsigned short
3306 channels[ListChannels];
3307
3308 /*
3309 Find the median value for each of the color.
3310 */
3311 for (channel=0; channel < 5; channel++)
3312 {
3313 list=pixel_list->lists+channel;
3314 color=65536L;
3315 count=0;
3316 do
3317 {
3318 color=list->nodes[color].next[0];
3319 count+=list->nodes[color].count;
3320 } while (count <= (ssize_t) (pixel_list->length >> 1));
3321 channels[channel]=(unsigned short) color;
3322 }
3323 GetMagickPixelPacket((const Image *) NULL,pixel);
3324 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3325 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3326 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3327 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3328 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3329}
3330
3331static void GetMinimumPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3332{
3333 SkipList
3334 *list;
3335
3336 ssize_t
3337 channel;
3338
3339 size_t
3340 color,
3341 minimum;
3342
3343 ssize_t
3344 count;
3345
3346 unsigned short
3347 channels[ListChannels];
3348
3349 /*
3350 Find the minimum value for each of the color.
3351 */
3352 for (channel=0; channel < 5; channel++)
3353 {
3354 list=pixel_list->lists+channel;
3355 count=0;
3356 color=65536UL;
3357 minimum=list->nodes[color].next[0];
3358 do
3359 {
3360 color=list->nodes[color].next[0];
3361 if (color < minimum)
3362 minimum=color;
3363 count+=list->nodes[color].count;
3364 } while (count < (ssize_t) pixel_list->length);
3365 channels[channel]=(unsigned short) minimum;
3366 }
3367 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3368 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3369 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3370 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3371 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3372}
3373
3374static void GetModePixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3375{
3376 SkipList
3377 *list;
3378
3379 ssize_t
3380 channel;
3381
3382 size_t
3383 color,
3384 max_count,
3385 mode;
3386
3387 ssize_t
3388 count;
3389
3390 unsigned short
3391 channels[5];
3392
3393 /*
3394 Make each pixel the 'predominant color' of the specified neighborhood.
3395 */
3396 for (channel=0; channel < 5; channel++)
3397 {
3398 list=pixel_list->lists+channel;
3399 color=65536L;
3400 mode=color;
3401 max_count=list->nodes[mode].count;
3402 count=0;
3403 do
3404 {
3405 color=list->nodes[color].next[0];
3406 if (list->nodes[color].count > max_count)
3407 {
3408 mode=color;
3409 max_count=list->nodes[mode].count;
3410 }
3411 count+=list->nodes[color].count;
3412 } while (count < (ssize_t) pixel_list->length);
3413 channels[channel]=(unsigned short) mode;
3414 }
3415 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3416 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3417 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3418 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3419 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3420}
3421
3422static void GetNonpeakPixelList(PixelList *pixel_list,MagickPixelPacket *pixel)
3423{
3424 SkipList
3425 *list;
3426
3427 ssize_t
3428 channel;
3429
3430 size_t
3431 color,
3432 next,
3433 previous;
3434
3435 ssize_t
3436 count;
3437
3438 unsigned short
3439 channels[5];
3440
3441 /*
3442 Finds the non peak value for each of the colors.
3443 */
3444 for (channel=0; channel < 5; channel++)
3445 {
3446 list=pixel_list->lists+channel;
3447 color=65536L;
3448 next=list->nodes[color].next[0];
3449 count=0;
3450 do
3451 {
3452 previous=color;
3453 color=next;
3454 next=list->nodes[color].next[0];
3455 count+=list->nodes[color].count;
3456 } while (count <= (ssize_t) (pixel_list->length >> 1));
3457 if ((previous == 65536UL) && (next != 65536UL))
3458 color=next;
3459 else
3460 if ((previous != 65536UL) && (next == 65536UL))
3461 color=previous;
3462 channels[channel]=(unsigned short) color;
3463 }
3464 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3465 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3466 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3467 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3468 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3469}
3470
3471static void GetRootMeanSquarePixelList(PixelList *pixel_list,
3472 MagickPixelPacket *pixel)
3473{
3474 MagickRealType
3475 sum;
3476
3477 SkipList
3478 *list;
3479
3480 ssize_t
3481 channel;
3482
3483 size_t
3484 color;
3485
3486 ssize_t
3487 count;
3488
3489 unsigned short
3490 channels[ListChannels];
3491
3492 /*
3493 Find the root mean square value for each of the color.
3494 */
3495 for (channel=0; channel < 5; channel++)
3496 {
3497 list=pixel_list->lists+channel;
3498 color=65536L;
3499 count=0;
3500 sum=0.0;
3501 do
3502 {
3503 color=list->nodes[color].next[0];
3504 sum+=(MagickRealType) (list->nodes[color].count*color*color);
3505 count+=list->nodes[color].count;
3506 } while (count < (ssize_t) pixel_list->length);
3507 sum/=pixel_list->length;
3508 channels[channel]=(unsigned short) sqrt(sum);
3509 }
3510 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3511 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3512 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3513 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3514 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3515}
3516
3517static void GetStandardDeviationPixelList(PixelList *pixel_list,
3518 MagickPixelPacket *pixel)
3519{
3520 MagickRealType
3521 sum,
3522 sum_squared;
3523
3524 SkipList
3525 *list;
3526
3527 size_t
3528 color;
3529
3530 ssize_t
3531 channel,
3532 count;
3533
3534 unsigned short
3535 channels[ListChannels];
3536
3537 /*
3538 Find the standard-deviation value for each of the color.
3539 */
3540 for (channel=0; channel < 5; channel++)
3541 {
3542 list=pixel_list->lists+channel;
3543 color=65536L;
3544 count=0;
3545 sum=0.0;
3546 sum_squared=0.0;
3547 do
3548 {
3549 ssize_t
3550 i;
3551
3552 color=list->nodes[color].next[0];
3553 sum+=(MagickRealType) list->nodes[color].count*color;
3554 for (i=0; i < (ssize_t) list->nodes[color].count; i++)
3555 sum_squared+=((MagickRealType) color)*((MagickRealType) color);
3556 count+=list->nodes[color].count;
3557 } while (count < (ssize_t) pixel_list->length);
3558 sum/=pixel_list->length;
3559 sum_squared/=pixel_list->length;
3560 channels[channel]=(unsigned short) sqrt(sum_squared-(sum*sum));
3561 }
3562 pixel->red=(MagickRealType) ScaleShortToQuantum(channels[0]);
3563 pixel->green=(MagickRealType) ScaleShortToQuantum(channels[1]);
3564 pixel->blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
3565 pixel->opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
3566 pixel->index=(MagickRealType) ScaleShortToQuantum(channels[4]);
3567}
3568
3569static inline void InsertPixelList(const Image *image,const PixelPacket *pixel,
3570 const IndexPacket *indexes,PixelList *pixel_list)
3571{
3572 size_t
3573 signature;
3574
3575 unsigned short
3576 index;
3577
3578 index=ScaleQuantumToShort(GetPixelRed(pixel));
3579 signature=pixel_list->lists[0].nodes[index].signature;
3580 if (signature == pixel_list->signature)
3581 pixel_list->lists[0].nodes[index].count++;
3582 else
3583 AddNodePixelList(pixel_list,0,index);
3584 index=ScaleQuantumToShort(GetPixelGreen(pixel));
3585 signature=pixel_list->lists[1].nodes[index].signature;
3586 if (signature == pixel_list->signature)
3587 pixel_list->lists[1].nodes[index].count++;
3588 else
3589 AddNodePixelList(pixel_list,1,index);
3590 index=ScaleQuantumToShort(GetPixelBlue(pixel));
3591 signature=pixel_list->lists[2].nodes[index].signature;
3592 if (signature == pixel_list->signature)
3593 pixel_list->lists[2].nodes[index].count++;
3594 else
3595 AddNodePixelList(pixel_list,2,index);
3596 index=ScaleQuantumToShort(GetPixelOpacity(pixel));
3597 signature=pixel_list->lists[3].nodes[index].signature;
3598 if (signature == pixel_list->signature)
3599 pixel_list->lists[3].nodes[index].count++;
3600 else
3601 AddNodePixelList(pixel_list,3,index);
3602 if (image->colorspace == CMYKColorspace)
3603 index=ScaleQuantumToShort(GetPixelIndex(indexes));
3604 signature=pixel_list->lists[4].nodes[index].signature;
3605 if (signature == pixel_list->signature)
3606 pixel_list->lists[4].nodes[index].count++;
3607 else
3608 AddNodePixelList(pixel_list,4,index);
3609}
3610
3611static void ResetPixelList(PixelList *pixel_list)
3612{
3613 int
3614 level;
3615
3616 ListNode
3617 *root;
3618
3619 SkipList
3620 *list;
3621
3622 ssize_t
3623 channel;
3624
3625 /*
3626 Reset the skip-list.
3627 */
3628 for (channel=0; channel < 5; channel++)
3629 {
3630 list=pixel_list->lists+channel;
3631 root=list->nodes+65536UL;
3632 list->level=0;
3633 for (level=0; level < 9; level++)
3634 root->next[level]=65536UL;
3635 }
3636 pixel_list->seed=pixel_list->signature++;
3637}
3638
3639MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
3640 const size_t width,const size_t height,ExceptionInfo *exception)
3641{
3642 Image
3643 *statistic_image;
3644
3645 statistic_image=StatisticImageChannel(image,DefaultChannels,type,width,
3646 height,exception);
3647 return(statistic_image);
3648}
3649
3650MagickExport Image *StatisticImageChannel(const Image *image,
3651 const ChannelType channel,const StatisticType type,const size_t width,
3652 const size_t height,ExceptionInfo *exception)
3653{
3654#define StatisticImageTag "Statistic/Image"
3655
3656 CacheView
3657 *image_view,
3658 *statistic_view;
3659
3660 Image
3661 *statistic_image;
3662
3663 MagickBooleanType
3664 status;
3665
3666 MagickOffsetType
3667 progress;
3668
3669 PixelList
3670 **magick_restrict pixel_list;
3671
3672 size_t
3673 neighbor_height,
3674 neighbor_width;
3675
3676 ssize_t
3677 y;
3678
3679 /*
3680 Initialize statistics image attributes.
3681 */
3682 assert(image != (Image *) NULL);
3683 assert(image->signature == MagickCoreSignature);
3684 if (IsEventLogging() != MagickFalse)
3685 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3686 assert(exception != (ExceptionInfo *) NULL);
3687 assert(exception->signature == MagickCoreSignature);
3688 statistic_image=CloneImage(image,0,0,MagickTrue,exception);
3689 if (statistic_image == (Image *) NULL)
3690 return((Image *) NULL);
3691 if (SetImageStorageClass(statistic_image,DirectClass) == MagickFalse)
3692 {
3693 InheritException(exception,&statistic_image->exception);
3694 statistic_image=DestroyImage(statistic_image);
3695 return((Image *) NULL);
3696 }
3697 neighbor_width=width == 0 ? GetOptimalKernelWidth2D((double) width,0.5) :
3698 width;
3699 neighbor_height=height == 0 ? GetOptimalKernelWidth2D((double) height,0.5) :
3700 height;
3701 pixel_list=AcquirePixelListTLS(neighbor_width,neighbor_height);
3702 if (pixel_list == (PixelList **) NULL)
3703 {
3704 statistic_image=DestroyImage(statistic_image);
3705 ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3706 }
3707 /*
3708 Make each pixel the min / max / median / mode / etc. of the neighborhood.
3709 */
3710 status=MagickTrue;
3711 progress=0;
3712 image_view=AcquireVirtualCacheView(image,exception);
3713 statistic_view=AcquireAuthenticCacheView(statistic_image,exception);
3714#if defined(MAGICKCORE_OPENMP_SUPPORT)
3715 #pragma omp parallel for schedule(static) shared(progress,status) \
3716 magick_number_threads(image,statistic_image,statistic_image->rows,1)
3717#endif
3718 for (y=0; y < (ssize_t) statistic_image->rows; y++)
3719 {
3720 const int
3721 id = GetOpenMPThreadId();
3722
3723 const IndexPacket
3724 *magick_restrict indexes;
3725
3726 const PixelPacket
3727 *magick_restrict p;
3728
3729 IndexPacket
3730 *magick_restrict statistic_indexes;
3731
3732 PixelPacket
3733 *magick_restrict q;
3734
3735 ssize_t
3736 x;
3737
3738 if (status == MagickFalse)
3739 continue;
3740 p=GetCacheViewVirtualPixels(image_view,-((ssize_t) neighbor_width/2L),y-
3741 (ssize_t) (neighbor_height/2L),image->columns+neighbor_width,
3742 neighbor_height,exception);
3743 q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns, 1,exception);
3744 if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
3745 {
3746 status=MagickFalse;
3747 continue;
3748 }
3749 indexes=GetCacheViewVirtualIndexQueue(image_view);
3750 statistic_indexes=GetCacheViewAuthenticIndexQueue(statistic_view);
3751 for (x=0; x < (ssize_t) statistic_image->columns; x++)
3752 {
3753 MagickPixelPacket
3754 pixel;
3755
3756 const IndexPacket
3757 *magick_restrict s;
3758
3759 const PixelPacket
3760 *magick_restrict r;
3761
3762 ssize_t
3763 u,
3764 v;
3765
3766 r=p;
3767 s=indexes+x;
3768 ResetPixelList(pixel_list[id]);
3769 for (v=0; v < (ssize_t) neighbor_height; v++)
3770 {
3771 for (u=0; u < (ssize_t) neighbor_width; u++)
3772 InsertPixelList(image,r+u,s+u,pixel_list[id]);
3773 r+=(ptrdiff_t) image->columns+neighbor_width;
3774 s+=(ptrdiff_t) image->columns+neighbor_width;
3775 }
3776 GetMagickPixelPacket(image,&pixel);
3777 SetMagickPixelPacket(image,p+neighbor_width*neighbor_height/2,indexes+x+
3778 neighbor_width*neighbor_height/2,&pixel);
3779 switch (type)
3780 {
3781 case GradientStatistic:
3782 {
3783 MagickPixelPacket
3784 maximum,
3785 minimum;
3786
3787 GetMinimumPixelList(pixel_list[id],&pixel);
3788 minimum=pixel;
3789 GetMaximumPixelList(pixel_list[id],&pixel);
3790 maximum=pixel;
3791 pixel.red=MagickAbsoluteValue(maximum.red-minimum.red);
3792 pixel.green=MagickAbsoluteValue(maximum.green-minimum.green);
3793 pixel.blue=MagickAbsoluteValue(maximum.blue-minimum.blue);
3794 pixel.opacity=MagickAbsoluteValue(maximum.opacity-minimum.opacity);
3795 if (image->colorspace == CMYKColorspace)
3796 pixel.index=MagickAbsoluteValue(maximum.index-minimum.index);
3797 break;
3798 }
3799 case MaximumStatistic:
3800 {
3801 GetMaximumPixelList(pixel_list[id],&pixel);
3802 break;
3803 }
3804 case MeanStatistic:
3805 {
3806 GetMeanPixelList(pixel_list[id],&pixel);
3807 break;
3808 }
3809 case MedianStatistic:
3810 default:
3811 {
3812 GetMedianPixelList(pixel_list[id],&pixel);
3813 break;
3814 }
3815 case MinimumStatistic:
3816 {
3817 GetMinimumPixelList(pixel_list[id],&pixel);
3818 break;
3819 }
3820 case ModeStatistic:
3821 {
3822 GetModePixelList(pixel_list[id],&pixel);
3823 break;
3824 }
3825 case NonpeakStatistic:
3826 {
3827 GetNonpeakPixelList(pixel_list[id],&pixel);
3828 break;
3829 }
3830 case RootMeanSquareStatistic:
3831 {
3832 GetRootMeanSquarePixelList(pixel_list[id],&pixel);
3833 break;
3834 }
3835 case StandardDeviationStatistic:
3836 {
3837 GetStandardDeviationPixelList(pixel_list[id],&pixel);
3838 break;
3839 }
3840 }
3841 if ((channel & RedChannel) != 0)
3842 SetPixelRed(q,ClampToQuantum(pixel.red));
3843 if ((channel & GreenChannel) != 0)
3844 SetPixelGreen(q,ClampToQuantum(pixel.green));
3845 if ((channel & BlueChannel) != 0)
3846 SetPixelBlue(q,ClampToQuantum(pixel.blue));
3847 if ((channel & OpacityChannel) != 0)
3848 SetPixelOpacity(q,ClampToQuantum(pixel.opacity));
3849 if (((channel & IndexChannel) != 0) &&
3850 (image->colorspace == CMYKColorspace))
3851 SetPixelIndex(statistic_indexes+x,ClampToQuantum(pixel.index));
3852 p++;
3853 q++;
3854 }
3855 if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
3856 status=MagickFalse;
3857 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3858 {
3859 MagickBooleanType
3860 proceed;
3861
3862 proceed=SetImageProgress(image,StatisticImageTag,progress++,
3863 image->rows);
3864 if (proceed == MagickFalse)
3865 status=MagickFalse;
3866 }
3867 }
3868 statistic_view=DestroyCacheView(statistic_view);
3869 image_view=DestroyCacheView(image_view);
3870 pixel_list=DestroyPixelListTLS(pixel_list);
3871 if (status == MagickFalse)
3872 statistic_image=DestroyImage(statistic_image);
3873 return(statistic_image);
3874}
_ListNode
Definition statistic.c:3034
_PixelList
Definition statistic.c:3051
_SkipList
Definition statistic.c:3042