]> granicus.if.org Git - imagemagick/blobdiff - MagickCore/compare.c
(no commit message)
[imagemagick] / MagickCore / compare.c
index 10f56cbbfd5a2a031e352a902113d4ca12c5bb4f..fc5702022de7a102b4c20795760a2c086998d75c 100644 (file)
@@ -17,7 +17,7 @@
 %                               December 2003                                 %
 %                                                                             %
 %                                                                             %
-%  Copyright 1999-2011 ImageMagick Studio LLC, a non-profit organization      %
+%  Copyright 1999-2013 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  %
@@ -43,6 +43,7 @@
 #include "MagickCore/studio.h"
 #include "MagickCore/artifact.h"
 #include "MagickCore/cache-view.h"
+#include "MagickCore/channel.h"
 #include "MagickCore/client.h"
 #include "MagickCore/color.h"
 #include "MagickCore/color-private.h"
@@ -64,6 +65,7 @@
 #include "MagickCore/resource_.h"
 #include "MagickCore/string_.h"
 #include "MagickCore/statistic.h"
+#include "MagickCore/thread-private.h"
 #include "MagickCore/transform.h"
 #include "MagickCore/utility.h"
 #include "MagickCore/version.h"
@@ -79,7 +81,7 @@
 %                                                                             %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-%  CompareImages() compares one or more pixel channels of an image to a 
+%  CompareImages() compares one or more pixel channels of an image to a
 %  reconstructed image and returns the difference image.
 %
 %  The format of the CompareImages method is:
@@ -161,32 +163,24 @@ MagickExport Image *CompareImages(Image *image,const Image *reconstruct_image,
       return((Image *) NULL);
     }
   (void) SetImageAlphaChannel(highlight_image,OpaqueAlphaChannel,exception);
-  (void) QueryMagickColorCompliance("#f1001ecc",AllCompliance,&highlight,
-    exception);
+  (void) QueryColorCompliance("#f1001ecc",AllCompliance,&highlight,exception);
   artifact=GetImageArtifact(image,"highlight-color");
   if (artifact != (const char *) NULL)
-    (void) QueryMagickColorCompliance(artifact,AllCompliance,&highlight,
-      exception);
-  (void) QueryMagickColorCompliance("#ffffffcc",AllCompliance,&lowlight,
-    exception);
+    (void) QueryColorCompliance(artifact,AllCompliance,&highlight,exception);
+  (void) QueryColorCompliance("#ffffffcc",AllCompliance,&lowlight,exception);
   artifact=GetImageArtifact(image,"lowlight-color");
   if (artifact != (const char *) NULL)
-    (void) QueryMagickColorCompliance(artifact,AllCompliance,&lowlight,
-      exception);
-  if (highlight_image->colorspace == CMYKColorspace)
-    {
-      ConvertRGBToCMYK(&highlight);
-      ConvertRGBToCMYK(&lowlight);
-    }
+    (void) QueryColorCompliance(artifact,AllCompliance,&lowlight,exception);
   /*
     Generate difference image.
   */
   status=MagickTrue;
