/src/imagemagick/MagickCore/compare.c

Source
/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%               CCCC   OOO   M   M  PPPP    AAA   RRRR    EEEEE               %
%              C      O   O  MM MM  P   P  A   A  R   R   E                   %
%              C      O   O  M M M  PPPP   AAAAA  RRRR    EEE                 %
%              C      O   O  M   M  P      A   A  R R     E                   %
%               CCCC   OOO   M   M  P      A   A  R  R    EEEEE               %
%                                                                             %
%                                                                             %
%                      MagickCore Image Comparison Methods                    %
%                                                                             %
%                              Software Design                                %
%                                  Cristy                                     %
%                               December 2003                                 %
%                                                                             %
%                                                                             %
%  Copyright @ 1999 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  %
%  obtain a copy of the License at                                            %
%                                                                             %
%    https://imagemagick.org/license/                                         %
%                                                                             %
%  Unless required by applicable law or agreed to in writing, software        %
%  distributed under the License is distributed on an "AS IS" BASIS,          %
%  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
%  See the License for the specific language governing permissions and        %
%  limitations under the License.                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%
*/

/*
  Include declarations.
*/
#include "MagickCore/studio.h"
#include "MagickCore/artifact.h"
#include "MagickCore/attribute.h"
#include "MagickCore/cache-view.h"
#include "MagickCore/channel.h"
#include "MagickCore/client.h"
#include "MagickCore/color.h"
#include "MagickCore/color-private.h"
#include "MagickCore/colorspace.h"
#include "MagickCore/colorspace-private.h"
#include "MagickCore/compare.h"
#include "MagickCore/compare-private.h"
#include "MagickCore/composite-private.h"
#include "MagickCore/constitute.h"
#include "MagickCore/distort.h"
#include "MagickCore/exception-private.h"
#include "MagickCore/enhance.h"
#include "MagickCore/fourier.h"
#include "MagickCore/geometry.h"
#include "MagickCore/image-private.h"
#include "MagickCore/list.h"
#include "MagickCore/log.h"
#include "MagickCore/memory_.h"
#include "MagickCore/monitor.h"
#include "MagickCore/monitor-private.h"
#include "MagickCore/option.h"
#include "MagickCore/pixel-accessor.h"
#include "MagickCore/property.h"
#include "MagickCore/registry.h"
#include "MagickCore/resource_.h"
#include "MagickCore/string_.h"
#include "MagickCore/statistic.h"
#include "MagickCore/statistic-private.h"
#include "MagickCore/string-private.h"
#include "MagickCore/thread-private.h"
#include "MagickCore/threshold.h"
#include "MagickCore/transform.h"
#include "MagickCore/utility.h"
#include "MagickCore/version.h"

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   C o m p a r e I m a g e s                                                 %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  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:
%
%      Image *CompareImages(const Image *image,const Image *reconstruct_image,
%        const MetricType metric,double *distortion,ExceptionInfo *exception)
%
%  A description of each parameter follows:
%
%    o image: the image.
%
%    o reconstruct_image: the reconstruction image.
%
%    o metric: the metric.
%
%    o distortion: the computed distortion between the images.
%
%    o exception: return any errors or warnings in this structure.
%
*/
MagickExport Image *CompareImages(Image *image,const Image *reconstruct_image,
  const MetricType metric,double *distortion,ExceptionInfo *exception)
{
  CacheView
    *highlight_view,
    *image_view,
    *reconstruct_view;

  const char
    *artifact;

  Image
    *clone_image,
    *difference_image,
    *highlight_image;

  MagickBooleanType
    status = MagickTrue;

  PixelInfo
    highlight,
    lowlight,
    masklight;

  RectangleInfo
    geometry;

  size_t
    columns,
    rows;

  ssize_t
    y;

  assert(image != (Image *) NULL);
  assert(image->signature == MagickCoreSignature);
  assert(reconstruct_image != (const Image *) NULL);
  assert(reconstruct_image->signature == MagickCoreSignature);
  assert(distortion != (double *) NULL);
  if (IsEventLogging() != MagickFalse)
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
  *distortion=0.0;
  status=GetImageDistortion(image,reconstruct_image,metric,distortion,
    exception);
  if (status == MagickFalse)
    return((Image *) NULL);
  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
  SetGeometry(image,&geometry);
  geometry.width=columns;
  geometry.height=rows;
  clone_image=CloneImage(image,0,0,MagickTrue,exception);
  if (clone_image == (Image *) NULL)
    return((Image *) NULL);
  (void) SetImageMask(clone_image,ReadPixelMask,(Image *) NULL,exception);
  difference_image=ExtentImage(clone_image,&geometry,exception);
  clone_image=DestroyImage(clone_image);
  if (difference_image == (Image *) NULL)
    return((Image *) NULL);
  (void) ResetImagePage(difference_image,"0x0+0+0");
  (void) SetImageAlphaChannel(difference_image,OpaqueAlphaChannel,exception);
  highlight_image=CloneImage(image,columns,rows,MagickTrue,exception);
  if (highlight_image == (Image *) NULL)
    {
      difference_image=DestroyImage(difference_image);
      return((Image *) NULL);
    }
  status=SetImageStorageClass(highlight_image,DirectClass,exception);
  if (status == MagickFalse)
    {
      difference_image=DestroyImage(difference_image);
      highlight_image=DestroyImage(highlight_image);
      return((Image *) NULL);
    }
  (void) SetImageMask(highlight_image,ReadPixelMask,(Image *) NULL,exception);
  (void) SetImageAlphaChannel(highlight_image,OpaqueAlphaChannel,exception);
  (void) QueryColorCompliance("#f1001ecc",AllCompliance,&highlight,exception);
  artifact=GetImageArtifact(image,"compare:highlight-color");
  if (artifact != (const char *) NULL)
    (void) QueryColorCompliance(artifact,AllCompliance,&highlight,exception);
  (void) QueryColorCompliance("#ffffffcc",AllCompliance,&lowlight,exception);
  artifact=GetImageArtifact(image,"compare:lowlight-color");
  if (artifact != (const char *) NULL)
    (void) QueryColorCompliance(artifact,AllCompliance,&lowlight,exception);
  (void) QueryColorCompliance("#888888cc",AllCompliance,&masklight,exception);
  artifact=GetImageArtifact(image,"compare:masklight-color");
  if (artifact != (const char *) NULL)
    (void) QueryColorCompliance(artifact,AllCompliance,&masklight,exception);
  /*
    Generate difference 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(static) shared(status) \
    magick_number_threads(image,highlight_image,rows,1)
#endif
  for (y=0; y < (ssize_t) rows; y++)
  {
    const Quantum
      *magick_restrict p,
      *magick_restrict q;

    MagickBooleanType
      sync;

    Quantum
      *magick_restrict r;

    ssize_t
      x;

    if (status == MagickFalse)
      continue;
    p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
    q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
    r=QueueCacheViewAuthenticPixels(highlight_view,0,y,columns,1,exception);
    if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL) ||
        (r == (Quantum *) NULL))
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (ssize_t) columns; x++)
    {
      if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
          (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
        {
          SetPixelViaPixelInfo(highlight_image,&masklight,r);
          p+=(ptrdiff_t) GetPixelChannels(image);
          q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
          r+=(ptrdiff_t) GetPixelChannels(highlight_image);
          continue;
        }
      if (IsFuzzyEquivalencePixel(image,p,reconstruct_image,q) == MagickFalse)
        SetPixelViaPixelInfo(highlight_image,&highlight,r);
      else
        SetPixelViaPixelInfo(highlight_image,&lowlight,r);
      p+=(ptrdiff_t) GetPixelChannels(image);
      q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
      r+=(ptrdiff_t) GetPixelChannels(highlight_image);
    }
    sync=SyncCacheViewAuthenticPixels(highlight_view,exception);
    if (sync == MagickFalse)
      status=MagickFalse;
  }
  highlight_view=DestroyCacheView(highlight_view);
  reconstruct_view=DestroyCacheView(reconstruct_view);
  image_view=DestroyCacheView(image_view);
  if ((status != MagickFalse) && (difference_image != (Image *) NULL))
    status=CompositeImage(difference_image,highlight_image,image->compose,
      MagickTrue,0,0,exception);
  highlight_image=DestroyImage(highlight_image);
  if (status == MagickFalse)
    difference_image=DestroyImage(difference_image);
  return(difference_image);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   G e t I m a g e D i s t o r t i o n                                       %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  GetImageDistortion() compares one or more pixel channels of an image to a
%  reconstructed image and returns the specified distortion metric.
%
%  The format of the GetImageDistortion method is:
%
%      MagickBooleanType GetImageDistortion(const Image *image,
%        const Image *reconstruct_image,const MetricType metric,
%        double *distortion,ExceptionInfo *exception)
%
%  A description of each parameter follows:
%
%    o image: the image.
%
%    o reconstruct_image: the reconstruction image.
%
%    o metric: the metric.
%
%    o distortion: the computed distortion between the images.
%
%    o exception: return any errors or warnings in this structure.
%
*/

static MagickBooleanType GetAESimilarity(const Image *image,
  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
{
  CacheView
    *image_view,
    *reconstruct_view;

  double
    area,
    fuzz;

  MagickBooleanType
    status = MagickTrue;

  size_t
    columns,
    rows;

  ssize_t
    channels = 0,
    k,
    y;

  /*
    Compute the absolute error similarity.
  */
  fuzz=GetFuzzyColorDistance(image,reconstruct_image);
  (void) memset(similarity,0,(MaxPixelChannels+1)*sizeof(*similarity));
  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
  image_view=AcquireVirtualCacheView(image,exception);
  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(similarity,status) \
    magick_number_threads(image,image,rows,1)
#endif
  for (y=0; y < (ssize_t) rows; y++)
  {
    const Quantum
      *magick_restrict p,
      *magick_restrict q;

    double
      channel_similarity[MaxPixelChannels+1] = { 0.0 };

    ssize_t
      x;

    if (status == MagickFalse)
      continue;
    p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
    q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
    if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (ssize_t) columns; x++)
    {
      double
        Da,
        Sa;

      ssize_t
        i;

      if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
          (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
        {
          p+=(ptrdiff_t) GetPixelChannels(image);
          q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
          continue;
        }
      Sa=QuantumScale*(double) GetPixelAlpha(image,p);
      Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
      {
        double
          error;

        PixelChannel channel = GetPixelChannelChannel(image,i);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
          channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((reconstruct_traits & UpdatePixelTrait) == 0))
          continue;
        if (channel == AlphaPixelChannel)
          error=(double) p[i]-(double) GetPixelChannel(reconstruct_image,
            channel,q);
        else
          error=Sa*(double) p[i]-Da*(double) GetPixelChannel(reconstruct_image,channel,
            q);
        if (MagickSafeSignificantError(error*error,fuzz) != MagickFalse)
          {
            double ae = fabs(QuantumScale*error);
            channel_similarity[i]+=ae;
            channel_similarity[CompositePixelChannel]+=ae;
          }
      }
      p+=(ptrdiff_t) GetPixelChannels(image);
      q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
    }
#if defined(MAGICKCORE_OPENMP_SUPPORT)
    #pragma omp critical (MagickCore_GetAESimilarity)
#endif
    {
      ssize_t
        j;

      for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
      {
        PixelChannel channel = GetPixelChannelChannel(image,j);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
          channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((reconstruct_traits & UpdatePixelTrait) == 0))
          continue;
        similarity[j]+=channel_similarity[j];
      }
      similarity[CompositePixelChannel]+=
        channel_similarity[CompositePixelChannel];
    }
  }
  reconstruct_view=DestroyCacheView(reconstruct_view);
  image_view=DestroyCacheView(image_view);
  area=MagickSafeReciprocal((double) columns*(double) rows);
  for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
  {
    PixelChannel channel = GetPixelChannelChannel(image,k);
    PixelTrait traits = GetPixelChannelTraits(image,channel);
    PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
      channel);
    if (((traits & UpdatePixelTrait) == 0) ||
        ((reconstruct_traits & UpdatePixelTrait) == 0))
      continue;
    similarity[k]*=area;
    channels++;
  }
  similarity[CompositePixelChannel]*=area;
  if (channels != 0)
    similarity[CompositePixelChannel]/=(double) channels;
  return(status);
}

static MagickBooleanType GetDPCSimilarity(const Image *image,
  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
{
#define SimilarityImageTag  "Similarity/Image"

  CacheView
    *image_view,
    *reconstruct_view;

  ChannelStatistics
    *image_statistics,
    *reconstruct_statistics;

  double
    norm[MaxPixelChannels+1] = { 0.0 },
    reconstruct_norm[MaxPixelChannels+1] = { 0.0 };

  MagickBooleanType
    status = MagickTrue;

  MagickOffsetType
    progress = 0;

  size_t
    columns,
    rows;

  ssize_t
    k,
    y;

  /*
    Compute the dot product correlation similarity.
  */
  image_statistics=GetImageStatistics(image,exception);
  reconstruct_statistics=GetImageStatistics(reconstruct_image,exception);
  if ((image_statistics == (ChannelStatistics *) NULL) ||
      (reconstruct_statistics == (ChannelStatistics *) NULL))
    {
      if (image_statistics != (ChannelStatistics *) NULL)
        image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
          image_statistics);
      if (reconstruct_statistics != (ChannelStatistics *) NULL)
        reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
          reconstruct_statistics);
      return(MagickFalse);
    }
  (void) memset(similarity,0,(MaxPixelChannels+1)*sizeof(*similarity));
  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
  image_view=AcquireVirtualCacheView(image,exception);
  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(norm,reconstruct_norm,similarity,status) \
    magick_number_threads(image,image,rows,1)
#endif
  for (y=0; y < (ssize_t) rows; y++)
  {
    const Quantum
      *magick_restrict p,
      *magick_restrict q;

    double
      channel_norm[MaxPixelChannels+1] = { 0.0 },
      channel_reconstruct_norm[MaxPixelChannels+1] = { 0.0 },
      channel_similarity[MaxPixelChannels+1] = { 0.0 };

    ssize_t
      x;

    p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
    q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
    if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (ssize_t) columns; x++)
    {
      double
        Da,
        Sa;

      ssize_t
        i;

      if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
          (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
        {
          p+=(ptrdiff_t) GetPixelChannels(image);
          q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
          continue;
        }
      Sa=QuantumScale*(double) GetPixelAlpha(image,p);
      Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
      {
        double
          alpha,
          beta;

        PixelChannel channel = GetPixelChannelChannel(image,i);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
          channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((reconstruct_traits & UpdatePixelTrait) == 0))
          continue;
        if (channel == AlphaPixelChannel)
          {
            alpha=QuantumScale*((double) p[i]-image_statistics[channel].mean);
            beta=QuantumScale*((double) GetPixelChannel(reconstruct_image,
              channel,q)-reconstruct_statistics[channel].mean);
          }
        else
          {
            alpha=QuantumScale*(Sa*(double) p[i]-
              image_statistics[channel].mean);
            beta=QuantumScale*(Da*(double) GetPixelChannel(reconstruct_image,channel,
              q)-reconstruct_statistics[channel].mean);
          }
        channel_similarity[i]+=alpha*beta;
        channel_norm[i]+=alpha*alpha;
        channel_reconstruct_norm[i]+=beta*beta;
      }
      p+=(ptrdiff_t) GetPixelChannels(image);
      q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
    }
#if defined(MAGICKCORE_OPENMP_SUPPORT)
    #pragma omp critical (MagickCore_GetDPCSimilarity)
#endif
    {
      ssize_t
        j;

      for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
      {
        PixelChannel channel = GetPixelChannelChannel(image,j);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
          channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((reconstruct_traits & UpdatePixelTrait) == 0))
          continue;
        similarity[j]+=channel_similarity[j];
        similarity[CompositePixelChannel]+=channel_similarity[j];
        norm[j]+=channel_norm[j];
        norm[CompositePixelChannel]+=channel_norm[j];
        reconstruct_norm[j]+=channel_reconstruct_norm[j];
        reconstruct_norm[CompositePixelChannel]+=channel_reconstruct_norm[j];
      }
    }
    if (image->progress_monitor != (MagickProgressMonitor) NULL)
      {
        MagickBooleanType
          proceed;

#if defined(MAGICKCORE_OPENMP_SUPPORT)
        #pragma omp atomic
#endif
        progress++;
        proceed=SetImageProgress(image,SimilarityImageTag,progress,rows);
        if (proceed == MagickFalse)
          {
            status=MagickFalse;
            continue;
          }
      }
  }
  reconstruct_view=DestroyCacheView(reconstruct_view);
  image_view=DestroyCacheView(image_view);
  /*
    Compute dot product correlation: divide by mean.
  */
  for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
  {
    PixelChannel channel = GetPixelChannelChannel(image,k);
    PixelTrait traits = GetPixelChannelTraits(image,channel);
    PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
      channel);
    if (((traits & UpdatePixelTrait) == 0) ||
        ((reconstruct_traits & UpdatePixelTrait) == 0))
      continue;
    similarity[k]*=MagickSafeReciprocal(sqrt(norm[k]*reconstruct_norm[k]));
  }
  similarity[CompositePixelChannel]*=MagickSafeReciprocal(sqrt(
    norm[CompositePixelChannel]*reconstruct_norm[CompositePixelChannel]));
  /*
    Free resources.
  */
  reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
    reconstruct_statistics);
  image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
    image_statistics);
  return(status);
}

