43#include "MagickCore/studio.h"
44#include "MagickCore/accelerate-private.h"
45#include "MagickCore/animate.h"
46#include "MagickCore/artifact.h"
47#include "MagickCore/blob.h"
48#include "MagickCore/blob-private.h"
49#include "MagickCore/cache.h"
50#include "MagickCore/cache-private.h"
51#include "MagickCore/cache-view.h"
52#include "MagickCore/client.h"
53#include "MagickCore/color.h"
54#include "MagickCore/color-private.h"
55#include "MagickCore/colorspace.h"
56#include "MagickCore/colorspace-private.h"
57#include "MagickCore/composite.h"
58#include "MagickCore/composite-private.h"
59#include "MagickCore/compress.h"
60#include "MagickCore/constitute.h"
61#include "MagickCore/display.h"
62#include "MagickCore/draw.h"
63#include "MagickCore/enhance.h"
64#include "MagickCore/exception.h"
65#include "MagickCore/exception-private.h"
66#include "MagickCore/gem.h"
67#include "MagickCore/gem-private.h"
68#include "MagickCore/geometry.h"
69#include "MagickCore/list.h"
70#include "MagickCore/image-private.h"
71#include "MagickCore/magic.h"
72#include "MagickCore/magick.h"
73#include "MagickCore/memory_.h"
74#include "MagickCore/module.h"
75#include "MagickCore/monitor.h"
76#include "MagickCore/monitor-private.h"
77#include "MagickCore/option.h"
78#include "MagickCore/paint.h"
79#include "MagickCore/pixel-accessor.h"
80#include "MagickCore/profile.h"
81#include "MagickCore/property.h"
82#include "MagickCore/quantize.h"
83#include "MagickCore/quantum-private.h"
84#include "MagickCore/random_.h"
85#include "MagickCore/random-private.h"
86#include "MagickCore/resource_.h"
87#include "MagickCore/segment.h"
88#include "MagickCore/semaphore.h"
89#include "MagickCore/signature-private.h"
90#include "MagickCore/statistic.h"
91#include "MagickCore/statistic-private.h"
92#include "MagickCore/string_.h"
93#include "MagickCore/thread-private.h"
94#include "MagickCore/timer.h"
95#include "MagickCore/utility.h"
96#include "MagickCore/version.h"
138 channel[MaxPixelChannels];
151 rows=MagickMax(GetImageListLength(images),(
size_t)
152 GetMagickResourceLimit(ThreadResource));
153 for (i=0; i < (ssize_t) rows; i++)
155 pixels[i]=(
PixelChannels *) RelinquishMagickMemory(pixels[i]);
176 number_images=GetImageListLength(images);
177 rows=MagickMax(number_images,(
size_t) GetMagickResourceLimit(ThreadResource));
178 pixels=(
PixelChannels **) AcquireQuantumMemory(rows,
sizeof(*pixels));
181 (void) memset(pixels,0,rows*
sizeof(*pixels));
182 columns=MagickMax(number_images,MaxPixelChannels);
183 for (next=images; next != (
Image *) NULL; next=next->next)
184 columns=MagickMax(next->columns,columns);
185 for (i=0; i < (ssize_t) rows; i++)
190 pixels[i]=(
PixelChannels *) AcquireQuantumMemory(columns,
sizeof(**pixels));
192 return(DestroyPixelTLS(images,pixels));
193 for (j=0; j < (ssize_t) columns; j++)
198 for (k=0; k < MaxPixelChannels; k++)
199 pixels[i][j].channel[k]=0.0;
205static inline double EvaluateMax(
const double x,
const double y)
212#if defined(__cplusplus) || defined(c_plusplus)
216static int IntensityCompare(
const void *x,
const void *y)
231 for (i=0; i < MaxPixelChannels; i++)
232 distance+=color_1->channel[i]-(
double) color_2->channel[i];
233 return(distance < 0.0 ? -1 : distance > 0.0 ? 1 : 0);
236#if defined(__cplusplus) || defined(c_plusplus)
240static double ApplyEvaluateOperator(
RandomInfo *random_info,
const Quantum pixel,
241 const MagickEvaluateOperator op,
const double value)
252 case UndefinedEvaluateOperator:
254 case AbsEvaluateOperator:
256 result=(double) fabs((
double) pixel+value);
259 case AddEvaluateOperator:
261 result=(double) pixel+value;
264 case AddModulusEvaluateOperator:
272 result=(double) pixel+value;
273 result-=((double) QuantumRange+1.0)*floor(result/((
double)
277 case AndEvaluateOperator:
279 result=(double) ((ssize_t) pixel & (ssize_t) (value+0.5));
282 case CosineEvaluateOperator:
284 result=(double) QuantumRange*(0.5*cos((
double) (2.0*MagickPI*
285 QuantumScale*(double) pixel*value))+0.5);
288 case DivideEvaluateOperator:
290 result=(double) pixel/(value == 0.0 ? 1.0 : value);
293 case ExponentialEvaluateOperator:
295 result=(double) QuantumRange*exp(value*QuantumScale*(
double) pixel);
298 case GaussianNoiseEvaluateOperator:
300 result=(double) GenerateDifferentialNoise(random_info,pixel,GaussianNoise,
304 case ImpulseNoiseEvaluateOperator:
306 result=(double) GenerateDifferentialNoise(random_info,pixel,ImpulseNoise,
310 case InverseLogEvaluateOperator:
312 result=(double) QuantumRange*pow((value+1.0),QuantumScale*(double)
313 pixel-1.0)*PerceptibleReciprocal(value);
316 case LaplacianNoiseEvaluateOperator:
318 result=(double) GenerateDifferentialNoise(random_info,pixel,
319 LaplacianNoise,value);
322 case LeftShiftEvaluateOperator:
324 result=(double) pixel;
325 for (i=0; i < (ssize_t) value; i++)
329 case LogEvaluateOperator:
331 if ((QuantumScale*(
double) pixel) >= MagickEpsilon)
332 result=(double) QuantumRange*log(QuantumScale*value*
333 (
double) pixel+1.0)/log((
double) (value+1.0));
336 case MaxEvaluateOperator:
338 result=(double) EvaluateMax((
double) pixel,value);
341 case MeanEvaluateOperator:
343 result=(double) pixel+value;
346 case MedianEvaluateOperator:
348 result=(double) pixel+value;
351 case MinEvaluateOperator:
353 result=MagickMin((
double) pixel,value);
356 case MultiplicativeNoiseEvaluateOperator:
358 result=(double) GenerateDifferentialNoise(random_info,pixel,
359 MultiplicativeGaussianNoise,value);
362 case MultiplyEvaluateOperator:
364 result=(double) pixel*value;
367 case OrEvaluateOperator:
369 result=(double) ((ssize_t) pixel | (ssize_t) (value+0.5));
372 case PoissonNoiseEvaluateOperator:
374 result=(double) GenerateDifferentialNoise(random_info,pixel,PoissonNoise,
378 case PowEvaluateOperator:
380 if (((
double) pixel < 0) && ((value-floor(value)) > MagickEpsilon))
381 result=(
double) -((double) QuantumRange*pow(-(QuantumScale*(
double)
382 pixel),(
double) value));
384 result=(double) QuantumRange*pow(QuantumScale*(
double) pixel,
388 case RightShiftEvaluateOperator:
390 result=(double) pixel;
391 for (i=0; i < (ssize_t) value; i++)
395 case RootMeanSquareEvaluateOperator:
397 result=((double) pixel*(
double) pixel+value);
400 case SetEvaluateOperator:
405 case SineEvaluateOperator:
407 result=(double) QuantumRange*(0.5*sin((
double) (2.0*MagickPI*
408 QuantumScale*(double) pixel*value))+0.5);
411 case SubtractEvaluateOperator:
413 result=(double) pixel-value;
416 case SumEvaluateOperator:
418 result=(double) pixel+value;
421 case ThresholdEvaluateOperator:
423 result=(double) (((
double) pixel <= value) ? 0 : QuantumRange);
426 case ThresholdBlackEvaluateOperator:
428 result=(double) (((
double) pixel <= value) ? 0 : pixel);
431 case ThresholdWhiteEvaluateOperator:
433 result=(double) (((
double) pixel > value) ? QuantumRange : pixel);
436 case UniformNoiseEvaluateOperator:
438 result=(double) GenerateDifferentialNoise(random_info,pixel,UniformNoise,
442 case XorEvaluateOperator:
444 result=(double) ((ssize_t) pixel ^ (ssize_t) (value+0.5));
462 columns=images->columns;
464 for (p=images; p != (
Image *) NULL; p=p->next)
466 if (p->number_channels > q->number_channels)
468 if (p->columns > columns)
473 return(CloneImage(q,columns,rows,MagickTrue,exception));
476MagickExport
Image *EvaluateImages(
const Image *images,
479#define EvaluateImageTag "Evaluate/Image"
498 **magick_restrict evaluate_pixels;
501 **magick_restrict random_info;
510#if defined(MAGICKCORE_OPENMP_SUPPORT)
515 assert(images != (
Image *) NULL);
516 assert(images->signature == MagickCoreSignature);
518 assert(exception->signature == MagickCoreSignature);
519 if (IsEventLogging() != MagickFalse)
520 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",images->filename);
521 image=AcquireImageCanvas(images,exception);
522 if (image == (
Image *) NULL)
523 return((
Image *) NULL);
524 if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
526 image=DestroyImage(image);
527 return((
Image *) NULL);
529 number_images=GetImageListLength(images);
530 evaluate_pixels=AcquirePixelTLS(images);
533 image=DestroyImage(image);
534 (void) ThrowMagickException(exception,GetMagickModule(),
535 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",images->filename);
536 return((
Image *) NULL);
538 image_view=(
CacheView **) AcquireQuantumMemory(number_images,
539 sizeof(*image_view));
542 image=DestroyImage(image);
543 evaluate_pixels=DestroyPixelTLS(images,evaluate_pixels);
544 (void) ThrowMagickException(exception,GetMagickModule(),
545 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",images->filename);
549 for (n=0; n < (ssize_t) number_images; n++)
551 image_view[n]=AcquireVirtualCacheView(view,exception);
552 view=GetNextImageInList(view);
559 random_info=AcquireRandomInfoTLS();
560 evaluate_view=AcquireAuthenticCacheView(image,exception);
561 if (op == MedianEvaluateOperator)
563#if defined(MAGICKCORE_OPENMP_SUPPORT)
564 key=GetRandomSecretKey(random_info[0]);
565 #pragma omp parallel for schedule(static) shared(progress,status) \
566 magick_number_threads(image,images,image->rows,key == ~0UL)
568 for (y=0; y < (ssize_t) image->rows; y++)
571 id = GetOpenMPThreadId();
588 if (status == MagickFalse)
590 p=(
const Quantum **) AcquireQuantumMemory(number_images,
sizeof(*p));
591 if (p == (
const Quantum **) NULL)
594 (void) ThrowMagickException(exception,GetMagickModule(),
595 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",
599 for (j=0; j < (ssize_t) number_images; j++)
601 p[j]=GetCacheViewVirtualPixels(image_view[j],0,y,image->columns,1,
603 if (p[j] == (
const Quantum *) NULL)
606 q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,image->columns,1,
608 if ((j < (ssize_t) number_images) || (q == (Quantum *) NULL))
613 evaluate_pixel=evaluate_pixels[id];
614 for (x=0; x < (ssize_t) image->columns; x++)
623 for (j=0; j < (ssize_t) number_images; j++)
625 for (i=0; i < MaxPixelChannels; i++)
626 evaluate_pixel[j].channel[i]=0.0;
627 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
629 PixelChannel channel = GetPixelChannelChannel(image,i);
630 PixelTrait traits = GetPixelChannelTraits(next,channel);
631 PixelTrait evaluate_traits = GetPixelChannelTraits(image,channel);
632 if ((traits == UndefinedPixelTrait) ||
633 (evaluate_traits == UndefinedPixelTrait) ||
634 ((traits & UpdatePixelTrait) == 0))
636 evaluate_pixel[j].channel[i]=ApplyEvaluateOperator(
637 random_info[
id],GetPixelChannel(next,channel,p[j]),op,
638 evaluate_pixel[j].channel[i]);
640 p[j]+=GetPixelChannels(next);
641 next=GetNextImageInList(next);
643 qsort((
void *) evaluate_pixel,number_images,
sizeof(*evaluate_pixel),
645 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
647 PixelChannel channel = GetPixelChannelChannel(image,i);
648 PixelTrait traits = GetPixelChannelTraits(image,channel);
649 if ((traits == UndefinedPixelTrait) ||
650 ((traits & UpdatePixelTrait) == 0))
652 q[i]=ClampToQuantum(evaluate_pixel[number_images/2].channel[i]);
654 q+=GetPixelChannels(image);
656 p=(
const Quantum **) RelinquishMagickMemory((
void *) p);
657 if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
659 if (images->progress_monitor != (MagickProgressMonitor) NULL)
664#if defined(MAGICKCORE_OPENMP_SUPPORT)
668 proceed=SetImageProgress(images,EvaluateImageTag,progress,
670 if (proceed == MagickFalse)
677#if defined(MAGICKCORE_OPENMP_SUPPORT)
678 key=GetRandomSecretKey(random_info[0]);
679 #pragma omp parallel for schedule(static) shared(progress,status) \
680 magick_number_threads(image,images,image->rows,key == ~0UL)
682 for (y=0; y < (ssize_t) image->rows; y++)
688 id = GetOpenMPThreadId();
706 if (status == MagickFalse)
708 p=(
const Quantum **) AcquireQuantumMemory(number_images,
sizeof(*p));
709 if (p == (
const Quantum **) NULL)
712 (void) ThrowMagickException(exception,GetMagickModule(),
713 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",
717 for (j=0; j < (ssize_t) number_images; j++)
719 p[j]=GetCacheViewVirtualPixels(image_view[j],0,y,image->columns,1,
721 if (p[j] == (
const Quantum *) NULL)
724 q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,image->columns,1,
726 if ((j < (ssize_t) number_images) || (q == (Quantum *) NULL))
731 evaluate_pixel=evaluate_pixels[id];
732 for (j=0; j < (ssize_t) image->columns; j++)
733 for (i=0; i < MaxPixelChannels; i++)
734 evaluate_pixel[j].channel[i]=0.0;
736 for (j=0; j < (ssize_t) number_images; j++)
738 for (x=0; x < (ssize_t) image->columns; x++)
740 for (i=0; i < (ssize_t) GetPixelChannels(next); i++)
742 PixelChannel channel = GetPixelChannelChannel(image,i);
743 PixelTrait traits = GetPixelChannelTraits(next,channel);
744 PixelTrait evaluate_traits = GetPixelChannelTraits(image,channel);
745 if ((traits == UndefinedPixelTrait) ||
746 (evaluate_traits == UndefinedPixelTrait))
748 if ((traits & UpdatePixelTrait) == 0)
750 evaluate_pixel[x].channel[i]=ApplyEvaluateOperator(
751 random_info[
id],GetPixelChannel(next,channel,p[j]),j == 0 ?
752 AddEvaluateOperator : op,evaluate_pixel[x].channel[i]);
754 p[j]+=GetPixelChannels(next);
756 next=GetNextImageInList(next);
758 for (x=0; x < (ssize_t) image->columns; x++)
762 case MeanEvaluateOperator:
764 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
765 evaluate_pixel[x].channel[i]/=(
double) number_images;
768 case MultiplyEvaluateOperator:
770 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
772 for (j=0; j < (ssize_t) (number_images-1); j++)
773 evaluate_pixel[x].channel[i]*=QuantumScale;
777 case RootMeanSquareEvaluateOperator:
779 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
780 evaluate_pixel[x].channel[i]=sqrt(evaluate_pixel[x].channel[i]/
788 for (x=0; x < (ssize_t) image->columns; x++)
790 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
792 PixelChannel channel = GetPixelChannelChannel(image,i);
793 PixelTrait traits = GetPixelChannelTraits(image,channel);
794 if ((traits == UndefinedPixelTrait) ||
795 ((traits & UpdatePixelTrait) == 0))
797 q[i]=ClampToQuantum(evaluate_pixel[x].channel[i]);
799 q+=GetPixelChannels(image);
801 p=(
const Quantum **) RelinquishMagickMemory((
void *) p);
802 if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
804 if (images->progress_monitor != (MagickProgressMonitor) NULL)
809#if defined(MAGICKCORE_OPENMP_SUPPORT)
813 proceed=SetImageProgress(images,EvaluateImageTag,progress,
815 if (proceed == MagickFalse)
820 for (n=0; n < (ssize_t) number_images; n++)
821 image_view[n]=DestroyCacheView(image_view[n]);
822 image_view=(
CacheView **) RelinquishMagickMemory(image_view);
823 evaluate_view=DestroyCacheView(evaluate_view);
824 evaluate_pixels=DestroyPixelTLS(images,evaluate_pixels);
825 random_info=DestroyRandomInfoTLS(random_info);
826 if (status == MagickFalse)
827 image=DestroyImage(image);
831MagickExport MagickBooleanType EvaluateImage(
Image *image,
832 const MagickEvaluateOperator op,
const double value,
ExceptionInfo *exception)
848 **magick_restrict random_info;
853#if defined(MAGICKCORE_OPENMP_SUPPORT)
858 assert(image != (
Image *) NULL);
859 assert(image->signature == MagickCoreSignature);
861 assert(exception->signature == MagickCoreSignature);
862 if (IsEventLogging() != MagickFalse)
863 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
864 if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
869 artifact=GetImageArtifact(image,
"evaluate:clamp");
870 if (artifact != (
const char *) NULL)
871 clamp=IsStringTrue(artifact);
872 random_info=AcquireRandomInfoTLS();
873 image_view=AcquireAuthenticCacheView(image,exception);
874#if defined(MAGICKCORE_OPENMP_SUPPORT)
875 key=GetRandomSecretKey(random_info[0]);
876 #pragma omp parallel for schedule(static) shared(progress,status) \
877 magick_number_threads(image,image,image->rows,key == ~0UL)
879 for (y=0; y < (ssize_t) image->rows; y++)
882 id = GetOpenMPThreadId();
890 if (status == MagickFalse)
892 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
893 if (q == (Quantum *) NULL)
898 for (x=0; x < (ssize_t) image->columns; x++)
906 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
908 PixelChannel channel = GetPixelChannelChannel(image,i);
909 PixelTrait traits = GetPixelChannelTraits(image,channel);
910 if (traits == UndefinedPixelTrait)
912 if ((traits & CopyPixelTrait) != 0)
914 if ((traits & UpdatePixelTrait) == 0)
916 result=ApplyEvaluateOperator(random_info[
id],q[i],op,value);
917 if (op == MeanEvaluateOperator)
919 q[i]=clamp != MagickFalse ? ClampPixel(result) : ClampToQuantum(result);
921 q+=GetPixelChannels(image);
923 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
925 if (image->progress_monitor != (MagickProgressMonitor) NULL)
930#if defined(MAGICKCORE_OPENMP_SUPPORT)
934 proceed=SetImageProgress(image,EvaluateImageTag,progress,image->rows);
935 if (proceed == MagickFalse)
939 image_view=DestroyCacheView(image_view);
940 random_info=DestroyRandomInfoTLS(random_info);
978static Quantum ApplyFunction(Quantum pixel,
const MagickFunction function,
979 const size_t number_parameters,
const double *parameters,
992 case PolynomialFunction:
999 for (i=0; i < (ssize_t) number_parameters; i++)
1000 result=result*QuantumScale*(
double) pixel+parameters[i];
1001 result*=(double) QuantumRange;
1004 case SinusoidFunction:
1015 frequency=(number_parameters >= 1) ? parameters[0] : 1.0;
1016 phase=(number_parameters >= 2) ? parameters[1] : 0.0;
1017 amplitude=(number_parameters >= 3) ? parameters[2] : 0.5;
1018 bias=(number_parameters >= 4) ? parameters[3] : 0.5;
1019 result=(double) QuantumRange*(amplitude*sin((
double) (2.0*
1020 MagickPI*(frequency*QuantumScale*(
double) pixel+phase/360.0)))+bias);
1023 case ArcsinFunction:
1035 width=(number_parameters >= 1) ? parameters[0] : 1.0;
1036 center=(number_parameters >= 2) ? parameters[1] : 0.5;
1037 range=(number_parameters >= 3) ? parameters[2] : 1.0;
1038 bias=(number_parameters >= 4) ? parameters[3] : 0.5;
1039 result=2.0*PerceptibleReciprocal(width)*(QuantumScale*(double) pixel-
1042 result=bias-range/2.0;
1045 result=bias+range/2.0;
1047 result=(double) (range/MagickPI*asin((
double) result)+bias);
1048 result*=(double) QuantumRange;
1051 case ArctanFunction:
1062 slope=(number_parameters >= 1) ? parameters[0] : 1.0;
1063 center=(number_parameters >= 2) ? parameters[1] : 0.5;
1064 range=(number_parameters >= 3) ? parameters[2] : 1.0;
1065 bias=(number_parameters >= 4) ? parameters[3] : 0.5;
1066 result=MagickPI*slope*(QuantumScale*(double) pixel-center);
1067 result=(double) QuantumRange*(range/MagickPI*atan((
double) result)+bias);
1070 case UndefinedFunction:
1073 return(ClampToQuantum(result));
1076MagickExport MagickBooleanType FunctionImage(
Image *image,
1077 const MagickFunction function,
const size_t number_parameters,
1080#define FunctionImageTag "Function/Image "
1094 assert(image != (
Image *) NULL);
1095 assert(image->signature == MagickCoreSignature);
1097 assert(exception->signature == MagickCoreSignature);
1098 if (IsEventLogging() != MagickFalse)
1099 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
1100#if defined(MAGICKCORE_OPENCL_SUPPORT)
1101 if (AccelerateFunctionImage(image,function,number_parameters,parameters,
1102 exception) != MagickFalse)
1105 if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
1106 return(MagickFalse);
1109 image_view=AcquireAuthenticCacheView(image,exception);
1110#if defined(MAGICKCORE_OPENMP_SUPPORT)
1111 #pragma omp parallel for schedule(static) shared(progress,status) \
1112 magick_number_threads(image,image,image->rows,1)
1114 for (y=0; y < (ssize_t) image->rows; y++)
1122 if (status == MagickFalse)
1124 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
1125 if (q == (Quantum *) NULL)
1130 for (x=0; x < (ssize_t) image->columns; x++)
1135 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1137 PixelChannel channel = GetPixelChannelChannel(image,i);
1138 PixelTrait traits = GetPixelChannelTraits(image,channel);
1139 if (traits == UndefinedPixelTrait)
1141 if ((traits & UpdatePixelTrait) == 0)
1143 q[i]=ApplyFunction(q[i],function,number_parameters,parameters,
1146 q+=GetPixelChannels(image);
1148 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1150 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1155#if defined(MAGICKCORE_OPENMP_SUPPORT)
1159 proceed=SetImageProgress(image,FunctionImageTag,progress,image->rows);
1160 if (proceed == MagickFalse)
1164 image_view=DestroyCacheView(image_view);
1195MagickExport MagickBooleanType GetImageEntropy(
const Image *image,
1199 *channel_statistics;
1201 assert(image != (
Image *) NULL);
1202 assert(image->signature == MagickCoreSignature);
1203 if (IsEventLogging() != MagickFalse)
1204 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
1205 channel_statistics=GetImageStatistics(image,exception);
1207 return(MagickFalse);
1208 *entropy=channel_statistics[CompositePixelChannel].entropy;
1210 channel_statistics);
1243MagickExport MagickBooleanType GetImageExtrema(
const Image *image,
1253 assert(image != (
Image *) NULL);
1254 assert(image->signature == MagickCoreSignature);
1255 if (IsEventLogging() != MagickFalse)
1256 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
1257 status=GetImageRange(image,&min,&max,exception);
1258 *minima=(size_t) ceil(min-0.5);
1259 *maxima=(size_t) floor(max+0.5);
1293MagickExport MagickBooleanType GetImageKurtosis(
const Image *image,
1297 *channel_statistics;
1299 assert(image != (
Image *) NULL);
1300 assert(image->signature == MagickCoreSignature);
1301 if (IsEventLogging() != MagickFalse)
1302 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
1303 channel_statistics=GetImageStatistics(image,exception);
1305 return(MagickFalse);
1306 *kurtosis=channel_statistics[CompositePixelChannel].kurtosis;
1307 *skewness=channel_statistics[CompositePixelChannel].skewness;
1309 channel_statistics);
1343MagickExport MagickBooleanType GetImageMean(
const Image *image,
double *mean,
1347 *channel_statistics;
1349 assert(image != (
Image *) NULL);
1350 assert(image->signature == MagickCoreSignature);
1351 if (IsEventLogging() != MagickFalse)
1352 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
1353 channel_statistics=GetImageStatistics(image,exception);
1355 return(MagickFalse);
1356 *mean=channel_statistics[CompositePixelChannel].mean;
1357 *standard_deviation=
1358 channel_statistics[CompositePixelChannel].standard_deviation;
1360 channel_statistics);
1391MagickExport MagickBooleanType GetImageMedian(
const Image *image,
double *median,
1395 *channel_statistics;
1397 assert(image != (
Image *) NULL);
1398 assert(image->signature == MagickCoreSignature);
1399 if (IsEventLogging() != MagickFalse)
1400 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
1401 channel_statistics=GetImageStatistics(image,exception);
1403 return(MagickFalse);
1404 *median=channel_statistics[CompositePixelChannel].median;
1406 channel_statistics);
1439#define MaxNumberImageMoments 8
1449 M00[2*MaxPixelChannels+1],
1450 M01[2*MaxPixelChannels+1],
1451 M02[2*MaxPixelChannels+1],
1452 M03[2*MaxPixelChannels+1],
1453 M10[2*MaxPixelChannels+1],
1454 M11[2*MaxPixelChannels+1],
1455 M12[2*MaxPixelChannels+1],
1456 M20[2*MaxPixelChannels+1],
1457 M21[2*MaxPixelChannels+1],
1458 M22[2*MaxPixelChannels+1],
1459 M30[2*MaxPixelChannels+1];
1462 centroid[2*MaxPixelChannels+1];
1468 assert(image != (
Image *) NULL);
1469 assert(image->signature == MagickCoreSignature);
1470 if (IsEventLogging() != MagickFalse)
1471 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
1472 channel_moments=(
ChannelMoments *) AcquireQuantumMemory(MaxPixelChannels+1,
1473 sizeof(*channel_moments));
1475 return(channel_moments);
1476 (void) memset(channel_moments,0,(MaxPixelChannels+1)*
1477 sizeof(*channel_moments));
1478 (void) memset(centroid,0,
sizeof(centroid));
1479 (void) memset(M00,0,
sizeof(M00));
1480 (void) memset(M01,0,
sizeof(M01));
1481 (void) memset(M02,0,
sizeof(M02));
1482 (void) memset(M03,0,
sizeof(M03));
1483 (void) memset(M10,0,
sizeof(M10));
1484 (void) memset(M11,0,
sizeof(M11));
1485 (void) memset(M12,0,
sizeof(M12));
1486 (void) memset(M20,0,
sizeof(M20));
1487 (void) memset(M21,0,
sizeof(M21));
1488 (void) memset(M22,0,
sizeof(M22));
1489 (void) memset(M30,0,
sizeof(M30));
1490 image_view=AcquireVirtualCacheView(image,exception);
1491 for (y=0; y < (ssize_t) image->rows; y++)
1502 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1503 if (p == (
const Quantum *) NULL)
1505 for (x=0; x < (ssize_t) image->columns; x++)
1510 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1512 PixelChannel channel = GetPixelChannelChannel(image,i);
1513 PixelTrait traits = GetPixelChannelTraits(image,channel);
1514 if (traits == UndefinedPixelTrait)
1516 if ((traits & UpdatePixelTrait) == 0)
1518 M00[channel]+=QuantumScale*(double) p[i];
1519 M00[MaxPixelChannels]+=QuantumScale*(double) p[i];
1520 M10[channel]+=x*QuantumScale*(double) p[i];
1521 M10[MaxPixelChannels]+=x*QuantumScale*(double) p[i];
1522 M01[channel]+=y*QuantumScale*(double) p[i];
1523 M01[MaxPixelChannels]+=y*QuantumScale*(double) p[i];
1525 p+=GetPixelChannels(image);
1528 for (c=0; c <= MaxPixelChannels; c++)
1533 centroid[c].x=M10[c]*PerceptibleReciprocal(M00[c]);
1534 centroid[c].y=M01[c]*PerceptibleReciprocal(M00[c]);
1536 for (y=0; y < (ssize_t) image->rows; y++)
1547 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1548 if (p == (
const Quantum *) NULL)
1550 for (x=0; x < (ssize_t) image->columns; x++)
1555 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1557 PixelChannel channel = GetPixelChannelChannel(image,i);
1558 PixelTrait traits = GetPixelChannelTraits(image,channel);
1559 if (traits == UndefinedPixelTrait)
1561 if ((traits & UpdatePixelTrait) == 0)
1563 M11[channel]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
1564 QuantumScale*(double) p[i];
1565 M11[MaxPixelChannels]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
1566 QuantumScale*(double) p[i];
1567 M20[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1568 QuantumScale*(double) p[i];
1569 M20[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1570 QuantumScale*(double) p[i];
1571 M02[channel]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
1572 QuantumScale*(double) p[i];
1573 M02[MaxPixelChannels]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
1574 QuantumScale*(double) p[i];
1575 M21[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1576 (y-centroid[channel].y)*QuantumScale*(
double) p[i];
1577 M21[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1578 (y-centroid[channel].y)*QuantumScale*(
double) p[i];
1579 M12[channel]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
1580 (y-centroid[channel].y)*QuantumScale*(
double) p[i];
1581 M12[MaxPixelChannels]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
1582 (y-centroid[channel].y)*QuantumScale*(
double) p[i];
1583 M22[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1584 (y-centroid[channel].y)*(y-centroid[channel].y)*QuantumScale*(double)
1586 M22[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1587 (y-centroid[channel].y)*(y-centroid[channel].y)*QuantumScale*(double)
1589 M30[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1590 (x-centroid[channel].x)*QuantumScale*(
double) p[i];
1591 M30[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1592 (x-centroid[channel].x)*QuantumScale*(
double) p[i];
1593 M03[channel]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
1594 (y-centroid[channel].y)*QuantumScale*(
double) p[i];
1595 M03[MaxPixelChannels]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
1596 (y-centroid[channel].y)*QuantumScale*(
double) p[i];
1598 p+=GetPixelChannels(image);
1601 channels=(double) GetImageChannels(image);
1602 M00[MaxPixelChannels]/=channels;
1603 M01[MaxPixelChannels]/=channels;
1604 M02[MaxPixelChannels]/=channels;
1605 M03[MaxPixelChannels]/=channels;
1606 M10[MaxPixelChannels]/=channels;
1607 M11[MaxPixelChannels]/=channels;
1608 M12[MaxPixelChannels]/=channels;
1609 M20[MaxPixelChannels]/=channels;
1610 M21[MaxPixelChannels]/=channels;
1611 M22[MaxPixelChannels]/=channels;
1612 M30[MaxPixelChannels]/=channels;
1613 for (c=0; c <= MaxPixelChannels; c++)
1618 channel_moments[c].centroid=centroid[c];
1619 channel_moments[c].ellipse_axis.x=sqrt((2.0*PerceptibleReciprocal(M00[c]))*
1620 ((M20[c]+M02[c])+sqrt(4.0*M11[c]*M11[c]+(M20[c]-M02[c])*(M20[c]-M02[c]))));
1621 channel_moments[c].ellipse_axis.y=sqrt((2.0*PerceptibleReciprocal(M00[c]))*
1622 ((M20[c]+M02[c])-sqrt(4.0*M11[c]*M11[c]+(M20[c]-M02[c])*(M20[c]-M02[c]))));
1623 channel_moments[c].ellipse_angle=RadiansToDegrees(1.0/2.0*atan(2.0*
1624 M11[c]*PerceptibleReciprocal(M20[c]-M02[c])));
1625 if (fabs(M11[c]) < 0.0)
1627 if ((fabs(M20[c]-M02[c]) >= 0.0) &&
1628 ((M20[c]-M02[c]) < 0.0))
1629 channel_moments[c].ellipse_angle+=90.0;
1634 if (fabs(M20[c]-M02[c]) >= 0.0)
1636 if ((M20[c]-M02[c]) < 0.0)
1637 channel_moments[c].ellipse_angle+=90.0;
1639 channel_moments[c].ellipse_angle+=180.0;
1643 if ((fabs(M20[c]-M02[c]) >= 0.0) && ((M20[c]-M02[c]) < 0.0))
1644 channel_moments[c].ellipse_angle+=90.0;
1645 channel_moments[c].ellipse_eccentricity=sqrt(1.0-(
1646 channel_moments[c].ellipse_axis.y*
1647 channel_moments[c].ellipse_axis.y*PerceptibleReciprocal(
1648 channel_moments[c].ellipse_axis.x*
1649 channel_moments[c].ellipse_axis.x)));
1650 channel_moments[c].ellipse_intensity=M00[c]*
1651 PerceptibleReciprocal(MagickPI*channel_moments[c].ellipse_axis.x*
1652 channel_moments[c].ellipse_axis.y+MagickEpsilon);
1654 for (c=0; c <= MaxPixelChannels; c++)
1661 M11[c]*=PerceptibleReciprocal(pow(M00[c],1.0+(1.0+1.0)/2.0));
1662 M20[c]*=PerceptibleReciprocal(pow(M00[c],1.0+(2.0+0.0)/2.0));
1663 M02[c]*=PerceptibleReciprocal(pow(M00[c],1.0+(0.0+2.0)/2.0));
1664 M21[c]*=PerceptibleReciprocal(pow(M00[c],1.0+(2.0+1.0)/2.0));
1665 M12[c]*=PerceptibleReciprocal(pow(M00[c],1.0+(1.0+2.0)/2.0));
1666 M22[c]*=PerceptibleReciprocal(pow(M00[c],1.0+(2.0+2.0)/2.0));
1667 M30[c]*=PerceptibleReciprocal(pow(M00[c],1.0+(3.0+0.0)/2.0));
1668 M03[c]*=PerceptibleReciprocal(pow(M00[c],1.0+(0.0+3.0)/2.0));
1671 image_view=DestroyCacheView(image_view);
1672 for (c=0; c <= MaxPixelChannels; c++)
1677 channel_moments[c].invariant[0]=M20[c]+M02[c];
1678 channel_moments[c].invariant[1]=(M20[c]-M02[c])*(M20[c]-M02[c])+4.0*M11[c]*
1680 channel_moments[c].invariant[2]=(M30[c]-3.0*M12[c])*(M30[c]-3.0*M12[c])+
1681 (3.0*M21[c]-M03[c])*(3.0*M21[c]-M03[c]);
1682 channel_moments[c].invariant[3]=(M30[c]+M12[c])*(M30[c]+M12[c])+
1683 (M21[c]+M03[c])*(M21[c]+M03[c]);
1684 channel_moments[c].invariant[4]=(M30[c]-3.0*M12[c])*(M30[c]+M12[c])*
1685 ((M30[c]+M12[c])*(M30[c]+M12[c])-3.0*(M21[c]+M03[c])*(M21[c]+M03[c]))+
1686 (3.0*M21[c]-M03[c])*(M21[c]+M03[c])*(3.0*(M30[c]+M12[c])*(M30[c]+M12[c])-
1687 (M21[c]+M03[c])*(M21[c]+M03[c]));
1688 channel_moments[c].invariant[5]=(M20[c]-M02[c])*((M30[c]+M12[c])*
1689 (M30[c]+M12[c])-(M21[c]+M03[c])*(M21[c]+M03[c]))+4.0*M11[c]*
1690 (M30[c]+M12[c])*(M21[c]+M03[c]);
1691 channel_moments[c].invariant[6]=(3.0*M21[c]-M03[c])*(M30[c]+M12[c])*
1692 ((M30[c]+M12[c])*(M30[c]+M12[c])-3.0*(M21[c]+M03[c])*(M21[c]+M03[c]))-
1693 (M30[c]-3*M12[c])*(M21[c]+M03[c])*(3.0*(M30[c]+M12[c])*(M30[c]+M12[c])-
1694 (M21[c]+M03[c])*(M21[c]+M03[c]));
1695 channel_moments[c].invariant[7]=M11[c]*((M30[c]+M12[c])*(M30[c]+M12[c])-
1696 (M03[c]+M21[c])*(M03[c]+M21[c]))-(M20[c]-M02[c])*(M30[c]+M12[c])*
1699 if (y < (ssize_t) image->rows)
1700 channel_moments=(
ChannelMoments *) RelinquishMagickMemory(channel_moments);
1701 return(channel_moments);
1751 MaxPixelChannels+1UL,
sizeof(*perceptual_hash));
1754 artifact=GetImageArtifact(image,
"phash:colorspaces");
1755 if (artifact != NULL)
1756 colorspaces=AcquireString(artifact);
1758 colorspaces=AcquireString(
"xyY,HSB");
1759 perceptual_hash[0].number_colorspaces=0;
1760 perceptual_hash[0].number_channels=0;
1762 for (i=0; (p=StringToken(
",",&q)) != (
char *) NULL; i++)
1777 if (i >= MaximumNumberOfPerceptualColorspaces)
1779 colorspace=ParseCommandOption(MagickColorspaceOptions,MagickFalse,p);
1782 perceptual_hash[0].colorspace[i]=(ColorspaceType) colorspace;
1783 hash_image=BlurImage(image,0.0,1.0,exception);
1784 if (hash_image == (
Image *) NULL)
1786 hash_image->depth=8;
1787 status=TransformImageColorspace(hash_image,(ColorspaceType) colorspace,
1789 if (status == MagickFalse)
1791 moments=GetImageMoments(hash_image,exception);
1792 perceptual_hash[0].number_colorspaces++;
1793 perceptual_hash[0].number_channels+=GetImageChannels(hash_image);
1794 hash_image=DestroyImage(hash_image);
1797 for (channel=0; channel <= MaxPixelChannels; channel++)
1798 for (j=0; j < MaximumNumberOfImageMoments; j++)
1799 perceptual_hash[channel].phash[i][j]=
1800 (-MagickLog10(moments[channel].invariant[j]));
1803 colorspaces=DestroyString(colorspaces);
1804 return(perceptual_hash);
1836MagickExport MagickBooleanType GetImageRange(
const Image *image,
double *minima,
1849 assert(image != (
Image *) NULL);
1850 assert(image->signature == MagickCoreSignature);
1851 if (IsEventLogging() != MagickFalse)
1852 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
1854 initialize=MagickTrue;
1857 image_view=AcquireVirtualCacheView(image,exception);
1858#if defined(MAGICKCORE_OPENMP_SUPPORT)
1859 #pragma omp parallel for schedule(static) shared(status,initialize) \
1860 magick_number_threads(image,image,image->rows,1)
1862 for (y=0; y < (ssize_t) image->rows; y++)
1877 if (status == MagickFalse)
1879 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1880 if (p == (
const Quantum *) NULL)
1885 row_initialize=MagickTrue;
1886 for (x=0; x < (ssize_t) image->columns; x++)
1891 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1893 PixelChannel channel = GetPixelChannelChannel(image,i);
1894 PixelTrait traits = GetPixelChannelTraits(image,channel);
1895 if (traits == UndefinedPixelTrait)
1897 if ((traits & UpdatePixelTrait) == 0)
1899 if (row_initialize != MagickFalse)
1901 row_minima=(double) p[i];
1902 row_maxima=(double) p[i];
1903 row_initialize=MagickFalse;
1907 if ((
double) p[i] < row_minima)
1908 row_minima=(
double) p[i];
1909 if ((
double) p[i] > row_maxima)
1910 row_maxima=(
double) p[i];
1913 p+=GetPixelChannels(image);
1915#if defined(MAGICKCORE_OPENMP_SUPPORT)
1916#pragma omp critical (MagickCore_GetImageRange)
1919 if (initialize != MagickFalse)
1923 initialize=MagickFalse;
1927 if (row_minima < *minima)
1929 if (row_maxima > *maxima)
1934 image_view=DestroyCacheView(image_view);
1972static ssize_t GetMedianPixel(Quantum *pixels,
const size_t n)
1974#define SwapPixels(alpha,beta) \
1976 Quantum gamma=(alpha); \
1977 (alpha)=(beta);(beta)=gamma; \
1982 high = (ssize_t) n-1,
1983 median = (low+high)/2;
1994 if (high == (low+1))
1996 if (pixels[low] > pixels[high])
1997 SwapPixels(pixels[low],pixels[high]);
2000 if (pixels[mid] > pixels[high])
2001 SwapPixels(pixels[mid],pixels[high]);
2002 if (pixels[low] > pixels[high])
2003 SwapPixels(pixels[low], pixels[high]);
2004 if (pixels[mid] > pixels[low])
2005 SwapPixels(pixels[mid],pixels[low]);
2006 SwapPixels(pixels[mid],pixels[low+1]);
2009 do l++;
while (pixels[low] > pixels[l]);
2010 do h--;
while (pixels[h] > pixels[low]);
2013 SwapPixels(pixels[l],pixels[h]);
2015 SwapPixels(pixels[low],pixels[h]);
2023static inline long double PerceptibleReciprocalLD(
const long double x)
2031 sign=x < 0.0 ? -1.0 : 1.0;
2032 if ((sign*x) >= MagickEpsilon)
2034 return(sign/MagickEpsilon);
2041 *channel_statistics;
2069 assert(image != (
Image *) NULL);
2070 assert(image->signature == MagickCoreSignature);
2071 if (IsEventLogging() != MagickFalse)
2072 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
2073 histogram=(
double *) AcquireQuantumMemory(MaxMap+1UL,
2074 MagickMax(GetPixelChannels(image),1)*
sizeof(*histogram));
2076 MaxPixelChannels+1,
sizeof(*channel_statistics));
2078 (histogram == (
double *) NULL))
2080 if (histogram != (
double *) NULL)
2081 histogram=(
double *) RelinquishMagickMemory(histogram);
2084 channel_statistics);
2085 (void) ThrowMagickException(exception,GetMagickModule(),
2086 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
2087 return(channel_statistics);
2089 (void) memset(channel_statistics,0,(MaxPixelChannels+1)*
2090 sizeof(*channel_statistics));
2091 for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
2096 cs->maxima=(-MagickMaximumValue);
2097 cs->minima=MagickMaximumValue;
2098 cs->sum=cs->sumLD=0;
2100 cs->standard_deviation = cs->variance = 0.0;
2101 cs->skewness = cs->kurtosis = 0.0;
2104 (void) memset(histogram,0,(MaxMap+1)*GetPixelChannels(image)*
2105 sizeof(*histogram));
2106 for (y=0; y < (ssize_t) image->rows; y++)
2117 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
2118 if (p == (
const Quantum *) NULL)
2120 for (x=0; x < (ssize_t) image->columns; x++)
2122 if (GetPixelReadMask(image,p) <= (QuantumRange/2))
2124 p+=GetPixelChannels(image);
2127 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2132 PixelChannel channel = GetPixelChannelChannel(image,i);
2133 PixelTrait traits = GetPixelChannelTraits(image,channel);
2134 if (traits == UndefinedPixelTrait)
2136 cs=channel_statistics+i;
2137 if (cs->depth != MAGICKCORE_QUANTUM_DEPTH)
2140 range=GetQuantumRange(depth);
2141 status=p[i] != ScaleAnyToQuantum(ScaleQuantumToAny(p[i],range),
2142 range) ? MagickTrue : MagickFalse;
2143 if (status != MagickFalse)
2146 if (cs->depth > channel_statistics[CompositePixelChannel].depth)
2147 channel_statistics[CompositePixelChannel].depth=cs->depth;
2155 if ((
double) p[i] < cs->minima)
2156 cs->minima=(double) p[i];
2157 if ((
double) p[i] > cs->maxima)
2158 cs->maxima=(
double) p[i];
2160 histogram[(ssize_t) GetPixelChannels(image)*ScaleQuantumToMap(
2161 ClampToQuantum((
double) p[i]))+i]++;
2163 cs->sumLD += (
long double) p[i];
2167 cs->sum_squared+=(double) p[i]*(
double) p[i];
2168 cs->sum_cubed+=(double) p[i]*(
double) p[i]*(double) p[i];
2169 cs->sum_fourth_power+=(double) p[i]*(
double) p[i]*(double) p[i]*
2175 long double delta, delta_n, delta_n2, term1;
2176 double n1 = cs->area-1;
2177 double n = cs->area;
2179 delta = (double) p[i] - cs->M1;
2180 delta_n = delta / n;
2181 delta_n2 = delta_n * delta_n;
2182 term1 = delta * delta_n * n1;
2183 cs->M4 += term1 * delta_n2 * (n*n - 3*n + 3) + 6 * delta_n2 *
2184 cs->M2 - 4 * delta_n * cs->M3;
2185 cs->M3 += term1 * delta_n * (n - 2) - 3 * delta_n * cs->M2;
2190 p+=GetPixelChannels(image);
2194 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2199 PixelChannel channel = GetPixelChannelChannel(image,i);
2200 PixelTrait traits = GetPixelChannelTraits(image,channel);
2201 double AdjArea = 1.0;
2202 if (traits == UndefinedPixelTrait)
2204 cs = &channel_statistics[i];
2207 cs->mean = cs->sumLD / (
long double) cs->area;
2209 AdjArea = cs->area / ( cs->area - 1.0);
2213 cs->sum = (double)cs->sum;
2215 cs->standard_deviation = 0.0;
2221 cs->standard_deviation = sqrtl(cs->M2/((
long double)cs->area-1.0));
2223 cs->standard_deviation = sqrtl(cs->M2/((
long double)cs->area));
2225 cs->variance = cs->standard_deviation * cs->standard_deviation;
2226 cs->skewness = sqrtl(cs->area) * cs->M3 / powl (cs->M2*AdjArea,1.5);
2227 cs->kurtosis = cs->area * cs->M4 / (cs->M2*cs->M2*AdjArea*AdjArea) - 3.0;
2231 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2239 PixelChannel channel = GetPixelChannelChannel(image,i);
2242 cs->sum /= cs->area;
2243 cs->sum_squared /= cs->area;
2244 cs->sum_cubed /= cs->area;
2245 cs->sum_fourth_power /= cs->area;
2252 for (j=0; j <= (ssize_t) MaxMap; j++)
2253 if (histogram[(ssize_t) GetPixelChannels(image)*j+i] > 0.0)
2255 area=PerceptibleReciprocalLD(channel_statistics[channel].area);
2256 for (j=0; j <= (ssize_t) MaxMap; j++)
2261 count=area*histogram[(ssize_t) GetPixelChannels(image)*j+i];
2262 channel_statistics[channel].entropy+=((
long double) -count*
2263 MagickLog10(count)*PerceptibleReciprocalLD((
long double)
2264 MagickLog10(number_bins)));
2265 channel_statistics[CompositePixelChannel].entropy+=((
long double) -count*
2266 MagickLog10(count)*PerceptibleReciprocalLD((
long double)
2267 MagickLog10(number_bins))/GetPixelChannels(image));
2270 histogram=(
double *) RelinquishMagickMemory(histogram);
2272 median_info=AcquireVirtualMemory(image->columns,image->rows*
sizeof(*median));
2274 (void) ThrowMagickException(exception,GetMagickModule(),
2275 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
2278 median=(Quantum *) GetVirtualMemoryBlob(median_info);
2279 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2287 PixelChannel channel = GetPixelChannelChannel(image,i);
2288 PixelTrait traits = GetPixelChannelTraits(image,channel);
2289 if (traits == UndefinedPixelTrait)
2291 if ((traits & UpdatePixelTrait) == 0)
2293 for (y=0; y < (ssize_t) image->rows; y++)
2301 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
2302 if (p == (
const Quantum *) NULL)
2304 for (x=0; x < (ssize_t) image->columns; x++)
2306 if (GetPixelReadMask(image,p) <= (QuantumRange/2))
2308 p+=GetPixelChannels(image);
2312 p+=GetPixelChannels(image);
2315 channel_statistics[channel].median=(double) median[
2316 GetMedianPixel(median,n)];
2318 median_info=RelinquishVirtualMemory(median_info);
2323 csComp->sum_squared=0.0;
2324 csComp->sum_cubed=0.0;
2325 csComp->sum_fourth_power=0.0;
2326 csComp->maxima=(-MagickMaximumValue);
2327 csComp->minima=MagickMaximumValue;
2331 csComp->variance=0.0;
2332 csComp->standard_deviation=0.0;
2333 csComp->entropy=0.0;
2334 csComp->skewness=0.0;
2335 csComp->kurtosis=0.0;
2336 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2341 PixelChannel channel = GetPixelChannelChannel(image,i);
2342 PixelTrait traits = GetPixelChannelTraits(image,channel);
2343 if (traits == UndefinedPixelTrait)
2345 if ((traits & UpdatePixelTrait) == 0)
2347 cs=channel_statistics+i;
2348 if (csComp->maxima < cs->maxima)
2349 csComp->maxima = cs->maxima;
2350 if (csComp->minima > cs->minima)
2351 csComp->minima = cs->minima;
2352 csComp->sum += cs->sum;
2353 csComp->sum_squared += cs->sum_squared;
2354 csComp->sum_cubed += cs->sum_cubed;
2355 csComp->sum_fourth_power += cs->sum_fourth_power;
2356 csComp->median += cs->median;
2357 csComp->area += cs->area;
2358 csComp->mean += cs->mean;
2359 csComp->variance += cs->variance;
2360 csComp->standard_deviation += cs->standard_deviation;
2361 csComp->skewness += cs->skewness;
2362 csComp->kurtosis += cs->kurtosis;
2363 csComp->entropy += cs->entropy;
2365 channels=(double) GetImageChannels(image);
2366 csComp->sum /= channels;
2367 csComp->sum_squared /= channels;
2368 csComp->sum_cubed /= channels;
2369 csComp->sum_fourth_power /= channels;
2370 csComp->median /= channels;
2371 csComp->area /= channels;
2372 csComp->mean /= channels;
2373 csComp->variance /= channels;
2374 csComp->standard_deviation /= channels;
2375 csComp->skewness /= channels;
2376 csComp->kurtosis /= channels;
2377 csComp->entropy /= channels;
2379 if (y < (ssize_t) image->rows)
2381 channel_statistics);
2382 return(channel_statistics);
2418MagickExport
Image *PolynomialImage(
const Image *images,
2419 const size_t number_terms,
const double *terms,
ExceptionInfo *exception)
2421#define PolynomialImageTag "Polynomial/Image"
2436 **magick_restrict polynomial_pixels;
2444 assert(images != (
Image *) NULL);
2445 assert(images->signature == MagickCoreSignature);
2446 if (IsEventLogging() != MagickFalse)
2447 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",images->filename);
2449 assert(exception->signature == MagickCoreSignature);
2450 image=AcquireImageCanvas(images,exception);
2451 if (image == (
Image *) NULL)
2452 return((
Image *) NULL);
2453 if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
2455 image=DestroyImage(image);
2456 return((
Image *) NULL);
2458 number_images=GetImageListLength(images);
2459 polynomial_pixels=AcquirePixelTLS(images);
2462 image=DestroyImage(image);
2463 (void) ThrowMagickException(exception,GetMagickModule(),
2464 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",images->filename);
2465 return((
Image *) NULL);
2472 polynomial_view=AcquireAuthenticCacheView(image,exception);
2473#if defined(MAGICKCORE_OPENMP_SUPPORT)
2474 #pragma omp parallel for schedule(static) shared(progress,status) \
2475 magick_number_threads(image,image,image->rows,1)
2477 for (y=0; y < (ssize_t) image->rows; y++)
2486 id = GetOpenMPThreadId();
2499 if (status == MagickFalse)
2501 q=QueueCacheViewAuthenticPixels(polynomial_view,0,y,image->columns,1,
2503 if (q == (Quantum *) NULL)
2508 polynomial_pixel=polynomial_pixels[id];
2509 for (j=0; j < (ssize_t) image->columns; j++)
2510 for (i=0; i < MaxPixelChannels; i++)
2511 polynomial_pixel[j].channel[i]=0.0;
2513 for (j=0; j < (ssize_t) number_images; j++)
2518 if (j >= (ssize_t) number_terms)
2520 image_view=AcquireVirtualCacheView(next,exception);
2521 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2522 if (p == (
const Quantum *) NULL)
2524 image_view=DestroyCacheView(image_view);
2527 for (x=0; x < (ssize_t) image->columns; x++)
2529 for (i=0; i < (ssize_t) GetPixelChannels(next); i++)
2535 PixelChannel channel = GetPixelChannelChannel(image,i);
2536 PixelTrait traits = GetPixelChannelTraits(next,channel);
2537 PixelTrait polynomial_traits=GetPixelChannelTraits(image,channel);
2538 if ((traits == UndefinedPixelTrait) ||
2539 (polynomial_traits == UndefinedPixelTrait))
2541 if ((traits & UpdatePixelTrait) == 0)
2543 coefficient=(MagickRealType) terms[2*j];
2544 degree=(MagickRealType) terms[(j << 1)+1];
2545 polynomial_pixel[x].channel[i]+=coefficient*
2546 pow(QuantumScale*(
double) GetPixelChannel(image,channel,p),degree);
2548 p+=GetPixelChannels(next);
2550 image_view=DestroyCacheView(image_view);
2551 next=GetNextImageInList(next);
2553 for (x=0; x < (ssize_t) image->columns; x++)
2555 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2557 PixelChannel channel = GetPixelChannelChannel(image,i);
2558 PixelTrait traits = GetPixelChannelTraits(image,channel);
2559 if (traits == UndefinedPixelTrait)
2561 if ((traits & UpdatePixelTrait) == 0)
2563 q[i]=ClampToQuantum((
double) QuantumRange*
2564 polynomial_pixel[x].channel[i]);
2566 q+=GetPixelChannels(image);
2568 if (SyncCacheViewAuthenticPixels(polynomial_view,exception) == MagickFalse)
2570 if (images->progress_monitor != (MagickProgressMonitor) NULL)
2575#if defined(MAGICKCORE_OPENMP_SUPPORT)
2579 proceed=SetImageProgress(images,PolynomialImageTag,progress,
2581 if (proceed == MagickFalse)
2585 polynomial_view=DestroyCacheView(polynomial_view);
2586 polynomial_pixels=DestroyPixelTLS(images,polynomial_pixels);
2587 if (status == MagickFalse)
2588 image=DestroyImage(image);
2659 if (pixel_list->skip_list.nodes != (
SkipNode *) NULL)
2660 pixel_list->skip_list.nodes=(
SkipNode *) RelinquishAlignedMemory(
2661 pixel_list->skip_list.nodes);
2662 pixel_list=(
PixelList *) RelinquishMagickMemory(pixel_list);
2671 assert(pixel_list != (
PixelList **) NULL);
2672 for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
2673 if (pixel_list[i] != (
PixelList *) NULL)
2674 pixel_list[i]=DestroyPixelList(pixel_list[i]);
2675 pixel_list=(
PixelList **) RelinquishMagickMemory(pixel_list);
2679static PixelList *AcquirePixelList(
const size_t width,
const size_t height)
2684 pixel_list=(
PixelList *) AcquireMagickMemory(
sizeof(*pixel_list));
2687 (void) memset((
void *) pixel_list,0,
sizeof(*pixel_list));
2688 pixel_list->length=width*height;
2689 pixel_list->skip_list.nodes=(
SkipNode *) AcquireAlignedMemory(65537UL,
2690 sizeof(*pixel_list->skip_list.nodes));
2691 if (pixel_list->skip_list.nodes == (
SkipNode *) NULL)
2692 return(DestroyPixelList(pixel_list));
2693 (void) memset(pixel_list->skip_list.nodes,0,65537UL*
2694 sizeof(*pixel_list->skip_list.nodes));
2695 pixel_list->signature=MagickCoreSignature;
2699static PixelList **AcquirePixelListTLS(
const size_t width,
2700 const size_t height)
2711 number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
2712 pixel_list=(
PixelList **) AcquireQuantumMemory(number_threads,
2713 sizeof(*pixel_list));
2716 (void) memset(pixel_list,0,number_threads*
sizeof(*pixel_list));
2717 for (i=0; i < (ssize_t) number_threads; i++)
2719 pixel_list[i]=AcquirePixelList(width,height);
2720 if (pixel_list[i] == (
PixelList *) NULL)
2721 return(DestroyPixelListTLS(pixel_list));
2726static void AddNodePixelList(
PixelList *pixel_list,
const size_t color)
2741 p=(&pixel_list->skip_list);
2742 p->nodes[color].signature=pixel_list->signature;
2743 p->nodes[color].count=1;
2748 (void) memset(update,0,
sizeof(update));
2749 for (level=p->level; level >= 0; level--)
2751 while (p->nodes[search].next[level] < color)
2752 search=p->nodes[search].next[level];
2753 update[level]=search;
2758 for (level=0; ; level++)
2760 pixel_list->seed=(pixel_list->seed*42893621L)+1L;
2761 if ((pixel_list->seed & 0x300) != 0x300)
2766 if (level > (p->level+2))
2771 while (level > p->level)
2774 update[p->level]=65536UL;
2781 p->nodes[color].next[level]=p->nodes[update[level]].next[level];
2782 p->nodes[update[level]].next[level]=color;
2783 }
while (level-- > 0);
2786static inline void GetMedianPixelList(
PixelList *pixel_list,Quantum *pixel)
2800 p=(&pixel_list->skip_list);
2805 color=p->nodes[color].next[0];
2806 count+=(ssize_t) p->nodes[color].count;
2807 }
while (count <= (ssize_t) (pixel_list->length >> 1));
2808 *pixel=ScaleShortToQuantum((
unsigned short) color);
2811static inline void GetModePixelList(
PixelList *pixel_list,Quantum *pixel)
2827 p=(&pixel_list->skip_list);
2830 max_count=p->nodes[mode].count;
2834 color=p->nodes[color].next[0];
2835 if (p->nodes[color].count > max_count)
2838 max_count=p->nodes[mode].count;
2840 count+=(ssize_t) p->nodes[color].count;
2841 }
while (count < (ssize_t) pixel_list->length);
2842 *pixel=ScaleShortToQuantum((
unsigned short) mode);
2845static inline void GetNonpeakPixelList(
PixelList *pixel_list,Quantum *pixel)
2861 p=(&pixel_list->skip_list);
2863 next=p->nodes[color].next[0];
2869 next=p->nodes[color].next[0];
2870 count+=(ssize_t) p->nodes[color].count;
2871 }
while (count <= (ssize_t) (pixel_list->length >> 1));
2872 if ((previous == 65536UL) && (next != 65536UL))
2875 if ((previous != 65536UL) && (next == 65536UL))
2877 *pixel=ScaleShortToQuantum((
unsigned short) color);
2880static inline void InsertPixelList(
const Quantum pixel,
PixelList *pixel_list)
2888 index=ScaleQuantumToShort(pixel);
2889 signature=pixel_list->skip_list.nodes[index].signature;
2890 if (signature == pixel_list->signature)
2892 pixel_list->skip_list.nodes[index].count++;
2895 AddNodePixelList(pixel_list,index);
2898static void ResetPixelList(
PixelList *pixel_list)
2912 p=(&pixel_list->skip_list);
2913 root=p->nodes+65536UL;
2915 for (level=0; level < 9; level++)
2916 root->next[level]=65536UL;
2917 pixel_list->seed=pixel_list->signature++;
2920MagickExport
Image *StatisticImage(
const Image *image,
const StatisticType type,
2921 const size_t width,
const size_t height,
ExceptionInfo *exception)
2923#define StatisticImageTag "Statistic/Image"
2939 **magick_restrict pixel_list;
2948 assert(image != (
Image *) NULL);
2949 assert(image->signature == MagickCoreSignature);
2950 if (IsEventLogging() != MagickFalse)
2951 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
2953 assert(exception->signature == MagickCoreSignature);
2954 statistic_image=CloneImage(image,0,0,MagickTrue,
2956 if (statistic_image == (
Image *) NULL)
2957 return((
Image *) NULL);
2958 status=SetImageStorageClass(statistic_image,DirectClass,exception);
2959 if (status == MagickFalse)
2961 statistic_image=DestroyImage(statistic_image);
2962 return((
Image *) NULL);
2964 pixel_list=AcquirePixelListTLS(MagickMax(width,1),MagickMax(height,1));
2967 statistic_image=DestroyImage(statistic_image);
2968 ThrowImageException(ResourceLimitError,
"MemoryAllocationFailed");
2973 center=(ssize_t) GetPixelChannels(image)*((ssize_t) image->columns+
2974 MagickMax((ssize_t) width,1L))*(MagickMax((ssize_t) height,1)/2L)+(ssize_t)
2975 GetPixelChannels(image)*(MagickMax((ssize_t) width,1L)/2L);
2978 image_view=AcquireVirtualCacheView(image,exception);
2979 statistic_view=AcquireAuthenticCacheView(statistic_image,exception);
2980#if defined(MAGICKCORE_OPENMP_SUPPORT)
2981 #pragma omp parallel for schedule(static) shared(progress,status) \
2982 magick_number_threads(image,statistic_image,statistic_image->rows,1)
2984 for (y=0; y < (ssize_t) statistic_image->rows; y++)
2987 id = GetOpenMPThreadId();
2998 if (status == MagickFalse)
3000 p=GetCacheViewVirtualPixels(image_view,-((ssize_t) MagickMax(width,1)/2L),y-
3001 (ssize_t) (MagickMax(height,1)/2L),image->columns+MagickMax(width,1),
3002 MagickMax(height,1),exception);
3003 q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns, 1,exception);
3004 if ((p == (
const Quantum *) NULL) || (q == (Quantum *) NULL))
3009 for (x=0; x < (ssize_t) statistic_image->columns; x++)
3014 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3027 *magick_restrict pixels;
3035 PixelChannel channel = GetPixelChannelChannel(image,i);
3036 PixelTrait traits = GetPixelChannelTraits(image,channel);
3037 PixelTrait statistic_traits=GetPixelChannelTraits(statistic_image,
3039 if ((traits == UndefinedPixelTrait) ||
3040 (statistic_traits == UndefinedPixelTrait))
3042 if (((statistic_traits & CopyPixelTrait) != 0) ||
3043 (GetPixelWriteMask(image,p) <= (QuantumRange/2)))
3045 SetPixelChannel(statistic_image,channel,p[center+i],q);
3048 if ((statistic_traits & UpdatePixelTrait) == 0)
3056 ResetPixelList(pixel_list[
id]);
3057 for (v=0; v < (ssize_t) MagickMax(height,1); v++)
3059 for (u=0; u < (ssize_t) MagickMax(width,1); u++)
3061 if ((type == MedianStatistic) || (type == ModeStatistic) ||
3062 (type == NonpeakStatistic))
3064 InsertPixelList(pixels[i],pixel_list[
id]);
3065 pixels+=GetPixelChannels(image);
3069 if ((
double) pixels[i] < minimum)
3070 minimum=(double) pixels[i];
3071 if ((
double) pixels[i] > maximum)
3072 maximum=(
double) pixels[i];
3073 sum+=(double) pixels[i];
3074 sum_squared+=(double) pixels[i]*(
double) pixels[i];
3075 pixels+=GetPixelChannels(image);
3077 pixels+=GetPixelChannels(image)*image->columns;
3081 case ContrastStatistic:
3083 pixel=ClampToQuantum(MagickAbsoluteValue((maximum-minimum)*
3084 PerceptibleReciprocal(maximum+minimum)));
3087 case GradientStatistic:
3089 pixel=ClampToQuantum(MagickAbsoluteValue(maximum-minimum));
3092 case MaximumStatistic:
3094 pixel=ClampToQuantum(maximum);
3100 pixel=ClampToQuantum(sum/area);
3103 case MedianStatistic:
3105 GetMedianPixelList(pixel_list[
id],&pixel);
3108 case MinimumStatistic:
3110 pixel=ClampToQuantum(minimum);
3115 GetModePixelList(pixel_list[
id],&pixel);
3118 case NonpeakStatistic:
3120 GetNonpeakPixelList(pixel_list[
id],&pixel);
3123 case RootMeanSquareStatistic:
3125 pixel=ClampToQuantum(sqrt(sum_squared/area));
3128 case StandardDeviationStatistic:
3130 pixel=ClampToQuantum(sqrt(sum_squared/area-(sum/area*sum/area)));
3134 SetPixelChannel(statistic_image,channel,pixel,q);
3136 p+=GetPixelChannels(image);
3137 q+=GetPixelChannels(statistic_image);
3139 if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
3141 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3146#if defined(MAGICKCORE_OPENMP_SUPPORT)
3150 proceed=SetImageProgress(image,StatisticImageTag,progress,image->rows);
3151 if (proceed == MagickFalse)
3155 statistic_view=DestroyCacheView(statistic_view);
3156 image_view=DestroyCacheView(image_view);
3157 pixel_list=DestroyPixelListTLS(pixel_list);
3158 if (status == MagickFalse)
3159 statistic_image=DestroyImage(statistic_image);
3160 return(statistic_image);