Support similarity and dissimilarity metrics

diff --git a/MagickCore/cache.c b/MagickCore/cache.c
index eb2ac1684b2bd0bd1cfd295afb18aea556ffb579..40f98e74f34b367071dca27187f7ef1896a77f5b 100644 (file)
--- a/MagickCore/cache.c
+++ b/MagickCore/cache.c
@@ -2293,9 +2293,9 @@ MagickPrivate ClassType GetPixelCacheStorageClass(const Cache cache)
 %
 %    o image: the image.
 %
-%    o width: the optimize cache tile width in pixels.
+%    o width: the optimized cache tile width in pixels.
 %
-%    o height: the optimize cache tile height in pixels.
+%    o height: the optimized cache tile height in pixels.
 %
 */
 MagickPrivate void GetPixelCacheTileSize(const Image *image,size_t *width,

diff --git a/MagickCore/compare.c b/MagickCore/compare.c
index 48a7425d934f8b8646977d256c402095c112a415..e53386f3631fd208de4f3db5ce1237954ec06d40 100644 (file)
--- a/MagickCore/compare.c
+++ b/MagickCore/compare.c
@@ -67,6 +67,7 @@
 #include "MagickCore/resource_.h"
 #include "MagickCore/string_.h"
 #include "MagickCore/statistic.h"
+#include "MagickCore/string-private.h"
 #include "MagickCore/thread-private.h"
 #include "MagickCore/transform.h"
 #include "MagickCore/utility.h"
@@ -1316,15 +1317,182 @@ static MagickBooleanType GetRootMeanSquaredDistortion(const Image *image,
 static MagickBooleanType GetStructuralSimilarityDistortion(const Image *image,
   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
 {
+#define SSIMBlocksize  8
+#define SSIMK1  0.01
+#define SSIMK2  0.03
+
+  CacheView
+    *image_view,
+    *reconstruct_view;
+
+  const char
+    *artifact;
+
+  double
+    c1,
+    c2;
+
   MagickBooleanType
     status;
 
   register ssize_t
     i;
 
-  status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
-  for (i=0; i <= MaxPixelChannels; i++)
-    distortion[i]=sqrt(distortion[i]);
+  size_t
+    blocksize,
+    n;
+
+  ssize_t
+    tile_y;
+
+  /*
+    Compute structural similarity index.
+  */
+  blocksize=SSIMBlocksize;
+  artifact=GetImageArtifact(image,"compare:blocksize");
+  if (artifact != (const char *) NULL)
+    if (StringToDouble(artifact,(char **) NULL) > 0.0)
+      blocksize=(size_t) StringToDouble(artifact,(char **) NULL);
+  c1=pow(SSIMK1*QuantumRange,2.0);
+  artifact=GetImageArtifact(image,"compare:k1");
+  if (artifact != (const char *) NULL)
+    c1=pow(StringToDouble(artifact,(char **) NULL)*QuantumRange,2.0);
+  c2=pow(SSIMK2*QuantumRange,2.0);
+  artifact=GetImageArtifact(image,"compare:k2");
+  if (artifact != (const char *) NULL)
+    c2=pow(StringToDouble(artifact,(char **) NULL)*QuantumRange,2.0);
+  status=MagickTrue;
+  n=0;
+  image_view=AcquireVirtualCacheView(image,exception);
+  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
+#if defined(MMAGICKCORE_OPENMP_SUPPORT)
+  #pragma omp parallel for schedule(static,4) shared(status) \
+    magick_threads(image,image,1,1)
+#endif
+  for (tile_y=0; tile_y < (ssize_t) image->rows; tile_y+=(ssize_t) blocksize)
+  {
+    double
+      channel_distortion[MaxPixelChannels+1];
+
+    register ssize_t
+      tile_x;
+
+    ssize_t
+      j;
+
+    if (status == MagickFalse)
+      continue;
+    (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
+    tile_x=0;
+    for ( ; tile_x < (ssize_t) image->columns; tile_x+=(ssize_t) blocksize)
+    {
+      double
+        image_reconstruct_sum[MaxPixelChannels+1],
+        image_sum[MaxPixelChannels+1],
+        image_sum_squared[MaxPixelChannels+1],
+        reconstruct_sum[MaxPixelChannels+1],
+        reconstruct_sum_squared[MaxPixelChannels+1];
+
+      register const Quantum
+        *magick_restrict p,
+        *magick_restrict q;
+
+      register ssize_t
+        y;
+
+      p=GetCacheViewVirtualPixels(image_view,tile_x,tile_y,blocksize,blocksize,
+        exception);
+      q=GetCacheViewVirtualPixels(reconstruct_view,tile_x,tile_y,blocksize,
+        blocksize,exception);
+      if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
+        {
+          status=MagickFalse;
+          break;
+        }
+      (void) ResetMagickMemory(image_sum,0,(MaxPixelChannels+1)*
+        sizeof(*image_sum));
+      (void) ResetMagickMemory(image_sum_squared,0,(MaxPixelChannels+1)*
+        sizeof(*image_sum_squared));
+      (void) ResetMagickMemory(reconstruct_sum,0,(MaxPixelChannels+1)*
+        sizeof(*reconstruct_sum));
+      (void) ResetMagickMemory(reconstruct_sum_squared,0,(MaxPixelChannels+1)*
+        sizeof(*reconstruct_sum_squared));
+      (void) ResetMagickMemory(image_reconstruct_sum,0,(MaxPixelChannels+1)*
+        sizeof(*image_reconstruct_sum));
+      for (y=0; y < (ssize_t) blocksize; y++)
+      {
+        register ssize_t
+          x;
+
+        for (x=0; x < (ssize_t) blocksize; x++)
+        {
+          register ssize_t
+            i;
+
+          for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
+          {
+            PixelChannel channel=GetPixelChannelChannel(image,i);
+            PixelTrait traits=GetPixelChannelTraits(image,channel);
+            PixelTrait reconstruct_traits=GetPixelChannelTraits(
+              reconstruct_image,channel);
+            if ((traits == UndefinedPixelTrait) ||
+                (reconstruct_traits == UndefinedPixelTrait) ||
+                ((reconstruct_traits & UpdatePixelTrait) == 0))
+              continue;
+            image_sum[i]+=p[i];
+            image_sum_squared[i]+=(p[i]*p[i]);
+            reconstruct_sum[i]+=GetPixelChannel(reconstruct_image,channel,q);
+            reconstruct_sum_squared[i]+=(GetPixelChannel(reconstruct_image,
+              channel,q)*GetPixelChannel(reconstruct_image,channel,q));
+            image_reconstruct_sum[i]+=(p[i]*GetPixelChannel(reconstruct_image,
+              channel,q));
+          }
+          p+=GetPixelChannels(image);
+          q+=GetPixelChannels(reconstruct_image);
+        }
+        for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
+        {
+          double
+            image_mean,
+            image_variance,
+            reconstruct_mean,
+            reconstruct_variance,
+            image_reconstruct_covariance;
+
+          image_mean=image_sum[i]/(blocksize*blocksize);
+          image_variance=image_sum_squared[i]/(blocksize*blocksize)-
+            (image_mean*image_mean);
+          reconstruct_mean=reconstruct_sum[i]/(blocksize*blocksize);
+          reconstruct_variance=reconstruct_sum_squared[i]/(blocksize*blocksize)-
+            (reconstruct_mean*reconstruct_mean);
+          image_reconstruct_covariance=image_reconstruct_sum[i]/
+            (blocksize*blocksize)-(image_mean*reconstruct_mean);
+          channel_distortion[i]+=((2.0*image_mean*reconstruct_mean+c1)*(2.0*
+            image_reconstruct_covariance+c2))/((image_mean*image_mean+
+            reconstruct_mean*reconstruct_mean+c1)*(image_variance+
+            reconstruct_variance+c2));
+        }
+        n++;
+      }
+    }
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+    #pragma omp critical (MagickCore_GetStructuralSimilarityDistortion)
+#endif
+    for (j=0; j <= MaxPixelChannels; j++)
+      distortion[j]+=channel_distortion[j];
+  }
+  image_view=DestroyCacheView(image_view);
+  reconstruct_view=DestroyCacheView(reconstruct_view);
+  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
+  {
+    PixelChannel channel=GetPixelChannelChannel(image,i);
+    PixelTrait traits=GetPixelChannelTraits(image,channel);
+    if ((traits == UndefinedPixelTrait) || ((traits & UpdatePixelTrait) == 0))
+      continue;
+    distortion[i]/=(double) n;
+    distortion[CompositePixelChannel]+=distortion[i];
+  }
+  distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
   return(status);
 }
 
@@ -1334,9 +1502,13 @@ static MagickBooleanType GetStructuralDisimilarityDistortion(const Image *image,
   MagickBooleanType
     status;
 
+  register ssize_t
+    i;
+
   status=GetStructuralSimilarityDistortion(image,reconstruct_image,
     distortion,exception);
-  *distortion=(1.0-(*distortion))/2.0;
+  for (i=0; i <= MaxPixelChannels; i++)
+    distortion[i]=(1.0-(distortion[i]))/2.0;
   return(status);
 }
 
@@ -1583,6 +1755,18 @@ MagickExport double *GetImageDistortions(Image *image,
         channel_distortion,exception);
       break;
     }
+    case StructuralSimilarityErrorMetric:
+    {
+      status=GetStructuralSimilarityDistortion(image,reconstruct_image,
+        channel_distortion,exception);
+      break;
+    }
+    case StructuralDissimilarityErrorMetric:
+    {
+      status=GetStructuralDisimilarityDistortion(image,reconstruct_image,
+        channel_distortion,exception);
+      break;
+    }
   }
   if (status == MagickFalse)
     {

diff --git a/MagickCore/option.c b/MagickCore/option.c
index 63d99946e21185ae60e90241779883052403709c..267c8a5785babd1b28a7adda27e3daba3b4e4ae7 100644 (file)
--- a/MagickCore/option.c
+++ b/MagickCore/option.c
@@ -1434,6 +1434,7 @@ static const OptionInfo
   {
     { "Undefined", UndefinedErrorMetric, UndefinedOptionFlag, MagickTrue },
     { "AE", AbsoluteErrorMetric, UndefinedOptionFlag, MagickFalse },
+    { "DSSIM", StructuralDissimilarityErrorMetric, UndefinedOptionFlag, MagickFalse },
     { "Fuzz", FuzzErrorMetric, UndefinedOptionFlag, MagickFalse },
     { "MAE", MeanAbsoluteErrorMetric, UndefinedOptionFlag, MagickFalse },
     { "MEPP", MeanErrorPerPixelErrorMetric, UndefinedOptionFlag, MagickFalse },

diff --git a/MagickWand/compare.c b/MagickWand/compare.c
index 3648809a6ba7530e2b232624ea8885657acb9941..06bf41be4141aab3f29873b9f99fcc703f8e2b85 100644 (file)
--- a/MagickWand/compare.c
+++ b/MagickWand/compare.c
@@ -1239,6 +1239,8 @@ WandExport MagickBooleanType CompareImagesCommand(ImageInfo *image_info,
             case NormalizedCrossCorrelationErrorMetric:
             case PeakSignalToNoiseRatioErrorMetric:
             case PerceptualHashErrorMetric:
+            case StructuralSimilarityErrorMetric:
+            case StructuralDissimilarityErrorMetric:
             {
               (void) FormatLocaleFile(stderr,"%g",distortion);
               break;
@@ -1340,6 +1342,8 @@ WandExport MagickBooleanType CompareImagesCommand(ImageInfo *image_info,
             case NormalizedCrossCorrelationErrorMetric:
             case PeakSignalToNoiseRatioErrorMetric:
             case PerceptualHashErrorMetric:
+            case StructuralSimilarityErrorMetric:
+            case StructuralDissimilarityErrorMetric:
             {
               switch (image->colorspace)
               {