static MagickBooleanType GetFUZZSimilarity(const Image *image,
  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
{
  CacheView
    *image_view,
    *reconstruct_view;

  double
    area = 0.0,
    fuzz = 0.0;

  MagickBooleanType
    status = MagickTrue;

  size_t
    columns,
    rows;

  ssize_t
    k,
    y;

  /*
    Compute the MSE similarity within tolerance (fuzz).
  */
  fuzz=GetFuzzyColorDistance(image,reconstruct_image);
  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
  image_view=AcquireVirtualCacheView(image,exception);
  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(area,similarity,status) \
    magick_number_threads(image,image,rows,1)
#endif
  for (y=0; y < (ssize_t) rows; y++)
  {
    const Quantum
      *magick_restrict p,
      *magick_restrict q;

    double
      channel_area = 0.0,
      channel_similarity[MaxPixelChannels+1] = { 0.0 };

    ssize_t
      x;

    if (status == MagickFalse)
      continue;
    p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
    q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
    if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (ssize_t) columns; x++)
    {
      double
        Da,
        Sa;

      ssize_t
        i;

      if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
          (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
        {
          p+=(ptrdiff_t) GetPixelChannels(image);
          q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
          continue;
        }
      Sa=QuantumScale*(double) GetPixelAlpha(image,p);
      Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
      {
        double
          error;

        PixelChannel channel = GetPixelChannelChannel(image,i);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
          channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((reconstruct_traits & UpdatePixelTrait) == 0))
          continue;
        if (channel == AlphaPixelChannel)
          error=(double) p[i]-(double) GetPixelChannel(reconstruct_image,
            channel,q);
        else
          error=Sa*(double) p[i]-Da*(double) GetPixelChannel(reconstruct_image,
            channel,q);
        if (MagickSafeSignificantError(error*error,fuzz) != MagickFalse)
          {
            channel_similarity[i]+=QuantumScale*error*QuantumScale*error;
            channel_similarity[CompositePixelChannel]+=QuantumScale*error*
              QuantumScale*error;
            channel_area++;
          }
      }
      p+=(ptrdiff_t) GetPixelChannels(image);
      q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
    }
#if defined(MAGICKCORE_OPENMP_SUPPORT)
    #pragma omp critical (MagickCore_GetFUZZSimilarity)
#endif
    {
      ssize_t
        j;

      area+=channel_area;
      for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
      {
        PixelChannel channel = GetPixelChannelChannel(image,j);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
          channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((reconstruct_traits & UpdatePixelTrait) == 0))
          continue;
        similarity[j]+=channel_similarity[j];
      }
      similarity[CompositePixelChannel]+=
        channel_similarity[CompositePixelChannel];
    }
  }
  reconstruct_view=DestroyCacheView(reconstruct_view);
  image_view=DestroyCacheView(image_view);
  area=MagickSafeReciprocal(area);
  for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
  {
    PixelChannel channel = GetPixelChannelChannel(image,k);
    PixelTrait traits = GetPixelChannelTraits(image,channel);
    PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
      channel);
    if (((traits & UpdatePixelTrait) == 0) ||
        ((reconstruct_traits & UpdatePixelTrait) == 0))
      continue;
    similarity[k]*=area;
  }
  similarity[CompositePixelChannel]*=area;
  return(status);
}

static MagickBooleanType GetMAESimilarity(const Image *image,
  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
{
  CacheView
    *image_view,
    *reconstruct_view;

  double
    area = 0.0;

  MagickBooleanType
    status = MagickTrue;

  size_t
    columns,
    rows;

  ssize_t
    channels = 0,
    k,
    y;

  /*
    Compute the mean absolute error similarity.
  */
  (void) memset(similarity,0,(MaxPixelChannels+1)*sizeof(*similarity));
  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
  image_view=AcquireVirtualCacheView(image,exception);
  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(area,similarity,status) \
    magick_number_threads(image,image,rows,1)
#endif
  for (y=0; y < (ssize_t) rows; y++)
  {
    const Quantum
      *magick_restrict p,
      *magick_restrict q;

    double
      channel_area = 0.0,
      channel_similarity[MaxPixelChannels+1] = { 0.0 };

    ssize_t
      x;

    if (status == MagickFalse)
      continue;
    p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
    q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
    if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (ssize_t) columns; x++)
    {
      double
        Da,
        Sa;

      ssize_t
        i;

      if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
          (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
        {
          p+=(ptrdiff_t) GetPixelChannels(image);
          q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
          continue;
        }
      Sa=QuantumScale*(double) GetPixelAlpha(image,p);
      Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
      {
        double
          error;

        PixelChannel channel = GetPixelChannelChannel(image,i);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
          channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((reconstruct_traits & UpdatePixelTrait) == 0))
          continue;
        if (channel == AlphaPixelChannel)
          error=QuantumScale*fabs((double) p[i]-(double) GetPixelChannel(
            reconstruct_image,channel,q));
        else
          error=QuantumScale*fabs(Sa*(double) p[i]-Da*(double)
            GetPixelChannel(reconstruct_image,channel,q));
        channel_similarity[i]+=error;
        channel_similarity[CompositePixelChannel]+=error;
      }
      channel_area++;
      p+=(ptrdiff_t) GetPixelChannels(image);
      q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
    }
#if defined(MAGICKCORE_OPENMP_SUPPORT)
    #pragma omp critical (MagickCore_GetMAESimilarity)
#endif
    {
      ssize_t
        j;

      area+=channel_area;
      for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
      {
        PixelChannel channel = GetPixelChannelChannel(image,j);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
          channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((reconstruct_traits & UpdatePixelTrait) == 0))
          continue;
        similarity[j]+=channel_similarity[j];
      }
      similarity[CompositePixelChannel]+=
        channel_similarity[CompositePixelChannel];
    }
  }
  reconstruct_view=DestroyCacheView(reconstruct_view);
  image_view=DestroyCacheView(image_view);
  area=MagickSafeReciprocal(area);
  for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
  {
    PixelChannel channel = GetPixelChannelChannel(image,k);
    PixelTrait traits = GetPixelChannelTraits(image,channel);
    PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
      channel);
    if (((traits & UpdatePixelTrait) == 0) ||
        ((reconstruct_traits & UpdatePixelTrait) == 0))
      continue;
    similarity[k]*=area;
    channels++;
  }
  similarity[CompositePixelChannel]*=area;
  if (channels != 0)
    similarity[CompositePixelChannel]/=(double) channels;
  return(status);
}

static MagickBooleanType GetMEPPSimilarity(Image *image,
  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
{
  CacheView
    *image_view,
    *reconstruct_view;

  double
    area = 0.0,
    maximum_error = -MagickMaximumValue,
    mean_error = 0.0;

  MagickBooleanType
    status = MagickTrue;

  size_t
    columns,
    rows;

  ssize_t
    channels = 0,
    k,
    y;

  /*
    Compute the mean error per pixel similarity.
  */
  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
  image_view=AcquireVirtualCacheView(image,exception);
  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(area,similarity,maximum_error,mean_error,status) \
    magick_number_threads(image,image,rows,1)
#endif
  for (y=0; y < (ssize_t) rows; y++)
  {
    const Quantum
      *magick_restrict p,
      *magick_restrict q;

    double
      channel_area = 0.0,
      channel_similarity[MaxPixelChannels+1] = { 0.0 },
      channel_maximum_error = maximum_error,
      channel_mean_error = 0.0;

    ssize_t
      x;

    if (status == MagickFalse)
      continue;
    p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
    q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
    if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (ssize_t) columns; x++)
    {
      double
        Da,
        Sa;

      ssize_t
        i;

      if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
          (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
        {
          p+=(ptrdiff_t) GetPixelChannels(image);
          q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
          continue;
        }
      Sa=QuantumScale*(double) GetPixelAlpha(image,p);
      Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
      {
        double
          error;

        PixelChannel channel = GetPixelChannelChannel(image,i);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
          channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((reconstruct_traits & UpdatePixelTrait) == 0))
          continue;
        if (channel == AlphaPixelChannel)
          error=QuantumScale*fabs((double) p[i]-(double) GetPixelChannel(
            reconstruct_image,channel,q));
        else
          error=QuantumScale*fabs(Sa*(double) p[i]-Da*(double)
            GetPixelChannel(reconstruct_image,channel,q));
        channel_similarity[i]+=error;
        channel_similarity[CompositePixelChannel]+=error;
        channel_mean_error+=error*error;
        if (error > channel_maximum_error)
          channel_maximum_error=error;
      }
      channel_area++;
      p+=(ptrdiff_t) GetPixelChannels(image);
      q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
    }
#if defined(MAGICKCORE_OPENMP_SUPPORT)
    #pragma omp critical (MagickCore_GetMEPPSimilarity)
#endif
    {
      ssize_t
        j;

      area+=channel_area;
      for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
      {
        PixelChannel channel = GetPixelChannelChannel(image,j);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
          channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((reconstruct_traits & UpdatePixelTrait) == 0))
          continue;
        similarity[j]+=channel_similarity[j];
      }
      similarity[CompositePixelChannel]+=
        channel_similarity[CompositePixelChannel];
      mean_error+=channel_mean_error;
      if (channel_maximum_error > maximum_error)
        maximum_error=channel_maximum_error;
    }
  }
  reconstruct_view=DestroyCacheView(reconstruct_view);
  image_view=DestroyCacheView(image_view);
  area=MagickSafeReciprocal(area);
  for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
  {
    PixelChannel channel = GetPixelChannelChannel(image,k);
    PixelTrait traits = GetPixelChannelTraits(image,channel);
    PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
      channel);
    if (((traits & UpdatePixelTrait) == 0) ||
        ((reconstruct_traits & UpdatePixelTrait) == 0))
      continue;
    similarity[k]*=area;
    channels++;
  }
  similarity[CompositePixelChannel]*=area;
  if (channels != 0)
    similarity[CompositePixelChannel]/=(double) channels;
  image->error.mean_error_per_pixel=(double) QuantumRange*
    similarity[CompositePixelChannel];
  image->error.normalized_mean_error=mean_error*area;
  image->error.normalized_maximum_error=maximum_error;
  return(status);
}

static MagickBooleanType GetMSESimilarity(const Image *image,
  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
{
  CacheView
    *image_view,
    *reconstruct_view;

  double
    area = 0.0;

  MagickBooleanType
    status = MagickTrue;

  size_t
    columns,
    rows;

  ssize_t
    channels = 0,
    k,
    y;

  /*
    Compute the mean sequared error similarity.
  */
  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
  image_view=AcquireVirtualCacheView(image,exception);
  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(area,similarity,status) \
    magick_number_threads(image,image,rows,1)
#endif
  for (y=0; y < (ssize_t) rows; y++)
  {
    const Quantum
      *magick_restrict p,
      *magick_restrict q;

    double
      channel_area = 0.0,
      channel_similarity[MaxPixelChannels+1] = { 0.0 };

    ssize_t
      x;

    if (status == MagickFalse)
      continue;
    p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
    q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
    if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (ssize_t) columns; x++)
    {
      double
        Da,
        Sa;

      ssize_t
        i;

      if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
          (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
        {
          p+=(ptrdiff_t) GetPixelChannels(image);
          q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
          continue;
        }
      Sa=QuantumScale*(double) GetPixelAlpha(image,p);
      Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
      {
        double
          error;

        PixelChannel channel = GetPixelChannelChannel(image,i);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
          channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((reconstruct_traits & UpdatePixelTrait) == 0))
          continue;
        if (channel == AlphaPixelChannel)
          error=QuantumScale*((double) p[i]-(double) GetPixelChannel(
            reconstruct_image,channel,q));
        else
          error=QuantumScale*(Sa*(double) p[i]-Da*(double) GetPixelChannel(
            reconstruct_image,channel,q));
        channel_similarity[i]+=error*error;
        channel_similarity[CompositePixelChannel]+=error*error;
      }
      channel_area++;
      p+=(ptrdiff_t) GetPixelChannels(image);
      q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
    }
#if defined(MAGICKCORE_OPENMP_SUPPORT)
    #pragma omp critical (MagickCore_GetMSESimilarity)
#endif
    {
      ssize_t
        j;

      area+=channel_area;
      for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
      {
        PixelChannel channel = GetPixelChannelChannel(image,j);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
          channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((reconstruct_traits & UpdatePixelTrait) == 0))
          continue;
        similarity[j]+=channel_similarity[j];
      }
      similarity[CompositePixelChannel]+=
        channel_similarity[CompositePixelChannel];
    }
  }
  reconstruct_view=DestroyCacheView(reconstruct_view);
  image_view=DestroyCacheView(image_view);
  area=MagickSafeReciprocal(area);
  for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
  {
    PixelChannel channel = GetPixelChannelChannel(image,k);
    PixelTrait traits = GetPixelChannelTraits(image,channel);
    PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
      channel);
    if (((traits & UpdatePixelTrait) == 0) ||
        ((reconstruct_traits & UpdatePixelTrait) == 0))
      continue;
    similarity[k]*=area;
    channels++;
  }
  similarity[CompositePixelChannel]*=area;
  if (channels != 0)
   similarity[CompositePixelChannel]/=(double) channels;
  return(status);
}

static MagickBooleanType GetNCCSimilarity(const Image *image,
  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
{
  CacheView
    *image_view,
    *reconstruct_view;

  ChannelStatistics
    *image_statistics,
    *reconstruct_statistics;

  double
    reconstruct_variance[MaxPixelChannels+1] = { 0.0 },
    variance[MaxPixelChannels+1] = { 0.0 };

  MagickBooleanType
    status = MagickTrue;

  MagickOffsetType
    progress = 0;

  size_t
    columns,
    rows;

  ssize_t
    k,
    y;

  /*
    Compute the normalized criss-correlation similarity.
  */
  image_statistics=GetImageStatistics(image,exception);
  reconstruct_statistics=GetImageStatistics(reconstruct_image,exception);
  if ((image_statistics == (ChannelStatistics *) NULL) ||
      (reconstruct_statistics == (ChannelStatistics *) NULL))
    {
      if (image_statistics != (ChannelStatistics *) NULL)
        image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
          image_statistics);
      if (reconstruct_statistics != (ChannelStatistics *) NULL)
        reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
          reconstruct_statistics);
      return(MagickFalse);
    }
  (void) memset(similarity,0,(MaxPixelChannels+1)*sizeof(*similarity));
  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
  image_view=AcquireVirtualCacheView(image,exception);
  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(variance,reconstruct_variance,similarity,status) \
    magick_number_threads(image,image,rows,1)
#endif
  for (y=0; y < (ssize_t) rows; y++)
  {
    const Quantum
      *magick_restrict p,
      *magick_restrict q;

    double
      channel_reconstruct_variance[MaxPixelChannels+1] = { 0.0 },
      channel_similarity[MaxPixelChannels+1] = { 0.0 },
      channel_variance[MaxPixelChannels+1] = { 0.0 };

    ssize_t
      x;

    p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
    q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
    if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (ssize_t) columns; x++)
    {
      double
        Da,
        Sa;

      ssize_t
        i;

      if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
          (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
        {
          p+=(ptrdiff_t) GetPixelChannels(image);
          q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
          continue;
        }
      Sa=QuantumScale*(double) GetPixelAlpha(image,p);
      Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
      {
        double
          alpha,
          beta;

        PixelChannel channel = GetPixelChannelChannel(image,i);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
          channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((reconstruct_traits & UpdatePixelTrait) == 0))
          continue;
        if (channel == AlphaPixelChannel)
          {
            alpha=QuantumScale*((double) p[i]-image_statistics[channel].mean);
            beta=QuantumScale*((double) GetPixelChannel(reconstruct_image,
              channel,q)-reconstruct_statistics[channel].mean);
          }
        else
          {
            alpha=QuantumScale*(Sa*(double) p[i]-
              image_statistics[channel].mean);
            beta=QuantumScale*(Da*(double) GetPixelChannel(reconstruct_image,channel,
              q)-reconstruct_statistics[channel].mean);
          }
        channel_similarity[i]+=alpha*beta;
        channel_variance[i]+=alpha*alpha;
        channel_reconstruct_variance[i]+=beta*beta;
      }
      p+=(ptrdiff_t) GetPixelChannels(image);
      q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
    }
#if defined(MAGICKCORE_OPENMP_SUPPORT)
    #pragma omp critical (MagickCore_GetNCCSimilarity)
#endif
    {
      ssize_t
        j;

      for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
      {
        PixelChannel channel = GetPixelChannelChannel(image,j);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
          channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((reconstruct_traits & UpdatePixelTrait) == 0))
          continue;
        similarity[j]+=channel_similarity[j];
        similarity[CompositePixelChannel]+=channel_similarity[j];
        variance[j]+=channel_variance[j];
        variance[CompositePixelChannel]+=channel_variance[j];
        reconstruct_variance[j]+=channel_reconstruct_variance[j];
        reconstruct_variance[CompositePixelChannel]+=
          channel_reconstruct_variance[j];
      }
    }
    if (image->progress_monitor != (MagickProgressMonitor) NULL)
      {
        MagickBooleanType
          proceed;

#if defined(MAGICKCORE_OPENMP_SUPPORT)
        #pragma omp atomic
#endif
        progress++;
        proceed=SetImageProgress(image,SimilarityImageTag,progress,rows);
        if (proceed == MagickFalse)
          {
            status=MagickFalse;
            continue;
          }
      }
  }
  reconstruct_view=DestroyCacheView(reconstruct_view);
  image_view=DestroyCacheView(image_view);
  /*
    Compute normalized cross correlation: divide by standard deviation.
  */
  for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
  {
    PixelChannel channel = GetPixelChannelChannel(image,k);
    PixelTrait traits = GetPixelChannelTraits(image,channel);
    PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
      channel);
    if (((traits & UpdatePixelTrait) == 0) ||
        ((reconstruct_traits & UpdatePixelTrait) == 0))
      continue;
    similarity[k]*=MagickSafeReciprocal(sqrt(variance[k])*
      sqrt(reconstruct_variance[k]));
  }
  similarity[CompositePixelChannel]*=MagickSafeReciprocal(sqrt(
    variance[CompositePixelChannel])*sqrt(
    reconstruct_variance[CompositePixelChannel]));
  /*
    Free resources.
  */
  reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
    reconstruct_statistics);
  image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
    image_statistics);
  return(status);
}

