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