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/string_.h"
92 #include "MagickCore/thread-private.h"
93 #include "MagickCore/timer.h"
94 #include "MagickCore/utility.h"
95 #include "MagickCore/version.h"
137 channel[MaxPixelChannels];
150 rows=MagickMax(GetImageListLength(images),(
size_t)
151 GetMagickResourceLimit(ThreadResource));
152 for (i=0; i < (ssize_t) rows; i++)
154 pixels[i]=(
PixelChannels *) RelinquishMagickMemory(pixels[i]);
175 number_images=GetImageListLength(images);
176 rows=MagickMax(number_images,(
size_t) GetMagickResourceLimit(ThreadResource));
177 pixels=(
PixelChannels **) AcquireQuantumMemory(rows,
sizeof(*pixels));
180 (void) memset(pixels,0,rows*
sizeof(*pixels));
181 columns=MagickMax(number_images,MaxPixelChannels);
182 for (next=images; next != (
Image *) NULL; next=next->next)
183 columns=MagickMax(next->columns,columns);
184 for (i=0; i < (ssize_t) rows; i++)
189 pixels[i]=(
PixelChannels *) AcquireQuantumMemory(columns,
sizeof(**pixels));
191 return(DestroyPixelTLS(images,pixels));
192 for (j=0; j < (ssize_t) columns; j++)
197 for (k=0; k < MaxPixelChannels; k++)
198 pixels[i][j].channel[k]=0.0;
204 static inline double EvaluateMax(
const double x,
const double y)
211 #if defined(__cplusplus) || defined(c_plusplus)
215 static int IntensityCompare(
const void *x,
const void *y)
230 for (i=0; i < MaxPixelChannels; i++)
231 distance+=color_1->channel[i]-(
double) color_2->channel[i];
232 return(distance < 0.0 ? -1 : distance > 0.0 ? 1 : 0);
235 #if defined(__cplusplus) || defined(c_plusplus)
239 static double ApplyEvaluateOperator(
RandomInfo *random_info,
const Quantum pixel,
240 const MagickEvaluateOperator op,
const double value)
251 case UndefinedEvaluateOperator:
253 case AbsEvaluateOperator:
255 result=(double) fabs((
double) (pixel+value));
258 case AddEvaluateOperator:
260 result=(double) (pixel+value);
263 case AddModulusEvaluateOperator:
272 result-=(QuantumRange+1.0)*floor((
double) result/(QuantumRange+1.0));
275 case AndEvaluateOperator:
277 result=(double) ((ssize_t) pixel & (ssize_t) (value+0.5));
280 case CosineEvaluateOperator:
282 result=(double) (QuantumRange*(0.5*cos((
double) (2.0*MagickPI*
283 QuantumScale*pixel*value))+0.5));
286 case DivideEvaluateOperator:
288 result=pixel/(value == 0.0 ? 1.0 : value);
291 case ExponentialEvaluateOperator:
293 result=(double) (QuantumRange*exp((
double) (value*QuantumScale*pixel)));
296 case GaussianNoiseEvaluateOperator:
298 result=(double) GenerateDifferentialNoise(random_info,pixel,GaussianNoise,
302 case ImpulseNoiseEvaluateOperator:
304 result=(double) GenerateDifferentialNoise(random_info,pixel,ImpulseNoise,
308 case InverseLogEvaluateOperator:
310 result=(QuantumRange*pow((value+1.0),QuantumScale*pixel)-1.0)*
311 PerceptibleReciprocal(value);
314 case LaplacianNoiseEvaluateOperator:
316 result=(double) GenerateDifferentialNoise(random_info,pixel,
317 LaplacianNoise,value);
320 case LeftShiftEvaluateOperator:
322 result=(double) pixel;
323 for (i=0; i < (ssize_t) value; i++)
327 case LogEvaluateOperator:
329 if ((QuantumScale*pixel) >= MagickEpsilon)
330 result=(
double) (QuantumRange*log((
double) (QuantumScale*value*pixel+
331 1.0))/log((
double) (value+1.0)));
334 case MaxEvaluateOperator:
336 result=(double) EvaluateMax((
double) pixel,value);
339 case MeanEvaluateOperator:
341 result=(double) (pixel+value);
344 case MedianEvaluateOperator:
346 result=(double) (pixel+value);
349 case MinEvaluateOperator:
351 result=(double) MagickMin((
double) pixel,value);
354 case MultiplicativeNoiseEvaluateOperator:
356 result=(double) GenerateDifferentialNoise(random_info,pixel,
357 MultiplicativeGaussianNoise,value);
360 case MultiplyEvaluateOperator:
362 result=(double) (value*pixel);
365 case OrEvaluateOperator:
367 result=(double) ((ssize_t) pixel | (ssize_t) (value+0.5));
370 case PoissonNoiseEvaluateOperator:
372 result=(double) GenerateDifferentialNoise(random_info,pixel,PoissonNoise,
376 case PowEvaluateOperator:
379 result=(double) -(QuantumRange*pow((
double) -(QuantumScale*pixel),
382 result=(double) (QuantumRange*pow((
double) (QuantumScale*pixel),
386 case RightShiftEvaluateOperator:
388 result=(double) pixel;
389 for (i=0; i < (ssize_t) value; i++)
393 case RootMeanSquareEvaluateOperator:
395 result=((double) pixel*pixel+value);
398 case SetEvaluateOperator:
403 case SineEvaluateOperator:
405 result=(double) (QuantumRange*(0.5*sin((
double) (2.0*MagickPI*
406 QuantumScale*pixel*value))+0.5));
409 case SubtractEvaluateOperator:
411 result=(double) (pixel-value);
414 case SumEvaluateOperator:
416 result=(double) (pixel+value);
419 case ThresholdEvaluateOperator:
421 result=(double) (((
double) pixel <= value) ? 0 : QuantumRange);
424 case ThresholdBlackEvaluateOperator:
426 result=(double) (((
double) pixel <= value) ? 0 : pixel);
429 case ThresholdWhiteEvaluateOperator:
431 result=(double) (((
double) pixel > value) ? QuantumRange : pixel);
434 case UniformNoiseEvaluateOperator:
436 result=(double) GenerateDifferentialNoise(random_info,pixel,UniformNoise,
440 case XorEvaluateOperator:
442 result=(double) ((ssize_t) pixel ^ (ssize_t) (value+0.5));
460 columns=images->columns;
462 for (p=images; p != (
Image *) NULL; p=p->next)
464 if (p->number_channels > q->number_channels)
466 if (p->columns > columns)
471 return(CloneImage(q,columns,rows,MagickTrue,exception));
474 MagickExport
Image *EvaluateImages(
const Image *images,
477 #define EvaluateImageTag "Evaluate/Image"
496 **magick_restrict evaluate_pixels;
499 **magick_restrict random_info;
508 #if defined(MAGICKCORE_OPENMP_SUPPORT)
513 assert(images != (
Image *) NULL);
514 assert(images->signature == MagickCoreSignature);
516 assert(exception->signature == MagickCoreSignature);
517 if (IsEventLogging() != MagickFalse)
518 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",images->filename);
519 image=AcquireImageCanvas(images,exception);
520 if (image == (
Image *) NULL)
521 return((
Image *) NULL);
522 if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
524 image=DestroyImage(image);
525 return((
Image *) NULL);
527 number_images=GetImageListLength(images);
528 evaluate_pixels=AcquirePixelTLS(images);
531 image=DestroyImage(image);
532 (void) ThrowMagickException(exception,GetMagickModule(),
533 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",images->filename);
534 return((
Image *) NULL);
536 image_view=(
CacheView **) AcquireQuantumMemory(number_images,
537 sizeof(*image_view));
540 image=DestroyImage(image);
541 evaluate_pixels=DestroyPixelTLS(images,evaluate_pixels);
542 (void) ThrowMagickException(exception,GetMagickModule(),
543 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",images->filename);
547 for (n=0; n < (ssize_t) number_images; n++)
549 image_view[n]=AcquireVirtualCacheView(view,exception);
550 view=GetNextImageInList(view);
557 random_info=AcquireRandomInfoTLS();
558 evaluate_view=AcquireAuthenticCacheView(image,exception);
559 if (op == MedianEvaluateOperator)
561 #if defined(MAGICKCORE_OPENMP_SUPPORT)
562 key=GetRandomSecretKey(random_info[0]);
563 #pragma omp parallel for schedule(static) shared(progress,status) \
564 magick_number_threads(image,images,image->rows,key == ~0UL)
566 for (y=0; y < (ssize_t) image->rows; y++)
569 id = GetOpenMPThreadId();
586 if (status == MagickFalse)
588 p=(
const Quantum **) AcquireQuantumMemory(number_images,
sizeof(*p));
589 if (p == (
const Quantum **) NULL)
592 (void) ThrowMagickException(exception,GetMagickModule(),
593 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",
597 for (j=0; j < (ssize_t) number_images; j++)
599 p[j]=GetCacheViewVirtualPixels(image_view[j],0,y,image->columns,1,
601 if (p[j] == (
const Quantum *) NULL)
604 q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,image->columns,1,
606 if ((j < (ssize_t) number_images) || (q == (Quantum *) NULL))
611 evaluate_pixel=evaluate_pixels[id];
612 for (x=0; x < (ssize_t) image->columns; x++)
621 for (j=0; j < (ssize_t) number_images; j++)
623 for (i=0; i < MaxPixelChannels; i++)
624 evaluate_pixel[j].channel[i]=0.0;
625 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
627 PixelChannel channel = GetPixelChannelChannel(image,i);
628 PixelTrait traits = GetPixelChannelTraits(next,channel);
629 PixelTrait evaluate_traits = GetPixelChannelTraits(image,channel);
630 if ((traits == UndefinedPixelTrait) ||
631 (evaluate_traits == UndefinedPixelTrait) ||
632 ((traits & UpdatePixelTrait) == 0))
634 evaluate_pixel[j].channel[i]=ApplyEvaluateOperator(
635 random_info[
id],GetPixelChannel(next,channel,p[j]),op,
636 evaluate_pixel[j].channel[i]);
638 p[j]+=GetPixelChannels(next);
639 next=GetNextImageInList(next);
641 qsort((
void *) evaluate_pixel,number_images,
sizeof(*evaluate_pixel),
643 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
645 PixelChannel channel = GetPixelChannelChannel(image,i);
646 PixelTrait traits = GetPixelChannelTraits(image,channel);
647 if ((traits == UndefinedPixelTrait) ||
648 ((traits & UpdatePixelTrait) == 0))
650 q[i]=ClampToQuantum(evaluate_pixel[number_images/2].channel[i]);
652 q+=GetPixelChannels(image);
654 p=(
const Quantum **) RelinquishMagickMemory((
void *) p);
655 if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
657 if (images->progress_monitor != (MagickProgressMonitor) NULL)
662 #if defined(MAGICKCORE_OPENMP_SUPPORT)
666 proceed=SetImageProgress(images,EvaluateImageTag,progress,
668 if (proceed == MagickFalse)
675 #if defined(MAGICKCORE_OPENMP_SUPPORT)
676 key=GetRandomSecretKey(random_info[0]);
677 #pragma omp parallel for schedule(static) shared(progress,status) \
678 magick_number_threads(image,images,image->rows,key == ~0UL)
680 for (y=0; y < (ssize_t) image->rows; y++)
686 id = GetOpenMPThreadId();
704 if (status == MagickFalse)
706 p=(
const Quantum **) AcquireQuantumMemory(number_images,
sizeof(*p));
707 if (p == (
const Quantum **) NULL)
710 (void) ThrowMagickException(exception,GetMagickModule(),
711 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",
715 for (j=0; j < (ssize_t) number_images; j++)
717 p[j]=GetCacheViewVirtualPixels(image_view[j],0,y,image->columns,1,
719 if (p[j] == (
const Quantum *) NULL)
722 q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,image->columns,1,
724 if ((j < (ssize_t) number_images) || (q == (Quantum *) NULL))
729 evaluate_pixel=evaluate_pixels[id];
730 for (j=0; j < (ssize_t) image->columns; j++)
731 for (i=0; i < MaxPixelChannels; i++)
732 evaluate_pixel[j].channel[i]=0.0;
734 for (j=0; j < (ssize_t) number_images; j++)
736 for (x=0; x < (ssize_t) image->columns; x++)
738 for (i=0; i < (ssize_t) GetPixelChannels(next); i++)
740 PixelChannel channel = GetPixelChannelChannel(image,i);
741 PixelTrait traits = GetPixelChannelTraits(next,channel);
742 PixelTrait evaluate_traits = GetPixelChannelTraits(image,channel);
743 if ((traits == UndefinedPixelTrait) ||
744 (evaluate_traits == UndefinedPixelTrait))
746 if ((traits & UpdatePixelTrait) == 0)
748 evaluate_pixel[x].channel[i]=ApplyEvaluateOperator(
749 random_info[
id],GetPixelChannel(next,channel,p[j]),j == 0 ?
750 AddEvaluateOperator : op,evaluate_pixel[x].channel[i]);
752 p[j]+=GetPixelChannels(next);
754 next=GetNextImageInList(next);
756 for (x=0; x < (ssize_t) image->columns; x++)
760 case MeanEvaluateOperator:
762 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
763 evaluate_pixel[x].channel[i]/=(
double) number_images;
766 case MultiplyEvaluateOperator:
768 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
770 for (j=0; j < (ssize_t) (number_images-1); j++)
771 evaluate_pixel[x].channel[i]*=QuantumScale;
775 case RootMeanSquareEvaluateOperator:
777 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
778 evaluate_pixel[x].channel[i]=sqrt(evaluate_pixel[x].channel[i]/
786 for (x=0; x < (ssize_t) image->columns; x++)
788 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
790 PixelChannel channel = GetPixelChannelChannel(image,i);
791 PixelTrait traits = GetPixelChannelTraits(image,channel);
792 if ((traits == UndefinedPixelTrait) ||
793 ((traits & UpdatePixelTrait) == 0))
795 q[i]=ClampToQuantum(evaluate_pixel[x].channel[i]);
797 q+=GetPixelChannels(image);
799 p=(
const Quantum **) RelinquishMagickMemory((
void *) p);
800 if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse)
802 if (images->progress_monitor != (MagickProgressMonitor) NULL)
807 #if defined(MAGICKCORE_OPENMP_SUPPORT)
811 proceed=SetImageProgress(images,EvaluateImageTag,progress,
813 if (proceed == MagickFalse)
818 for (n=0; n < (ssize_t) number_images; n++)
819 image_view[n]=DestroyCacheView(image_view[n]);
820 image_view=(
CacheView **) RelinquishMagickMemory(image_view);
821 evaluate_view=DestroyCacheView(evaluate_view);
822 evaluate_pixels=DestroyPixelTLS(images,evaluate_pixels);
823 random_info=DestroyRandomInfoTLS(random_info);
824 if (status == MagickFalse)
825 image=DestroyImage(image);
829 MagickExport MagickBooleanType EvaluateImage(
Image *image,
830 const MagickEvaluateOperator op,
const double value,
ExceptionInfo *exception)
846 **magick_restrict random_info;
851 #if defined(MAGICKCORE_OPENMP_SUPPORT)
856 assert(image != (
Image *) NULL);
857 assert(image->signature == MagickCoreSignature);
859 assert(exception->signature == MagickCoreSignature);
860 if (IsEventLogging() != MagickFalse)
861 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
862 if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
867 artifact=GetImageArtifact(image,
"evaluate:clamp");
868 if (artifact != (
const char *) NULL)
869 clamp=IsStringTrue(artifact);
870 random_info=AcquireRandomInfoTLS();
871 image_view=AcquireAuthenticCacheView(image,exception);
872 #if defined(MAGICKCORE_OPENMP_SUPPORT)
873 key=GetRandomSecretKey(random_info[0]);
874 #pragma omp parallel for schedule(static) shared(progress,status) \
875 magick_number_threads(image,image,image->rows,key == ~0UL)
877 for (y=0; y < (ssize_t) image->rows; y++)
880 id = GetOpenMPThreadId();
888 if (status == MagickFalse)
890 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
891 if (q == (Quantum *) NULL)
896 for (x=0; x < (ssize_t) image->columns; x++)
904 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
906 PixelChannel channel = GetPixelChannelChannel(image,i);
907 PixelTrait traits = GetPixelChannelTraits(image,channel);
908 if (traits == UndefinedPixelTrait)
910 if ((traits & CopyPixelTrait) != 0)
912 if ((traits & UpdatePixelTrait) == 0)
914 result=ApplyEvaluateOperator(random_info[
id],q[i],op,value);
915 if (op == MeanEvaluateOperator)
917 q[i]=clamp != MagickFalse ? ClampPixel(result) : ClampToQuantum(result);
919 q+=GetPixelChannels(image);
921 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
923 if (image->progress_monitor != (MagickProgressMonitor) NULL)
928 #if defined(MAGICKCORE_OPENMP_SUPPORT)
932 proceed=SetImageProgress(image,EvaluateImageTag,progress,image->rows);
933 if (proceed == MagickFalse)
937 image_view=DestroyCacheView(image_view);
938 random_info=DestroyRandomInfoTLS(random_info);
976 static Quantum ApplyFunction(Quantum pixel,
const MagickFunction
function,
977 const size_t number_parameters,
const double *parameters,
990 case PolynomialFunction:
997 for (i=0; i < (ssize_t) number_parameters; i++)
998 result=result*QuantumScale*pixel+parameters[i];
999 result*=QuantumRange;
1002 case SinusoidFunction:
1013 frequency=(number_parameters >= 1) ? parameters[0] : 1.0;
1014 phase=(number_parameters >= 2) ? parameters[1] : 0.0;
1015 amplitude=(number_parameters >= 3) ? parameters[2] : 0.5;
1016 bias=(number_parameters >= 4) ? parameters[3] : 0.5;
1017 result=(double) (QuantumRange*(amplitude*sin((
double) (2.0*
1018 MagickPI*(frequency*QuantumScale*pixel+phase/360.0)))+bias));
1021 case ArcsinFunction:
1033 width=(number_parameters >= 1) ? parameters[0] : 1.0;
1034 center=(number_parameters >= 2) ? parameters[1] : 0.5;
1035 range=(number_parameters >= 3) ? parameters[2] : 1.0;
1036 bias=(number_parameters >= 4) ? parameters[3] : 0.5;
1037 result=2.0*PerceptibleReciprocal(width)*(QuantumScale*pixel-center);
1039 result=bias-range/2.0;
1042 result=bias+range/2.0;
1044 result=(double) (range/MagickPI*asin((
double) result)+bias);
1045 result*=QuantumRange;
1048 case ArctanFunction:
1059 slope=(number_parameters >= 1) ? parameters[0] : 1.0;
1060 center=(number_parameters >= 2) ? parameters[1] : 0.5;
1061 range=(number_parameters >= 3) ? parameters[2] : 1.0;
1062 bias=(number_parameters >= 4) ? parameters[3] : 0.5;
1063 result=(double) (MagickPI*slope*(QuantumScale*pixel-center));
1064 result=(double) (QuantumRange*(range/MagickPI*atan((
double)
1068 case UndefinedFunction:
1071 return(ClampToQuantum(result));
1074 MagickExport MagickBooleanType FunctionImage(
Image *image,
1075 const MagickFunction
function,
const size_t number_parameters,
1078 #define FunctionImageTag "Function/Image "
1092 assert(image != (
Image *) NULL);
1093 assert(image->signature == MagickCoreSignature);
1095 assert(exception->signature == MagickCoreSignature);
1096 if (IsEventLogging() != MagickFalse)
1097 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
1098 #if defined(MAGICKCORE_OPENCL_SUPPORT)
1099 if (AccelerateFunctionImage(image,
function,number_parameters,parameters,
1100 exception) != MagickFalse)
1103 if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
1104 return(MagickFalse);
1107 image_view=AcquireAuthenticCacheView(image,exception);
1108 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1109 #pragma omp parallel for schedule(static) shared(progress,status) \
1110 magick_number_threads(image,image,image->rows,1)
1112 for (y=0; y < (ssize_t) image->rows; y++)
1120 if (status == MagickFalse)
1122 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
1123 if (q == (Quantum *) NULL)
1128 for (x=0; x < (ssize_t) image->columns; x++)
1133 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1135 PixelChannel channel = GetPixelChannelChannel(image,i);
1136 PixelTrait traits = GetPixelChannelTraits(image,channel);
1137 if (traits == UndefinedPixelTrait)
1139 if ((traits & UpdatePixelTrait) == 0)
1141 q[i]=ApplyFunction(q[i],
function,number_parameters,parameters,
1144 q+=GetPixelChannels(image);
1146 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1148 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1153 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1157 proceed=SetImageProgress(image,FunctionImageTag,progress,image->rows);
1158 if (proceed == MagickFalse)
1162 image_view=DestroyCacheView(image_view);
1193 MagickExport MagickBooleanType GetImageEntropy(
const Image *image,
1197 *channel_statistics;
1199 assert(image != (
Image *) NULL);
1200 assert(image->signature == MagickCoreSignature);
1201 if (IsEventLogging() != MagickFalse)
1202 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
1203 channel_statistics=GetImageStatistics(image,exception);
1205 return(MagickFalse);
1206 *entropy=channel_statistics[CompositePixelChannel].entropy;
1208 channel_statistics);
1241 MagickExport MagickBooleanType GetImageExtrema(
const Image *image,
1251 assert(image != (
Image *) NULL);
1252 assert(image->signature == MagickCoreSignature);
1253 if (IsEventLogging() != MagickFalse)
1254 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
1255 status=GetImageRange(image,&min,&max,exception);
1256 *minima=(size_t) ceil(min-0.5);
1257 *maxima=(size_t) floor(max+0.5);
1291 MagickExport MagickBooleanType GetImageKurtosis(
const Image *image,
1295 *channel_statistics;
1297 assert(image != (
Image *) NULL);
1298 assert(image->signature == MagickCoreSignature);
1299 if (IsEventLogging() != MagickFalse)
1300 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
1301 channel_statistics=GetImageStatistics(image,exception);
1303 return(MagickFalse);
1304 *kurtosis=channel_statistics[CompositePixelChannel].kurtosis;
1305 *skewness=channel_statistics[CompositePixelChannel].skewness;
1307 channel_statistics);
1341 MagickExport MagickBooleanType GetImageMean(
const Image *image,
double *mean,
1345 *channel_statistics;
1347 assert(image != (
Image *) NULL);
1348 assert(image->signature == MagickCoreSignature);
1349 if (IsEventLogging() != MagickFalse)
1350 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
1351 channel_statistics=GetImageStatistics(image,exception);
1353 return(MagickFalse);
1354 *mean=channel_statistics[CompositePixelChannel].mean;
1355 *standard_deviation=
1356 channel_statistics[CompositePixelChannel].standard_deviation;
1358 channel_statistics);
1389 MagickExport MagickBooleanType GetImageMedian(
const Image *image,
double *median,
1393 *channel_statistics;
1395 assert(image != (
Image *) NULL);
1396 assert(image->signature == MagickCoreSignature);
1397 if (IsEventLogging() != MagickFalse)
1398 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
1399 channel_statistics=GetImageStatistics(image,exception);
1401 return(MagickFalse);
1402 *median=channel_statistics[CompositePixelChannel].median;
1404 channel_statistics);
1437 #define MaxNumberImageMoments 8
1447 M00[2*MaxPixelChannels+1],
1448 M01[2*MaxPixelChannels+1],
1449 M02[2*MaxPixelChannels+1],
1450 M03[2*MaxPixelChannels+1],
1451 M10[2*MaxPixelChannels+1],
1452 M11[2*MaxPixelChannels+1],
1453 M12[2*MaxPixelChannels+1],
1454 M20[2*MaxPixelChannels+1],
1455 M21[2*MaxPixelChannels+1],
1456 M22[2*MaxPixelChannels+1],
1457 M30[2*MaxPixelChannels+1];
1460 centroid[2*MaxPixelChannels+1];
1466 assert(image != (
Image *) NULL);
1467 assert(image->signature == MagickCoreSignature);
1468 if (IsEventLogging() != MagickFalse)
1469 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
1470 channel_moments=(
ChannelMoments *) AcquireQuantumMemory(MaxPixelChannels+1,
1471 sizeof(*channel_moments));
1473 return(channel_moments);
1474 (void) memset(channel_moments,0,(MaxPixelChannels+1)*
1475 sizeof(*channel_moments));
1476 (void) memset(centroid,0,
sizeof(centroid));
1477 (void) memset(M00,0,
sizeof(M00));
1478 (void) memset(M01,0,
sizeof(M01));
1479 (void) memset(M02,0,
sizeof(M02));
1480 (void) memset(M03,0,
sizeof(M03));
1481 (void) memset(M10,0,
sizeof(M10));
1482 (void) memset(M11,0,
sizeof(M11));
1483 (void) memset(M12,0,
sizeof(M12));
1484 (void) memset(M20,0,
sizeof(M20));
1485 (void) memset(M21,0,
sizeof(M21));
1486 (void) memset(M22,0,
sizeof(M22));
1487 (void) memset(M30,0,
sizeof(M30));
1488 image_view=AcquireVirtualCacheView(image,exception);
1489 for (y=0; y < (ssize_t) image->rows; y++)
1500 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1501 if (p == (
const Quantum *) NULL)
1503 for (x=0; x < (ssize_t) image->columns; x++)
1508 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1510 PixelChannel channel = GetPixelChannelChannel(image,i);
1511 PixelTrait traits = GetPixelChannelTraits(image,channel);
1512 if (traits == UndefinedPixelTrait)
1514 if ((traits & UpdatePixelTrait) == 0)
1516 M00[channel]+=QuantumScale*p[i];
1517 M00[MaxPixelChannels]+=QuantumScale*p[i];
1518 M10[channel]+=x*QuantumScale*p[i];
1519 M10[MaxPixelChannels]+=x*QuantumScale*p[i];
1520 M01[channel]+=y*QuantumScale*p[i];
1521 M01[MaxPixelChannels]+=y*QuantumScale*p[i];
1523 p+=GetPixelChannels(image);
1526 for (c=0; c <= MaxPixelChannels; c++)
1531 centroid[c].x=M10[c]*PerceptibleReciprocal(M00[c]);
1532 centroid[c].y=M01[c]*PerceptibleReciprocal(M00[c]);
1534 for (y=0; y < (ssize_t) image->rows; y++)
1545 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1546 if (p == (
const Quantum *) NULL)
1548 for (x=0; x < (ssize_t) image->columns; x++)
1553 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1555 PixelChannel channel = GetPixelChannelChannel(image,i);
1556 PixelTrait traits = GetPixelChannelTraits(image,channel);
1557 if (traits == UndefinedPixelTrait)
1559 if ((traits & UpdatePixelTrait) == 0)
1561 M11[channel]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
1563 M11[MaxPixelChannels]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
1565 M20[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1567 M20[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1569 M02[channel]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
1571 M02[MaxPixelChannels]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
1573 M21[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1574 (y-centroid[channel].y)*QuantumScale*p[i];
1575 M21[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1576 (y-centroid[channel].y)*QuantumScale*p[i];
1577 M12[channel]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
1578 (y-centroid[channel].y)*QuantumScale*p[i];
1579 M12[MaxPixelChannels]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
1580 (y-centroid[channel].y)*QuantumScale*p[i];
1581 M22[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1582 (y-centroid[channel].y)*(y-centroid[channel].y)*QuantumScale*p[i];
1583 M22[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1584 (y-centroid[channel].y)*(y-centroid[channel].y)*QuantumScale*p[i];
1585 M30[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1586 (x-centroid[channel].x)*QuantumScale*p[i];
1587 M30[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
1588 (x-centroid[channel].x)*QuantumScale*p[i];
1589 M03[channel]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
1590 (y-centroid[channel].y)*QuantumScale*p[i];
1591 M03[MaxPixelChannels]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
1592 (y-centroid[channel].y)*QuantumScale*p[i];
1594 p+=GetPixelChannels(image);
1597 channels=(double) GetImageChannels(image);
1598 M00[MaxPixelChannels]/=channels;
1599 M01[MaxPixelChannels]/=channels;
1600 M02[MaxPixelChannels]/=channels;
1601 M03[MaxPixelChannels]/=channels;
1602 M10[MaxPixelChannels]/=channels;
1603 M11[MaxPixelChannels]/=channels;
1604 M12[MaxPixelChannels]/=channels;
1605 M20[MaxPixelChannels]/=channels;
1606 M21[MaxPixelChannels]/=channels;
1607 M22[MaxPixelChannels]/=channels;
1608 M30[MaxPixelChannels]/=channels;
1609 for (c=0; c <= MaxPixelChannels; c++)
1614 channel_moments[c].centroid=centroid[c];
1615 channel_moments[c].ellipse_axis.x=sqrt((2.0*PerceptibleReciprocal(M00[c]))*
1616 ((M20[c]+M02[c])+sqrt(4.0*M11[c]*M11[c]+(M20[c]-M02[c])*(M20[c]-M02[c]))));
1617 channel_moments[c].ellipse_axis.y=sqrt((2.0*PerceptibleReciprocal(M00[c]))*
1618 ((M20[c]+M02[c])-sqrt(4.0*M11[c]*M11[c]+(M20[c]-M02[c])*(M20[c]-M02[c]))));
1619 channel_moments[c].ellipse_angle=RadiansToDegrees(1.0/2.0*atan(2.0*
1620 M11[c]*PerceptibleReciprocal(M20[c]-M02[c])));
1621 if (fabs(M11[c]) < 0.0)
1623 if ((fabs(M20[c]-M02[c]) >= 0.0) &&
1624 ((M20[c]-M02[c]) < 0.0))
1625 channel_moments[c].ellipse_angle+=90.0;
1630 if (fabs(M20[c]-M02[c]) >= 0.0)
1632 if ((M20[c]-M02[c]) < 0.0)
1633 channel_moments[c].ellipse_angle+=90.0;
1635 channel_moments[c].ellipse_angle+=180.0;
1639 if ((fabs(M20[c]-M02[c]) >= 0.0) && ((M20[c]-M02[c]) < 0.0))
1640 channel_moments[c].ellipse_angle+=90.0;
1641 channel_moments[c].ellipse_eccentricity=sqrt(1.0-(
1642 channel_moments[c].ellipse_axis.y*
1643 channel_moments[c].ellipse_axis.y*PerceptibleReciprocal(
1644 channel_moments[c].ellipse_axis.x*
1645 channel_moments[c].ellipse_axis.x)));
1646 channel_moments[c].ellipse_intensity=M00[c]*
1647 PerceptibleReciprocal(MagickPI*channel_moments[c].ellipse_axis.x*
1648 channel_moments[c].ellipse_axis.y+MagickEpsilon);
1650 for (c=0; c <= MaxPixelChannels; c++)
1657 M11[c]*=PerceptibleReciprocal(pow(M00[c],1.0+(1.0+1.0)/2.0));
1658 M20[c]*=PerceptibleReciprocal(pow(M00[c],1.0+(2.0+0.0)/2.0));
1659 M02[c]*=PerceptibleReciprocal(pow(M00[c],1.0+(0.0+2.0)/2.0));
1660 M21[c]*=PerceptibleReciprocal(pow(M00[c],1.0+(2.0+1.0)/2.0));
1661 M12[c]*=PerceptibleReciprocal(pow(M00[c],1.0+(1.0+2.0)/2.0));
1662 M22[c]*=PerceptibleReciprocal(pow(M00[c],1.0+(2.0+2.0)/2.0));
1663 M30[c]*=PerceptibleReciprocal(pow(M00[c],1.0+(3.0+0.0)/2.0));
1664 M03[c]*=PerceptibleReciprocal(pow(M00[c],1.0+(0.0+3.0)/2.0));
1667 image_view=DestroyCacheView(image_view);
1668 for (c=0; c <= MaxPixelChannels; c++)
1673 channel_moments[c].invariant[0]=M20[c]+M02[c];
1674 channel_moments[c].invariant[1]=(M20[c]-M02[c])*(M20[c]-M02[c])+4.0*M11[c]*
1676 channel_moments[c].invariant[2]=(M30[c]-3.0*M12[c])*(M30[c]-3.0*M12[c])+
1677 (3.0*M21[c]-M03[c])*(3.0*M21[c]-M03[c]);
1678 channel_moments[c].invariant[3]=(M30[c]+M12[c])*(M30[c]+M12[c])+
1679 (M21[c]+M03[c])*(M21[c]+M03[c]);
1680 channel_moments[c].invariant[4]=(M30[c]-3.0*M12[c])*(M30[c]+M12[c])*
1681 ((M30[c]+M12[c])*(M30[c]+M12[c])-3.0*(M21[c]+M03[c])*(M21[c]+M03[c]))+
1682 (3.0*M21[c]-M03[c])*(M21[c]+M03[c])*(3.0*(M30[c]+M12[c])*(M30[c]+M12[c])-
1683 (M21[c]+M03[c])*(M21[c]+M03[c]));
1684 channel_moments[c].invariant[5]=(M20[c]-M02[c])*((M30[c]+M12[c])*
1685 (M30[c]+M12[c])-(M21[c]+M03[c])*(M21[c]+M03[c]))+4.0*M11[c]*
1686 (M30[c]+M12[c])*(M21[c]+M03[c]);
1687 channel_moments[c].invariant[6]=(3.0*M21[c]-M03[c])*(M30[c]+M12[c])*
1688 ((M30[c]+M12[c])*(M30[c]+M12[c])-3.0*(M21[c]+M03[c])*(M21[c]+M03[c]))-
1689 (M30[c]-3*M12[c])*(M21[c]+M03[c])*(3.0*(M30[c]+M12[c])*(M30[c]+M12[c])-
1690 (M21[c]+M03[c])*(M21[c]+M03[c]));
1691 channel_moments[c].invariant[7]=M11[c]*((M30[c]+M12[c])*(M30[c]+M12[c])-
1692 (M03[c]+M21[c])*(M03[c]+M21[c]))-(M20[c]-M02[c])*(M30[c]+M12[c])*
1695 if (y < (ssize_t) image->rows)
1696 channel_moments=(
ChannelMoments *) RelinquishMagickMemory(channel_moments);
1697 return(channel_moments);
1727 static inline double MagickLog10(
const double x)
1729 #define Log10Epsilon (1.0e-11)
1731 if (fabs(x) < Log10Epsilon)
1732 return(log10(Log10Epsilon));
1733 return(log10(fabs(x)));
1757 MaxPixelChannels+1UL,
sizeof(*perceptual_hash));
1760 artifact=GetImageArtifact(image,
"phash:colorspaces");
1761 if (artifact != NULL)
1762 colorspaces=AcquireString(artifact);
1764 colorspaces=AcquireString(
"sRGB,HCLp");
1765 perceptual_hash[0].number_colorspaces=0;
1766 perceptual_hash[0].number_channels=0;
1768 for (i=0; (p=StringToken(
",",&q)) != (
char *) NULL; i++)
1783 if (i >= MaximumNumberOfPerceptualColorspaces)
1785 colorspace=ParseCommandOption(MagickColorspaceOptions,MagickFalse,p);
1788 perceptual_hash[0].colorspace[i]=(ColorspaceType) colorspace;
1789 hash_image=BlurImage(image,0.0,1.0,exception);
1790 if (hash_image == (
Image *) NULL)
1792 hash_image->depth=8;
1793 status=TransformImageColorspace(hash_image,(ColorspaceType) colorspace,
1795 if (status == MagickFalse)
1797 moments=GetImageMoments(hash_image,exception);
1798 perceptual_hash[0].number_colorspaces++;
1799 perceptual_hash[0].number_channels+=GetImageChannels(hash_image);
1800 hash_image=DestroyImage(hash_image);
1803 for (channel=0; channel <= MaxPixelChannels; channel++)
1804 for (j=0; j < MaximumNumberOfImageMoments; j++)
1805 perceptual_hash[channel].phash[i][j]=
1806 (-MagickLog10(moments[channel].invariant[j]));
1809 colorspaces=DestroyString(colorspaces);
1810 return(perceptual_hash);
1842 MagickExport MagickBooleanType GetImageRange(
const Image *image,
double *minima,
1855 assert(image != (
Image *) NULL);
1856 assert(image->signature == MagickCoreSignature);
1857 if (IsEventLogging() != MagickFalse)
1858 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
1860 initialize=MagickTrue;
1863 image_view=AcquireVirtualCacheView(image,exception);
1864 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1865 #pragma omp parallel for schedule(static) shared(status,initialize) \
1866 magick_number_threads(image,image,image->rows,1)
1868 for (y=0; y < (ssize_t) image->rows; y++)
1883 if (status == MagickFalse)
1885 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1886 if (p == (
const Quantum *) NULL)
1891 row_initialize=MagickTrue;
1892 for (x=0; x < (ssize_t) image->columns; x++)
1897 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1899 PixelChannel channel = GetPixelChannelChannel(image,i);
1900 PixelTrait traits = GetPixelChannelTraits(image,channel);
1901 if (traits == UndefinedPixelTrait)
1903 if ((traits & UpdatePixelTrait) == 0)
1905 if (row_initialize != MagickFalse)
1907 row_minima=(double) p[i];
1908 row_maxima=(double) p[i];
1909 row_initialize=MagickFalse;
1913 if ((
double) p[i] < row_minima)
1914 row_minima=(
double) p[i];
1915 if ((
double) p[i] > row_maxima)
1916 row_maxima=(
double) p[i];
1919 p+=GetPixelChannels(image);
1921 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1922 #pragma omp critical (MagickCore_GetImageRange)
1925 if (initialize != MagickFalse)
1929 initialize=MagickFalse;
1933 if (row_minima < *minima)
1935 if (row_maxima > *maxima)
1940 image_view=DestroyCacheView(image_view);
1978 static ssize_t GetMedianPixel(Quantum *pixels,
const size_t n)
1980 #define SwapPixels(alpha,beta) \
1982 Quantum gamma=(alpha); \
1983 (alpha)=(beta);(beta)=gamma; \
1988 high = (ssize_t) n-1,
1989 median = (low+high)/2;
2000 if (high == (low+1))
2002 if (pixels[low] > pixels[high])
2003 SwapPixels(pixels[low],pixels[high]);
2006 if (pixels[mid] > pixels[high])
2007 SwapPixels(pixels[mid],pixels[high]);
2008 if (pixels[low] > pixels[high])
2009 SwapPixels(pixels[low], pixels[high]);
2010 if (pixels[mid] > pixels[low])
2011 SwapPixels(pixels[mid],pixels[low]);
2012 SwapPixels(pixels[mid],pixels[low+1]);
2015 do l++;
while (pixels[low] > pixels[l]);
2016 do h--;
while (pixels[h] > pixels[low]);
2019 SwapPixels(pixels[l],pixels[h]);
2021 SwapPixels(pixels[low],pixels[h]);
2033 *channel_statistics;
2060 assert(image != (
Image *) NULL);
2061 assert(image->signature == MagickCoreSignature);
2062 if (IsEventLogging() != MagickFalse)
2063 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
2064 histogram=(
double *) AcquireQuantumMemory(MaxMap+1UL,GetPixelChannels(image)*
2065 sizeof(*histogram));
2067 MaxPixelChannels+1,
sizeof(*channel_statistics));
2069 (histogram == (
double *) NULL))
2071 if (histogram != (
double *) NULL)
2072 histogram=(
double *) RelinquishMagickMemory(histogram);
2075 channel_statistics);
2076 return(channel_statistics);
2078 (void) memset(channel_statistics,0,(MaxPixelChannels+1)*
2079 sizeof(*channel_statistics));
2080 for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
2082 channel_statistics[i].depth=1;
2083 channel_statistics[i].maxima=(-MagickMaximumValue);
2084 channel_statistics[i].minima=MagickMaximumValue;
2086 (void) memset(histogram,0,(MaxMap+1)*GetPixelChannels(image)*
2087 sizeof(*histogram));
2088 for (y=0; y < (ssize_t) image->rows; y++)
2099 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
2100 if (p == (
const Quantum *) NULL)
2102 for (x=0; x < (ssize_t) image->columns; x++)
2104 if (GetPixelReadMask(image,p) <= (QuantumRange/2))
2106 p+=GetPixelChannels(image);
2109 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2111 PixelChannel channel = GetPixelChannelChannel(image,i);
2112 PixelTrait traits = GetPixelChannelTraits(image,channel);
2113 if (traits == UndefinedPixelTrait)
2115 if ((traits & UpdatePixelTrait) == 0)
2117 if (channel_statistics[channel].depth != MAGICKCORE_QUANTUM_DEPTH)
2119 depth=channel_statistics[channel].depth;
2120 range=GetQuantumRange(depth);
2121 status=p[i] != ScaleAnyToQuantum(ScaleQuantumToAny(p[i],range),
2122 range) ? MagickTrue : MagickFalse;
2123 if (status != MagickFalse)
2125 channel_statistics[channel].depth++;
2126 if (channel_statistics[channel].depth >
2127 channel_statistics[CompositePixelChannel].depth)
2128 channel_statistics[CompositePixelChannel].depth=
2129 channel_statistics[channel].depth;
2134 if ((
double) p[i] < channel_statistics[channel].minima)
2135 channel_statistics[channel].minima=(double) p[i];
2136 if ((
double) p[i] > channel_statistics[channel].maxima)
2137 channel_statistics[channel].maxima=(
double) p[i];
2138 channel_statistics[channel].sum+=p[i];
2139 channel_statistics[channel].sum_squared+=(double) p[i]*p[i];
2140 channel_statistics[channel].sum_cubed+=(double) p[i]*p[i]*p[i];
2141 channel_statistics[channel].sum_fourth_power+=(double) p[i]*p[i]*p[i]*
2143 channel_statistics[channel].area++;
2144 if ((
double) p[i] < channel_statistics[CompositePixelChannel].minima)
2145 channel_statistics[CompositePixelChannel].minima=(
double) p[i];
2146 if ((
double) p[i] > channel_statistics[CompositePixelChannel].maxima)
2147 channel_statistics[CompositePixelChannel].maxima=(
double) p[i];
2148 histogram[GetPixelChannels(image)*ScaleQuantumToMap(
2149 ClampToQuantum((
double) p[i]))+i]++;
2150 channel_statistics[CompositePixelChannel].sum+=(double) p[i];
2151 channel_statistics[CompositePixelChannel].sum_squared+=(double)
2153 channel_statistics[CompositePixelChannel].sum_cubed+=(double)
2155 channel_statistics[CompositePixelChannel].sum_fourth_power+=(double)
2156 p[i]*p[i]*p[i]*p[i];
2157 channel_statistics[CompositePixelChannel].area++;
2159 p+=GetPixelChannels(image);
2162 for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
2167 area=PerceptibleReciprocal(channel_statistics[i].area);
2168 channel_statistics[i].sum*=area;
2169 channel_statistics[i].sum_squared*=area;
2170 channel_statistics[i].sum_cubed*=area;
2171 channel_statistics[i].sum_fourth_power*=area;
2172 channel_statistics[i].mean=channel_statistics[i].sum;
2173 channel_statistics[i].variance=channel_statistics[i].sum_squared;
2174 standard_deviation=sqrt(channel_statistics[i].variance-
2175 (channel_statistics[i].mean*channel_statistics[i].mean));
2176 standard_deviation=sqrt(PerceptibleReciprocal(channel_statistics[i].area-
2177 1.0)*channel_statistics[i].area*standard_deviation*standard_deviation);
2178 channel_statistics[i].standard_deviation=standard_deviation;
2180 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2191 PixelChannel channel = GetPixelChannelChannel(image,i);
2193 for (j=0; j <= (ssize_t) MaxMap; j++)
2194 if (histogram[GetPixelChannels(image)*j+i] > 0.0)
2196 area=PerceptibleReciprocal(channel_statistics[channel].area);
2197 for (j=0; j <= (ssize_t) MaxMap; j++)
2202 count=area*histogram[GetPixelChannels(image)*j+i];
2203 channel_statistics[channel].entropy+=-count*MagickLog10(count)*
2204 PerceptibleReciprocal(MagickLog10(number_bins));
2205 channel_statistics[CompositePixelChannel].entropy+=-count*
2206 MagickLog10(count)*PerceptibleReciprocal(MagickLog10(number_bins))/
2207 GetPixelChannels(image);
2210 histogram=(
double *) RelinquishMagickMemory(histogram);
2211 for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
2216 standard_deviation=PerceptibleReciprocal(
2217 channel_statistics[i].standard_deviation);
2218 channel_statistics[i].skewness=(channel_statistics[i].sum_cubed-3.0*
2219 channel_statistics[i].mean*channel_statistics[i].sum_squared+2.0*
2220 channel_statistics[i].mean*channel_statistics[i].mean*
2221 channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2222 standard_deviation);
2223 channel_statistics[i].kurtosis=(channel_statistics[i].sum_fourth_power-4.0*
2224 channel_statistics[i].mean*channel_statistics[i].sum_cubed+6.0*
2225 channel_statistics[i].mean*channel_statistics[i].mean*
2226 channel_statistics[i].sum_squared-3.0*channel_statistics[i].mean*
2227 channel_statistics[i].mean*1.0*channel_statistics[i].mean*
2228 channel_statistics[i].mean)*(standard_deviation*standard_deviation*
2229 standard_deviation*standard_deviation)-3.0;
2231 median_info=AcquireVirtualMemory(image->columns,image->rows*
sizeof(*median));
2233 (void) ThrowMagickException(exception,GetMagickModule(),
2234 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
2237 median=(Quantum *) GetVirtualMemoryBlob(median_info);
2238 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2246 PixelChannel channel = GetPixelChannelChannel(image,i);
2247 PixelTrait traits = GetPixelChannelTraits(image,channel);
2248 if (traits == UndefinedPixelTrait)
2250 if ((traits & UpdatePixelTrait) == 0)
2252 for (y=0; y < (ssize_t) image->rows; y++)
2260 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
2261 if (p == (
const Quantum *) NULL)
2263 for (x=0; x < (ssize_t) image->columns; x++)
2265 if (GetPixelReadMask(image,p) <= (QuantumRange/2))
2267 p+=GetPixelChannels(image);
2272 p+=GetPixelChannels(image);
2274 channel_statistics[channel].median=(double) median[
2275 GetMedianPixel(median,n)];
2277 median_info=RelinquishVirtualMemory(median_info);
2279 channel_statistics[CompositePixelChannel].mean=0.0;
2280 channel_statistics[CompositePixelChannel].median=0.0;
2281 channel_statistics[CompositePixelChannel].standard_deviation=0.0;
2282 channel_statistics[CompositePixelChannel].entropy=0.0;
2283 for (i=0; i < (ssize_t) MaxPixelChannels; i++)
2285 channel_statistics[CompositePixelChannel].mean+=
2286 channel_statistics[i].mean;
2287 channel_statistics[CompositePixelChannel].median+=
2288 channel_statistics[i].median;
2289 channel_statistics[CompositePixelChannel].standard_deviation+=
2290 channel_statistics[i].standard_deviation;
2291 channel_statistics[CompositePixelChannel].entropy+=
2292 channel_statistics[i].entropy;
2294 channels=(double) GetImageChannels(image);
2295 channel_statistics[CompositePixelChannel].mean/=channels;
2296 channel_statistics[CompositePixelChannel].median/=channels;
2297 channel_statistics[CompositePixelChannel].standard_deviation/=channels;
2298 channel_statistics[CompositePixelChannel].entropy/=channels;
2299 if (y < (ssize_t) image->rows)
2301 channel_statistics);
2302 return(channel_statistics);
2338 MagickExport
Image *PolynomialImage(
const Image *images,
2339 const size_t number_terms,
const double *terms,
ExceptionInfo *exception)
2341 #define PolynomialImageTag "Polynomial/Image"
2356 **magick_restrict polynomial_pixels;
2364 assert(images != (
Image *) NULL);
2365 assert(images->signature == MagickCoreSignature);
2366 if (IsEventLogging() != MagickFalse)
2367 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",images->filename);
2369 assert(exception->signature == MagickCoreSignature);
2370 image=AcquireImageCanvas(images,exception);
2371 if (image == (
Image *) NULL)
2372 return((
Image *) NULL);
2373 if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
2375 image=DestroyImage(image);
2376 return((
Image *) NULL);
2378 number_images=GetImageListLength(images);
2379 polynomial_pixels=AcquirePixelTLS(images);
2382 image=DestroyImage(image);
2383 (void) ThrowMagickException(exception,GetMagickModule(),
2384 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",images->filename);
2385 return((
Image *) NULL);
2392 polynomial_view=AcquireAuthenticCacheView(image,exception);
2393 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2394 #pragma omp parallel for schedule(static) shared(progress,status) \
2395 magick_number_threads(image,image,image->rows,1)
2397 for (y=0; y < (ssize_t) image->rows; y++)
2406 id = GetOpenMPThreadId();
2419 if (status == MagickFalse)
2421 q=QueueCacheViewAuthenticPixels(polynomial_view,0,y,image->columns,1,
2423 if (q == (Quantum *) NULL)
2428 polynomial_pixel=polynomial_pixels[id];
2429 for (j=0; j < (ssize_t) image->columns; j++)
2430 for (i=0; i < MaxPixelChannels; i++)
2431 polynomial_pixel[j].channel[i]=0.0;
2433 for (j=0; j < (ssize_t) number_images; j++)
2438 if (j >= (ssize_t) number_terms)
2440 image_view=AcquireVirtualCacheView(next,exception);
2441 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2442 if (p == (
const Quantum *) NULL)
2444 image_view=DestroyCacheView(image_view);
2447 for (x=0; x < (ssize_t) image->columns; x++)
2449 for (i=0; i < (ssize_t) GetPixelChannels(next); i++)
2455 PixelChannel channel = GetPixelChannelChannel(image,i);
2456 PixelTrait traits = GetPixelChannelTraits(next,channel);
2457 PixelTrait polynomial_traits=GetPixelChannelTraits(image,channel);
2458 if ((traits == UndefinedPixelTrait) ||
2459 (polynomial_traits == UndefinedPixelTrait))
2461 if ((traits & UpdatePixelTrait) == 0)
2463 coefficient=(MagickRealType) terms[2*j];
2464 degree=(MagickRealType) terms[(j << 1)+1];
2465 polynomial_pixel[x].channel[i]+=coefficient*
2466 pow(QuantumScale*GetPixelChannel(image,channel,p),degree);
2468 p+=GetPixelChannels(next);
2470 image_view=DestroyCacheView(image_view);
2471 next=GetNextImageInList(next);
2473 for (x=0; x < (ssize_t) image->columns; x++)
2475 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2477 PixelChannel channel = GetPixelChannelChannel(image,i);
2478 PixelTrait traits = GetPixelChannelTraits(image,channel);
2479 if (traits == UndefinedPixelTrait)
2481 if ((traits & UpdatePixelTrait) == 0)
2483 q[i]=ClampToQuantum(QuantumRange*polynomial_pixel[x].channel[i]);
2485 q+=GetPixelChannels(image);
2487 if (SyncCacheViewAuthenticPixels(polynomial_view,exception) == MagickFalse)
2489 if (images->progress_monitor != (MagickProgressMonitor) NULL)
2494 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2498 proceed=SetImageProgress(images,PolynomialImageTag,progress,
2500 if (proceed == MagickFalse)
2504 polynomial_view=DestroyCacheView(polynomial_view);
2505 polynomial_pixels=DestroyPixelTLS(images,polynomial_pixels);
2506 if (status == MagickFalse)
2507 image=DestroyImage(image);
2578 if (pixel_list->skip_list.nodes != (
SkipNode *) NULL)
2579 pixel_list->skip_list.nodes=(
SkipNode *) RelinquishAlignedMemory(
2580 pixel_list->skip_list.nodes);
2581 pixel_list=(
PixelList *) RelinquishMagickMemory(pixel_list);
2590 assert(pixel_list != (
PixelList **) NULL);
2591 for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
2592 if (pixel_list[i] != (
PixelList *) NULL)
2593 pixel_list[i]=DestroyPixelList(pixel_list[i]);
2594 pixel_list=(
PixelList **) RelinquishMagickMemory(pixel_list);
2598 static PixelList *AcquirePixelList(
const size_t width,
const size_t height)
2603 pixel_list=(
PixelList *) AcquireMagickMemory(
sizeof(*pixel_list));
2606 (void) memset((
void *) pixel_list,0,
sizeof(*pixel_list));
2607 pixel_list->length=width*height;
2608 pixel_list->skip_list.nodes=(
SkipNode *) AcquireAlignedMemory(65537UL,
2609 sizeof(*pixel_list->skip_list.nodes));
2610 if (pixel_list->skip_list.nodes == (
SkipNode *) NULL)
2611 return(DestroyPixelList(pixel_list));
2612 (void) memset(pixel_list->skip_list.nodes,0,65537UL*
2613 sizeof(*pixel_list->skip_list.nodes));
2614 pixel_list->signature=MagickCoreSignature;
2618 static PixelList **AcquirePixelListTLS(
const size_t width,
2619 const size_t height)
2630 number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
2631 pixel_list=(
PixelList **) AcquireQuantumMemory(number_threads,
2632 sizeof(*pixel_list));
2635 (void) memset(pixel_list,0,number_threads*
sizeof(*pixel_list));
2636 for (i=0; i < (ssize_t) number_threads; i++)
2638 pixel_list[i]=AcquirePixelList(width,height);
2639 if (pixel_list[i] == (
PixelList *) NULL)
2640 return(DestroyPixelListTLS(pixel_list));
2645 static void AddNodePixelList(
PixelList *pixel_list,
const size_t color)
2660 p=(&pixel_list->skip_list);
2661 p->nodes[color].signature=pixel_list->signature;
2662 p->nodes[color].count=1;
2667 for (level=p->level; level >= 0; level--)
2669 while (p->nodes[search].next[level] < color)
2670 search=p->nodes[search].next[level];
2671 update[level]=search;
2676 for (level=0; ; level++)
2678 pixel_list->seed=(pixel_list->seed*42893621L)+1L;
2679 if ((pixel_list->seed & 0x300) != 0x300)
2684 if (level > (p->level+2))
2689 while (level > p->level)
2692 update[p->level]=65536UL;
2699 p->nodes[color].next[level]=p->nodes[update[level]].next[level];
2700 p->nodes[update[level]].next[level]=color;
2701 }
while (level-- > 0);
2704 static inline void GetMedianPixelList(
PixelList *pixel_list,Quantum *pixel)
2718 p=(&pixel_list->skip_list);
2723 color=p->nodes[color].next[0];
2724 count+=p->nodes[color].count;
2725 }
while (count <= (ssize_t) (pixel_list->length >> 1));
2726 *pixel=ScaleShortToQuantum((
unsigned short) color);
2729 static inline void GetModePixelList(
PixelList *pixel_list,Quantum *pixel)
2745 p=(&pixel_list->skip_list);
2748 max_count=p->nodes[mode].count;
2752 color=p->nodes[color].next[0];
2753 if (p->nodes[color].count > max_count)
2756 max_count=p->nodes[mode].count;
2758 count+=p->nodes[color].count;
2759 }
while (count < (ssize_t) pixel_list->length);
2760 *pixel=ScaleShortToQuantum((
unsigned short) mode);
2763 static inline void GetNonpeakPixelList(
PixelList *pixel_list,Quantum *pixel)
2779 p=(&pixel_list->skip_list);
2781 next=p->nodes[color].next[0];
2787 next=p->nodes[color].next[0];
2788 count+=p->nodes[color].count;
2789 }
while (count <= (ssize_t) (pixel_list->length >> 1));
2790 if ((previous == 65536UL) && (next != 65536UL))
2793 if ((previous != 65536UL) && (next == 65536UL))
2795 *pixel=ScaleShortToQuantum((
unsigned short) color);
2798 static inline void InsertPixelList(
const Quantum pixel,
PixelList *pixel_list)
2806 index=ScaleQuantumToShort(pixel);
2807 signature=pixel_list->skip_list.nodes[index].signature;
2808 if (signature == pixel_list->signature)
2810 pixel_list->skip_list.nodes[index].count++;
2813 AddNodePixelList(pixel_list,index);
2816 static void ResetPixelList(
PixelList *pixel_list)
2830 p=(&pixel_list->skip_list);
2831 root=p->nodes+65536UL;
2833 for (level=0; level < 9; level++)
2834 root->next[level]=65536UL;
2835 pixel_list->seed=pixel_list->signature++;
2838 MagickExport
Image *StatisticImage(
const Image *image,
const StatisticType type,
2839 const size_t width,
const size_t height,
ExceptionInfo *exception)
2841 #define StatisticImageTag "Statistic/Image"
2857 **magick_restrict pixel_list;
2866 assert(image != (
Image *) NULL);
2867 assert(image->signature == MagickCoreSignature);
2868 if (IsEventLogging() != MagickFalse)
2869 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
2871 assert(exception->signature == MagickCoreSignature);
2872 statistic_image=CloneImage(image,0,0,MagickTrue,
2874 if (statistic_image == (
Image *) NULL)
2875 return((
Image *) NULL);
2876 status=SetImageStorageClass(statistic_image,DirectClass,exception);
2877 if (status == MagickFalse)
2879 statistic_image=DestroyImage(statistic_image);
2880 return((
Image *) NULL);
2882 pixel_list=AcquirePixelListTLS(MagickMax(width,1),MagickMax(height,1));
2885 statistic_image=DestroyImage(statistic_image);
2886 ThrowImageException(ResourceLimitError,
"MemoryAllocationFailed");
2891 center=(ssize_t) GetPixelChannels(image)*(image->columns+MagickMax(width,1))*
2892 (MagickMax(height,1)/2L)+GetPixelChannels(image)*(MagickMax(width,1)/2L);
2895 image_view=AcquireVirtualCacheView(image,exception);
2896 statistic_view=AcquireAuthenticCacheView(statistic_image,exception);
2897 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2898 #pragma omp parallel for schedule(static) shared(progress,status) \
2899 magick_number_threads(image,statistic_image,statistic_image->rows,1)
2901 for (y=0; y < (ssize_t) statistic_image->rows; y++)
2904 id = GetOpenMPThreadId();
2915 if (status == MagickFalse)
2917 p=GetCacheViewVirtualPixels(image_view,-((ssize_t) MagickMax(width,1)/2L),y-
2918 (ssize_t) (MagickMax(height,1)/2L),image->columns+MagickMax(width,1),
2919 MagickMax(height,1),exception);
2920 q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns, 1,exception);
2921 if ((p == (
const Quantum *) NULL) || (q == (Quantum *) NULL))
2926 for (x=0; x < (ssize_t) statistic_image->columns; x++)
2931 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2944 *magick_restrict pixels;
2952 PixelChannel channel = GetPixelChannelChannel(image,i);
2953 PixelTrait traits = GetPixelChannelTraits(image,channel);
2954 PixelTrait statistic_traits=GetPixelChannelTraits(statistic_image,
2956 if ((traits == UndefinedPixelTrait) ||
2957 (statistic_traits == UndefinedPixelTrait))
2959 if (((statistic_traits & CopyPixelTrait) != 0) ||
2960 (GetPixelWriteMask(image,p) <= (QuantumRange/2)))
2962 SetPixelChannel(statistic_image,channel,p[center+i],q);
2965 if ((statistic_traits & UpdatePixelTrait) == 0)
2973 ResetPixelList(pixel_list[
id]);
2974 for (v=0; v < (ssize_t) MagickMax(height,1); v++)
2976 for (u=0; u < (ssize_t) MagickMax(width,1); u++)
2978 if ((type == MedianStatistic) || (type == ModeStatistic) ||
2979 (type == NonpeakStatistic))
2981 InsertPixelList(pixels[i],pixel_list[
id]);
2982 pixels+=GetPixelChannels(image);
2986 if (pixels[i] < minimum)
2987 minimum=(double) pixels[i];
2988 if (pixels[i] > maximum)
2989 maximum=(double) pixels[i];
2990 sum+=(double) pixels[i];
2991 sum_squared+=(double) pixels[i]*pixels[i];
2992 pixels+=GetPixelChannels(image);
2994 pixels+=GetPixelChannels(image)*image->columns;
2998 case ContrastStatistic:
3000 pixel=ClampToQuantum(MagickAbsoluteValue((maximum-minimum)*
3001 PerceptibleReciprocal(maximum+minimum)));
3004 case GradientStatistic:
3006 pixel=ClampToQuantum(MagickAbsoluteValue(maximum-minimum));
3009 case MaximumStatistic:
3011 pixel=ClampToQuantum(maximum);
3017 pixel=ClampToQuantum(sum/area);
3020 case MedianStatistic:
3022 GetMedianPixelList(pixel_list[
id],&pixel);
3025 case MinimumStatistic:
3027 pixel=ClampToQuantum(minimum);
3032 GetModePixelList(pixel_list[
id],&pixel);
3035 case NonpeakStatistic:
3037 GetNonpeakPixelList(pixel_list[
id],&pixel);
3040 case RootMeanSquareStatistic:
3042 pixel=ClampToQuantum(sqrt(sum_squared/area));
3045 case StandardDeviationStatistic:
3047 pixel=ClampToQuantum(sqrt(sum_squared/area-(sum/area*sum/area)));
3051 SetPixelChannel(statistic_image,channel,pixel,q);
3053 p+=GetPixelChannels(image);
3054 q+=GetPixelChannels(statistic_image);
3056 if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
3058 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3063 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3067 proceed=SetImageProgress(image,StatisticImageTag,progress,image->rows);
3068 if (proceed == MagickFalse)
3072 statistic_view=DestroyCacheView(statistic_view);
3073 image_view=DestroyCacheView(image_view);
3074 pixel_list=DestroyPixelListTLS(pixel_list);
3075 if (status == MagickFalse)
3076 statistic_image=DestroyImage(statistic_image);
3077 return(statistic_image);