static MagickBooleanType GetPASimilarity(const Image *image,
  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
{
  CacheView
    *image_view,
    *reconstruct_view;

  MagickBooleanType
    status = MagickTrue;

  size_t
    columns,
    rows;

  ssize_t
    y;

  /*
    Compute the peak absolute similarity.
  */
  (void) memset(similarity,0,(MaxPixelChannels+1)*sizeof(*similarity));
  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
  image_view=AcquireVirtualCacheView(image,exception);
  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(similarity,status) \
    magick_number_threads(image,image,rows,1)
#endif
  for (y=0; y < (ssize_t) rows; y++)
  {
    const Quantum
      *magick_restrict p,
      *magick_restrict q;

    double
      channel_similarity[MaxPixelChannels+1] = { 0.0 };

    ssize_t
      x;

    if (status == MagickFalse)
      continue;
    p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
    q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
    if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (ssize_t) columns; x++)
    {
      double
        Da,
        Sa;

      ssize_t
        i;

      if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
          (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
        {
          p+=(ptrdiff_t) GetPixelChannels(image);
          q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
          continue;
        }
      Sa=QuantumScale*(double) GetPixelAlpha(image,p);
      Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
      {
        double
          distance;

        PixelChannel channel = GetPixelChannelChannel(image,i);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
          channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((reconstruct_traits & UpdatePixelTrait) == 0))
          continue;
        if (channel == AlphaPixelChannel)
          distance=QuantumScale*fabs((double) p[i]-(double)
            GetPixelChannel(reconstruct_image,channel,q));
        else
          distance=QuantumScale*fabs(Sa*(double) p[i]-Da*(double) GetPixelChannel(
            reconstruct_image,channel,q));
        if (distance > channel_similarity[i])
          channel_similarity[i]=distance;
        if (distance > channel_similarity[CompositePixelChannel])
          channel_similarity[CompositePixelChannel]=distance;
      }
      p+=(ptrdiff_t) GetPixelChannels(image);
      q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
    }
#if defined(MAGICKCORE_OPENMP_SUPPORT)
    #pragma omp critical (MagickCore_GetPASimilarity)
#endif
    {
      ssize_t
        j;

      for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
      {
        PixelChannel channel = GetPixelChannelChannel(image,j);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
          channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((reconstruct_traits & UpdatePixelTrait) == 0))
          continue;
        if (channel_similarity[j] > similarity[j])
          similarity[j]=channel_similarity[j];
      }
      if (channel_similarity[CompositePixelChannel] > similarity[CompositePixelChannel])
        similarity[CompositePixelChannel]=
          channel_similarity[CompositePixelChannel];
    }
  }
  reconstruct_view=DestroyCacheView(reconstruct_view);
  image_view=DestroyCacheView(image_view);
  return(status);
}

static MagickBooleanType GetPDCSimilarity(const Image *image,
  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
{
  CacheView
    *image_view,
    *reconstruct_view;

  double
    fuzz;

  MagickBooleanType
    status = MagickTrue;

  size_t
    columns,
    rows;

  ssize_t
    y;

  /*
    Compute the pixel difference count similarity.
  */
  fuzz=GetFuzzyColorDistance(image,reconstruct_image);
  (void) memset(similarity,0,(MaxPixelChannels+1)*sizeof(*similarity));
  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
  image_view=AcquireVirtualCacheView(image,exception);
  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(similarity,status) \
    magick_number_threads(image,image,rows,1)
#endif
  for (y=0; y < (ssize_t) rows; y++)
  {
    const Quantum
      *magick_restrict p,
      *magick_restrict q;

    double
      channel_similarity[MaxPixelChannels+1] = { 0.0 };

    ssize_t
      x;

    if (status == MagickFalse)
      continue;
    p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
    q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
    if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (ssize_t) columns; x++)
    {
      double
        Da,
        Sa;

      size_t
        count = 0;

      ssize_t
        i;

      if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
          (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
        {
          p+=(ptrdiff_t) GetPixelChannels(image);
          q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
          continue;
        }
      Sa=QuantumScale*(double) GetPixelAlpha(image,p);
      Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
      {
        double
          error;

        PixelChannel channel = GetPixelChannelChannel(image,i);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
          channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((reconstruct_traits & UpdatePixelTrait) == 0))
          continue;
        if (channel == AlphaPixelChannel)
          error=(double) p[i]-(double) GetPixelChannel(reconstruct_image,
            channel,q);
        else
          error=Sa*(double) p[i]-Da*(double) GetPixelChannel(reconstruct_image,
            channel,q);
        if (MagickSafeSignificantError(error*error,fuzz) != MagickFalse)
          {
            channel_similarity[i]++;
            count++;
          }
      }
      if (count != 0)
        channel_similarity[CompositePixelChannel]++;
      p+=(ptrdiff_t) GetPixelChannels(image);
      q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
    }
#if defined(MAGICKCORE_OPENMP_SUPPORT)
    #pragma omp critical (MagickCore_GetPDCSimilarity)
#endif
    {
      ssize_t
        j;

      for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
      {
        PixelChannel channel = GetPixelChannelChannel(image,j);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
          channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((reconstruct_traits & UpdatePixelTrait) == 0))
          continue;
        similarity[j]+=channel_similarity[j];
      }
      similarity[CompositePixelChannel]+=
        channel_similarity[CompositePixelChannel];
    }
  }
  reconstruct_view=DestroyCacheView(reconstruct_view);
  image_view=DestroyCacheView(image_view);
  return(status);
}

static Image *GetEdgeCorrelationSurface(const Image *image,
  ExceptionInfo *exception)
{
  Image
    *surface;

  KernelInfo
    *kernel;

  /*
    Build a spatial-domain edge/correlation surface using a LoG convolution.
  */
  kernel=AcquireKernelInfo("LoG:0,1",exception);
  if (kernel == (KernelInfo *) NULL)
    return((Image *) NULL);
  surface=MorphologyImage(image,ConvolveMorphology,1,kernel,exception);
  kernel=DestroyKernelInfo(kernel);
  return(surface);
}

static MagickBooleanType GetPHASESimilarity(const Image *image,
  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
{
  CacheView
    *edge_reconstruct_view,
    *edge_view,
    *image_view,
    *reconstruct_view;

  double
    area = 0,
    correlation[MaxPixelChannels+1] = { 0.0 },
    image_sum[MaxPixelChannels+1] = { 0.0 },
    image_sum_squared[MaxPixelChannels+1] = { 0.0 },
    reconstruct_sum[MaxPixelChannels+1] = { 0.0 },
    reconstruct_sum_squared[MaxPixelChannels+1] = { 0.0 };

  Image
    *edge_image,
    *edge_reconstruct;

  MagickBooleanType
    status = MagickTrue;

  size_t
    columns = 0,
    rows = 0;

  ssize_t
    channels = 0,
    j,
    y;

  /*
    Compute a Pearson-correlation similarity between two Laplacian-of-
    Gaussian edge surfaces.
  */
  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
  edge_image=GetEdgeCorrelationSurface(image,exception);
  edge_reconstruct=GetEdgeCorrelationSurface(reconstruct_image,exception);
  if ((edge_image == (Image *) NULL) || (edge_reconstruct == (Image *) NULL))
    {
      if (edge_image != (Image *) NULL)
        edge_image=DestroyImage(edge_image);
      if (edge_reconstruct != (Image *) NULL)
        edge_reconstruct=DestroyImage(edge_reconstruct);
      return(MagickFalse);
    }
  image_view=AcquireVirtualCacheView(image,exception);
  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
  edge_view=AcquireVirtualCacheView(edge_image,exception);
  edge_reconstruct_view=AcquireVirtualCacheView(edge_reconstruct,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(status) \
    magick_number_threads(edge_image,edge_reconstruct,rows,1)
#endif
  for (y=0; y < (ssize_t) rows; y++)
  {
    const Quantum
      *magick_restrict p,
      *magick_restrict pm,
      *magick_restrict q,
      *magick_restrict qm;

    double
      channel_area = 0,
      channel_correlation[MaxPixelChannels+1] = { 0.0 },
      channel_image_sum[MaxPixelChannels+1] = { 0.0 },
      channel_image_sum_squared[MaxPixelChannels+1] = { 0.0 },
      channel_reconstruct_sum[MaxPixelChannels+1] = { 0.0 },
      channel_reconstruct_sum_squared[MaxPixelChannels+1] = { 0.0 };

    ssize_t
      i,
      x;

    if (status == MagickFalse)
      continue;
    pm=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
    qm=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
    p=GetCacheViewVirtualPixels(edge_view,0,y,columns,1,exception);
    q=GetCacheViewVirtualPixels(edge_reconstruct_view,0,y,columns,1,exception);
    if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL) ||
        (pm == (const Quantum *) NULL) || (qm == (const Quantum *) NULL))
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (ssize_t) columns; x++)
    {
      if ((GetPixelReadMask(image,pm) <= (QuantumRange/2)) ||
          (GetPixelReadMask(reconstruct_image,qm) <= (QuantumRange/2)))
        {
          p+=(ptrdiff_t) GetPixelChannels(edge_image);
          q+=(ptrdiff_t) GetPixelChannels(edge_reconstruct);
          pm+=(ptrdiff_t) GetPixelChannels(image);
          qm+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
          continue;
        }
      for (i=0; i < (ssize_t) GetPixelChannels(edge_image); i++)
      {
        double
          alpha,
          beta;

        PixelChannel
          channel;

        PixelTrait
          reconstruct_traits,
          traits;

        ssize_t
          offset;

        channel=GetPixelChannelChannel(edge_image,i);
        traits=GetPixelChannelTraits(edge_image,channel);
        reconstruct_traits=GetPixelChannelTraits(edge_reconstruct,channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((reconstruct_traits & UpdatePixelTrait) == 0))
          continue;
        offset=GetPixelChannelOffset(edge_reconstruct,channel);
        if (offset < 0)
          continue;
        alpha=QuantumScale*(double) p[i];
        beta=QuantumScale*(double) q[offset];
        channel_image_sum[i]+=alpha;
        channel_image_sum_squared[i]+=alpha*alpha;
        channel_reconstruct_sum[i]+=beta;
        channel_reconstruct_sum_squared[i]+=beta*beta;
        channel_correlation[i]+=alpha*beta;
      }
      channel_area++;
      p+=(ptrdiff_t) GetPixelChannels(edge_image);
      q+=(ptrdiff_t) GetPixelChannels(edge_reconstruct);
      pm+=(ptrdiff_t) GetPixelChannels(image);
      qm+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
    }
#if defined(MAGICKCORE_OPENMP_SUPPORT)
    #pragma omp critical (MagickCore_GetPHASESimilarity)
#endif
    {
      area+=channel_area;
      for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
      {
        correlation[i]+=channel_correlation[i];
        image_sum[i]+=channel_image_sum[i];
        image_sum_squared[i]+=channel_image_sum_squared[i];
        reconstruct_sum[i]+=channel_reconstruct_sum[i];
        reconstruct_sum_squared[i]+=channel_reconstruct_sum_squared[i];
      }
    }
  }
  edge_reconstruct_view=DestroyCacheView(edge_reconstruct_view);
  edge_view=DestroyCacheView(edge_view);
  reconstruct_view=DestroyCacheView(reconstruct_view);
  image_view=DestroyCacheView(image_view);
  edge_image=DestroyImage(edge_image);
  edge_reconstruct=DestroyImage(edge_reconstruct);
  if (status == MagickFalse)
    return(MagickFalse);
  if (area < 1.0)
    {
      (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
        "InsufficientImageDataInRaster","`%s'",image->filename);
      return(MagickFalse);
    }
  /*
    Reduce to a per-channel Pearson coefficient and average.
  */
  similarity[CompositePixelChannel]=0.0;
  for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
  {
    double
      denominator,
      image_variance,
      numerator,
      pearson,
      reconstruct_variance;

    PixelChannel
      channel;

    PixelTrait
      reconstruct_traits,
      traits;

    channel=GetPixelChannelChannel(image,j);
    traits=GetPixelChannelTraits(image,channel);
    reconstruct_traits=GetPixelChannelTraits(reconstruct_image,channel);
    if (((traits & UpdatePixelTrait) == 0) ||
        ((reconstruct_traits & UpdatePixelTrait) == 0))
      continue;
    numerator=area*correlation[j]-image_sum[j]*reconstruct_sum[j];
    image_variance=area*image_sum_squared[j]-image_sum[j]*image_sum[j];
    reconstruct_variance=area*reconstruct_sum_squared[j]-reconstruct_sum[j]*
      reconstruct_sum[j];
    if ((image_variance < MagickEpsilon) &&
        (reconstruct_variance < MagickEpsilon))
      pearson=(fabs(image_sum[j]-reconstruct_sum[j]) < MagickEpsilon) ?
        1.0 : 0.0;
    else
      {
        denominator=sqrt(image_variance)*sqrt(reconstruct_variance);
        pearson=denominator < MagickEpsilon ? 0.0 : numerator/denominator;
      }
    if (pearson < -1.0)
      pearson=(-1.0);
    if (pearson > 1.0)
      pearson=1.0;
    similarity[j]=pearson;
    similarity[CompositePixelChannel]+=pearson;
    channels++;
  }
  if (channels != 0)
    similarity[CompositePixelChannel]/=(double) channels;
  return(MagickTrue);
}

static MagickBooleanType GetPHASHSimilarity(const Image *image,
  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
{
  ChannelPerceptualHash
    *channel_phash,
    *reconstruct_phash;

  const char
    *artifact;

  ssize_t
    channels = 0,
    i;

  /*
    Compute the perceptual hash similarity.
  */
  channel_phash=GetImagePerceptualHash(image,exception);
  if (channel_phash == (ChannelPerceptualHash *) NULL)
    return(MagickFalse);
  reconstruct_phash=GetImagePerceptualHash(reconstruct_image,exception);
  if (reconstruct_phash == (ChannelPerceptualHash *) NULL)
    {
      channel_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
        channel_phash);
      return(MagickFalse);
    }
  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
  {
    double
      difference = 0.0;

    ssize_t
      j;

    PixelChannel channel = GetPixelChannelChannel(image,i);
    PixelTrait traits = GetPixelChannelTraits(image,channel);
    PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
      channel);
    if (((traits & UpdatePixelTrait) == 0) ||
        ((reconstruct_traits & UpdatePixelTrait) == 0))
      continue;
    for (j=0; j < (ssize_t) channel_phash[0].number_colorspaces; j++)
    {
      double
        alpha,
        beta;

      ssize_t
        k;

      for (k=0; k < MaximumNumberOfPerceptualHashes; k++)
      {
        double
          error;

        alpha=channel_phash[i].phash[j][k];
        beta=reconstruct_phash[i].phash[j][k];
        error=beta-alpha;
        if (IsNaN(error) != 0)
          error=0.0;
        difference+=error*error;
      }
    }
    similarity[i]+=difference;
    similarity[CompositePixelChannel]+=difference;
    channels++;
  }
  if (channels != 0)
    similarity[CompositePixelChannel]/=(double) channels;
  artifact=GetImageArtifact(image,"phash:normalize");
  if (IsStringTrue(artifact) != MagickFalse)
    {
      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 & UpdatePixelTrait) == 0) ||
            ((reconstruct_traits & UpdatePixelTrait) == 0))
          continue;
        similarity[i]=sqrt(similarity[i]/(double)
          channel_phash[0].number_colorspaces);
      }
      similarity[CompositePixelChannel]=sqrt(similarity[CompositePixelChannel]/
        (double) channel_phash[0].number_colorspaces);
    }
  /*
    Free resources.
  */
  reconstruct_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
    reconstruct_phash);
  channel_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(channel_phash);
  return(MagickTrue);
}

static MagickBooleanType GetPSNRSimilarity(const Image *image,
  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
{
  MagickBooleanType
    status = MagickTrue;

  ssize_t
    i;

  /*
    Compute the peak signal-to-noise ratio similarity.
  */
  status=GetMSESimilarity(image,reconstruct_image,similarity,exception);
  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 & UpdatePixelTrait) == 0) ||
        ((reconstruct_traits & UpdatePixelTrait) == 0))
      continue;
    similarity[i]=10.0*MagickSafeLog10(MagickSafeReciprocal(
      similarity[i]))/MagickSafePSNRRecipicol(10.0);
  }
  similarity[CompositePixelChannel]=10.0*MagickSafeLog10(
    MagickSafeReciprocal(similarity[CompositePixelChannel]))/
    MagickSafePSNRRecipicol(10.0);
  return(status);
}

static MagickBooleanType GetRMSESimilarity(const Image *image,
  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
{
#define RMSESquareRoot(x)  sqrt((x) < 0.0 ? 0.0 : (x))

  MagickBooleanType
    status = MagickTrue;

  ssize_t
    i;

  /*
    Compute the root mean-squared error similarity.
  */
  status=GetMSESimilarity(image,reconstruct_image,similarity,exception);
  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 & UpdatePixelTrait) == 0) ||
        ((reconstruct_traits & UpdatePixelTrait) == 0))
      continue;
    similarity[i]=RMSESquareRoot(similarity[i]);
  }
  similarity[CompositePixelChannel]=RMSESquareRoot(
    similarity[CompositePixelChannel]);
  return(status);
}