-  image_view=AcquireCacheView(image);
-  reconstruct_view=AcquireCacheView(reconstruct_image);
-  highlight_view=AcquireCacheView(highlight_image);
+  image_view=AcquireVirtualCacheView(image,exception);
+  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
+  highlight_view=AcquireAuthenticCacheView(highlight_image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(dynamic,4) shared(status)
+  #pragma omp parallel for schedule(static,4) shared(status) \
+    magick_threads(image,highlight_image,image->rows,1)
 #endif
   for (y=0; y < (ssize_t) image->rows; y++)
   {
@@ -218,42 +212,48 @@ MagickExport Image *CompareImages(Image *image,const Image *reconstruct_image,
       }
     for (x=0; x < (ssize_t) image->columns; x++)
     {
+      double
+        Da,
+        Sa;
+
       MagickStatusType
         difference;
 
       register ssize_t
         i;
 
+      if (GetPixelReadMask(image,p) == 0)
+        {
+          SetPixelInfoPixel(highlight_image,&lowlight,r);
+          p+=GetPixelChannels(image);
+          q+=GetPixelChannels(reconstruct_image);
+          r+=GetPixelChannels(highlight_image);
+          continue;
+        }
       difference=MagickFalse;
+      Sa=QuantumScale*GetPixelAlpha(image,p);
+      Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
       {
-        MagickRealType
+        double
           distance;
 
-        PixelChannel
-          channel;
-
-        PixelTrait
-          reconstruct_traits,
-          traits;
-
-        traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
-        channel=GetPixelChannelMapChannel(image,(PixelChannel) i);
-        reconstruct_traits=GetPixelChannelMapTraits(reconstruct_image,channel);
+        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 == UndefinedPixelTrait) ||
+            ((reconstruct_traits & UpdatePixelTrait) == 0))
           continue;
-        if ((reconstruct_traits & UpdatePixelTrait) == 0)
-          continue;
-        distance=p[i]-(MagickRealType)
-          GetPixelChannel(reconstruct_image,channel,q);
+        distance=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
         if (fabs((double) distance) >= MagickEpsilon)
           difference=MagickTrue;
       }
       if (difference == MagickFalse)
-        SetPixelPixelInfo(highlight_image,&lowlight,r);
+        SetPixelInfoPixel(highlight_image,&lowlight,r);
       else
-        SetPixelPixelInfo(highlight_image,&highlight,r);
+        SetPixelInfoPixel(highlight_image,&highlight,r);
       p+=GetPixelChannels(image);
       q+=GetPixelChannels(reconstruct_image);
       r+=GetPixelChannels(highlight_image);
@@ -265,7 +265,8 @@ MagickExport Image *CompareImages(Image *image,const Image *reconstruct_image,
   highlight_view=DestroyCacheView(highlight_view);
   reconstruct_view=DestroyCacheView(reconstruct_view);
   image_view=DestroyCacheView(image_view);
-  (void) CompositeImage(difference_image,image->compose,highlight_image,0,0);
+  (void) CompositeImage(difference_image,highlight_image,image->compose,
+    MagickTrue,0,0,exception);
   highlight_image=DestroyImage(highlight_image);
   if (status == MagickFalse)
     difference_image=DestroyImage(difference_image);
@@ -283,10 +284,10 @@ MagickExport Image *CompareImages(Image *image,const Image *reconstruct_image,
 %                                                                             %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-%  GetImageDistortion() compares one or more pixel channels of an image to a 
+%  GetImageDistortion() compares one or more pixel channels of an image to a
 %  reconstructed image and returns the specified distortion metric.
 %
-%  The format of the CompareImages method is:
+%  The format of the GetImageDistortion method is:
 %
 %      MagickBooleanType GetImageDistortion(const Image *image,
 %        const Image *reconstruct_image,const MetricType metric,
@@ -306,6 +307,13 @@ MagickExport Image *CompareImages(Image *image,const Image *reconstruct_image,
 %
 */
 
+static inline double MagickMax(const double x,const double y)
+{
+  if (x > y)
+    return(x);
+  return(y);
+}
+
 static MagickBooleanType GetAbsoluteDistortion(const Image *image,
   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
 {
@@ -313,6 +321,9 @@ static MagickBooleanType GetAbsoluteDistortion(const Image *image,
     *image_view,
     *reconstruct_view;
 
+  double
+    fuzz;
+
   MagickBooleanType
     status;
 
@@ -323,10 +334,21 @@ static MagickBooleanType GetAbsoluteDistortion(const Image *image,
     Compute the absolute difference in pixels between two images.
   */
   status=MagickTrue;
-  image_view=AcquireCacheView(image);
-  reconstruct_view=AcquireCacheView(reconstruct_image);
+  if (image->fuzz == 0.0)
+    fuzz=MagickMax(reconstruct_image->fuzz,MagickSQ1_2)*
+      MagickMax(reconstruct_image->fuzz,MagickSQ1_2);
+  else
+    if (reconstruct_image->fuzz == 0.0)
+      fuzz=MagickMax(image->fuzz,MagickSQ1_2)*
+        MagickMax(image->fuzz,MagickSQ1_2);
+    else
+      fuzz=MagickMax(image->fuzz,MagickSQ1_2)*
+        MagickMax(reconstruct_image->fuzz,MagickSQ1_2);
+  image_view=AcquireVirtualCacheView(image,exception);
+  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(dynamic,4) shared(status)
+  #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++)
   {
@@ -354,43 +376,55 @@ static MagickBooleanType GetAbsoluteDistortion(const Image *image,
     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
     for (x=0; x < (ssize_t) image->columns; x++)
     {
+      double
+        Da,
+        Sa;
+
       MagickBooleanType
         difference;
 
       register ssize_t
         i;
 
+      if (GetPixelReadMask(image,p) == 0)
+        {
+          p+=GetPixelChannels(image);
+          q+=GetPixelChannels(reconstruct_image);
+          continue;
+        }
       difference=MagickFalse;
+      Sa=QuantumScale*GetPixelAlpha(image,p);
+      Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
       {
-        PixelChannel
-          channel;
-
-        PixelTrait
-          reconstruct_traits,
-          traits;
+        double
+          distance;
 
-        traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
-        channel=GetPixelChannelMapChannel(image,(PixelChannel) i);
-        reconstruct_traits=GetPixelChannelMapTraits(reconstruct_image,channel);
+        PixelChannel channel=GetPixelChannelChannel(image,i);
+        PixelTrait traits=GetPixelChannelTraits(image,channel);
+        PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
+          channel);
         if ((traits == UndefinedPixelTrait) ||
-            (reconstruct_traits == UndefinedPixelTrait))
-          continue;
-        if ((reconstruct_traits & UpdatePixelTrait) == 0)
+            (reconstruct_traits == UndefinedPixelTrait) ||
+            ((reconstruct_traits & UpdatePixelTrait) == 0))
           continue;
-        if (p[i] != GetPixelChannel(reconstruct_image,channel,q))
-          difference=MagickTrue;
+        distance=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
+        if ((distance*distance) > fuzz)
+          {
+            difference=MagickTrue;
+            break;
+          }
       }
       if (difference != MagickFalse)
         {
           channel_distortion[i]++;
-          channel_distortion[MaxPixelChannels]++;
+          channel_distortion[CompositePixelChannel]++;
         }
       p+=GetPixelChannels(image);
       q+=GetPixelChannels(reconstruct_image);
     }
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_GetAbsoluteError)
+    #pragma omp critical (MagickCore_GetAbsoluteError)
 #endif
     for (i=0; i <= MaxPixelChannels; i++)
       distortion[i]+=channel_distortion[i];
@@ -411,10 +445,8 @@ static size_t GetImageChannels(const Image *image)
   channels=0;
   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
   {
-    PixelTrait
-      traits;
-
-    traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
+    PixelChannel channel=GetPixelChannelChannel(image,i);
+    PixelTrait traits=GetPixelChannelTraits(image,channel);
     if ((traits & UpdatePixelTrait) != 0)
       channels++;
   }
@@ -438,10 +470,11 @@ static MagickBooleanType GetFuzzDistortion(const Image *image,
     y;
 
   status=MagickTrue;
-  image_view=AcquireCacheView(image);
-  reconstruct_view=AcquireCacheView(reconstruct_image);
+  image_view=AcquireVirtualCacheView(image,exception);
+  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(dynamic,4) shared(status)
+  #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++)
   {
@@ -469,40 +502,44 @@ static MagickBooleanType GetFuzzDistortion(const Image *image,
     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
     for (x=0; x < (ssize_t) image->columns; x++)
     {
+      double
+        Da,
+        Sa;
+
       register ssize_t
         i;
 
+      if (GetPixelReadMask(image,p) == 0)
+        {
+          p+=GetPixelChannels(image);
+          q+=GetPixelChannels(reconstruct_image);
+          continue;
+        }
+      Sa=QuantumScale*GetPixelAlpha(image,p);
+      Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
       {
-        MagickRealType
+        double
           distance;
 
-        PixelChannel
-          channel;
-
-        PixelTrait
-          reconstruct_traits,
-          traits;
-
-        traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
-        channel=GetPixelChannelMapChannel(image,(PixelChannel) i);
-        reconstruct_traits=GetPixelChannelMapTraits(reconstruct_image,channel);
+        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 == UndefinedPixelTrait) ||
+            ((reconstruct_traits & UpdatePixelTrait) == 0))
           continue;
-        if ((reconstruct_traits & UpdatePixelTrait) == 0)
-          continue;
-        distance=QuantumScale*(p[i]-(MagickRealType) GetPixelChannel(
-          reconstruct_image,channel,q));
-        distance*=distance;
-        channel_distortion[i]+=distance;
-        channel_distortion[MaxPixelChannels]+=distance;
+        distance=QuantumScale*(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
+          channel,q));
+        channel_distortion[i]+=distance*distance;
+        channel_distortion[CompositePixelChannel]+=distance*distance;
       }
       p+=GetPixelChannels(image);
       q+=GetPixelChannels(reconstruct_image);
     }
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_GetMeanSquaredError)
+    #pragma omp critical (MagickCore_GetFuzzDistortion)
 #endif
     for (i=0; i <= MaxPixelChannels; i++)
       distortion[i]+=channel_distortion[i];
@@ -511,8 +548,8 @@ static MagickBooleanType GetFuzzDistortion(const Image *image,
   image_view=DestroyCacheView(image_view);
   for (i=0; i <= MaxPixelChannels; i++)
     distortion[i]/=((double) image->columns*image->rows);
-  distortion[MaxPixelChannels]/=(double) GetImageChannels(image);
-  distortion[MaxPixelChannels]=sqrt(distortion[MaxPixelChannels]);
+  distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
+  distortion[CompositePixelChannel]=sqrt(distortion[CompositePixelChannel]);
   return(status);
 }
 
@@ -533,10 +570,11 @@ static MagickBooleanType GetMeanAbsoluteDistortion(const Image *image,
     y;
 
   status=MagickTrue;
-  image_view=AcquireCacheView(image);
-  reconstruct_view=AcquireCacheView(reconstruct_image);
+  image_view=AcquireVirtualCacheView(image,exception);
+  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(dynamic,4) shared(status)
+  #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++)
   {
@@ -564,39 +602,44 @@ static MagickBooleanType GetMeanAbsoluteDistortion(const Image *image,
     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
     for (x=0; x < (ssize_t) image->columns; x++)
     {
+      double
+        Da,
+        Sa;
+
       register ssize_t
         i;
 
+      if (GetPixelReadMask(image,p) == 0)
+        {
+          p+=GetPixelChannels(image);
+          q+=GetPixelChannels(reconstruct_image);
+          continue;
+        }
+      Sa=QuantumScale*GetPixelAlpha(image,p);
+      Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
       {
-        MagickRealType
+        double
           distance;
 
-        PixelChannel
-          channel;
-
-        PixelTrait
-          reconstruct_traits,
-          traits;
-
-        traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
-        channel=GetPixelChannelMapChannel(image,(PixelChannel) i);
-        reconstruct_traits=GetPixelChannelMapTraits(reconstruct_image,channel);
+        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 == UndefinedPixelTrait) ||
+            ((reconstruct_traits & UpdatePixelTrait) == 0))
           continue;
-        if ((reconstruct_traits & UpdatePixelTrait) == 0)
-          continue;
-        distance=QuantumScale*fabs(p[i]-(MagickRealType) GetPixelChannel(
-          reconstruct_image,channel,q));
+        distance=QuantumScale*fabs(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
+          channel,q));
         channel_distortion[i]+=distance;
-        channel_distortion[MaxPixelChannels]+=distance;
+        channel_distortion[CompositePixelChannel]+=distance;
       }
       p+=GetPixelChannels(image);
       q+=GetPixelChannels(reconstruct_image);
     }
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_GetMeanAbsoluteError)
+    #pragma omp critical (MagickCore_GetMeanAbsoluteError)
 #endif
     for (i=0; i <= MaxPixelChannels; i++)
       distortion[i]+=channel_distortion[i];
@@ -605,7 +648,7 @@ static MagickBooleanType GetMeanAbsoluteDistortion(const Image *image,
   image_view=DestroyCacheView(image_view);
   for (i=0; i <= MaxPixelChannels; i++)
     distortion[i]/=((double) image->columns*image->rows);
-  distortion[MaxPixelChannels]/=(double) GetImageChannels(image);
+  distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
   return(status);
 }
 
@@ -619,10 +662,8 @@ static MagickBooleanType GetMeanErrorPerPixel(Image *image,
   MagickBooleanType
     status;
 
-  MagickRealType
-    alpha,
+  double
     area,
-    beta,
     maximum_error,
     mean_error;
 
@@ -630,13 +671,11 @@ static MagickBooleanType GetMeanErrorPerPixel(Image *image,
     y;
 
   status=MagickTrue;
-  alpha=1.0;
-  beta=1.0;
   area=0.0;
   maximum_error=0.0;
   mean_error=0.0;
-  image_view=AcquireCacheView(image);
-  reconstruct_view=AcquireCacheView(reconstruct_image);
+  image_view=AcquireVirtualCacheView(image,exception);
+  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
   for (y=0; y < (ssize_t) image->rows; y++)
   {
     register const Quantum
@@ -656,33 +695,38 @@ static MagickBooleanType GetMeanErrorPerPixel(Image *image,
       }
     for (x=0; x < (ssize_t) image->columns; x++)
     {
+      double
+        Da,
+        Sa;
+
       register ssize_t
         i;
 
+      if (GetPixelReadMask(image,p) == 0)
+        {
+          p+=GetPixelChannels(image);
+          q+=GetPixelChannels(reconstruct_image);
+          continue;
+        }
+      Sa=QuantumScale*GetPixelAlpha(image,p);
+      Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
       {
-        MagickRealType
+        double
           distance;
 
-        PixelChannel
-          channel;
-
-        PixelTrait
-          reconstruct_traits,
-          traits;
-
-        traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
-        channel=GetPixelChannelMapChannel(image,(PixelChannel) i);
-        reconstruct_traits=GetPixelChannelMapTraits(reconstruct_image,channel);
+        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 == UndefinedPixelTrait) ||
+            ((reconstruct_traits & UpdatePixelTrait) == 0))
           continue;
-        if ((reconstruct_traits & UpdatePixelTrait) == 0)
-          continue;
-        distance=fabs((double) (alpha*p[i]-beta*GetPixelChannel(
-          reconstruct_image,channel,q)));
+        distance=fabs((double) (Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
+          channel,q)));
         distortion[i]+=distance;
-        distortion[MaxPixelChannels]+=distance;
+        distortion[CompositePixelChannel]+=distance;
         mean_error+=distance*distance;
         if (distance > maximum_error)
           maximum_error=distance;
@@ -694,7 +738,7 @@ static MagickBooleanType GetMeanErrorPerPixel(Image *image,
   }
   reconstruct_view=DestroyCacheView(reconstruct_view);
   image_view=DestroyCacheView(image_view);
-  image->error.mean_error_per_pixel=distortion[MaxPixelChannels]/area;
+  image->error.mean_error_per_pixel=distortion[CompositePixelChannel]/area;
   image->error.normalized_mean_error=QuantumScale*QuantumScale*mean_error/area;
   image->error.normalized_maximum_error=QuantumScale*maximum_error;
   return(status);
@@ -717,10 +761,11 @@ static MagickBooleanType GetMeanSquaredDistortion(const Image *image,
     y;
 
   status=MagickTrue;
-  image_view=AcquireCacheView(image);
-  reconstruct_view=AcquireCacheView(reconstruct_image);
+  image_view=AcquireVirtualCacheView(image,exception);
+  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(dynamic,4) shared(status)
+  #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++)
   {
@@ -748,40 +793,44 @@ static MagickBooleanType GetMeanSquaredDistortion(const Image *image,
     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
     for (x=0; x < (ssize_t) image->columns; x++)
     {
+      double
+        Da,
+        Sa;
+
       register ssize_t
         i;
 
+      if (GetPixelReadMask(image,p) == 0)
+        {
+          p+=GetPixelChannels(image);
+          q+=GetPixelChannels(reconstruct_image);
+          continue;
+        }
+      Sa=QuantumScale*GetPixelAlpha(image,p);
+      Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
       {
-        MagickRealType
+        double
           distance;
 
-        PixelChannel
-          channel;
-
-        PixelTrait
-          reconstruct_traits,
-          traits;
-
-        traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
-        channel=GetPixelChannelMapChannel(image,(PixelChannel) i);
-        reconstruct_traits=GetPixelChannelMapTraits(reconstruct_image,channel);
+        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 == UndefinedPixelTrait) ||
+            ((reconstruct_traits & UpdatePixelTrait) == 0))
           continue;
-        if ((reconstruct_traits & UpdatePixelTrait) == 0)
-          continue;
-        distance=QuantumScale*(p[i]-(MagickRealType) GetPixelChannel(
-          reconstruct_image,channel,q));
-        distance*=distance;
-        channel_distortion[i]+=distance;
-        channel_distortion[MaxPixelChannels]+=distance;
+        distance=QuantumScale*(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
+          channel,q));
+        channel_distortion[i]+=distance*distance;
+        channel_distortion[CompositePixelChannel]+=distance*distance;
       }
       p+=GetPixelChannels(image);
       q+=GetPixelChannels(reconstruct_image);
     }
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_GetMeanSquaredError)
+    #pragma omp critical (MagickCore_GetMeanSquaredError)
 #endif
     for (i=0; i <= MaxPixelChannels; i++)
       distortion[i]+=channel_distortion[i];
@@ -790,7 +839,7 @@ static MagickBooleanType GetMeanSquaredDistortion(const Image *image,
   image_view=DestroyCacheView(image_view);
   for (i=0; i <= MaxPixelChannels; i++)
     distortion[i]/=((double) image->columns*image->rows);
-  distortion[MaxPixelChannels]/=GetImageChannels(image);
+  distortion[CompositePixelChannel]/=GetImageChannels(image);
   return(status);
 }
 
@@ -808,15 +857,15 @@ static MagickBooleanType GetNormalizedCrossCorrelationDistortion(
     *image_statistics,
     *reconstruct_statistics;
 
+  double
+    area;
+
   MagickBooleanType
     status;
 
   MagickOffsetType
     progress;
 
-  MagickRealType
-    area;
-
   register ssize_t
     i;
 
@@ -832,9 +881,9 @@ static MagickBooleanType GetNormalizedCrossCorrelationDistortion(
   progress=0;
   for (i=0; i <= MaxPixelChannels; i++)
     distortion[i]=0.0;
-  area=1.0/((MagickRealType) image->columns*image->rows);
-  image_view=AcquireCacheView(image);
-  reconstruct_view=AcquireCacheView(reconstruct_image);
+  area=1.0/((double) image->columns*image->rows-1);
+  image_view=AcquireVirtualCacheView(image,exception);
+  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
   for (y=0; y < (ssize_t) image->rows; y++)
   {
     register const Quantum
@@ -856,32 +905,37 @@ static MagickBooleanType GetNormalizedCrossCorrelationDistortion(
       }
     for (x=0; x < (ssize_t) image->columns; x++)
     {
+      double
+        Da,
+        Sa;
+
       register ssize_t
         i;
 
+      if (GetPixelReadMask(image,p) == 0)
+        {
+          p+=GetPixelChannels(image);
+          q+=GetPixelChannels(reconstruct_image);
+          continue;
+        }
+      Sa=QuantumScale*GetPixelAlpha(image,p);
+      Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
       {
-        PixelChannel
-          channel;
-
-        PixelTrait
-          reconstruct_traits,
-          traits;
-
-        traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
-        channel=GetPixelChannelMapChannel(image,(PixelChannel) i);
-        reconstruct_traits=GetPixelChannelMapTraits(reconstruct_image,channel);
+        PixelChannel channel=GetPixelChannelChannel(image,i);
+        PixelTrait traits=GetPixelChannelTraits(image,channel);
+        PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
+          channel);
         if ((traits == UndefinedPixelTrait) ||
-            (reconstruct_traits == UndefinedPixelTrait))
-          continue;
-        if ((reconstruct_traits & UpdatePixelTrait) == 0)
+            (reconstruct_traits == UndefinedPixelTrait) ||
+            ((reconstruct_traits & UpdatePixelTrait) == 0))
           continue;
-        distortion[i]+=area*QuantumScale*(p[i]-image_statistics[i].mean)*
-          (GetPixelChannel(reconstruct_image,channel,q)-
+        distortion[i]+=area*QuantumScale*(Sa*p[i]-image_statistics[i].mean)*
+          (Da*GetPixelChannel(reconstruct_image,channel,q)-
           reconstruct_statistics[channel].mean);
       }
       p+=GetPixelChannels(image);
-      q+=GetPixelChannels(image);
+      q+=GetPixelChannels(reconstruct_image);
     }
     if (image->progress_monitor != (MagickProgressMonitor) NULL)
       {
@@ -899,23 +953,20 @@ static MagickBooleanType GetNormalizedCrossCorrelationDistortion(
   /*
     Divide by the standard deviation.
   */
-  distortion[MaxPixelChannels]=0.0;
+  distortion[CompositePixelChannel]=0.0;
   for (i=0; i < MaxPixelChannels; i++)
   {
-    MagickRealType
+    double
       gamma;
 
-    PixelChannel
-      channel;
-
-    channel=GetPixelChannelMapChannel(image,(PixelChannel) i);
+    PixelChannel channel=GetPixelChannelChannel(image,i);
     gamma=image_statistics[i].standard_deviation*
       reconstruct_statistics[channel].standard_deviation;
-    gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
+    gamma=PerceptibleReciprocal(gamma);
     distortion[i]=QuantumRange*gamma*distortion[i];
-    distortion[MaxPixelChannels]+=distortion[i]*distortion[i];
+    distortion[CompositePixelChannel]+=distortion[i]*distortion[i];
   }
-  distortion[MaxPixelChannels]=sqrt(distortion[MaxPixelChannels]/
+  distortion[CompositePixelChannel]=sqrt(distortion[CompositePixelChannel]/
     GetImageChannels(image));
   /*
     Free resources.
@@ -941,10 +992,11 @@ static MagickBooleanType GetPeakAbsoluteDistortion(const Image *image,
     y;
 
   status=MagickTrue;
-  image_view=AcquireCacheView(image);
-  reconstruct_view=AcquireCacheView(reconstruct_image);
+  image_view=AcquireVirtualCacheView(image,exception);
+  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(dynamic,4) shared(status)
+  #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++)
   {
@@ -962,8 +1014,8 @@ static MagickBooleanType GetPeakAbsoluteDistortion(const Image *image,
     if (status == MagickFalse)
       continue;
     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
-    q=GetCacheViewVirtualPixels(reconstruct_view,0,y,
-      reconstruct_image->columns,1,exception);
+    q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
+      1,exception);
     if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
       {
         status=MagickFalse;
@@ -972,41 +1024,46 @@ static MagickBooleanType GetPeakAbsoluteDistortion(const Image *image,
     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
     for (x=0; x < (ssize_t) image->columns; x++)
     {
+      double
+        Da,
+        Sa;
+
       register ssize_t
         i;
 
+      if (GetPixelReadMask(image,p) == 0)
+        {
+          p+=GetPixelChannels(image);
+          q+=GetPixelChannels(reconstruct_image);
+          continue;
+        }
+      Sa=QuantumScale*GetPixelAlpha(image,p);
+      Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
       {
-        MagickRealType
+        double
           distance;
 
-        PixelChannel
-          channel;
-
-        PixelTrait
-          reconstruct_traits,
-          traits;
-
-        traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
-        channel=GetPixelChannelMapChannel(image,(PixelChannel) i);
-        reconstruct_traits=GetPixelChannelMapTraits(reconstruct_image,channel);
+        PixelChannel channel=GetPixelChannelChannel(image,i);
+        PixelTrait traits=GetPixelChannelTraits(image,channel);
+        PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
+          channel);
         if ((traits == UndefinedPixelTrait) ||
-            (reconstruct_traits == UndefinedPixelTrait))
-          continue;
-        if ((reconstruct_traits & UpdatePixelTrait) == 0)
+            (reconstruct_traits == UndefinedPixelTrait) ||
+            ((reconstruct_traits & UpdatePixelTrait) == 0))
           continue;
-        distance=QuantumScale*fabs(p[i]-(MagickRealType) GetPixelChannel(
-          reconstruct_image,channel,q));
+        distance=QuantumScale*fabs(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
+          channel,q));
         if (distance > channel_distortion[i])
           channel_distortion[i]=distance;
-        if (distance > channel_distortion[MaxPixelChannels])
-          channel_distortion[MaxPixelChannels]=distance;
+        if (distance > channel_distortion[CompositePixelChannel])
+          channel_distortion[CompositePixelChannel]=distance;
       }
       p+=GetPixelChannels(image);
-      q+=GetPixelChannels(image);
+      q+=GetPixelChannels(reconstruct_image);
     }
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_GetPeakAbsoluteError)
+    #pragma omp critical (MagickCore_GetPeakAbsoluteError)
 #endif
     for (i=0; i <= MaxPixelChannels; i++)
       if (channel_distortion[i] > distortion[i])
@@ -1141,7 +1198,7 @@ MagickExport MagickBooleanType GetImageDistortion(Image *image,
       break;
     }
   }
-  *distortion=channel_distortion[MaxPixelChannels];
+  *distortion=channel_distortion[CompositePixelChannel];
   channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
   return(status);
 }
@@ -1157,11 +1214,11 @@ MagickExport MagickBooleanType GetImageDistortion(Image *image,
 %                                                                             %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-%  GetImageDistrortion() compares the pixel channels of an image to a
+%  GetImageDistortions() compares the pixel channels of an image to a
 %  reconstructed image and returns the specified distortion metric for each
 %  channel.
 %
-%  The format of the CompareImages method is:
+%  The format of the GetImageDistortions method is:
 %
 %      double *GetImageDistortions(const Image *image,
 %        const Image *reconstruct_image,const MetricType metric,
@@ -1202,8 +1259,8 @@ MagickExport double *GetImageDistortions(Image *image,
   if ((reconstruct_image->columns != image->columns) ||
       (reconstruct_image->rows != image->rows))
     {
-      (void) ThrowMagickException(&image->exception,GetMagickModule(),
-        ImageError,"ImageSizeDiffers","`%s'",image->filename);
+      (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
+        "ImageSizeDiffers","`%s'",image->filename);
       return((double *) NULL);
     }
   /*
@@ -1342,7 +1399,7 @@ MagickExport MagickBooleanType IsImagesEqual(Image *image,
   MagickBooleanType
     status;
 
-  MagickRealType
+  double
     area,
     maximum_error,
     mean_error,
@@ -1362,8 +1419,8 @@ MagickExport MagickBooleanType IsImagesEqual(Image *image,
   maximum_error=0.0;
   mean_error_per_pixel=0.0;
   mean_error=0.0;
-  image_view=AcquireCacheView(image);
-  reconstruct_view=AcquireCacheView(reconstruct_image);
+  image_view=AcquireVirtualCacheView(image,exception);
+  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
   for (y=0; y < (ssize_t) image->rows; y++)
   {
     register const Quantum
@@ -1383,27 +1440,26 @@ MagickExport MagickBooleanType IsImagesEqual(Image *image,
       register ssize_t
         i;
 
+      if (GetPixelReadMask(image,p) == 0)
+        {
+          p+=GetPixelChannels(image);
+          q+=GetPixelChannels(reconstruct_image);
+          continue;
+        }
       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
       {
-        MagickRealType
+        double
           distance;
 
-        PixelChannel
-          channel;
-
-        PixelTrait
-          reconstruct_traits,
-          traits;
-
-        traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
-        channel=GetPixelChannelMapChannel(image,(PixelChannel) i);
-        reconstruct_traits=GetPixelChannelMapTraits(reconstruct_image,channel);
+        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 == UndefinedPixelTrait) ||
+            ((reconstruct_traits & UpdatePixelTrait) == 0))
           continue;
-        if ((reconstruct_traits & UpdatePixelTrait) == 0)
-          continue;
-        distance=fabs(p[i]-(MagickRealType) GetPixelChannel(reconstruct_image,
+        distance=fabs(p[i]-(double) GetPixelChannel(reconstruct_image,
           channel,q));
         mean_error_per_pixel+=distance;
         mean_error+=distance*distance;
@@ -1444,6 +1500,7 @@ MagickExport MagickBooleanType IsImagesEqual(Image *image,
 %  The format of the SimilarityImageImage method is:
 %
 %      Image *SimilarityImage(const Image *image,const Image *reference,
+%        const MetricType metric,const double similarity_threshold,
 %        RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
 %
 %  A description of each parameter follows:
@@ -1452,7 +1509,11 @@ MagickExport MagickBooleanType IsImagesEqual(Image *image,
 %
 %    o reference: find an area of the image that closely resembles this image.
 %
-%    o the best match offset of the reference image within the image.
+%    o metric: the metric.
+%
+%    o similarity_threshold: minimum distortion for (sub)image match.
+%
+%    o offset: the best match offset of the reference image within the image.
 %
 %    o similarity: the computed similarity between the images.
 %
@@ -1460,111 +1521,9 @@ MagickExport MagickBooleanType IsImagesEqual(Image *image,
 %
 */
 
-static double GetNCCDistortion(const Image *image,
-  const Image *reconstruct_image,
-  const ChannelStatistics *reconstruct_statistics,ExceptionInfo *exception)
-{
-#define SimilarityImageTag  "Similarity/Image"
-
-  CacheView
-    *image_view,
-    *reconstruct_view;
-
-  ChannelStatistics
-    *image_statistics;
-
-  double
-    distortion;
-
-  MagickBooleanType
-    status;
-
-  MagickRealType
-    area,
-    gamma;
-
-  ssize_t
-    y;
-
-  /*
-    Normalize to account for variation due to lighting and exposure condition.
-  */
-  image_statistics=GetImageStatistics(image,exception);
-  status=MagickTrue;
-  distortion=0.0;
-  area=1.0/((MagickRealType) image->columns*image->rows);
-  image_view=AcquireCacheView(image);
-  reconstruct_view=AcquireCacheView(reconstruct_image);
-  for (y=0; y < (ssize_t) image->rows; y++)
-  {
-    register const Quantum
-      *restrict p,
-      *restrict q;
-
-    register ssize_t
-      x;
-
-    if (status == MagickFalse)
-      continue;
-    p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
-    q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
-      1,exception);
-    if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
-      {
-        status=MagickFalse;
-        continue;
-      }
-    for (x=0; x < (ssize_t) image->columns; x++)
-    {
-     register ssize_t
-        i;
-
-      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
-      {
-        PixelChannel
-          channel;
-
-        PixelTrait
-          reconstruct_traits,
-          traits;
-
-        traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
-        channel=GetPixelChannelMapChannel(image,(PixelChannel) i);
-        reconstruct_traits=GetPixelChannelMapTraits(reconstruct_image,channel);
-        if ((traits == UndefinedPixelTrait) ||
-            (reconstruct_traits == UndefinedPixelTrait))
-          continue;
-        if ((reconstruct_traits & UpdatePixelTrait) == 0)
-          continue;
-        distortion+=area*QuantumScale*(p[i]-image_statistics[i].mean)*
-          (GetPixelChannel(reconstruct_image,channel,q)-
-          reconstruct_statistics[channel].mean);
-      }
-      p+=GetPixelChannels(image);
-      q+=GetPixelChannels(reconstruct_image);
-    }
-  }
-  reconstruct_view=DestroyCacheView(reconstruct_view);
-  image_view=DestroyCacheView(image_view);
-  /*
-    Divide by the standard deviation.
-  */
-  gamma=image_statistics[MaxPixelChannels].standard_deviation*
-    reconstruct_statistics[MaxPixelChannels].standard_deviation;
-  gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
-  distortion=QuantumRange*gamma*distortion;
-  distortion=sqrt(distortion/GetImageChannels(image));
-  /*
-    Free resources.
-  */
-  image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
-    image_statistics);
-  return(1.0-distortion);
-}
-
 static double GetSimilarityMetric(const Image *image,const Image *reference,
-  const ChannelStatistics *reference_statistics,const ssize_t x_offset,
-  const ssize_t y_offset,ExceptionInfo *exception)
+  const MetricType metric,const ssize_t x_offset,const ssize_t y_offset,
+  ExceptionInfo *exception)
 {
   double
     distortion;
@@ -1572,6 +1531,9 @@ static double GetSimilarityMetric(const Image *image,const Image *reference,
   Image
     *similarity_image;
 
+  MagickBooleanType
+    status;
+
   RectangleInfo
     geometry;
 
@@ -1581,13 +1543,17 @@ static double GetSimilarityMetric(const Image *image,const Image *reference,
   similarity_image=CropImage(image,&geometry,exception);
   if (similarity_image == (Image *) NULL)
     return(0.0);
-  distortion=GetNCCDistortion(reference,similarity_image,reference_statistics,
+  distortion=0.0;
+  status=GetImageDistortion(similarity_image,reference,metric,&distortion,
     exception);
   similarity_image=DestroyImage(similarity_image);
+  if (status == MagickFalse)
+    return(0.0);
   return(distortion);
 }
 
 MagickExport Image *SimilarityImage(Image *image,const Image *reference,
+  const MetricType metric,const double similarity_threshold,
   RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
 {
 #define SimilarityImageTag  "Similarity/Image"
@@ -1595,9 +1561,6 @@ MagickExport Image *SimilarityImage(Image *image,const Image *reference,
   CacheView
     *similarity_view;
 
-  ChannelStatistics
-    *reference_statistics;
-
   Image
     *similarity_image;
 
@@ -1631,15 +1594,18 @@ MagickExport Image *SimilarityImage(Image *image,const Image *reference,
       similarity_image=DestroyImage(similarity_image);
       return((Image *) NULL);
     }
+  (void) SetImageAlphaChannel(similarity_image,DeactivateAlphaChannel,
+    exception);
   /*
     Measure similarity of reference image against image.
   */
   status=MagickTrue;
   progress=0;
-  reference_statistics=GetImageStatistics(reference,exception);
-  similarity_view=AcquireCacheView(similarity_image);
+  similarity_view=AcquireAuthenticCacheView(similarity_image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
+  #pragma omp parallel for schedule(static,4) \
+    shared(progress,status,similarity_metric) \
+    magick_threads(image,image,image->rows,1)
 #endif
   for (y=0; y < (ssize_t) (image->rows-reference->rows+1); y++)
   {
@@ -1654,6 +1620,11 @@ MagickExport Image *SimilarityImage(Image *image,const Image *reference,
 
     if (status == MagickFalse)
       continue;
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+    #pragma omp flush(similarity_metric)
+#endif
+    if (*similarity_metric <= similarity_threshold)
+      continue;
     q=GetCacheViewAuthenticPixels(similarity_view,0,y,similarity_image->columns,
       1,exception);
     if (q == (Quantum *) NULL)
@@ -1666,33 +1637,35 @@ MagickExport Image *SimilarityImage(Image *image,const Image *reference,
       register ssize_t
         i;
 
-      similarity=GetSimilarityMetric(image,reference,reference_statistics,x,y,
-        exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_SimilarityImage)
+      #pragma omp flush(similarity_metric)
+#endif
+      if (*similarity_metric <= similarity_threshold)
+        break;
+      similarity=GetSimilarityMetric(image,reference,metric,x,y,exception);
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+      #pragma omp critical (MagickCore_SimilarityImage)
 #endif
       if (similarity < *similarity_metric)
         {
-          *similarity_metric=similarity;
           offset->x=x;
           offset->y=y;
+          *similarity_metric=similarity;
         }
-      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
+      if (GetPixelReadMask(similarity_image,q) == 0)
+        {
+          q+=GetPixelChannels(similarity_image);
+          continue;
+        }
+      for (i=0; i < (ssize_t) GetPixelChannels(similarity_image); i++)
       {
-        PixelChannel
-          channel;
-
-        PixelTrait
-          similarity_traits,
-          traits;
-
-        traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
-        channel=GetPixelChannelMapChannel(image,(PixelChannel) i);
-        similarity_traits=GetPixelChannelMapTraits(similarity_image,channel);
+        PixelChannel channel=GetPixelChannelChannel(image,i);
+        PixelTrait traits=GetPixelChannelTraits(image,channel);
+        PixelTrait similarity_traits=GetPixelChannelTraits(similarity_image,
+          channel);
         if ((traits == UndefinedPixelTrait) ||
-            (similarity_traits == UndefinedPixelTrait))
-          continue;
-        if ((similarity_traits & UpdatePixelTrait) == 0)
+            (similarity_traits == UndefinedPixelTrait) ||
+            ((similarity_traits & UpdatePixelTrait) == 0))
           continue;
         SetPixelChannel(similarity_image,channel,ClampToQuantum(QuantumRange-
           QuantumRange*similarity),q);
@@ -1707,7 +1680,7 @@ MagickExport Image *SimilarityImage(Image *image,const Image *reference,
           proceed;
 
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_SimilarityImage)
+        #pragma omp critical (MagickCore_SimilarityImage)
 #endif
         proceed=SetImageProgress(image,SimilarityImageTag,progress++,
           image->rows);
@@ -1716,7 +1689,7 @@ MagickExport Image *SimilarityImage(Image *image,const Image *reference,
       }
   }
   similarity_view=DestroyCacheView(similarity_view);
-  reference_statistics=(ChannelStatistics *) RelinquishMagickMemory(
-    reference_statistics);
+  if (status == MagickFalse)
+    similarity_image=DestroyImage(similarity_image);
   return(similarity_image);
 }