]> granicus.if.org Git - imagemagick/blobdiff - MagickCore/statistic.c
(no commit message)
[imagemagick] / MagickCore / statistic.c
index 7c570f64288fcc1d3a57d901ead5d203852693cb..42ea5facbb4b210791e6c34bee856980f04092fb 100644 (file)
 %                     MagickCore Image Statistical Methods                    %
 %                                                                             %
 %                              Software Design                                %
-%                                John Cristy                                  %
+%                                   Cristy                                    %
 %                                 July 1992                                   %
 %                                                                             %
 %                                                                             %
-%  Copyright 1999-2013 ImageMagick Studio LLC, a non-profit organization      %
+%  Copyright 1999-2014 ImageMagick Studio LLC, a non-profit organization      %
 %  dedicated to making software imaging solutions freely available.           %
 %                                                                             %
 %  You may not use this file except in compliance with the License.  You may  %
@@ -650,7 +650,7 @@ MagickExport Image *EvaluateImages(const Image *images,
             register ssize_t
               i;
 
-            if (GetPixelMask(next,p) != 0)
+            if (GetPixelReadMask(next,p) == 0)
               {
                 p+=GetPixelChannels(next);
                 continue;
@@ -708,7 +708,7 @@ MagickExport Image *EvaluateImages(const Image *images,
           register ssize_t
             i;
 
-          if (GetPixelMask(image,q) != 0)
+          if (GetPixelReadMask(image,q) == 0)
             {
               q+=GetPixelChannels(image);
               continue;
@@ -823,7 +823,7 @@ MagickExport MagickBooleanType EvaluateImage(Image *image,
         if (traits == UndefinedPixelTrait)
           continue;
         if (((traits & CopyPixelTrait) != 0) ||
-            (GetPixelMask(image,q) != 0))
+            (GetPixelReadMask(image,q) == 0))
           continue;
         q[i]=ClampToQuantum(ApplyEvaluateOperator(random_info[id],q[i],op,
           value));
@@ -1036,7 +1036,7 @@ MagickExport MagickBooleanType FunctionImage(Image *image,
       register ssize_t
         i;
 
-      if (GetPixelMask(image,q) != 0)
+      if (GetPixelReadMask(image,q) == 0)
         {
           q+=GetPixelChannels(image);
           continue;
@@ -1127,6 +1127,144 @@ MagickExport MagickBooleanType GetImageExtrema(const Image *image,
 %                                                                             %
 %                                                                             %
 %                                                                             %
+%   G e t I m a g e K u r t o s i s                                           %
+%                                                                             %
+%                                                                             %
+%                                                                             %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  GetImageKurtosis() returns the kurtosis and skewness of one or more image
+%  channels.
+%
+%  The format of the GetImageKurtosis method is:
+%
+%      MagickBooleanType GetImageKurtosis(const Image *image,double *kurtosis,
+%        double *skewness,ExceptionInfo *exception)
+%
+%  A description of each parameter follows:
+%
+%    o image: the image.
+%
+%    o kurtosis: the kurtosis of the channel.
+%
+%    o skewness: the skewness of the channel.
+%
+%    o exception: return any errors or warnings in this structure.
+%
+*/
+MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
+  double *kurtosis,double *skewness,ExceptionInfo *exception)
+{
+  CacheView
+    *image_view;
+
+  double
+    area,
+    mean,
+    standard_deviation,
+    sum_squares,
+    sum_cubes,
+    sum_fourth_power;
+
+  MagickBooleanType
+    status;
+
+  ssize_t
+    y;
+
+  assert(image != (Image *) NULL);
+  assert(image->signature == MagickSignature);
+  if (image->debug != MagickFalse)
+    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
+  status=MagickTrue;
+  *kurtosis=0.0;
+  *skewness=0.0;
+  area=0.0;
+  mean=0.0;
+  standard_deviation=0.0;
+  sum_squares=0.0;
+  sum_cubes=0.0;
+  sum_fourth_power=0.0;
+  image_view=AcquireVirtualCacheView(image,exception);
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+  #pragma omp parallel for schedule(static,4) shared(status) \
+    magick_threads(image,image,image->rows,1)
+#endif
+  for (y=0; y < (ssize_t) image->rows; y++)
+  {
+    register const Quantum
+      *restrict p;
+
+    register ssize_t
+      x;
+
+    if (status == MagickFalse)
+      continue;
+    p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
+    if (p == (const Quantum *) NULL)
+      {
+        status=MagickFalse;
+        continue;
+      }
+    for (x=0; x < (ssize_t) image->columns; x++)
+    {
+      register ssize_t
+        i;
+
+      if (GetPixelReadMask(image,p) == 0)
+        {
+          p+=GetPixelChannels(image);
+          continue;
+        }
+      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
+      {
+        PixelChannel channel=GetPixelChannelChannel(image,i);
+        PixelTrait traits=GetPixelChannelTraits(image,channel);
+        if (traits == UndefinedPixelTrait)
+          continue;
+        if ((traits & UpdatePixelTrait) == 0)
+          continue;
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+        #pragma omp critical (MagickCore_GetImageKurtosis)
+#endif
+        {
+          mean+=p[i];
+          sum_squares+=(double) p[i]*p[i];
+          sum_cubes+=(double) p[i]*p[i]*p[i];
+          sum_fourth_power+=(double) p[i]*p[i]*p[i]*p[i];
+          area++;
+        }
+      }
+      p+=GetPixelChannels(image);
+    }
+  }
+  image_view=DestroyCacheView(image_view);
+  if (area != 0.0)
+    {
+      mean/=area;
+      sum_squares/=area;
+      sum_cubes/=area;
+      sum_fourth_power/=area;
+    }
+  standard_deviation=sqrt(sum_squares-(mean*mean));
+  if (standard_deviation != 0.0)
+    {
+      *kurtosis=sum_fourth_power-4.0*mean*sum_cubes+6.0*mean*mean*sum_squares-
+        3.0*mean*mean*mean*mean;
+      *kurtosis/=standard_deviation*standard_deviation*standard_deviation*
+        standard_deviation;
+      *kurtosis-=3.0;
+      *skewness=sum_cubes-3.0*mean*sum_squares+2.0*mean*mean*mean;
+      *skewness/=standard_deviation*standard_deviation*standard_deviation;
+    }
+  return(status);
+}
+\f
+/*
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                             %
+%                                                                             %
+%                                                                             %
 %   G e t I m a g e M e a n                                                   %
 %                                                                             %
 %                                                                             %
@@ -1204,69 +1342,81 @@ MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean,
 %                                                                             %
 %                                                                             %
 %                                                                             %
-%   G e t I m a g e K u r t o s i s                                           %
+%   G e t I m a g e M o m e n t s                                             %
 %                                                                             %
 %                                                                             %
 %                                                                             %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-%  GetImageKurtosis() returns the kurtosis and skewness of one or more
-%  image channels.
+%  GetImageMoments() returns the normalized moments of one or more image
+%  channels.
 %
-%  The format of the GetImageKurtosis method is:
+%  The format of the GetImageMoments method is:
 %
-%      MagickBooleanType GetImageKurtosis(const Image *image,double *kurtosis,
-%        double *skewness,ExceptionInfo *exception)
+%      ChannelMoments *GetImageMoments(const Image *image,
+%        ExceptionInfo *exception)
 %
 %  A description of each parameter follows:
 %
 %    o image: the image.
 %
-%    o kurtosis: the kurtosis of the channel.
-%
-%    o skewness: the skewness of the channel.
-%
 %    o exception: return any errors or warnings in this structure.
 %
 */
-MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
-  double *kurtosis,double *skewness,ExceptionInfo *exception)
+MagickExport ChannelMoments *GetImageMoments(const Image *image,
+  ExceptionInfo *exception)
 {
+#define MaxNumberImageMoments  8
+
   CacheView
     *image_view;
 
-  double
-    area,
-    mean,
-    standard_deviation,
-    sum_squares,
-    sum_cubes,
-    sum_fourth_power;
+  ChannelMoments
+    *channel_moments;
 
-  MagickBooleanType
-    status;
+  double
+    M00[MaxPixelChannels+1],
+    M01[MaxPixelChannels+1],
+    M02[MaxPixelChannels+1],
+    M03[MaxPixelChannels+1],
+    M10[MaxPixelChannels+1],
+    M11[MaxPixelChannels+1],
+    M12[MaxPixelChannels+1],
+    M20[MaxPixelChannels+1],
+    M21[MaxPixelChannels+1],
+    M22[MaxPixelChannels+1],
+    M30[MaxPixelChannels+1];
+
+  PointInfo
+    centroid[MaxPixelChannels+1];
 
   ssize_t
+    channel,
     y;
 
   assert(image != (Image *) NULL);
   assert(image->signature == MagickSignature);
   if (image->debug != MagickFalse)
     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
-  status=MagickTrue;
-  *kurtosis=0.0;
-  *skewness=0.0;
-  area=0.0;
-  mean=0.0;
-  standard_deviation=0.0;
-  sum_squares=0.0;
-  sum_cubes=0.0;
-  sum_fourth_power=0.0;
+  channel_moments=(ChannelMoments *) AcquireQuantumMemory(MaxPixelChannels+1,
+    sizeof(*channel_moments));
+  if (channel_moments == (ChannelMoments *) NULL)
+    return(channel_moments);
+  (void) ResetMagickMemory(channel_moments,0,(MaxPixelChannels+1)*
+    sizeof(*channel_moments));
+  (void) ResetMagickMemory(centroid,0,sizeof(centroid));
+  (void) ResetMagickMemory(M00,0,sizeof(M00));
+  (void) ResetMagickMemory(M01,0,sizeof(M01));
+  (void) ResetMagickMemory(M02,0,sizeof(M02));
+  (void) ResetMagickMemory(M03,0,sizeof(M03));
+  (void) ResetMagickMemory(M10,0,sizeof(M10));
+  (void) ResetMagickMemory(M11,0,sizeof(M11));
+  (void) ResetMagickMemory(M12,0,sizeof(M12));
+  (void) ResetMagickMemory(M20,0,sizeof(M20));
+  (void) ResetMagickMemory(M21,0,sizeof(M21));
+  (void) ResetMagickMemory(M22,0,sizeof(M22));
+  (void) ResetMagickMemory(M30,0,sizeof(M30));
   image_view=AcquireVirtualCacheView(image,exception);
-#if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(static,4) shared(status) \
-    magick_threads(image,image,image->rows,1)
-#endif
   for (y=0; y < (ssize_t) image->rows; y++)
   {
     register const Quantum
@@ -1275,62 +1425,326 @@ MagickExport MagickBooleanType GetImageKurtosis(const Image *image,
     register ssize_t
       x;
 
-    if (status == MagickFalse)
-      continue;
+    /*
+      Compute center of mass (centroid).
+    */
     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
     if (p == (const Quantum *) NULL)
+      break;
+    for (x=0; x < (ssize_t) image->columns; x++)
+    {
+      register ssize_t
+        i;
+
+      if (GetPixelReadMask(image,p) == 0)
+        {
+          p+=GetPixelChannels(image);
+          continue;
+        }
+      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
       {
-        status=MagickFalse;
+        PixelChannel channel=GetPixelChannelChannel(image,i);
+        PixelTrait traits=GetPixelChannelTraits(image,channel);
+        if (traits == UndefinedPixelTrait)
+          continue;
+        if ((traits & UpdatePixelTrait) == 0)
+          continue;
+        M00[channel]+=QuantumScale*p[i];
+        M00[MaxPixelChannels]+=QuantumScale*p[i];
+        M10[channel]+=x*QuantumScale*p[i];
+        M10[MaxPixelChannels]+=QuantumScale*p[i];
+        M01[channel]+=y*QuantumScale*p[i];
+        M01[MaxPixelChannels]+=QuantumScale*p[i];
+      }
+      p+=GetPixelChannels(image);
+    }
+  }
+  for (channel=0; channel <= MaxPixelChannels; channel++)
+  {
+    /*
+       Compute center of mass (centroid).
+    */
+    if (M00[channel] < MagickEpsilon)
+      {
+        M00[channel]+=MagickEpsilon;
+        centroid[channel].x=image->columns/2.0;
+        centroid[channel].y=image->rows/2.0;
         continue;
       }
+    M00[channel]+=MagickEpsilon;
+    centroid[channel].x=M10[channel]/M00[channel];
+    centroid[channel].y=M01[channel]/M00[channel];
+  }
+  for (y=0; y < (ssize_t) image->rows; y++)
+  {
+    register const Quantum
+      *restrict p;
+
+    register ssize_t
+      x;
+
+    /*
+      Compute the image moments.
+    */
+    p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
+    if (p == (const Quantum *) NULL)
+      break;
     for (x=0; x < (ssize_t) image->columns; x++)
     {
       register ssize_t
         i;
 
+      if (GetPixelReadMask(image,p) == 0)
+        {
+          p+=GetPixelChannels(image);
+          continue;
+        }
       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
       {
         PixelChannel channel=GetPixelChannelChannel(image,i);
         PixelTrait traits=GetPixelChannelTraits(image,channel);
         if (traits == UndefinedPixelTrait)
           continue;
-        if (((traits & UpdatePixelTrait) == 0) ||
-            (GetPixelMask(image,p) != 0))
+        if ((traits & UpdatePixelTrait) == 0)
           continue;
-#if defined(MAGICKCORE_OPENMP_SUPPORT)
-        #pragma omp critical (MagickCore_GetImageKurtosis)
-#endif
-        {
-          mean+=p[i];
-          sum_squares+=(double) p[i]*p[i];
-          sum_cubes+=(double) p[i]*p[i]*p[i];
-          sum_fourth_power+=(double) p[i]*p[i]*p[i]*p[i];
-          area++;
-        }
+        M11[channel]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
+          QuantumScale*p[i];
+        M11[MaxPixelChannels]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
+          QuantumScale*p[i];
+        M20[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
+          QuantumScale*p[i];
+        M20[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
+          QuantumScale*p[i];
+        M02[channel]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
+          QuantumScale*p[i];
+        M02[MaxPixelChannels]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
+          QuantumScale*p[i];
+        M21[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
+          (y-centroid[channel].y)*QuantumScale*p[i];
+        M21[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
+          (y-centroid[channel].y)*QuantumScale*p[i];
+        M12[channel]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
+          (y-centroid[channel].y)*QuantumScale*p[i];
+        M12[MaxPixelChannels]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
+          (y-centroid[channel].y)*QuantumScale*p[i];
+        M22[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
+          (y-centroid[channel].y)*(y-centroid[channel].y)*QuantumScale*p[i];
+        M22[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
+          (y-centroid[channel].y)*(y-centroid[channel].y)*QuantumScale*p[i];
+        M30[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
+          (x-centroid[channel].x)*QuantumScale*p[i];
+        M30[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)*
+          (x-centroid[channel].x)*QuantumScale*p[i];
+        M03[channel]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
+          (y-centroid[channel].y)*QuantumScale*p[i];
+        M03[MaxPixelChannels]+=(y-centroid[channel].y)*(y-centroid[channel].y)*
+          (y-centroid[channel].y)*QuantumScale*p[i];
       }
       p+=GetPixelChannels(image);
     }
   }
+  M00[MaxPixelChannels]/=GetPixelChannels(image);
+  M01[MaxPixelChannels]/=GetPixelChannels(image);
+  M02[MaxPixelChannels]/=GetPixelChannels(image);
+  M03[MaxPixelChannels]/=GetPixelChannels(image);
+  M10[MaxPixelChannels]/=GetPixelChannels(image);
+  M11[MaxPixelChannels]/=GetPixelChannels(image);
+  M12[MaxPixelChannels]/=GetPixelChannels(image);
+  M20[MaxPixelChannels]/=GetPixelChannels(image);
+  M21[MaxPixelChannels]/=GetPixelChannels(image);
+  M22[MaxPixelChannels]/=GetPixelChannels(image);
+  M30[MaxPixelChannels]/=GetPixelChannels(image);
+  for (channel=0; channel <= MaxPixelChannels; channel++)
+  {
+    /*
+      Compute elliptical angle, major and minor axes, eccentricity, & intensity.
+    */
+    channel_moments[channel].centroid=centroid[channel];
+    channel_moments[channel].ellipse_axis.x=sqrt((2.0/M00[channel])*
+      ((M20[channel]+M02[channel])+sqrt(4.0*M11[channel]*M11[channel]+
+      (M20[channel]-M02[channel])*(M20[channel]-M02[channel]))));
+    channel_moments[channel].ellipse_axis.y=sqrt((2.0/M00[channel])*
+      ((M20[channel]+M02[channel])-sqrt(4.0*M11[channel]*M11[channel]+
+      (M20[channel]-M02[channel])*(M20[channel]-M02[channel]))));
+    channel_moments[channel].ellipse_angle=RadiansToDegrees(0.5*atan(2.0*
+      M11[channel]/(M20[channel]-M02[channel]+MagickEpsilon)));
+    channel_moments[channel].ellipse_eccentricity=sqrt(1.0-(
+      channel_moments[channel].ellipse_axis.y/
+      (channel_moments[channel].ellipse_axis.x+MagickEpsilon)));
+    channel_moments[channel].ellipse_intensity=M00[channel]/
+      (MagickPI*channel_moments[channel].ellipse_axis.x*
+      channel_moments[channel].ellipse_axis.y+MagickEpsilon);
+  }
+  for (channel=0; channel <= MaxPixelChannels; channel++)
+  {
+    /*
+      Normalize image moments.
+    */
+    M10[channel]=0.0;
+    M01[channel]=0.0;
+    M11[channel]/=pow(M00[channel],1.0+(1.0+1.0)/2.0);
+    M20[channel]/=pow(M00[channel],1.0+(2.0+0.0)/2.0);
+    M02[channel]/=pow(M00[channel],1.0+(0.0+2.0)/2.0);
+    M21[channel]/=pow(M00[channel],1.0+(2.0+1.0)/2.0);
+    M12[channel]/=pow(M00[channel],1.0+(1.0+2.0)/2.0);
+    M22[channel]/=pow(M00[channel],1.0+(2.0+2.0)/2.0);
+    M30[channel]/=pow(M00[channel],1.0+(3.0+0.0)/2.0);
+    M03[channel]/=pow(M00[channel],1.0+(0.0+3.0)/2.0);
+    M00[channel]=1.0;
+  }
   image_view=DestroyCacheView(image_view);
-  if (area != 0.0)
+  for (channel=0; channel <= MaxPixelChannels; channel++)
+  {
+    /*
+      Compute Hu invariant moments.
+    */
+    channel_moments[channel].I[0]=M20[channel]+M02[channel];
+    channel_moments[channel].I[1]=(M20[channel]-M02[channel])*
+      (M20[channel]-M02[channel])+4.0*M11[channel]*M11[channel];
+    channel_moments[channel].I[2]=(M30[channel]-3.0*M12[channel])*
+      (M30[channel]-3.0*M12[channel])+(3.0*M21[channel]-M03[channel])*
+      (3.0*M21[channel]-M03[channel]);
+    channel_moments[channel].I[3]=(M30[channel]+M12[channel])*
+      (M30[channel]+M12[channel])+(M21[channel]+M03[channel])*
+      (M21[channel]+M03[channel]);
+    channel_moments[channel].I[4]=(M30[channel]-3.0*M12[channel])*
+      (M30[channel]+M12[channel])*((M30[channel]+M12[channel])*
+      (M30[channel]+M12[channel])-3.0*(M21[channel]+M03[channel])*
+      (M21[channel]+M03[channel]))+(3.0*M21[channel]-M03[channel])*
+      (M21[channel]+M03[channel])*(3.0*(M30[channel]+M12[channel])*
+      (M30[channel]+M12[channel])-(M21[channel]+M03[channel])*
+      (M21[channel]+M03[channel]));
+    channel_moments[channel].I[5]=(M20[channel]-M02[channel])*
+      ((M30[channel]+M12[channel])*(M30[channel]+M12[channel])-
+      (M21[channel]+M03[channel])*(M21[channel]+M03[channel]))+
+      4.0*M11[channel]*(M30[channel]+M12[channel])*(M21[channel]+M03[channel]);
+    channel_moments[channel].I[6]=(3.0*M21[channel]-M03[channel])*
+      (M30[channel]+M12[channel])*((M30[channel]+M12[channel])*
+      (M30[channel]+M12[channel])-3.0*(M21[channel]+M03[channel])*
+      (M21[channel]+M03[channel]))-(M30[channel]-3*M12[channel])*
+      (M21[channel]+M03[channel])*(3.0*(M30[channel]+M12[channel])*
+      (M30[channel]+M12[channel])-(M21[channel]+M03[channel])*
+      (M21[channel]+M03[channel]));
+    channel_moments[channel].I[7]=M11[channel]*((M30[channel]+M12[channel])*
+      (M30[channel]+M12[channel])-(M03[channel]+M21[channel])*
+      (M03[channel]+M21[channel]))-(M20[channel]-M02[channel])*
+      (M30[channel]+M12[channel])*(M03[channel]+M21[channel]);
+  }
+  if (y < (ssize_t) image->rows)
+    channel_moments=(ChannelMoments *) RelinquishMagickMemory(channel_moments);
+  return(channel_moments);
+}
+\f
+/*
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                             %
+%                                                                             %
+%                                                                             %
+%   G e t I m a g e C h a n n e l P e r c e p t u a l H a s h                 %
+%                                                                             %
+%                                                                             %
+%                                                                             %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  GetImagePerceptualHash() returns the perceptual hash of one or more
+%  image channels.
+%
+%  The format of the GetImagePerceptualHash method is:
+%
+%      ChannelPerceptualHash *GetImagePerceptualHash(const Image *image,
+%        ExceptionInfo *exception)
+%
+%  A description of each parameter follows:
+%
+%    o image: the image.
+%
+%    o exception: return any errors or warnings in this structure.
+%
+*/
+
+static inline double MagickLog10(const double x)
+{
+#define Log10Epsilon  (1.0e-11)
+
+ if (fabs(x) < Log10Epsilon)
+   return(log10(Log10Epsilon));
+ return(log10(fabs(x)));
+}
+
+MagickExport ChannelPerceptualHash *GetImagePerceptualHash(
+  const Image *image,ExceptionInfo *exception)
+{
+  ChannelMoments
+    *moments;
+
+  ChannelPerceptualHash
+    *perceptual_hash;
+
+  Image
+    *hash_image;
+
+  MagickBooleanType
+    status;
+
+  register ssize_t
+    i;
+
+  ssize_t
+    channel;
+
+  /*
+    Blur then transform to sRGB colorspace.
+  */
+  hash_image=BlurImage(image,0.0,1.0,exception);
+  if (hash_image == (Image *) NULL)
+    return((ChannelPerceptualHash *) NULL);
+  hash_image->depth=8;
+  status=TransformImageColorspace(hash_image,sRGBColorspace,exception);
+  if (status == MagickFalse)
+    return((ChannelPerceptualHash *) NULL);
+  moments=GetImageMoments(hash_image,exception);
+  hash_image=DestroyImage(hash_image);
+  if (moments == (ChannelMoments *) NULL)
+    return((ChannelPerceptualHash *) NULL);
+  perceptual_hash=(ChannelPerceptualHash *) AcquireQuantumMemory(
+    CompositeChannels+1UL,sizeof(*perceptual_hash));
+  if (perceptual_hash == (ChannelPerceptualHash *) NULL)
+    return((ChannelPerceptualHash *) NULL);
+  for (channel=0; channel <= MaxPixelChannels; channel++)
+    for (i=0; i < 7; i++)
+      perceptual_hash[channel].P[i]=(-MagickLog10(moments[channel].I[i]));
+  moments=(ChannelMoments *) RelinquishMagickMemory(moments);
+  /*
+    Blur then transform to HCLp colorspace.
+  */
+  hash_image=BlurImage(image,0.0,1.0,exception);
+  if (hash_image == (Image *) NULL)
     {
-      mean/=area;
-      sum_squares/=area;
-      sum_cubes/=area;
-      sum_fourth_power/=area;
+      perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
+        perceptual_hash);
+      return((ChannelPerceptualHash *) NULL);
     }
-  standard_deviation=sqrt(sum_squares-(mean*mean));
-  if (standard_deviation != 0.0)
+  hash_image->depth=8;
+  status=TransformImageColorspace(hash_image,HCLpColorspace,exception);
+  if (status == MagickFalse)
     {
-      *kurtosis=sum_fourth_power-4.0*mean*sum_cubes+6.0*mean*mean*sum_squares-
-        3.0*mean*mean*mean*mean;
-      *kurtosis/=standard_deviation*standard_deviation*standard_deviation*
-        standard_deviation;
-      *kurtosis-=3.0;
-      *skewness=sum_cubes-3.0*mean*sum_squares+2.0*mean*mean*mean;
-      *skewness/=standard_deviation*standard_deviation*standard_deviation;
+      perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
+        perceptual_hash);
+      return((ChannelPerceptualHash *) NULL);
     }
-  return(status);
+  moments=GetImageMoments(hash_image,exception);
+  hash_image=DestroyImage(hash_image);
+  if (moments == (ChannelMoments *) NULL)
+    {
+      perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
+        perceptual_hash);
+      return((ChannelPerceptualHash *) NULL);
+    }
+  for (channel=0; channel <= MaxPixelChannels; channel++)
+    for (i=0; i < 7; i++)
+      perceptual_hash[channel].Q[i]=(-MagickLog10(moments[channel].I[i]));
+  moments=(ChannelMoments *) RelinquishMagickMemory(moments);
+  return(perceptual_hash);
 }
 \f
 /*
@@ -1409,7 +1823,7 @@ MagickExport MagickBooleanType GetImageRange(const Image *image,double *minima,
       register ssize_t
         i;
 
-      if (GetPixelMask(image,p) != 0)
+      if (GetPixelReadMask(image,p) == 0)
         {
           p+=GetPixelChannels(image);
           continue;
@@ -1495,7 +1909,7 @@ static size_t GetImageChannels(const Image *image)
   {
     PixelChannel channel=GetPixelChannelChannel(image,i);
     PixelTrait traits=GetPixelChannelTraits(image,channel);
-    if ((traits & UpdatePixelTrait) != 0)
+    if (traits != UndefinedPixelTrait)
       channels++;
   }
   return(channels);
@@ -1530,14 +1944,14 @@ MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
   channel_statistics=(ChannelStatistics *) AcquireQuantumMemory(
     MaxPixelChannels+1,sizeof(*channel_statistics));
   if (channel_statistics == (ChannelStatistics *) NULL)
-    ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
+    return(channel_statistics);
   (void) ResetMagickMemory(channel_statistics,0,(MaxPixelChannels+1)*
     sizeof(*channel_statistics));
   for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
   {
     channel_statistics[i].depth=1;
-    channel_statistics[i].maxima=(-MagickHuge);
-    channel_statistics[i].minima=MagickHuge;
+    channel_statistics[i].maxima=(-MagickMaximumValue);
+    channel_statistics[i].minima=MagickMaximumValue;
   }
   for (y=0; y < (ssize_t) image->rows; y++)
   {
@@ -1555,7 +1969,7 @@ MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
       register ssize_t
         i;
 
-      if (GetPixelMask(image,p) != 0)
+      if (GetPixelReadMask(image,p) == 0)
         {
           p+=GetPixelChannels(image);
           continue;
@@ -1611,9 +2025,7 @@ MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
   }
   for (i=0; i < (ssize_t) MaxPixelChannels; i++)
   {
-    channel_statistics[CompositePixelChannel].depth=(size_t) EvaluateMax(
-      (double) channel_statistics[CompositePixelChannel].depth,(double)
-      channel_statistics[i].depth);
+    channel_statistics[CompositePixelChannel].area+=channel_statistics[i].area;
     channel_statistics[CompositePixelChannel].minima=MagickMin(
       channel_statistics[CompositePixelChannel].minima,
       channel_statistics[i].minima);
@@ -1636,6 +2048,7 @@ MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
       channel_statistics[i].mean;
   }
   channels=GetImageChannels(image);
+  channel_statistics[CompositePixelChannel].area/=channels;
   channel_statistics[CompositePixelChannel].sum/=channels;
   channel_statistics[CompositePixelChannel].sum_squared/=channels;
   channel_statistics[CompositePixelChannel].sum_cubed/=channels;
@@ -1651,6 +2064,8 @@ MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
     double
       standard_deviation;
 
+    if (channel_statistics[i].standard_deviation == 0.0)
+      continue;
     standard_deviation=PerceptibleReciprocal(
       channel_statistics[i].standard_deviation);
     channel_statistics[i].skewness=(channel_statistics[i].sum_cubed-3.0*
@@ -1666,6 +2081,9 @@ MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
       channel_statistics[i].mean)*(standard_deviation*standard_deviation*
       standard_deviation*standard_deviation)-3.0;
   }
+  if (y < (ssize_t) image->rows)
+    channel_statistics=(ChannelStatistics *) RelinquishMagickMemory(
+      channel_statistics);
   return(channel_statistics);
 }
 \f
@@ -1820,7 +2238,7 @@ MagickExport Image *PolynomialImage(const Image *images,
         register ssize_t
           i;
 
-        if (GetPixelMask(next,p) != 0)
+        if (GetPixelReadMask(next,p) == 0)
           {
             p+=GetPixelChannels(next);
             continue;
@@ -1854,7 +2272,7 @@ MagickExport Image *PolynomialImage(const Image *images,
       register ssize_t
         i;
 
-      if (GetPixelMask(image,q) != 0)
+      if (GetPixelReadMask(image,q) == 0)
         {
           q+=GetPixelChannels(image);
           continue;
@@ -2310,8 +2728,7 @@ static inline void GetStandardDeviationPixelList(PixelList *pixel_list,
   *pixel=ScaleShortToQuantum((unsigned short) sqrt(sum_squared-(sum*sum)));
 }
 
-static inline void InsertPixelList(const Image *image,const Quantum pixel,
-  PixelList *pixel_list)
+static inline void InsertPixelList(const Quantum pixel,PixelList *pixel_list)
 {
   size_t
     signature;
@@ -2480,7 +2897,7 @@ MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
             (statistic_traits == UndefinedPixelTrait))
           continue;
         if (((statistic_traits & CopyPixelTrait) != 0) ||
-            (GetPixelMask(image,p) != 0))
+            (GetPixelReadMask(image,p) == 0))
           {
             SetPixelChannel(statistic_image,channel,p[center+i],q);
             continue;
@@ -2491,10 +2908,10 @@ MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
         {
           for (u=0; u < (ssize_t) MagickMax(width,1); u++)
           {
-            InsertPixelList(image,pixels[i],pixel_list[id]);
+            InsertPixelList(pixels[i],pixel_list[id]);
             pixels+=GetPixelChannels(image);
           }
-          pixels+=image->columns*GetPixelChannels(image);
+          pixels+=(image->columns-1)*GetPixelChannels(image);
         }
         switch (type)
         {