static MagickBooleanType GetSSIMSimularity(const Image *image,
  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
{
#define SSIMRadius  5.0
#define SSIMSigma  1.5
#define SSIMK1  0.01
#define SSIMK2  0.03
#define SSIML  1.0

  CacheView
    *image_view,
    *reconstruct_view;

  char
    geometry[MagickPathExtent];

  const char
    *artifact;

  double
    area = 0.0,
    c1,
    c2,
    radius,
    sigma;

  KernelInfo
    *kernel_info;

  MagickBooleanType
    status = MagickTrue;

  size_t
    columns,
    rows;

  ssize_t
    channels = 0,
    l,
    y;

  /*
    Compute the structual similarity index similarity.
  */
  radius=SSIMRadius;
  artifact=GetImageArtifact(image,"compare:ssim-radius");
  if (artifact != (const char *) NULL)
    radius=StringToDouble(artifact,(char **) NULL);
  sigma=SSIMSigma;
  artifact=GetImageArtifact(image,"compare:ssim-sigma");
  if (artifact != (const char *) NULL)
    sigma=StringToDouble(artifact,(char **) NULL);
  (void) FormatLocaleString(geometry,MagickPathExtent,"gaussian:%.17gx%.17g",
    radius,sigma);
  kernel_info=AcquireKernelInfo(geometry,exception);
  if (kernel_info == (KernelInfo *) NULL)
    ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
      image->filename);
  c1=pow(SSIMK1*SSIML,2.0);
  artifact=GetImageArtifact(image,"compare:ssim-k1");
  if (artifact != (const char *) NULL)
    c1=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
  c2=pow(SSIMK2*SSIML,2.0);
  artifact=GetImageArtifact(image,"compare:ssim-k2");
  if (artifact != (const char *) NULL)
    c2=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
  image_view=AcquireVirtualCacheView(image,exception);
  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(area,similarity,status) \
    magick_number_threads(image,reconstruct_image,rows,1)
#endif
  for (y=0; y < (ssize_t) rows; y++)
  {
    const Quantum
      *magick_restrict p,
      *magick_restrict q;

    double
      channel_area = 0.0,
      channel_similarity[MaxPixelChannels+1] = { 0.0 };

    ssize_t
      i,
      x;

    if (status == MagickFalse)
      continue;
    p=GetCacheViewVirtualPixels(image_view,-((ssize_t) kernel_info->width/2L),y-
      ((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
      kernel_info->height,exception);
    q=GetCacheViewVirtualPixels(reconstruct_view,-((ssize_t) kernel_info->width/
      2L),y-((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
      kernel_info->height,exception);
    if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (ssize_t) columns; x++)
    {
      const Quantum
        *magick_restrict reconstruct,
        *magick_restrict test;

      double
        x_pixel_mu[MaxPixelChannels+1] = { 0.0 },
        x_pixel_sigma_squared[MaxPixelChannels+1] = { 0.0 },
        xy_sigma[MaxPixelChannels+1] = { 0.0 },
        y_pixel_mu[MaxPixelChannels+1] = { 0.0 },
        y_pixel_sigma_squared[MaxPixelChannels+1] = { 0.0 };

      MagickRealType
        *k;

      ssize_t
        v;

      if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
          (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
        {
          p+=(ptrdiff_t) GetPixelChannels(image);
          q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
          continue;
        }
      k=kernel_info->values;
      test=p;
      reconstruct=q;
      for (v=0; v < (ssize_t) kernel_info->height; v++)
      {
        ssize_t
          u;

        for (u=0; u < (ssize_t) kernel_info->width; u++)
        {
          for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
          {
            double
              x_pixel,
              y_pixel;

            PixelChannel channel = GetPixelChannelChannel(image,i);
            PixelTrait traits = GetPixelChannelTraits(image,channel);
            PixelTrait reconstruct_traits = GetPixelChannelTraits(
              reconstruct_image,channel);
            if (((traits & UpdatePixelTrait) == 0) ||
                ((reconstruct_traits & UpdatePixelTrait) == 0))
              continue;
            x_pixel=QuantumScale*(double) test[i];
            x_pixel_mu[i]+=(*k)*x_pixel;
            x_pixel_sigma_squared[i]+=(*k)*x_pixel*x_pixel;
            y_pixel=QuantumScale*(double)
              GetPixelChannel(reconstruct_image,channel,reconstruct);
            y_pixel_mu[i]+=(*k)*y_pixel;
            y_pixel_sigma_squared[i]+=(*k)*y_pixel*y_pixel;
            xy_sigma[i]+=(*k)*x_pixel*y_pixel;
          }
          k++;
          test+=(ptrdiff_t) GetPixelChannels(image);
          reconstruct+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
        }
        test+=(ptrdiff_t) ((double) GetPixelChannels(image)*(double) columns);
        reconstruct+=(ptrdiff_t) ((double) GetPixelChannels(reconstruct_image)*
          (double) columns);
      }
      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
      {
        double
          ssim,
          x_pixel_mu_squared,
          x_pixel_sigmas_squared,
          xy_mu,
          xy_sigmas,
          y_pixel_mu_squared,
          y_pixel_sigmas_squared;

        PixelChannel channel = GetPixelChannelChannel(image,i);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        PixelTrait reconstruct_traits = GetPixelChannelTraits(
          reconstruct_image,channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((reconstruct_traits & UpdatePixelTrait) == 0))
          continue;
        x_pixel_mu_squared=x_pixel_mu[i]*x_pixel_mu[i];
        y_pixel_mu_squared=y_pixel_mu[i]*y_pixel_mu[i];
        xy_mu=x_pixel_mu[i]*y_pixel_mu[i];
        xy_sigmas=xy_sigma[i]-xy_mu;
        x_pixel_sigmas_squared=x_pixel_sigma_squared[i]-x_pixel_mu_squared;
        y_pixel_sigmas_squared=y_pixel_sigma_squared[i]-y_pixel_mu_squared;
        ssim=((2.0*xy_mu+c1)*(2.0*xy_sigmas+c2))*
          MagickSafeReciprocal((x_pixel_mu_squared+y_pixel_mu_squared+c1)*
           (x_pixel_sigmas_squared+y_pixel_sigmas_squared+c2));
        channel_similarity[i]+=ssim;
        channel_similarity[CompositePixelChannel]+=ssim;
      }
      p+=(ptrdiff_t) GetPixelChannels(image);
      q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
      channel_area++;
    }
#if defined(MAGICKCORE_OPENMP_SUPPORT)
    #pragma omp critical (MagickCore_GetSSIMSimularity)
#endif
    {
      ssize_t
        j;

      area+=channel_area;
      for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
      {
        PixelChannel channel = GetPixelChannelChannel(image,j);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
          channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((reconstruct_traits & UpdatePixelTrait) == 0))
          continue;
        similarity[j]+=channel_similarity[j];
      }
      similarity[CompositePixelChannel]+=
        channel_similarity[CompositePixelChannel];
    }
  }
  image_view=DestroyCacheView(image_view);
  reconstruct_view=DestroyCacheView(reconstruct_view);
  kernel_info=DestroyKernelInfo(kernel_info);
  area=MagickSafeReciprocal(area);
  for (l=0; l < (ssize_t) GetPixelChannels(image); l++)
  {
    PixelChannel channel = GetPixelChannelChannel(image,l);
    PixelTrait traits = GetPixelChannelTraits(image,channel);
    PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
      channel);
    if (((traits & UpdatePixelTrait) == 0) ||
        ((reconstruct_traits & UpdatePixelTrait) == 0))
      continue;
    similarity[l]*=area;
    channels++;
  }
  similarity[CompositePixelChannel]*=area;
  if (channels != 0)
   similarity[CompositePixelChannel]/=(double) channels;
  return(status);
}

static MagickBooleanType GetDSSIMSimilarity(const Image *image,
  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
{
  MagickBooleanType
    status = MagickTrue;

  ssize_t
    i;

  /*
    Compute the structual dissimilarity index similarity.
  */
  status=GetSSIMSimularity(image,reconstruct_image,similarity,exception);
  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 & UpdatePixelTrait) == 0) ||
        ((reconstruct_traits & UpdatePixelTrait) == 0))
      continue;
    similarity[i]=(1.0-similarity[i])/2.0;
  }
  similarity[CompositePixelChannel]=(1.0-similarity[CompositePixelChannel])/2.0;
  return(status);
}

