]> granicus.if.org Git - imagemagick/blobdiff - MagickCore/statistic.c
(no commit message)
[imagemagick] / MagickCore / statistic.c
index ccc3ca13514e714d7dd36d9e11662f6f45ad45cc..42ea5facbb4b210791e6c34bee856980f04092fb 100644 (file)
@@ -1348,7 +1348,8 @@ MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean,
 %                                                                             %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-%  GetImageMoments() returns the moments of one or more image channels.
+%  GetImageMoments() returns the normalized moments of one or more image
+%  channels.
 %
 %  The format of the GetImageMoments method is:
 %
@@ -1384,8 +1385,7 @@ MagickExport ChannelMoments *GetImageMoments(const Image *image,
     M20[MaxPixelChannels+1],
     M21[MaxPixelChannels+1],
     M22[MaxPixelChannels+1],
-    M30[MaxPixelChannels+1],
-    scale;
+    M30[MaxPixelChannels+1];
 
   PointInfo
     centroid[MaxPixelChannels+1];
@@ -1416,7 +1416,6 @@ MagickExport ChannelMoments *GetImageMoments(const Image *image,
   (void) ResetMagickMemory(M21,0,sizeof(M21));
   (void) ResetMagickMemory(M22,0,sizeof(M22));
   (void) ResetMagickMemory(M30,0,sizeof(M30));
-  scale=(double) ((1UL << image->depth)-1)/QuantumRange;
   image_view=AcquireVirtualCacheView(image,exception);
   for (y=0; y < (ssize_t) image->rows; y++)
   {
@@ -1450,9 +1449,12 @@ MagickExport ChannelMoments *GetImageMoments(const Image *image,
           continue;
         if ((traits & UpdatePixelTrait) == 0)
           continue;
-        M00[channel]+=scale*p[i];
-        M10[channel]+=x*scale*p[i];
-        M01[channel]+=y*scale*p[i];
+        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);
     }
@@ -1462,10 +1464,14 @@ MagickExport ChannelMoments *GetImageMoments(const Image *image,
     /*
        Compute center of mass (centroid).
     */
-    centroid[channel].x=M10[channel]/image->columns/2.0;
-    centroid[channel].y=M01[channel]/image->rows/2.0;
-    if (fabs(M00[channel]) < MagickEpsilon)
-      continue;
+    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];
   }
@@ -1502,32 +1508,57 @@ MagickExport ChannelMoments *GetImageMoments(const Image *image,
         if ((traits & UpdatePixelTrait) == 0)
           continue;
         M11[channel]+=(x-centroid[channel].x)*(y-centroid[channel].y)*
-          scale*p[i];
+          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)*
-          scale*p[i];
+          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)*
-          scale*p[i];
+          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)*scale*p[i];
+          (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)*scale*p[i];
+          (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)*scale*p[i];
+          (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)*scale*p[i];
+          (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)*scale*p[i];
+          (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.
     */
-    if (fabs(M00[channel]) < MagickEpsilon)
-      continue;
     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]+
@@ -1536,21 +1567,19 @@ MagickExport ChannelMoments *GetImageMoments(const Image *image,
       ((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])));
+      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));
+      (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);
+      channel_moments[channel].ellipse_axis.y+MagickEpsilon);
   }
   for (channel=0; channel <= MaxPixelChannels; channel++)
   {
     /*
       Normalize image moments.
     */
-    if (fabs(M00[channel]) < MagickEpsilon)
-      continue;
     M10[channel]=0.0;
     M01[channel]=0.0;
     M11[channel]/=pow(M00[channel],1.0+(1.0+1.0)/2.0);
@@ -1569,8 +1598,6 @@ MagickExport ChannelMoments *GetImageMoments(const Image *image,
     /*
       Compute Hu invariant moments.
     */
-    if (fabs(M00[channel]) < MagickEpsilon)
-      continue;
     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];
@@ -1613,6 +1640,118 @@ MagickExport ChannelMoments *GetImageMoments(const Image *image,
 %                                                                             %
 %                                                                             %
 %                                                                             %
+%   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)
+    {
+      perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
+        perceptual_hash);
+      return((ChannelPerceptualHash *) NULL);
+    }
+  hash_image->depth=8;
+  status=TransformImageColorspace(hash_image,HCLpColorspace,exception);
+  if (status == MagickFalse)
+    {
+      perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory(
+        perceptual_hash);
+      return((ChannelPerceptualHash *) NULL);
+    }
+  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
+/*
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                             %
+%                                                                             %
+%                                                                             %
 %   G e t I m a g e R a n g e                                                 %
 %                                                                             %
 %                                                                             %
@@ -1811,8 +1950,8 @@ MagickExport ChannelStatistics *GetImageStatistics(const Image *image,
   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++)
   {