MagickExport MagickBooleanType GetImageDistortion(Image *image,
  const Image *reconstruct_image,const MetricType metric,double *distortion,
  ExceptionInfo *exception)
{
#define CompareMetricNotSupportedException  "metric not supported"

  double
    *channel_similarity;

  MagickBooleanType
    status = MagickTrue;

  size_t
    length;

  assert(image != (Image *) NULL);
  assert(image->signature == MagickCoreSignature);
  assert(reconstruct_image != (const Image *) NULL);
  assert(reconstruct_image->signature == MagickCoreSignature);
  assert(distortion != (double *) NULL);
  if (IsEventLogging() != MagickFalse)
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
  /*
    Get image distortion.
  */
  *distortion=0.0;
  length=MaxPixelChannels+1UL;
  channel_similarity=(double *) AcquireQuantumMemory(length,
    sizeof(*channel_similarity));
  if (channel_similarity == (double *) NULL)
    ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
  (void) memset(channel_similarity,0,length*sizeof(*channel_similarity));
  switch (metric)
  {
    case AbsoluteErrorMetric:
    {
      status=GetAESimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case DotProductCorrelationErrorMetric:
    {
      status=GetDPCSimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case FuzzErrorMetric:
    {
      status=GetFUZZSimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case MeanAbsoluteErrorMetric:
    {
      status=GetMAESimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case MeanErrorPerPixelErrorMetric:
    {
      status=GetMEPPSimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case MeanSquaredErrorMetric:
    {
      status=GetMSESimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case NormalizedCrossCorrelationErrorMetric:
    {
      status=GetNCCSimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case PeakAbsoluteErrorMetric:
    {
      status=GetPASimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case PeakSignalToNoiseRatioErrorMetric:
    {
      status=GetPSNRSimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case PerceptualHashErrorMetric:
    {
      status=GetPHASHSimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case PhaseCorrelationErrorMetric:
    {
      status=GetPHASESimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case PixelDifferenceCountErrorMetric:
    {
      status=GetPDCSimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case RootMeanSquaredErrorMetric:
    case UndefinedErrorMetric:
    default:
    {
      status=GetRMSESimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case StructuralDissimilarityErrorMetric:
    {
      status=GetDSSIMSimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case StructuralSimilarityErrorMetric:
    {
      status=GetSSIMSimularity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
  }
  *distortion=channel_similarity[CompositePixelChannel];
  switch (metric)
  {
    case DotProductCorrelationErrorMetric:
    case NormalizedCrossCorrelationErrorMetric:
    case PhaseCorrelationErrorMetric:
    case StructuralSimilarityErrorMetric:
    {
      *distortion=(1.0-(*distortion))/2.0;
      break;
    }
    default: break;
  }
  channel_similarity=(double *) RelinquishMagickMemory(channel_similarity);
  if ((IsNaN(*distortion) != 0) || (fabs(*distortion) < MagickEpsilon))
    *distortion=0.0;
  (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
    *distortion);
  return(status);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   G e t I m a g e D i s t o r t i o n s                                     %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  GetImageDistortions() compares the pixel channels of an image to a
%  reconstructed image and returns the specified metric for each channel.
%
%  The format of the GetImageDistortions method is:
%
%      double *GetImageDistortions(const Image *image,
%        const Image *reconstruct_image,const MetricType metric,
%        ExceptionInfo *exception)
%
%  A description of each parameter follows:
%
%    o image: the image.
%
%    o reconstruct_image: the reconstruction image.
%
%    o metric: the metric.
%
%    o exception: return any errors or warnings in this structure.
%
*/
MagickExport double *GetImageDistortions(Image *image,
  const Image *reconstruct_image,const MetricType metric,
  ExceptionInfo *exception)
{
  double
    *distortion,
    *channel_similarity;

  MagickBooleanType
    status = MagickTrue;

  size_t
    length;

  ssize_t
    i;

  assert(image != (Image *) NULL);
  assert(image->signature == MagickCoreSignature);
  assert(reconstruct_image != (const Image *) NULL);
  assert(reconstruct_image->signature == MagickCoreSignature);
  if (IsEventLogging() != MagickFalse)
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
  /*
    Get image distortion.
  */
  length=MaxPixelChannels+1UL;
  channel_similarity=(double *) AcquireQuantumMemory(length,
    sizeof(*channel_similarity));
  if (channel_similarity == (double *) NULL)
    ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
  (void) memset(channel_similarity,0,length*sizeof(*channel_similarity));
  switch (metric)
  {
    case AbsoluteErrorMetric:
    {
      status=GetAESimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case DotProductCorrelationErrorMetric:
    {
      status=GetDPCSimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case FuzzErrorMetric:
    {
      status=GetFUZZSimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case MeanAbsoluteErrorMetric:
    {
      status=GetMAESimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case MeanErrorPerPixelErrorMetric:
    {
      status=GetMEPPSimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case MeanSquaredErrorMetric:
    {
      status=GetMSESimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case NormalizedCrossCorrelationErrorMetric:
    {
      status=GetNCCSimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case PeakAbsoluteErrorMetric:
    {
      status=GetPASimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case PeakSignalToNoiseRatioErrorMetric:
    {
      status=GetPSNRSimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case PerceptualHashErrorMetric:
    {
      status=GetPHASHSimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case PhaseCorrelationErrorMetric:
    {
      status=GetPHASESimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case PixelDifferenceCountErrorMetric:
    {
      status=GetPDCSimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case RootMeanSquaredErrorMetric:
    case UndefinedErrorMetric:
    default:
    {
      status=GetRMSESimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case StructuralDissimilarityErrorMetric:
    {
      status=GetDSSIMSimilarity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
    case StructuralSimilarityErrorMetric:
    {
      status=GetSSIMSimularity(image,reconstruct_image,channel_similarity,
        exception);
      break;
    }
  }
  if (status == MagickFalse)
    {
      channel_similarity=(double *) RelinquishMagickMemory(channel_similarity);
      return((double *) NULL);
    }
  distortion=channel_similarity;
  switch (metric)
  {
    case DotProductCorrelationErrorMetric:
    case NormalizedCrossCorrelationErrorMetric:
    case PhaseCorrelationErrorMetric:
    case StructuralSimilarityErrorMetric:
    {
      for (i=0; i <= MaxPixelChannels; i++)
        distortion[i]=(1.0-distortion[i])/2.0;
      break;
    }
    default: break;
  }
  for (i=0; i <= MaxPixelChannels; i++)
    if (fabs(distortion[i]) < MagickEpsilon)
      distortion[i]=0.0;
  (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
    distortion[CompositePixelChannel]);
  return(distortion);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%  I s I m a g e s E q u a l                                                  %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  IsImagesEqual() compare the pixels of two images and returns immediately
%  if any pixel is not identical.
%
%  The format of the IsImagesEqual method is:
%
%      MagickBooleanType IsImagesEqual(const Image *image,
%        const Image *reconstruct_image,ExceptionInfo *exception)
%
%  A description of each parameter follows.
%
%    o image: the image.
%
%    o reconstruct_image: the reconstruction image.
%
%    o exception: return any errors or warnings in this structure.
%
*/
MagickExport MagickBooleanType IsImagesEqual(const Image *image,
  const Image *reconstruct_image,ExceptionInfo *exception)
{
  CacheView
    *image_view,
    *reconstruct_view;

  size_t
    columns,
    rows;

  ssize_t
    y;

  assert(image != (Image *) NULL);
  assert(image->signature == MagickCoreSignature);
  assert(reconstruct_image != (const Image *) NULL);
  assert(reconstruct_image->signature == MagickCoreSignature);
  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
  image_view=AcquireVirtualCacheView(image,exception);
  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
  for (y=0; y < (ssize_t) rows; y++)
  {
    const Quantum
      *magick_restrict p,
      *magick_restrict q;

    ssize_t
      x;

    p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
    q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
    if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
      break;
    for (x=0; x < (ssize_t) columns; x++)
    {
      ssize_t
        i;

      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
      {
        double
          distance;

        PixelChannel channel = GetPixelChannelChannel(image,i);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
          channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((reconstruct_traits & UpdatePixelTrait) == 0))
          continue;
        distance=fabs((double) p[i]-(double) GetPixelChannel(reconstruct_image,
          channel,q));
        if (distance >= MagickEpsilon)
          break;
      }
      if (i < (ssize_t) GetPixelChannels(image))
        break;
      p+=(ptrdiff_t) GetPixelChannels(image);
      q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
    }
    if (x < (ssize_t) columns)
      break;
  }
  reconstruct_view=DestroyCacheView(reconstruct_view);
  image_view=DestroyCacheView(image_view);
  return(y < (ssize_t) rows ? MagickFalse : MagickTrue);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%  S e t I m a g e C o l o r M e t r i c                                      %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  SetImageColorMetric() measures the difference between colors at each pixel
%  location of two images.  A value other than 0 means the colors match
%  exactly.  Otherwise an error measure is computed by summing over all
%  pixels in an image the distance squared in RGB space between each image
%  pixel and its corresponding pixel in the reconstruction image.  The error
%  measure is assigned to these image members:
%
%    o mean_error_per_pixel:  The mean error for any single pixel in
%      the image.
%
%    o normalized_mean_error:  The normalized mean quantization error for
%      any single pixel in the image.  This distance measure is normalized to
%      a range between 0 and 1.  It is independent of the range of red, green,
%      and blue values in the image.
%
%    o normalized_maximum_error:  The normalized maximum quantization
%      error for any single pixel in the image.  This distance measure is
%      normalized to a range between 0 and 1.  It is independent of the range
%      of red, green, and blue values in your image.
%
%  A small normalized mean square error, accessed as
%  image->normalized_mean_error, suggests the images are very similar in
%  spatial layout and color.
%
%  The format of the SetImageColorMetric method is:
%
%      MagickBooleanType SetImageColorMetric(Image *image,
%        const Image *reconstruct_image,ExceptionInfo *exception)
%
%  A description of each parameter follows.
%
%    o image: the image.
%
%    o reconstruct_image: the reconstruction image.
%
%    o exception: return any errors or warnings in this structure.
%
*/
MagickExport MagickBooleanType SetImageColorMetric(Image *image,
  const Image *reconstruct_image,ExceptionInfo *exception)
{
  double
    channel_similarity[MaxPixelChannels+1] = { 0.0 };

  MagickBooleanType
    status;

  status=GetMEPPSimilarity(image,reconstruct_image,channel_similarity,
    exception);
  if (status == MagickFalse)
    return(MagickFalse);
  status=fabs(image->error.mean_error_per_pixel) < MagickEpsilon ?
    MagickTrue : MagickFalse;
  return(status);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   S i m i l a r i t y I m a g e                                             %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  SimilarityImage() compares the reconstruction of the image and returns the
%  best match offset.  In addition, it returns a similarity image such that an
%  exact match location is completely white and if none of the pixels match,
%  black, otherwise some gray level in-between.
%
%  Contributed by Fred Weinhaus.
%
%  The format of the SimilarityImageImage method is:
%
%      Image *SimilarityImage(const Image *image,const Image *reconstruct,
%        const MetricType metric,const double similarity_threshold,
%        RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
%
%  A description of each parameter follows:
%
%    o image: the image.
%
%    o reconstruct: find an area of the image that closely resembles this image.
%
%    o metric: the metric.
%
%    o similarity_threshold: minimum similarity for (sub)image match.
%
%    o offset: the best match offset of the reconstruction image within the
%      image.
%
%    o similarity: the computed similarity between the images.
%
%    o exception: return any errors or warnings in this structure.
%
*/

#if defined(MAGICKCORE_HDRI_SUPPORT) && defined(MAGICKCORE_FFTW_DELEGATE)
static Image *SIMCrossCorrelationImage(const Image *alpha_image,
  const Image *beta_image,ExceptionInfo *exception)
{
  Image
    *alpha_fft = (Image *) NULL,
    *beta_fft = (Image *) NULL,
    *complex_conjugate = (Image *) NULL,
    *complex_multiplication = (Image *) NULL,
    *cross_correlation = (Image *) NULL,
    *temp_image = (Image *) NULL;

  /*
    Take the FFT of beta (reconstruction) image.
  */
  temp_image=CloneImage(beta_image,0,0,MagickTrue,exception);
  if (temp_image == (Image *) NULL)
    return((Image *) NULL);
  (void) SetImageArtifact(temp_image,"fourier:normalize","inverse");
  beta_fft=ForwardFourierTransformImage(temp_image,MagickFalse,exception);
  temp_image=DestroyImageList(temp_image);
  if (beta_fft == (Image *) NULL)
    return((Image *) NULL);
  /*
    Take the complex conjugate of beta_fft.
  */
  complex_conjugate=ComplexImages(beta_fft,ConjugateComplexOperator,exception);
  beta_fft=DestroyImageList(beta_fft);
  if (complex_conjugate == (Image *) NULL)
    return((Image *) NULL);
  /*
    Take the FFT of the alpha (test) image.
  */
  temp_image=CloneImage(alpha_image,0,0,MagickTrue,exception);
  if (temp_image == (Image *) NULL)
    {
      complex_conjugate=DestroyImageList(complex_conjugate);
      return((Image *) NULL);
    }
  (void) SetImageArtifact(temp_image,"fourier:normalize","inverse");
  alpha_fft=ForwardFourierTransformImage(temp_image,MagickFalse,exception);
  temp_image=DestroyImageList(temp_image);
  if (alpha_fft == (Image *) NULL)
    {
      complex_conjugate=DestroyImageList(complex_conjugate);
      return((Image *) NULL);
    }
  /*
    Do complex multiplication.
  */
  DisableCompositeClampUnlessSpecified(complex_conjugate);
  DisableCompositeClampUnlessSpecified(complex_conjugate->next);
  AppendImageToList(&complex_conjugate,alpha_fft);
  complex_multiplication=ComplexImages(complex_conjugate,
    MultiplyComplexOperator,exception);
  complex_conjugate=DestroyImageList(complex_conjugate);
  if (complex_multiplication == (Image *) NULL)
    return((Image *) NULL);
  /*
    Do the IFT and return the cross-correlation result.
  */
  cross_correlation=InverseFourierTransformImage(complex_multiplication,
    complex_multiplication->next,MagickFalse,exception);
  complex_multiplication=DestroyImageList(complex_multiplication);
  return(cross_correlation);
}

static Image *SIMDerivativeImage(const Image *image,const char *kernel,
  ExceptionInfo *exception)
{
  Image
    *derivative_image;

  KernelInfo
    *kernel_info;

  kernel_info=AcquireKernelInfo(kernel,exception);
  if (kernel_info == (KernelInfo *) NULL)
    return((Image *) NULL);
  derivative_image=MorphologyImage(image,ConvolveMorphology,1,kernel_info,
    exception);
  kernel_info=DestroyKernelInfo(kernel_info);
  return(derivative_image);
}

static Image *SIMDivideImage(const Image *numerator_image,
  const Image *denominator_image,ExceptionInfo *exception)
{
  CacheView
    *denominator_view,
    *numerator_view;

  Image
    *divide_image;

  MagickBooleanType
    status = MagickTrue;

  ssize_t
    y;

  /*
    Divide one image into another.
  */
  divide_image=CloneImage(numerator_image,0,0,MagickTrue,exception);
  if (divide_image == (Image *) NULL)
    return(divide_image);
  numerator_view=AcquireAuthenticCacheView(divide_image,exception);
  denominator_view=AcquireVirtualCacheView(denominator_image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(status) \
    magick_number_threads(denominator_image,divide_image,divide_image->rows,1)
#endif
  for (y=0; y < (ssize_t) divide_image->rows; y++)
  {
    const Quantum
      *magick_restrict p;

    Quantum
      *magick_restrict q;

    ssize_t
      x;

    if (status == MagickFalse)
      continue;
    p=GetCacheViewVirtualPixels(denominator_view,0,y,
      denominator_image->columns,1,exception);
    q=GetCacheViewAuthenticPixels(numerator_view,0,y,divide_image->columns,1,
      exception);
    if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (ssize_t) divide_image->columns; x++)
    {
      ssize_t
        i;

      for (i=0; i < (ssize_t) GetPixelChannels(divide_image); i++)
      {
        PixelChannel channel = GetPixelChannelChannel(divide_image,i);
        PixelTrait traits = GetPixelChannelTraits(divide_image,channel);
        PixelTrait denominator_traits = GetPixelChannelTraits(denominator_image,
          channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((denominator_traits & UpdatePixelTrait) == 0))
          continue;
        q[i]=(Quantum) ((double) q[i]*MagickSafeReciprocal(QuantumScale*
          (double) GetPixelChannel(denominator_image,channel,p)));
      }
      p+=(ptrdiff_t) GetPixelChannels(denominator_image);
      q+=(ptrdiff_t) GetPixelChannels(divide_image);
    }
    if (SyncCacheViewAuthenticPixels(numerator_view,exception) == MagickFalse)
      status=MagickFalse;
  }
  denominator_view=DestroyCacheView(denominator_view);
  numerator_view=DestroyCacheView(numerator_view);
  if (status == MagickFalse)
    divide_image=DestroyImage(divide_image);
  return(divide_image);
}

static Image *SIMDivideByMagnitude(Image *image,Image *magnitude_image,
  const Image *source_image,ExceptionInfo *exception)
{
  Image
    *divide_image,
    *result_image;

  RectangleInfo
    geometry;

  divide_image=SIMDivideImage(image,magnitude_image,exception);
  if (divide_image == (Image *) NULL)
    return((Image *) NULL);
  GetPixelInfoRGBA((Quantum) 0,(Quantum) 0,(Quantum) 0,(Quantum) 0,
    &divide_image->background_color);
  SetGeometry(source_image,&geometry);
  geometry.width=MagickMax(source_image->columns,divide_image->columns);
  geometry.height=MagickMax(source_image->rows,divide_image->rows);
  result_image=ExtentImage(divide_image,&geometry,exception);
  divide_image=DestroyImage(divide_image);
  return(result_image);
}

static MagickBooleanType SIMFilterImageNaNs(Image *image,
  ExceptionInfo *exception)
{
  CacheView
    *image_view;

  MagickBooleanType
    status = MagickTrue;

  ssize_t
    y;

  /*
    Square each pixel in the image.
  */
  image_view=AcquireAuthenticCacheView(image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(status) \
    magick_number_threads(image,image,image->rows,1)
#endif
  for (y=0; y < (ssize_t) image->rows; y++)
  {
    Quantum
      *magick_restrict q;

    ssize_t
      x;

    if (status == MagickFalse)
      continue;
    q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
    if (q == (Quantum *) NULL)
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (ssize_t) image->columns; x++)
    {
      ssize_t
        i;

      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
      {
        PixelChannel channel = GetPixelChannelChannel(image,i);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        if ((traits & UpdatePixelTrait) == 0)
          continue;
        if (IsNaN((double) q[i]) != 0)
          q[i]=(Quantum) 0;
      }
      q+=(ptrdiff_t) GetPixelChannels(image);
    }
    if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
      status=MagickFalse;
  }
  image_view=DestroyCacheView(image_view);
  return(status);
}

static Image *SIMSquareImage(const Image *image,ExceptionInfo *exception)
{
  CacheView
    *image_view;

  Image
    *square_image;

  MagickBooleanType
    status = MagickTrue;

  ssize_t
    y;

  /*
    Square each pixel in the image.
  */
  square_image=CloneImage(image,0,0,MagickTrue,exception);
  if (square_image == (Image *) NULL)
    return(square_image);
  image_view=AcquireAuthenticCacheView(square_image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(status) \
    magick_number_threads(square_image,square_image,square_image->rows,1)
#endif
  for (y=0; y < (ssize_t) square_image->rows; y++)
  {
    Quantum
      *magick_restrict q;

    ssize_t
      x;

    if (status == MagickFalse)
      continue;
    q=GetCacheViewAuthenticPixels(image_view,0,y,square_image->columns,1,
      exception);
    if (q == (Quantum *) NULL)
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (ssize_t) square_image->columns; x++)
    {
      ssize_t
        i;

      for (i=0; i < (ssize_t) GetPixelChannels(square_image); i++)
      {
        PixelChannel channel = GetPixelChannelChannel(square_image,i);
        PixelTrait traits = GetPixelChannelTraits(square_image,channel);
        if ((traits & UpdatePixelTrait) == 0)
          continue;
        q[i]=(Quantum) (QuantumScale*(double) q[i]*(double) q[i]);
      }
      q+=(ptrdiff_t) GetPixelChannels(square_image);
    }
    if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
      status=MagickFalse;
  }
  image_view=DestroyCacheView(image_view);
  if (status == MagickFalse)
    square_image=DestroyImage(square_image);
  return(square_image);
}

static Image *SIMMagnitudeImage(Image *alpha_image,Image *beta_image,
  ExceptionInfo *exception)
{
  Image
    *magnitude_image,
    *xsq_image,
    *ysq_image;

  MagickBooleanType
    status = MagickTrue;

  (void) SetImageArtifact(alpha_image,"compose:clamp","False");
  xsq_image=SIMSquareImage(alpha_image,exception);
  if (xsq_image == (Image *) NULL)
    return((Image *) NULL);
  (void) SetImageArtifact(beta_image,"compose:clamp","False");
  ysq_image=SIMSquareImage(beta_image,exception);
  if (ysq_image == (Image *) NULL)
    {
      xsq_image=DestroyImage(xsq_image);
      return((Image *) NULL);
    }
  status=CompositeImage(xsq_image,ysq_image,PlusCompositeOp,MagickTrue,0,0,
    exception);
  magnitude_image=xsq_image;
  ysq_image=DestroyImage(ysq_image);
  if (status == MagickFalse)
    {
      magnitude_image=DestroyImage(magnitude_image);
      return((Image *) NULL);
    }
  status=EvaluateImage(magnitude_image,PowEvaluateOperator,0.5,exception);
  if (status == MagickFalse)
    {
      magnitude_image=DestroyImage(magnitude_image);
      return (Image *) NULL;
    }
  return(magnitude_image);
}

static MagickBooleanType SIMMaximaImage(const Image *image,double *maxima,
  RectangleInfo *offset,ExceptionInfo *exception)
{
  typedef struct
  {
    double
      maxima;

    ssize_t
      x,
      y;
  } MaximaInfo;

  CacheView
    *image_view;

  const Quantum
    *magick_restrict q;

  MagickBooleanType
    status = MagickTrue;

  MaximaInfo
    maxima_info = { -MagickMaximumValue, 0, 0 };

  ssize_t
    y;

  /*
    Identify the maxima value in the image and its location.
  */
  image_view=AcquireVirtualCacheView(image,exception);
  q=GetCacheViewVirtualPixels(image_view,maxima_info.x,maxima_info.y,1,1,
    exception);
  if (q != (const Quantum *) NULL)
    maxima_info.maxima=IsNaN((double) q[0]) != 0 ? 0.0 : (double) q[0];
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(maxima_info,status) \
    magick_number_threads(image,image,image->rows,1)
#endif
  for (y=0; y < (ssize_t) image->rows; y++)
  {
    const Quantum
      *magick_restrict p;

    MaximaInfo
      channel_maxima = { -MagickMaximumValue, 0, 0 };

    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;
      }
    channel_maxima=maxima_info;
    for (x=0; x < (ssize_t) image->columns; x++)
    {
      ssize_t
        i;

      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
      {
        double
          pixel;

        PixelChannel channel = GetPixelChannelChannel(image,i);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        if ((traits & UpdatePixelTrait) == 0)
          continue;
        pixel=(double) p[i];
        if (IsNaN(pixel) != 0)
          pixel=0.0;
        if (pixel > channel_maxima.maxima)
          {
            channel_maxima.maxima=(double) p[i];
            channel_maxima.x=x;
            channel_maxima.y=y;
          }
      }
      p+=(ptrdiff_t) GetPixelChannels(image);
    }
#if defined(MAGICKCORE_OPENMP_SUPPORT)
    #pragma omp critical (MagickCore_SIMMaximaImage)
#endif
    if (channel_maxima.maxima > maxima_info.maxima)
      maxima_info=channel_maxima;
  }
  image_view=DestroyCacheView(image_view);
  *maxima=maxima_info.maxima;
  offset->x=maxima_info.x;
  offset->y=maxima_info.y;
  return(status);
}

static MagickBooleanType SIMMinimaImage(const Image *image,double *minima,
  RectangleInfo *offset,ExceptionInfo *exception)
{
  typedef struct
  {
    double
      minima;

    ssize_t
      x,
      y;
  } MinimaInfo;

  CacheView
    *image_view;

  const Quantum
    *magick_restrict q;

  MagickBooleanType
    status = MagickTrue;

  MinimaInfo
    minima_info = { MagickMaximumValue, 0, 0 };

  ssize_t
    y;

  /*
    Identify the minima value in the image and its location.
  */
  image_view=AcquireVirtualCacheView(image,exception);
  q=GetCacheViewVirtualPixels(image_view,minima_info.x,minima_info.y,1,1,
    exception);
  if (q != (const Quantum *) NULL)
    minima_info.minima=IsNaN((double) q[0]) != 0 ? 0.0 : (double) q[0];
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(minima_info,status) \
    magick_number_threads(image,image,image->rows,1)
#endif
  for (y=0; y < (ssize_t) image->rows; y++)
  {
    const Quantum
      *magick_restrict p;

    MinimaInfo
      channel_minima = { MagickMaximumValue, 0, 0 };

    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;
      }
    channel_minima=minima_info;
    for (x=0; x < (ssize_t) image->columns; x++)
    {
      ssize_t
        i;

      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
      {
        double
          pixel;

        PixelChannel channel = GetPixelChannelChannel(image,i);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        if ((traits & UpdatePixelTrait) == 0)
          continue;
        pixel=(double) p[i];
        if (IsNaN(pixel) != 0)
          pixel=0.0;
        if (pixel < channel_minima.minima)
          {
            channel_minima.minima=pixel;
            channel_minima.x=x;
            channel_minima.y=y;
          }
      }
      p+=(ptrdiff_t) GetPixelChannels(image);
    }
#if defined(MAGICKCORE_OPENMP_SUPPORT)
    #pragma omp critical (MagickCore_SIMMinimaImage)
#endif
    if (channel_minima.minima < minima_info.minima)
      minima_info=channel_minima;
  }
  image_view=DestroyCacheView(image_view);
  *minima=minima_info.minima;
  offset->x=minima_info.x;
  offset->y=minima_info.y;
  return(status);
}

static MagickBooleanType SIMMultiplyImage(Image *image,const double factor,
  const ChannelStatistics *channel_statistics,ExceptionInfo *exception)
{
  CacheView
    *image_view;

  MagickBooleanType
    status = MagickTrue;

  ssize_t
    y;

  /*
    Multiply each pixel by a factor.
  */
  image_view=AcquireAuthenticCacheView(image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(status) \
    magick_number_threads(image,image,image->rows,1)
#endif
  for (y=0; y < (ssize_t) image->rows; y++)
  {
    Quantum
      *magick_restrict q;

    ssize_t
      x;

    if (status == MagickFalse)
      continue;
    q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
    if (q == (Quantum *) NULL)
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (ssize_t) image->columns; x++)
    {
      ssize_t
        i;

      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
      {
        PixelChannel channel = GetPixelChannelChannel(image,i);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        if ((traits & UpdatePixelTrait) == 0)
          continue;
        if (channel_statistics != (const ChannelStatistics *) NULL)
          q[i]=(Quantum) (factor*(double) q[i]*QuantumScale*
            channel_statistics[channel].standard_deviation);
        else
          q[i]=(Quantum) (factor*(double) q[i]);
      }
      q+=(ptrdiff_t) GetPixelChannels(image);
    }
    if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
      status=MagickFalse;
  }
  image_view=DestroyCacheView(image_view);
  return(status);
}

static Image *SIMPhaseCorrelationImage(const Image *target_image,
  const Image *reconstruct_image,const Image *magnitude_image,
  ExceptionInfo *exception)
{
  Image
    *target_fft = (Image *) NULL,
    *reconstruct_fft = (Image *) NULL,
    *complex_multiplication = (Image *) NULL,
    *cross_correlation = (Image *) NULL;

  /*
    Take the FFT of the reconstruction image.
  */
  reconstruct_fft=CloneImage(reconstruct_image,0,0,MagickTrue,exception);
  if (reconstruct_fft == NULL)
    return((Image *) NULL);
  (void) SetImageArtifact(reconstruct_fft,"fourier:normalize","inverse");
  reconstruct_fft=ForwardFourierTransformImage(reconstruct_fft,MagickFalse,
    exception);
  if (reconstruct_fft == NULL)
    return((Image *) NULL);
  /*
    Take the FFT of the target image.
  */
  target_fft=CloneImage(target_image,0,0,MagickTrue,exception);
  if (target_fft == (Image *) NULL)
    {
      reconstruct_fft=DestroyImageList(reconstruct_fft);
      return((Image *) NULL);
    }
  (void) SetImageArtifact(target_fft,"fourier:normalize","inverse");
  target_fft=ForwardFourierTransformImage(target_fft,MagickFalse,exception);
  if (target_fft == (Image *) NULL)
    {
      reconstruct_fft=DestroyImageList(reconstruct_fft);
      return((Image *) NULL);
    }
  /*
    Take the complex conjugate of the reconstruction FFT.
  */
  reconstruct_fft=ComplexImages(reconstruct_fft,ConjugateComplexOperator,
    exception);
  if (reconstruct_fft == (Image *) NULL)
    {
      target_fft=DestroyImageList(target_fft);
      return((Image *) NULL);
    }
  /*
    Do complex multiplication.
  */
  AppendImageToList(&reconstruct_fft,target_fft);
  DisableCompositeClampUnlessSpecified(reconstruct_fft);
  DisableCompositeClampUnlessSpecified(reconstruct_fft->next);
  complex_multiplication=ComplexImages(reconstruct_fft,MultiplyComplexOperator,
    exception);
  reconstruct_fft=DestroyImageList(reconstruct_fft);
  if (complex_multiplication == (Image *) NULL)
    return((Image *) NULL);
  if (complex_multiplication->next != (Image *) NULL)
    {
      /*
        Normalize the cross-power spectrum by the product magnitude.
      */
      DisableCompositeClampUnlessSpecified(complex_multiplication);
      DisableCompositeClampUnlessSpecified(complex_multiplication->next);
      (void) CompositeImage(complex_multiplication,magnitude_image,
        DivideSrcCompositeOp,MagickTrue,0,0,exception);
      (void) CompositeImage(complex_multiplication->next,magnitude_image,
        DivideSrcCompositeOp,MagickTrue,0,0,exception);
    }
  /*
    Do the IFT and return the phase-correlation result.
  */
  (void) SetImageArtifact(complex_multiplication,"fourier:normalize","inverse");
  cross_correlation=InverseFourierTransformImage(complex_multiplication,
    complex_multiplication->next,MagickFalse,exception);
  complex_multiplication=DestroyImageList(complex_multiplication);
  return(cross_correlation);
}

static MagickBooleanType SIMSetImageMean(Image *image,
  const ChannelStatistics *channel_statistics,ExceptionInfo *exception)
{
  CacheView
    *image_view;

  MagickBooleanType
    status = MagickTrue;

  ssize_t
    y;

  /*
    Set image mean.
  */
  image_view=AcquireAuthenticCacheView(image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(status) \
    magick_number_threads(image,image,image->rows,1)
#endif
  for (y=0; y < (ssize_t) image->rows; y++)
  {
    Quantum
      *magick_restrict q;

    ssize_t
      x;

    if (status == MagickFalse)
      continue;
    q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
    if (q == (Quantum *) NULL)
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (ssize_t) image->columns; x++)
    {
      ssize_t
        i;

      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
      {
        PixelChannel channel = GetPixelChannelChannel(image,i);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        if ((traits & UpdatePixelTrait) == 0)
          continue;
        q[i]=(Quantum) channel_statistics[channel].mean;
      }
      q+=(ptrdiff_t) GetPixelChannels(image);
    }
    if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
      status=MagickFalse;
  }
  image_view=DestroyCacheView(image_view);
  return(status);
}

static Image *SIMSubtractImageMean(const Image *alpha_image,
  const Image *beta_image,const ChannelStatistics *channel_statistics,
  ExceptionInfo *exception)
{
  CacheView
    *beta_view,
    *image_view;

  Image
    *subtract_image;

  MagickBooleanType
    status = MagickTrue;

  ssize_t
    y;

  /*
    Subtract the image mean and pad.
  */
  subtract_image=CloneImage(beta_image,alpha_image->columns,alpha_image->rows,
    MagickTrue,exception);
  if (subtract_image == (Image *) NULL)
    return(subtract_image);
  image_view=AcquireAuthenticCacheView(subtract_image,exception);
  beta_view=AcquireVirtualCacheView(beta_image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(status) \
    magick_number_threads(beta_image,subtract_image,subtract_image->rows,1)
#endif
  for (y=0; y < (ssize_t) subtract_image->rows; y++)
  {
    const Quantum
      *magick_restrict p;

    Quantum
      *magick_restrict q;

    ssize_t
      x;

    if (status == MagickFalse)
      continue;
    p=GetCacheViewVirtualPixels(beta_view,0,y,beta_image->columns,1,exception);
    q=GetCacheViewAuthenticPixels(image_view,0,y,subtract_image->columns,1,
      exception);
    if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (ssize_t) subtract_image->columns; x++)
    {
      ssize_t
        i;

      for (i=0; i < (ssize_t) GetPixelChannels(subtract_image); i++)
      {
        PixelChannel channel = GetPixelChannelChannel(subtract_image,i);
        PixelTrait traits = GetPixelChannelTraits(subtract_image,channel);
        PixelTrait beta_traits = GetPixelChannelTraits(beta_image,channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((beta_traits & UpdatePixelTrait) == 0))
          continue;
        if ((x >= (ssize_t) beta_image->columns) ||
            (y >= (ssize_t) beta_image->rows))
          q[i]=(Quantum) 0;
        else
          q[i]=(Quantum) ((double) GetPixelChannel(beta_image,channel,p)-
            channel_statistics[channel].mean);
      }
      p+=(ptrdiff_t) GetPixelChannels(beta_image);
      q+=(ptrdiff_t) GetPixelChannels(subtract_image);
    }
    if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
      status=MagickFalse;
  }
  beta_view=DestroyCacheView(beta_view);
  image_view=DestroyCacheView(image_view);
  if (status == MagickFalse)
    subtract_image=DestroyImage(subtract_image);
  return(subtract_image);
}

static Image *SIMUnityImage(const Image *alpha_image,const Image *beta_image,
  ExceptionInfo *exception)
{
  CacheView
    *image_view;

  Image
    *unity_image;

  MagickBooleanType
    status = MagickTrue;

  ssize_t
    y;

  /*
    Create a padded unity image.
  */
  unity_image=CloneImage(alpha_image,alpha_image->columns,alpha_image->rows,
    MagickTrue,exception);
  if (unity_image == (Image *) NULL)
    return(unity_image);
  if (SetImageStorageClass(unity_image,DirectClass,exception) == MagickFalse)
    return(DestroyImage(unity_image));
  image_view=AcquireAuthenticCacheView(unity_image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(status) \
    magick_number_threads(unity_image,unity_image,unity_image->rows,1)
#endif
  for (y=0; y < (ssize_t) unity_image->rows; y++)
  {
    Quantum
      *magick_restrict q;

    ssize_t
      x;

    if (status == MagickFalse)
      continue;
    q=GetCacheViewAuthenticPixels(image_view,0,y,unity_image->columns,1,
      exception);
    if (q == (Quantum *) NULL)
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (ssize_t) unity_image->columns; x++)
    {
      ssize_t
        i;

      for (i=0; i < (ssize_t) GetPixelChannels(unity_image); i++)
      {
        PixelChannel channel = GetPixelChannelChannel(unity_image,i);
        PixelTrait traits = GetPixelChannelTraits(unity_image,channel);
        if ((traits & UpdatePixelTrait) == 0)
          continue;
        if ((x >= (ssize_t) beta_image->columns) ||
            (y >= (ssize_t) beta_image->rows))
          q[i]=(Quantum) 0;
        else
          q[i]=QuantumRange;
      }
      q+=(ptrdiff_t) GetPixelChannels(unity_image);
    }
    if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
      status=MagickFalse;
  }
  image_view=DestroyCacheView(image_view);
  if (status == MagickFalse)
    unity_image=DestroyImage(unity_image);
  return(unity_image);
}

static Image *SIMVarianceImage(Image *alpha_image,const Image *beta_image,
  ExceptionInfo *exception)
{
  CacheView
    *beta_view,
    *image_view;

  Image
    *variance_image;

  MagickBooleanType
    status = MagickTrue;

  ssize_t
    y;

  /*
    Compute the variance of the two images.
  */
  variance_image=CloneImage(alpha_image,0,0,MagickTrue,exception);
  if (variance_image == (Image *) NULL)
    return(variance_image);
  image_view=AcquireAuthenticCacheView(variance_image,exception);
  beta_view=AcquireVirtualCacheView(beta_image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(status) \
    magick_number_threads(beta_image,variance_image,variance_image->rows,1)
#endif
  for (y=0; y < (ssize_t) variance_image->rows; y++)
  {
    const Quantum
      *magick_restrict p;

    Quantum
      *magick_restrict q;

    ssize_t
      x;

    if (status == MagickFalse)
      continue;
    p=GetCacheViewVirtualPixels(beta_view,0,y,beta_image->columns,1,
      exception);
    q=GetCacheViewAuthenticPixels(image_view,0,y,variance_image->columns,1,
      exception);
    if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (ssize_t) variance_image->columns; x++)
    {
      ssize_t
        i;

      for (i=0; i < (ssize_t) GetPixelChannels(variance_image); i++)
      {
        double
          error;

        PixelChannel channel = GetPixelChannelChannel(variance_image,i);
        PixelTrait traits = GetPixelChannelTraits(variance_image,channel);
        PixelTrait beta_traits = GetPixelChannelTraits(beta_image,channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((beta_traits & UpdatePixelTrait) == 0))
          continue;
        error=(double) q[i]-(double) GetPixelChannel(beta_image,channel,p);
        q[i]=(Quantum) ((double) ClampToQuantum((double) QuantumRange*
          (sqrt(fabs(QuantumScale*error))/sqrt((double) QuantumRange))));
      }
      p+=(ptrdiff_t) GetPixelChannels(beta_image);
      q+=(ptrdiff_t) GetPixelChannels(variance_image);
    }
    if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
      status=MagickFalse;
  }
  beta_view=DestroyCacheView(beta_view);
  image_view=DestroyCacheView(image_view);
  if (status == MagickFalse)
    variance_image=DestroyImage(variance_image);
  return(variance_image);
}

static Image *DPCSimilarityImage(const Image *image,const Image *reconstruct,
  RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
{
#define ThrowDPCSimilarityException() \
{ \
  if (dot_product_image != (Image *) NULL) \
    dot_product_image=DestroyImage(dot_product_image); \
  if (magnitude_image != (Image *) NULL) \
    magnitude_image=DestroyImage(magnitude_image); \
  if (reconstruct_image != (Image *) NULL) \
    reconstruct_image=DestroyImage(reconstruct_image); \
  if (rx_image != (Image *) NULL) \
    rx_image=DestroyImage(rx_image); \
  if (ry_image != (Image *) NULL) \
    ry_image=DestroyImage(ry_image); \
  if (target_image != (Image *) NULL) \
    target_image=DestroyImage(target_image); \
  if (threshold_image != (Image *) NULL) \
    threshold_image=DestroyImage(threshold_image); \
  if (trx_image != (Image *) NULL) \
    trx_image=DestroyImage(trx_image); \
  if (try_image != (Image *) NULL) \
    try_image=DestroyImage(try_image); \
  if (tx_image != (Image *) NULL) \
    tx_image=DestroyImage(tx_image); \
  if (ty_image != (Image *) NULL) \
    ty_image=DestroyImage(ty_image); \
  return((Image *) NULL); \
}

  double
    edge_factor = 0.0,
    maxima = 0.0,
    mean = 0.0,
    standard_deviation = 0.0;

  Image
    *dot_product_image = (Image *) NULL,
    *magnitude_image = (Image *) NULL,
    *reconstruct_image = (Image *) NULL,
    *rx_image = (Image *) NULL,
    *ry_image = (Image *) NULL,
    *trx_image = (Image *) NULL,
    *target_image = (Image *) NULL,
    *threshold_image = (Image *) NULL,
    *try_image = (Image *) NULL,
    *tx_image = (Image *) NULL,
    *ty_image = (Image *) NULL;

  MagickBooleanType
    status = MagickTrue;

  RectangleInfo
    geometry;

  /*
    Dot product correlation-based image similarity using FFT local statistics.
  */
  target_image=CloneImage(image,0,0,MagickTrue,exception);
  if (target_image == (Image *) NULL)
    return((Image *) NULL);
  /*
    Compute the cross correlation of the test and reconstruct magnitudes.
  */
  reconstruct_image=CloneImage(reconstruct,0,0,MagickTrue,exception);
  if (reconstruct_image == (Image *) NULL)
    ThrowDPCSimilarityException();
  /*
    Compute X and Y derivatives of reference image.
  */
  (void) SetImageVirtualPixelMethod(reconstruct_image,EdgeVirtualPixelMethod,
    exception);
  rx_image=SIMDerivativeImage(reconstruct_image,"Sobel",exception);
  if (rx_image == (Image *) NULL)
    ThrowDPCSimilarityException();
  ry_image=SIMDerivativeImage(reconstruct_image,"Sobel:90",exception);
  reconstruct_image=DestroyImage(reconstruct_image);
  if (ry_image == (Image *) NULL)
    ThrowDPCSimilarityException();
  /*
    Compute magnitude of derivatives.
  */
  magnitude_image=SIMMagnitudeImage(rx_image,ry_image,exception);
  if (magnitude_image == (Image *) NULL)
    ThrowDPCSimilarityException();
  /*
    Compute an edge normalization correction.
  */
  threshold_image=CloneImage(magnitude_image,0,0,MagickTrue,exception);
  if (threshold_image == (Image *) NULL)
    ThrowDPCSimilarityException();
  status=BilevelImage(threshold_image,0.0,exception);
  if (status == MagickFalse)
    ThrowDPCSimilarityException();
  status=GetImageMean(threshold_image,&mean,&standard_deviation,exception);
  threshold_image=DestroyImage(threshold_image);
  if (status == MagickFalse)
    ThrowDPCSimilarityException();
  edge_factor=MagickSafeReciprocal(QuantumScale*mean*(double)
    reconstruct->columns*(double) reconstruct->rows)+QuantumScale;
  /*
    Divide X and Y derivitives of reference image by magnitude.
  */
  trx_image=SIMDivideByMagnitude(rx_image,magnitude_image,image,exception);
  rx_image=DestroyImage(rx_image);
  if (trx_image == (Image *) NULL)
    ThrowDPCSimilarityException();
  rx_image=trx_image;
  trx_image=(Image *) NULL;
  try_image=SIMDivideByMagnitude(ry_image,magnitude_image,image,exception);
  magnitude_image=DestroyImage(magnitude_image);
  ry_image=DestroyImage(ry_image);
  if (try_image == (Image *) NULL)
    ThrowDPCSimilarityException();
  ry_image=try_image;
  try_image=(Image *) NULL;
  /*
    Compute X and Y derivatives of image.
  */
  (void) SetImageVirtualPixelMethod(target_image,EdgeVirtualPixelMethod,
    exception);
  tx_image=SIMDerivativeImage(target_image,"Sobel",exception);
  if (tx_image == (Image *) NULL)
    ThrowDPCSimilarityException();
  ty_image=SIMDerivativeImage(target_image,"Sobel:90",exception);
  target_image=DestroyImage(target_image);
  if (ty_image == (Image *) NULL)
    ThrowDPCSimilarityException();
  /*
    Compute magnitude of derivatives.
  */
  magnitude_image=SIMMagnitudeImage(tx_image,ty_image,exception);
  if (magnitude_image == (Image *) NULL)
    ThrowDPCSimilarityException();
  /*
    Divide Lx and Ly by magnitude.
  */
  trx_image=SIMDivideByMagnitude(tx_image,magnitude_image,image,exception);
  tx_image=DestroyImage(tx_image);
  if (trx_image == (Image *) NULL)
    ThrowDPCSimilarityException();
  tx_image=trx_image;
  trx_image=(Image *) NULL;
  try_image=SIMDivideByMagnitude(ty_image,magnitude_image,image,exception);
  ty_image=DestroyImage(ty_image);
  magnitude_image=DestroyImage(magnitude_image);
  if (try_image == (Image *) NULL)
    ThrowDPCSimilarityException();
  ty_image=try_image;
  try_image=(Image *) NULL;
  /*
    Compute the cross correlation of the test and reference images.
  */
  trx_image=SIMCrossCorrelationImage(tx_image,rx_image,exception);
  rx_image=DestroyImage(rx_image);
  tx_image=DestroyImage(tx_image);
  if (trx_image == (Image *) NULL)
    ThrowDPCSimilarityException();
  try_image=SIMCrossCorrelationImage(ty_image,ry_image,exception);
  ry_image=DestroyImage(ry_image);
  ty_image=DestroyImage(ty_image);
  if (try_image == (Image *) NULL)
    ThrowDPCSimilarityException();
  /*
    Evaluate dot product correlation image.
  */
  (void) SetImageArtifact(try_image,"compose:clamp","false");
  status=CompositeImage(trx_image,try_image,PlusCompositeOp,MagickTrue,0,0,
    exception);
  try_image=DestroyImage(try_image);
  if (status == MagickFalse)
    ThrowDPCSimilarityException();
  status=SIMMultiplyImage(trx_image,edge_factor,
    (const ChannelStatistics *) NULL,exception);
  if (status == MagickFalse)
    ThrowDPCSimilarityException();
  /*
    Crop results.
  */
  SetGeometry(image,&geometry);
  geometry.width=image->columns;
  geometry.height=image->rows;
  (void) ResetImagePage(trx_image,"0x0+0+0");
  dot_product_image=CropImage(trx_image,&geometry,exception);
  trx_image=DestroyImage(trx_image);
  if (dot_product_image == (Image *) NULL)
    ThrowDPCSimilarityException();
  (void) ResetImagePage(dot_product_image,"0x0+0+0");
  /*
    Identify the maxima value in the image and its location.
  */
  status=GrayscaleImage(dot_product_image,AveragePixelIntensityMethod,
    exception);
  if (status == MagickFalse)
    ThrowDPCSimilarityException();
  dot_product_image->depth=32;
  dot_product_image->colorspace=GRAYColorspace;
  dot_product_image->alpha_trait=UndefinedPixelTrait;
  status=SIMFilterImageNaNs(dot_product_image,exception);
  if (status == MagickFalse)
    ThrowDPCSimilarityException();
  status=SIMMaximaImage(dot_product_image,&maxima,offset,exception);
  if (status == MagickFalse)
    ThrowDPCSimilarityException();
  if ((QuantumScale*maxima) > 1.0)
    {
      status=SIMMultiplyImage(dot_product_image,1.0/(QuantumScale*maxima),
        (const ChannelStatistics *) NULL,exception);
      maxima=(double) QuantumRange;
    }
  *similarity_metric=QuantumScale*maxima;
  return(dot_product_image);
}

static Image *MSESimilarityImage(const Image *image,const Image *reconstruct,
  RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
{
#define ThrowMSESimilarityException() \
{ \
  if (alpha_image != (Image *) NULL) \
    alpha_image=DestroyImage(alpha_image); \
  if (beta_image != (Image *) NULL) \
    beta_image=DestroyImage(beta_image); \
  if (channel_statistics != (ChannelStatistics *) NULL) \
    channel_statistics=(ChannelStatistics *) \
      RelinquishMagickMemory(channel_statistics); \
  if (mean_image != (Image *) NULL) \
    mean_image=DestroyImage(mean_image); \
  if (mse_image != (Image *) NULL) \
    mse_image=DestroyImage(mse_image); \
  if (reconstruct_image != (Image *) NULL) \
    reconstruct_image=DestroyImage(reconstruct_image); \
  if (sum_image != (Image *) NULL) \
    sum_image=DestroyImage(sum_image); \
  if (alpha_image != (Image *) NULL) \
    alpha_image=DestroyImage(alpha_image); \
  return((Image *) NULL); \
}

  ChannelStatistics
    *channel_statistics = (ChannelStatistics *) NULL;

  double
    minima = 0.0;

  Image
    *alpha_image = (Image *) NULL,
    *beta_image = (Image *) NULL,
    *mean_image = (Image *) NULL,
    *mse_image = (Image *) NULL,
    *reconstruct_image = (Image *) NULL,
    *sum_image = (Image *) NULL,
    *target_image = (Image *) NULL;

  MagickBooleanType
    status = MagickTrue;

  RectangleInfo
    geometry;

  /*
    MSE correlation-based image similarity using FFT local statistics.
  */
  target_image=SIMSquareImage(image,exception);
  if (target_image == (Image *) NULL)
    ThrowMSESimilarityException();
  reconstruct_image=SIMUnityImage(image,reconstruct,exception);
  if (reconstruct_image == (Image *) NULL)
    ThrowMSESimilarityException();
  /*
    Create (U * test)/# pixels.
  */
  alpha_image=SIMCrossCorrelationImage(target_image,reconstruct_image,
    exception);
  target_image=DestroyImage(target_image);
  if (alpha_image == (Image *) NULL)
    ThrowMSESimilarityException();
  status=SIMMultiplyImage(alpha_image,1.0/(double) reconstruct->columns/(double)
    reconstruct->rows,(const ChannelStatistics *) NULL,exception);
  if (status == MagickFalse)
    ThrowMSESimilarityException();
  /*
    Create 2*(test * reconstruction)# pixels.
  */
  (void) CompositeImage(reconstruct_image,reconstruct,CopyCompositeOp,
    MagickTrue,0,0,exception);
  beta_image=SIMCrossCorrelationImage(image,reconstruct_image,exception);
  if (beta_image == (Image *) NULL)
    {
      reconstruct_image=DestroyImage(reconstruct_image);
      ThrowMSESimilarityException();
    }
  status=SIMMultiplyImage(beta_image,-2.0/(double) reconstruct->columns/(double)
    reconstruct->rows,(const ChannelStatistics *) NULL,exception);
  reconstruct_image=DestroyImage(reconstruct_image);
  if (status == MagickFalse)
    ThrowMSESimilarityException();
  /*
    Mean of reconstruction squared.
  */
  sum_image=SIMSquareImage(reconstruct,exception);
  if (sum_image == (Image *) NULL)
    ThrowMSESimilarityException();
  channel_statistics=GetImageStatistics(sum_image,exception);
  if (channel_statistics == (ChannelStatistics *) NULL)
    ThrowMSESimilarityException();
  status=SetImageStorageClass(sum_image,DirectClass,exception);
  if (status == MagickFalse)
    ThrowMSESimilarityException();
  status=SIMSetImageMean(sum_image,channel_statistics,exception);
  channel_statistics=(ChannelStatistics *)
    RelinquishMagickMemory(channel_statistics);
  if (status == MagickFalse)
    ThrowMSESimilarityException();
  /*
    Create mean image.
  */
  AppendImageToList(&sum_image,alpha_image);
  AppendImageToList(&sum_image,beta_image);
  mean_image=EvaluateImages(sum_image,SumEvaluateOperator,exception);
  if (mean_image == (Image *) NULL)
    ThrowMSESimilarityException();
  sum_image=DestroyImage(sum_image);
  status=GrayscaleImage(mean_image,AveragePixelIntensityMethod,exception);
  if (status == MagickFalse)
    ThrowMSESimilarityException();
  /*
    Crop to difference of reconstruction and test images.
  */
  SetGeometry(image,&geometry);
  geometry.width=image->columns;
  geometry.height=image->rows;
  (void) ResetImagePage(mean_image,"0x0+0+0");
  mse_image=CropImage(mean_image,&geometry,exception);
  mean_image=DestroyImage(mean_image);
  if (mse_image == (Image *) NULL)
    ThrowMSESimilarityException();
  /*
    Identify the minima value in the correlation image and its location.
  */
  (void) ResetImagePage(mse_image,"0x0+0+0");
  (void) ClampImage(mse_image,exception);
  mse_image->depth=32;
  mse_image->colorspace=GRAYColorspace;
  mse_image->alpha_trait=UndefinedPixelTrait;
  status=SIMMinimaImage(mse_image,&minima,offset,exception);
  if (status == MagickFalse)
    ThrowMSESimilarityException();
  status=NegateImage(mse_image,MagickFalse,exception);
  if (status == MagickFalse)
    ThrowMSESimilarityException();
  alpha_image=DestroyImage(alpha_image);
  beta_image=DestroyImage(beta_image);
  if ((QuantumScale*minima) < (double) FLT_EPSILON)
    minima=0.0;
  *similarity_metric=QuantumScale*minima;
  return(mse_image);
}

static Image *NCCSimilarityImage(const Image *image,const Image *reconstruct,
  RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
{
#define ThrowNCCSimilarityException() \
{ \
  if (alpha_image != (Image *) NULL) \
    alpha_image=DestroyImage(alpha_image); \
  if (beta_image != (Image *) NULL) \
    beta_image=DestroyImage(beta_image); \
  if (channel_statistics != (ChannelStatistics *) NULL) \
    channel_statistics=(ChannelStatistics *) \
      RelinquishMagickMemory(channel_statistics); \
  if (correlation_image != (Image *) NULL) \
    correlation_image=DestroyImage(correlation_image); \
  if (divide_image != (Image *) NULL) \
    divide_image=DestroyImage(divide_image); \
  if (ncc_image != (Image *) NULL) \
    ncc_image=DestroyImage(ncc_image); \
  if (normalize_image != (Image *) NULL) \
    normalize_image=DestroyImage(normalize_image); \
  if (reconstruct_image != (Image *) NULL) \
    reconstruct_image=DestroyImage(reconstruct_image); \
  if (target_image != (Image *) NULL) \
    target_image=DestroyImage(target_image); \
  if (variance_image != (Image *) NULL) \
    variance_image=DestroyImage(variance_image); \
  return((Image *) NULL); \
}

  ChannelStatistics
    *channel_statistics = (ChannelStatistics *) NULL;

  double
    maxima = 0.0;

  Image
    *alpha_image = (Image *) NULL,
    *beta_image = (Image *) NULL,
    *correlation_image = (Image *) NULL,
    *divide_image = (Image *) NULL,
    *ncc_image = (Image *) NULL,
    *normalize_image = (Image *) NULL,
    *reconstruct_image = (Image *) NULL,
    *target_image = (Image *) NULL,
    *variance_image = (Image *) NULL;

  MagickBooleanType
    status = MagickTrue;

  RectangleInfo
    geometry;

  /*
    NCC correlation-based image similarity with FFT local statistics.
  */
  target_image=SIMSquareImage(image,exception);
  if (target_image == (Image *) NULL)
    ThrowNCCSimilarityException();
  reconstruct_image=SIMUnityImage(image,reconstruct,exception);
  if (reconstruct_image == (Image *) NULL)
    ThrowNCCSimilarityException();
  /*
    Compute the cross correlation of the test and reconstruction images.
  */
  alpha_image=SIMCrossCorrelationImage(target_image,reconstruct_image,
    exception);
  target_image=DestroyImage(target_image);
  if (alpha_image == (Image *) NULL)
    ThrowNCCSimilarityException();
  status=SIMMultiplyImage(alpha_image,(double) QuantumRange*
    (double) reconstruct->columns*(double) reconstruct->rows,
    (const ChannelStatistics *) NULL,exception);
  if (status == MagickFalse)
    ThrowNCCSimilarityException();
  /*
    Compute the cross correlation of the source and reconstruction images.
  */
  beta_image=SIMCrossCorrelationImage(image,reconstruct_image,exception);
  reconstruct_image=DestroyImage(reconstruct_image);
  if (beta_image == (Image *) NULL)
    ThrowNCCSimilarityException();
  target_image=SIMSquareImage(beta_image,exception);
  beta_image=DestroyImage(beta_image);
  if (target_image == (Image *) NULL)
    ThrowNCCSimilarityException();
  status=SIMMultiplyImage(target_image,(double) QuantumRange,
    (const ChannelStatistics *) NULL,exception);
  if (status == MagickFalse)
    ThrowNCCSimilarityException();
  /*
    Compute the variance of the two images.
  */
  variance_image=SIMVarianceImage(alpha_image,target_image,exception);
  target_image=DestroyImage(target_image);
  alpha_image=DestroyImage(alpha_image);
  if (variance_image == (Image *) NULL)
    ThrowNCCSimilarityException();
  /*
    Subtract the image mean.
  */
  channel_statistics=GetImageStatistics(reconstruct,exception);
  if (channel_statistics == (ChannelStatistics *) NULL)
    ThrowNCCSimilarityException();
  status=SIMMultiplyImage(variance_image,1.0,channel_statistics,exception);
  if (status == MagickFalse)
    ThrowNCCSimilarityException();
  normalize_image=SIMSubtractImageMean(image,reconstruct,channel_statistics,
    exception);
  channel_statistics=(ChannelStatistics *)
    RelinquishMagickMemory(channel_statistics);
  if (normalize_image == (Image *) NULL)
    ThrowNCCSimilarityException();
  correlation_image=SIMCrossCorrelationImage(image,normalize_image,exception);
  normalize_image=DestroyImage(normalize_image);
  if (correlation_image == (Image *) NULL)
    ThrowNCCSimilarityException();
  /*
    Divide the two images.
  */
  divide_image=SIMDivideImage(correlation_image,variance_image,exception);
  correlation_image=DestroyImage(correlation_image);
  variance_image=DestroyImage(variance_image);
  if (divide_image == (Image *) NULL)
    ThrowNCCSimilarityException();
  /*
    Crop padding.
  */
  SetGeometry(image,&geometry);
  geometry.width=image->columns;
  geometry.height=image->rows;
  (void) ResetImagePage(divide_image,"0x0+0+0");
  ncc_image=CropImage(divide_image,&geometry,exception);
  divide_image=DestroyImage(divide_image);
  if (ncc_image == (Image *) NULL)
    ThrowNCCSimilarityException();
  /*
    Identify the maxima value in the image and its location.
  */
  (void) ResetImagePage(ncc_image,"0x0+0+0");
  status=GrayscaleImage(ncc_image,AveragePixelIntensityMethod,exception);
  if (status == MagickFalse)
    ThrowNCCSimilarityException();
  ncc_image->depth=32;
  ncc_image->colorspace=GRAYColorspace;
  ncc_image->alpha_trait=UndefinedPixelTrait;
  status=SIMMaximaImage(ncc_image,&maxima,offset,exception);
  if (status == MagickFalse)
    ThrowNCCSimilarityException();
  if ((QuantumScale*maxima) > 1.0)
    {
      status=SIMMultiplyImage(ncc_image,1.0/(QuantumScale*maxima),
        (const ChannelStatistics *) NULL,exception);
      maxima=(double) QuantumRange;
    }
  *similarity_metric=QuantumScale*maxima;
  return(ncc_image);
}

static Image *PhaseSimilarityImage(const Image *image,const Image *reconstruct,
  RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
{
#define ThrowPhaseSimilarityException() \
{ \
  if (phase_image != (Image *) NULL) \
    phase_image=DestroyImage(phase_image); \
  if (gamma_image != (Image *) NULL) \
    gamma_image=DestroyImage(gamma_image); \
  if (test_magnitude != (Image *) NULL) \
    test_magnitude=DestroyImage(test_magnitude); \
  if (magnitude_image != (Image *) NULL) \
    magnitude_image=DestroyImage(magnitude_image); \
  if (reconstruct_magnitude != (Image *) NULL) \
    reconstruct_magnitude=DestroyImage(reconstruct_magnitude); \
  if (correlation_image != (Image *) NULL) \
    correlation_image=DestroyImage(correlation_image); \
  if (fft_images != (Image *) NULL) \
    fft_images=DestroyImageList(fft_images); \
  if (reconstruct_image != (Image *) NULL) \
    reconstruct_image=DestroyImage(reconstruct_image); \
  if (target_image != (Image *) NULL) \
    target_image=DestroyImage(target_image); \
  return((Image *) NULL); \
}

  double
    maxima = 0.0;

  Image
    *correlation_image = (Image *) NULL,
    *fft_images = (Image *) NULL,
    *gamma_image = (Image *) NULL,
    *magnitude_image = (Image *) NULL,
    *phase_image = (Image *) NULL,
    *reconstruct_image = (Image *) NULL,
    *reconstruct_magnitude = (Image *) NULL,
    *target_image = (Image *) NULL,
    *test_magnitude = (Image *) NULL;

  MagickBooleanType
    status = MagickTrue;

  RectangleInfo
    geometry;

  /*
    Phase correlation-based image similarity using FFT local statistics.
  */
  target_image=CloneImage(image,0,0,MagickTrue,exception);
  if (target_image == (Image *) NULL)
    ThrowPhaseSimilarityException();
  (void) ResetImagePage(target_image,"0x0+0+0");
  GetPixelInfoRGBA((Quantum) 0,(Quantum) 0,(Quantum) 0,(Quantum) 0,
    &target_image->background_color);
  status=SetImageExtent(target_image,2*CastDoubleToSizeT(ceil((double)
    image->columns/2.0)),2*CastDoubleToSizeT(ceil((double) image->rows/2.0)),
    exception);
  if (status == MagickFalse)
    ThrowPhaseSimilarityException();
  /*
    Compute the cross correlation of the target and reconstruct magnitudes.
  */
  reconstruct_image=CloneImage(reconstruct,0,0,MagickTrue,exception);
  if (reconstruct_image == (Image *) NULL)
    ThrowPhaseSimilarityException();
  (void) ResetImagePage(reconstruct_image,"0x0+0+0");
  GetPixelInfoRGBA((Quantum) 0,(Quantum) 0,(Quantum) 0,(Quantum) 0,
    &reconstruct_image->background_color);
  status=SetImageExtent(reconstruct_image,2*CastDoubleToSizeT(ceil((double)
    image->columns/2.0)),2*CastDoubleToSizeT(ceil((double) image->rows/2.0)),
    exception);
  if (status == MagickFalse)
    ThrowPhaseSimilarityException();
  /*
    Evaluate phase coorelation image and divide by the product magnitude.
  */
  (void) SetImageArtifact(target_image,"fourier:normalize","inverse");
  fft_images=ForwardFourierTransformImage(target_image,MagickTrue,exception);
  if (fft_images == (Image *) NULL)
    ThrowPhaseSimilarityException();
  test_magnitude=CloneImage(fft_images,0,0,MagickTrue,exception);
  fft_images=DestroyImageList(fft_images);
  if (test_magnitude == (Image *) NULL)
    ThrowPhaseSimilarityException();
  (void) SetImageArtifact(reconstruct_image,"fourier:normalize","inverse");
  fft_images=ForwardFourierTransformImage(reconstruct_image,MagickTrue,
    exception);
  if (fft_images == (Image *) NULL)
    ThrowPhaseSimilarityException();
  reconstruct_magnitude=CloneImage(fft_images,0,0,MagickTrue,exception);
  fft_images=DestroyImageList(fft_images);
  if (reconstruct_magnitude == (Image *) NULL)
    ThrowPhaseSimilarityException();
  magnitude_image=CloneImage(reconstruct_magnitude,0,0,MagickTrue,exception);
  if (magnitude_image == (Image *) NULL)
    ThrowPhaseSimilarityException();
  DisableCompositeClampUnlessSpecified(magnitude_image);
  (void) CompositeImage(magnitude_image,test_magnitude,MultiplyCompositeOp,
    MagickTrue,0,0,exception);
  /*
    Compute the cross correlation of the target and reconstruction images.
  */
  correlation_image=SIMPhaseCorrelationImage(target_image,reconstruct_image,
    magnitude_image,exception);
  target_image=DestroyImage(target_image);
  reconstruct_image=DestroyImage(reconstruct_image);
  test_magnitude=DestroyImage(test_magnitude);
  reconstruct_magnitude=DestroyImage(reconstruct_magnitude);
  if (correlation_image == (Image *) NULL)
    ThrowPhaseSimilarityException();
  gamma_image=CloneImage(correlation_image,0,0,MagickTrue,exception);
  correlation_image=DestroyImage(correlation_image);
  if (gamma_image == (Image *) NULL)
    ThrowPhaseSimilarityException();
  /*
    Crop padding.
  */
  SetGeometry(image,&geometry);
  geometry.width=image->columns;
  geometry.height=image->rows;
  (void) ResetImagePage(gamma_image,"0x0+0+0");
  phase_image=CropImage(gamma_image,&geometry,exception);
  gamma_image=DestroyImage(gamma_image);
  if (phase_image == (Image *) NULL)
    ThrowPhaseSimilarityException();
  (void) ResetImagePage(phase_image,"0x0+0+0");
  /*
    Identify the maxima value in the correlation image and its location.
  */
  status=GrayscaleImage(phase_image,AveragePixelIntensityMethod,exception);
  if (status == MagickFalse)
    ThrowPhaseSimilarityException();
  phase_image->depth=32;
  phase_image->colorspace=GRAYColorspace;
  phase_image->alpha_trait=UndefinedPixelTrait;
  status=SIMFilterImageNaNs(phase_image,exception);
  if (status == MagickFalse)
    ThrowPhaseSimilarityException();
  status=SIMMaximaImage(phase_image,&maxima,offset,exception);
  if (status == MagickFalse)
    ThrowPhaseSimilarityException();
  magnitude_image=DestroyImage(magnitude_image);
  *similarity_metric=QuantumScale*maxima;
  return(phase_image);
}

static Image *PSNRSimilarityImage(const Image *image,const Image *reconstruct,
  RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
{
  Image
    *psnr_image = (Image *) NULL;

  psnr_image=MSESimilarityImage(image,reconstruct,offset,similarity_metric,
    exception);
  if (psnr_image == (Image *) NULL)
    return(psnr_image);
  *similarity_metric=10.0*MagickSafeLog10(MagickSafeReciprocal(
    *similarity_metric))/MagickSafePSNRRecipicol(10.0);
  return(psnr_image);
}

static Image *RMSESimilarityImage(const Image *image,const Image *reconstruct,
  RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
{
  Image
    *rmse_image = (Image *) NULL;

  rmse_image=MSESimilarityImage(image,reconstruct,offset,similarity_metric,
    exception);
  if (rmse_image == (Image *) NULL)
    return(rmse_image);
  *similarity_metric=sqrt(*similarity_metric);
  return(rmse_image);
}
#endif

static double GetSimilarityMetric(const Image *image,
  const Image *reconstruct_image,const MetricType metric,
  const ssize_t x_offset,const ssize_t y_offset,ExceptionInfo *exception)
{
  double
    *channel_similarity,
    similarity = 0.0;

  ExceptionInfo
    *sans_exception = AcquireExceptionInfo();

  Image
    *similarity_image;

  MagickBooleanType
    status = MagickTrue;

  RectangleInfo
    geometry;

  size_t
    length = MaxPixelChannels+1UL;

  SetGeometry(reconstruct_image,&geometry);
  geometry.x=x_offset;
  geometry.y=y_offset;
  similarity_image=CropImage(image,&geometry,sans_exception);
  sans_exception=DestroyExceptionInfo(sans_exception);
  if (similarity_image == (Image *) NULL)
    return(NAN);
  /*
    Get image distortion.
  */
  channel_similarity=(double *) AcquireQuantumMemory(length,
    sizeof(*channel_similarity));
  if (channel_similarity == (double *) NULL)
    return(NAN);
  (void) memset(channel_similarity,0,length*sizeof(*channel_similarity));
  switch (metric)
  {
    case AbsoluteErrorMetric:
    {
      status=GetAESimilarity(similarity_image,reconstruct_image,
        channel_similarity,exception);
      break;
    }
    case DotProductCorrelationErrorMetric:
    {
      status=GetDPCSimilarity(similarity_image,reconstruct_image,
        channel_similarity,exception);
      break;
    }
    case FuzzErrorMetric:
    {
      status=GetFUZZSimilarity(similarity_image,reconstruct_image,
        channel_similarity,exception);
      break;
    }
    case MeanAbsoluteErrorMetric:
    {
      status=GetMAESimilarity(similarity_image,reconstruct_image,
        channel_similarity,exception);
      break;
    }
    case MeanErrorPerPixelErrorMetric:
    {
      status=GetMEPPSimilarity(similarity_image,reconstruct_image,
        channel_similarity,exception);
      break;
    }
    case MeanSquaredErrorMetric:
    {
      status=GetMSESimilarity(similarity_image,reconstruct_image,
        channel_similarity,exception);
      break;
    }
    case NormalizedCrossCorrelationErrorMetric:
    {
      status=GetNCCSimilarity(similarity_image,reconstruct_image,
        channel_similarity,exception);
      break;
    }
    case PeakAbsoluteErrorMetric:
    {
      status=GetPASimilarity(similarity_image,reconstruct_image,
        channel_similarity,exception);
      break;
    }
    case PeakSignalToNoiseRatioErrorMetric:
    {
      status=GetPSNRSimilarity(similarity_image,reconstruct_image,
        channel_similarity,exception);
      break;
    }
    case PerceptualHashErrorMetric:
    {
      status=GetPHASHSimilarity(similarity_image,reconstruct_image,
        channel_similarity,exception);
      break;
    }
    case PhaseCorrelationErrorMetric:
    {
      status=GetPHASESimilarity(similarity_image,reconstruct_image,
        channel_similarity,exception);
      break;
    }
    case PixelDifferenceCountErrorMetric:
    {
      status=GetPDCSimilarity(similarity_image,reconstruct_image,
        channel_similarity,exception);
      break;
    }
    case RootMeanSquaredErrorMetric:
    case UndefinedErrorMetric:
    default:
    {
      status=GetRMSESimilarity(similarity_image,reconstruct_image,
        channel_similarity,exception);
      break;
    }
    case StructuralDissimilarityErrorMetric:
    {
      status=GetDSSIMSimilarity(similarity_image,reconstruct_image,
        channel_similarity,exception);
      break;
    }
    case StructuralSimilarityErrorMetric:
    {
      status=GetSSIMSimularity(similarity_image,reconstruct_image,
        channel_similarity,exception);
      break;
    }
  }
  similarity_image=DestroyImage(similarity_image);
  similarity=channel_similarity[CompositePixelChannel];
  channel_similarity=(double *) RelinquishMagickMemory(channel_similarity);
  if (status == MagickFalse)
    return(NAN);
  return(similarity);
}

MagickExport Image *SimilarityImage(const Image *image,const Image *reconstruct,
  const MetricType metric,const double similarity_threshold,
  RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
{
#define SimilarityImageTag  "Similarity/Image"

  typedef struct
  {
    double
      similarity;

    ssize_t
      x,
      y;
  } SimilarityInfo;

  CacheView
    *similarity_view;

  Image
    *similarity_image = (Image *) NULL;

  MagickBooleanType
    status = MagickTrue,
    threshold_trigger = MagickFalse;

  MagickOffsetType
    progress = 0;

  SimilarityInfo
    similarity_info = { 0.0, 0, 0 };

  size_t
    columns,
    rows;

  ssize_t
    y;

  assert(image != (const Image *) NULL);
  assert(image->signature == MagickCoreSignature);
  assert(exception != (ExceptionInfo *) NULL);
  assert(exception->signature == MagickCoreSignature);
  assert(offset != (RectangleInfo *) NULL);
  if (IsEventLogging() != MagickFalse)
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
  SetGeometry(reconstruct,offset);
  *similarity_metric=0.0;
  offset->x=0;
  offset->y=0;
#if defined(MAGICKCORE_HDRI_SUPPORT) && defined(MAGICKCORE_FFTW_DELEGATE)
{
  const char *artifact = GetImageArtifact(image,"compare:frequency-domain");
  if (artifact == (const char *) NULL)
    artifact=GetImageArtifact(image,"compare:accelerate-ncc");
  if (((artifact == (const char *) NULL) ||
       (IsStringTrue(artifact) != MagickFalse)) &&
      ((image->channels & ReadMaskChannel) == 0))
    switch (metric)
    {
      case DotProductCorrelationErrorMetric:
      {
        similarity_image=DPCSimilarityImage(image,reconstruct,offset,
          similarity_metric,exception);
        return(similarity_image);
      }
      case MeanSquaredErrorMetric:
      {
        similarity_image=MSESimilarityImage(image,reconstruct,offset,
          similarity_metric,exception);
        return(similarity_image);
      }
      case NormalizedCrossCorrelationErrorMetric:
      {
        similarity_image=NCCSimilarityImage(image,reconstruct,offset,
          similarity_metric,exception);
        return(similarity_image);
      }
      case PeakSignalToNoiseRatioErrorMetric:
      {
        similarity_image=PSNRSimilarityImage(image,reconstruct,offset,
          similarity_metric,exception);
        return(similarity_image);
      }
      case PhaseCorrelationErrorMetric:
      {
        similarity_image=PhaseSimilarityImage(image,reconstruct,offset,
          similarity_metric,exception);
        return(similarity_image);
      }
      case RootMeanSquaredErrorMetric:
      case UndefinedErrorMetric:
      {
        similarity_image=RMSESimilarityImage(image,reconstruct,offset,
          similarity_metric,exception);
        return(similarity_image);
      }
      default:
        break;
    }
}
#endif
  if ((image->columns < reconstruct->columns) ||
      (image->rows < reconstruct->rows))
    {
      (void) ThrowMagickException(exception,GetMagickModule(),OptionWarning,
        "GeometryDoesNotContainImage","`%s'",image->filename);
      return((Image *) NULL);
    }
  SetImageCompareBounds(image,reconstruct,&columns,&rows);
  similarity_image=CloneImage(image,columns,rows,MagickTrue,exception);
  if (similarity_image == (Image *) NULL)
    return((Image *) NULL);
  similarity_image->depth=32;
  similarity_image->colorspace=GRAYColorspace;
  similarity_image->alpha_trait=UndefinedPixelTrait;
  status=SetImageStorageClass(similarity_image,DirectClass,exception);
  if (status == MagickFalse)
    return(DestroyImage(similarity_image));
  /*
    Measure similarity of reconstruction image against image.
  */
  similarity_info.similarity=GetSimilarityMetric(image,reconstruct,metric,
    similarity_info.x,similarity_info.y,exception);
  similarity_view=AcquireAuthenticCacheView(similarity_image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(similarity_info,status,threshold_trigger) \
    magick_number_threads(image,reconstruct,similarity_image->rows,1)
#endif
  for (y=0; y < (ssize_t) similarity_image->rows; y++)
  {
    double
      similarity;

    Quantum
      *magick_restrict q;

    SimilarityInfo
      channel_info = similarity_info;

    ssize_t
      x;

    if (status == MagickFalse)
      continue;
    if (threshold_trigger != MagickFalse)
      continue;
    q=QueueCacheViewAuthenticPixels(similarity_view,0,y,
      similarity_image->columns,1,exception);
    if (q == (Quantum *) NULL)
      {
        status=MagickFalse;
        continue;
      }
    for (x=0; x < (ssize_t) similarity_image->columns; x++)
    {
      ssize_t
        i;

      similarity=GetSimilarityMetric((Image *) image,reconstruct,metric,x,y,
        exception);
      switch (metric)
      {
        case DotProductCorrelationErrorMetric:
        case NormalizedCrossCorrelationErrorMetric:
        case PeakSignalToNoiseRatioErrorMetric:
        case PhaseCorrelationErrorMetric:
        case StructuralSimilarityErrorMetric:
        {
          if (similarity <= channel_info.similarity)
            break;
          channel_info.similarity=similarity;
          channel_info.x=x;
          channel_info.y=y;
          break;
        }
        default:
        {
          if (similarity >= channel_info.similarity)
            break;
          channel_info.similarity=similarity;
          channel_info.x=x;
          channel_info.y=y;
          break;
        }
      }
      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
      {
        PixelChannel channel = GetPixelChannelChannel(image,i);
        PixelTrait traits = GetPixelChannelTraits(image,channel);
        PixelTrait similarity_traits = GetPixelChannelTraits(similarity_image,
          channel);
        if (((traits & UpdatePixelTrait) == 0) ||
            ((similarity_traits & UpdatePixelTrait) == 0))
          continue;
        switch (metric)
        {
          case DotProductCorrelationErrorMetric:
          case NormalizedCrossCorrelationErrorMetric:
          case PeakSignalToNoiseRatioErrorMetric:
          case PhaseCorrelationErrorMetric:
          case StructuralSimilarityErrorMetric:
          {
            SetPixelChannel(similarity_image,channel,ClampToQuantum((double)
              QuantumRange*similarity),q);
            break;
          }
          default:
          {
            SetPixelChannel(similarity_image,channel,ClampToQuantum((double)
              QuantumRange*(1.0-similarity)),q);
            break;
          }
        }
      }
      q+=(ptrdiff_t) GetPixelChannels(similarity_image);
    }
#if defined(MAGICKCORE_OPENMP_SUPPORT)
    #pragma omp critical (MagickCore_GetSimilarityMetric)
#endif
    switch (metric)
    {
      case DotProductCorrelationErrorMetric:
      case NormalizedCrossCorrelationErrorMetric:
      case PeakSignalToNoiseRatioErrorMetric:
      case PhaseCorrelationErrorMetric:
      case StructuralSimilarityErrorMetric:
      {
        if (similarity_threshold != DefaultSimilarityThreshold)
          if (channel_info.similarity >= similarity_threshold)
            threshold_trigger=MagickTrue;
        if (channel_info.similarity >= similarity_info.similarity)
          similarity_info=channel_info;
        break;
      }
      default:
      {
        if (similarity_threshold != DefaultSimilarityThreshold)
          if (channel_info.similarity < similarity_threshold)
            threshold_trigger=MagickTrue;
        if (channel_info.similarity < similarity_info.similarity)
          similarity_info=channel_info;
        break;
      }
    }
    if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
      status=MagickFalse;
    if (image->progress_monitor != (MagickProgressMonitor) NULL)
      {
        MagickBooleanType
          proceed;

        progress++;
        proceed=SetImageProgress(image,SimilarityImageTag,progress,image->rows);
        if (proceed == MagickFalse)
          status=MagickFalse;
      }
  }
  similarity_view=DestroyCacheView(similarity_view);
  if (status == MagickFalse)
    similarity_image=DestroyImage(similarity_image);
  *similarity_metric=similarity_info.similarity;
  if (fabs(*similarity_metric) < MagickEpsilon)
    *similarity_metric=0.0;
  offset->x=similarity_info.x;
  offset->y=similarity_info.y;
  (void) FormatImageProperty((Image *) image,"similarity","%.*g",
    GetMagickPrecision(),*similarity_metric);
  (void) FormatImageProperty((Image *) image,"similarity.offset.x","%.*g",
    GetMagickPrecision(),(double) offset->x);
  (void) FormatImageProperty((Image *) image,"similarity.offset.y","%.*g",
    GetMagickPrecision(),(double) offset->y);
  return(similarity_image);
}

Coverage Report

Created: 2026-11-12 07:35