% December 2003 %
% %
% %
-% Copyright 1999-2009 ImageMagick Studio LLC, a non-profit organization %
+% Copyright 1999-2011 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 %
#include "magick/pixel-private.h"
#include "magick/resource_.h"
#include "magick/string_.h"
+#include "magick/statistic.h"
+#include "magick/transform.h"
#include "magick/utility.h"
#include "magick/version.h"
\f
Image
*highlight_image;
- highlight_image=CompareImageChannels(image,reconstruct_image,AllChannels,
+ highlight_image=CompareImageChannels(image,reconstruct_image,CompositeChannels,
metric,distortion,exception);
return(highlight_image);
}
const Image *reconstruct_image,const ChannelType channel,
const MetricType metric,double *distortion,ExceptionInfo *exception)
{
+ CacheView
+ *highlight_view,
+ *image_view,
+ *reconstruct_view;
+
const char
*artifact;
*difference_image,
*highlight_image;
- long
+ ssize_t
y;
MagickBooleanType
lowlight,
zero;
- CacheView
- *highlight_view,
- *image_view,
- *reconstruct_view;
-
assert(image != (Image *) NULL);
assert(image->signature == MagickSignature);
if (image->debug != MagickFalse)
#if defined(MAGICKCORE_OPENMP_SUPPORT)
#pragma omp parallel for schedule(dynamic,4) shared(status)
#endif
- for (y=0; y < (long) image->rows; y++)
+ for (y=0; y < (ssize_t) image->rows; y++)
{
MagickBooleanType
sync;
register IndexPacket
*restrict highlight_indexes;
- register long
+ register ssize_t
x;
register PixelPacket
highlight_indexes=GetCacheViewAuthenticIndexQueue(highlight_view);
pixel=zero;
reconstruct_pixel=zero;
- for (x=0; x < (long) image->columns; x++)
+ for (x=0; x < (ssize_t) image->columns; x++)
{
MagickStatusType
difference;
SetMagickPixelPacket(reconstruct_image,q,reconstruct_indexes+x,
&reconstruct_pixel);
difference=MagickFalse;
- if (channel == AllChannels)
+ if (channel == CompositeChannels)
{
if (IsMagickColorSimilar(&pixel,&reconstruct_pixel) == MagickFalse)
difference=MagickTrue;
}
else
{
- if (((channel & RedChannel) != 0) && (p->red != q->red))
+ if (((channel & RedChannel) != 0) &&
+ (GetRedPixelComponent(p) != GetRedPixelComponent(q)))
difference=MagickTrue;
- if (((channel & GreenChannel) != 0) && (p->green != q->green))
+ if (((channel & GreenChannel) != 0) &&
+ (GetGreenPixelComponent(p) != GetGreenPixelComponent(q)))
difference=MagickTrue;
- if (((channel & BlueChannel) != 0) && (p->blue != q->blue))
+ if (((channel & BlueChannel) != 0) &&
+ (GetBluePixelComponent(p) != GetBluePixelComponent(q)))
difference=MagickTrue;
if (((channel & OpacityChannel) != 0) &&
- (image->matte != MagickFalse) && (p->opacity != q->opacity))
+ (image->matte != MagickFalse) &&
+ (GetOpacityPixelComponent(p) != GetOpacityPixelComponent(q)))
difference=MagickTrue;
if ((((channel & IndexChannel) != 0) &&
(image->colorspace == CMYKColorspace) &&
(reconstruct_image->colorspace == CMYKColorspace)) &&
- (indexes[x] != reconstruct_indexes[x]))
+ (GetIndexPixelComponent(indexes+x) !=
+ GetIndexPixelComponent(reconstruct_indexes+x)))
difference=MagickTrue;
}
if (difference != MagickFalse)
MagickBooleanType
status;
- status=GetImageChannelDistortion(image,reconstruct_image,AllChannels,
+ status=GetImageChannelDistortion(image,reconstruct_image,CompositeChannels,
metric,distortion,exception);
return(status);
}
-static MagickBooleanType GetAbsoluteError(const Image *image,
+static MagickBooleanType GetAbsoluteDistortion(const Image *image,
const Image *reconstruct_image,const ChannelType channel,double *distortion,
ExceptionInfo *exception)
{
- long
- y;
+ CacheView
+ *image_view,
+ *reconstruct_view;
MagickBooleanType
status;
MagickPixelPacket
zero;
- CacheView
- *image_view,
- *reconstruct_view;
+ ssize_t
+ y;
/*
Compute the absolute difference in pixels between two images.
#if defined(MAGICKCORE_OPENMP_SUPPORT)
#pragma omp parallel for schedule(dynamic,4) shared(status)
#endif
- for (y=0; y < (long) image->rows; y++)
+ for (y=0; y < (ssize_t) image->rows; y++)
{
double
- channel_distortion[AllChannels+1];
+ channel_distortion[CompositeChannels+1];
MagickPixelPacket
pixel,
*restrict p,
*restrict q;
- register long
+ register ssize_t
i,
x;
pixel=zero;
reconstruct_pixel=pixel;
(void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
- for (x=0; x < (long) image->columns; x++)
+ for (x=0; x < (ssize_t) image->columns; x++)
{
SetMagickPixelPacket(image,p,indexes+x,&pixel);
SetMagickPixelPacket(reconstruct_image,q,reconstruct_indexes+x,
if (((channel & IndexChannel) != 0) &&
(image->colorspace == CMYKColorspace))
channel_distortion[BlackChannel]++;
- channel_distortion[AllChannels]++;
+ channel_distortion[CompositeChannels]++;
}
p++;
q++;
#if defined(MAGICKCORE_OPENMP_SUPPORT)
#pragma omp critical (MagickCore_GetAbsoluteError)
#endif
- for (i=0; i <= (long) AllChannels; i++)
+ for (i=0; i <= (ssize_t) CompositeChannels; i++)
distortion[i]+=channel_distortion[i];
}
reconstruct_view=DestroyCacheView(reconstruct_view);
return(status);
}
-static unsigned long GetNumberChannels(const Image *image,
+static size_t GetNumberChannels(const Image *image,
const ChannelType channel)
{
- unsigned long
+ size_t
channels;
channels=0;
return(channels);
}
-static MagickBooleanType GetMeanAbsoluteError(const Image *image,
+static MagickBooleanType GetFuzzDistortion(const Image *image,
const Image *reconstruct_image,const ChannelType channel,
double *distortion,ExceptionInfo *exception)
{
- long
- y;
+ CacheView
+ *image_view,
+ *reconstruct_view;
MagickBooleanType
status;
- register long
+ register ssize_t
i;
+ ssize_t
+ y;
+
+ status=MagickTrue;
+ image_view=AcquireCacheView(image);
+ reconstruct_view=AcquireCacheView(reconstruct_image);
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+ #pragma omp parallel for schedule(dynamic,4) shared(status)
+#endif
+ for (y=0; y < (ssize_t) image->rows; y++)
+ {
+ double
+ channel_distortion[CompositeChannels+1];
+
+ register const IndexPacket
+ *restrict indexes,
+ *restrict reconstruct_indexes;
+
+ register const PixelPacket
+ *restrict p,
+ *restrict q;
+
+ register ssize_t
+ i,
+ x;
+
+ if (status == MagickFalse)
+ continue;
+ p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
+ q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
+ 1,exception);
+ if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
+ {
+ status=MagickFalse;
+ continue;
+ }
+ indexes=GetCacheViewVirtualIndexQueue(image_view);
+ reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
+ (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
+ for (x=0; x < (ssize_t) image->columns; x++)
+ {
+ MagickRealType
+ distance;
+
+ if ((channel & RedChannel) != 0)
+ {
+ distance=QuantumScale*(GetRedPixelComponent(p)-(MagickRealType)
+ GetRedPixelComponent(q));
+ channel_distortion[RedChannel]+=distance*distance;
+ channel_distortion[CompositeChannels]+=distance*distance;
+ }
+ if ((channel & GreenChannel) != 0)
+ {
+ distance=QuantumScale*(GetGreenPixelComponent(p)-(MagickRealType)
+ GetGreenPixelComponent(q));
+ channel_distortion[GreenChannel]+=distance*distance;
+ channel_distortion[CompositeChannels]+=distance*distance;
+ }
+ if ((channel & BlueChannel) != 0)
+ {
+ distance=QuantumScale*(GetBluePixelComponent(p)-(MagickRealType)
+ GetBluePixelComponent(q));
+ channel_distortion[BlueChannel]+=distance*distance;
+ channel_distortion[CompositeChannels]+=distance*distance;
+ }
+ if (((channel & OpacityChannel) != 0) && ((image->matte != MagickFalse) ||
+ (reconstruct_image->matte != MagickFalse)))
+ {
+ distance=QuantumScale*((image->matte != MagickFalse ?
+ GetOpacityPixelComponent(p) : OpaqueOpacity)-
+ (reconstruct_image->matte != MagickFalse ?
+ GetOpacityPixelComponent(q): OpaqueOpacity));
+ channel_distortion[OpacityChannel]+=distance*distance;
+ channel_distortion[CompositeChannels]+=distance*distance;
+ }
+ if (((channel & IndexChannel) != 0) &&
+ (image->colorspace == CMYKColorspace) &&
+ (reconstruct_image->colorspace == CMYKColorspace))
+ {
+ distance=QuantumScale*(GetIndexPixelComponent(indexes+x)-
+ (MagickRealType) GetIndexPixelComponent(reconstruct_indexes+x));
+ channel_distortion[BlackChannel]+=distance*distance;
+ channel_distortion[CompositeChannels]+=distance*distance;
+ }
+ p++;
+ q++;
+ }
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+ #pragma omp critical (MagickCore_GetMeanSquaredError)
+#endif
+ for (i=0; i <= (ssize_t) CompositeChannels; i++)
+ distortion[i]+=channel_distortion[i];
+ }
+ reconstruct_view=DestroyCacheView(reconstruct_view);
+ image_view=DestroyCacheView(image_view);
+ for (i=0; i <= (ssize_t) CompositeChannels; i++)
+ distortion[i]/=((double) image->columns*image->rows);
+ if (((channel & OpacityChannel) != 0) && ((image->matte != MagickFalse) ||
+ (reconstruct_image->matte != MagickFalse)))
+ distortion[CompositeChannels]/=(double) (GetNumberChannels(image,channel)-1);
+ else
+ distortion[CompositeChannels]/=(double) GetNumberChannels(image,channel);
+ distortion[CompositeChannels]=sqrt(distortion[CompositeChannels]);
+ return(status);
+}
+
+static MagickBooleanType GetMeanAbsoluteDistortion(const Image *image,
+ const Image *reconstruct_image,const ChannelType channel,
+ double *distortion,ExceptionInfo *exception)
+{
CacheView
*image_view,
*reconstruct_view;
+ MagickBooleanType
+ status;
+
+ register ssize_t
+ i;
+
+ ssize_t
+ y;
+
status=MagickTrue;
image_view=AcquireCacheView(image);
reconstruct_view=AcquireCacheView(reconstruct_image);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
#pragma omp parallel for schedule(dynamic,4) shared(status)
#endif
- for (y=0; y < (long) image->rows; y++)
+ for (y=0; y < (ssize_t) image->rows; y++)
{
double
- channel_distortion[AllChannels+1];
+ channel_distortion[CompositeChannels+1];
register const IndexPacket
*restrict indexes,
*restrict p,
*restrict q;
- register long
+ register ssize_t
i,
x;
indexes=GetCacheViewVirtualIndexQueue(image_view);
reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
(void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
- for (x=0; x < (long) image->columns; x++)
+ for (x=0; x < (ssize_t) image->columns; x++)
{
MagickRealType
distance;
if ((channel & RedChannel) != 0)
{
- distance=QuantumScale*fabs(p->red-(double) q->red);
+ distance=QuantumScale*fabs(GetRedPixelComponent(p)-(double)
+ GetRedPixelComponent(q));
channel_distortion[RedChannel]+=distance;
- channel_distortion[AllChannels]+=distance;
+ channel_distortion[CompositeChannels]+=distance;
}
if ((channel & GreenChannel) != 0)
{
- distance=QuantumScale*fabs(p->green-(double) q->green);
+ distance=QuantumScale*fabs(GetGreenPixelComponent(p)-(double)
+ GetGreenPixelComponent(q));
channel_distortion[GreenChannel]+=distance;
- channel_distortion[AllChannels]+=distance;
+ channel_distortion[CompositeChannels]+=distance;
}
if ((channel & BlueChannel) != 0)
{
- distance=QuantumScale*fabs(p->blue-(double) q->blue);
+ distance=QuantumScale*fabs(GetBluePixelComponent(p)-(double)
+ GetBluePixelComponent(q));
channel_distortion[BlueChannel]+=distance;
- channel_distortion[AllChannels]+=distance;
+ channel_distortion[CompositeChannels]+=distance;
}
if (((channel & OpacityChannel) != 0) &&
(image->matte != MagickFalse))
{
- distance=QuantumScale*fabs(p->opacity-(double) q->opacity);
+ distance=QuantumScale*fabs(GetOpacityPixelComponent(p)-(double)
+ GetOpacityPixelComponent(q));
channel_distortion[OpacityChannel]+=distance;
- channel_distortion[AllChannels]+=distance;
+ channel_distortion[CompositeChannels]+=distance;
}
if (((channel & IndexChannel) != 0) &&
(image->colorspace == CMYKColorspace))
{
- distance=QuantumScale*fabs(indexes[x]-(double)
- reconstruct_indexes[x]);
+ distance=QuantumScale*fabs(GetIndexPixelComponent(indexes+x)-(double)
+ GetIndexPixelComponent(reconstruct_indexes+x));
channel_distortion[BlackChannel]+=distance;
- channel_distortion[AllChannels]+=distance;
+ channel_distortion[CompositeChannels]+=distance;
}
p++;
q++;
#if defined(MAGICKCORE_OPENMP_SUPPORT)
#pragma omp critical (MagickCore_GetMeanAbsoluteError)
#endif
- for (i=0; i <= (long) AllChannels; i++)
+ for (i=0; i <= (ssize_t) CompositeChannels; i++)
distortion[i]+=channel_distortion[i];
}
reconstruct_view=DestroyCacheView(reconstruct_view);
image_view=DestroyCacheView(image_view);
- for (i=0; i <= (long) AllChannels; i++)
+ for (i=0; i <= (ssize_t) CompositeChannels; i++)
distortion[i]/=((double) image->columns*image->rows);
- distortion[AllChannels]/=(double) GetNumberChannels(image,channel);
+ distortion[CompositeChannels]/=(double) GetNumberChannels(image,channel);
return(status);
}
const Image *reconstruct_image,const ChannelType channel,double *distortion,
ExceptionInfo *exception)
{
- long
- y;
+ CacheView
+ *image_view,
+ *reconstruct_view;
MagickBooleanType
status;
maximum_error,
mean_error;
- CacheView
- *image_view,
- *reconstruct_view;
+ ssize_t
+ y;
status=MagickTrue;
alpha=1.0;
mean_error=0.0;
image_view=AcquireCacheView(image);
reconstruct_view=AcquireCacheView(reconstruct_image);
- for (y=0; y < (long) image->rows; y++)
+ for (y=0; y < (ssize_t) image->rows; y++)
{
register const IndexPacket
*restrict indexes,
*restrict p,
*restrict q;
- register long
+ register ssize_t
x;
p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
}
indexes=GetCacheViewVirtualIndexQueue(image_view);
reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
- for (x=0; x < (long) image->columns; x++)
+ for (x=0; x < (ssize_t) image->columns; x++)
{
MagickRealType
distance;
if ((channel & OpacityChannel) != 0)
{
if (image->matte != MagickFalse)
- alpha=(MagickRealType) (QuantumScale*(QuantumRange-p->opacity));
+ alpha=(MagickRealType) (QuantumScale*(GetAlphaPixelComponent(p)));
if (reconstruct_image->matte != MagickFalse)
- beta=(MagickRealType) (QuantumScale*(QuantumRange-q->opacity));
+ beta=(MagickRealType) (QuantumScale*GetAlphaPixelComponent(q));
}
if ((channel & RedChannel) != 0)
{
- distance=fabs(alpha*p->red-beta*q->red);
+ distance=fabs(alpha*GetRedPixelComponent(p)-beta*
+ GetRedPixelComponent(q));
distortion[RedChannel]+=distance;
- distortion[AllChannels]+=distance;
+ distortion[CompositeChannels]+=distance;
mean_error+=distance*distance;
if (distance > maximum_error)
maximum_error=distance;
}
if ((channel & GreenChannel) != 0)
{
- distance=fabs(alpha*p->green-beta*q->green);
+ distance=fabs(alpha*GetGreenPixelComponent(p)-beta*
+ GetGreenPixelComponent(q));
distortion[GreenChannel]+=distance;
- distortion[AllChannels]+=distance;
+ distortion[CompositeChannels]+=distance;
mean_error+=distance*distance;
if (distance > maximum_error)
maximum_error=distance;
}
if ((channel & BlueChannel) != 0)
{
- distance=fabs(alpha*p->blue-beta*q->blue);
+ distance=fabs(alpha*GetBluePixelComponent(p)-beta*
+ GetBluePixelComponent(q));
distortion[BlueChannel]+=distance;
- distortion[AllChannels]+=distance;
+ distortion[CompositeChannels]+=distance;
mean_error+=distance*distance;
if (distance > maximum_error)
maximum_error=distance;
if (((channel & OpacityChannel) != 0) &&
(image->matte != MagickFalse))
{
- distance=fabs((double) p->opacity-q->opacity);
+ distance=fabs((double) GetOpacityPixelComponent(p)-
+ GetOpacityPixelComponent(q));
distortion[OpacityChannel]+=distance;
- distortion[AllChannels]+=distance;
+ distortion[CompositeChannels]+=distance;
mean_error+=distance*distance;
if (distance > maximum_error)
maximum_error=distance;
(image->colorspace == CMYKColorspace) &&
(reconstruct_image->colorspace == CMYKColorspace))
{
- distance=fabs(alpha*indexes[x]-beta*reconstruct_indexes[x]);
+ distance=fabs(alpha*GetIndexPixelComponent(indexes+x)-beta*
+ GetIndexPixelComponent(reconstruct_indexes+x));
distortion[BlackChannel]+=distance;
- distortion[AllChannels]+=distance;
+ distortion[CompositeChannels]+=distance;
mean_error+=distance*distance;
if (distance > maximum_error)
maximum_error=distance;
}
reconstruct_view=DestroyCacheView(reconstruct_view);
image_view=DestroyCacheView(image_view);
- image->error.mean_error_per_pixel=distortion[AllChannels]/area;
+ image->error.mean_error_per_pixel=distortion[CompositeChannels]/area;
image->error.normalized_mean_error=QuantumScale*QuantumScale*mean_error/area;
image->error.normalized_maximum_error=QuantumScale*maximum_error;
return(status);
}
-static MagickBooleanType GetMeanSquaredError(const Image *image,
+static MagickBooleanType GetMeanSquaredDistortion(const Image *image,
const Image *reconstruct_image,const ChannelType channel,
double *distortion,ExceptionInfo *exception)
{
- long
- y;
+ CacheView
+ *image_view,
+ *reconstruct_view;
MagickBooleanType
status;
- register long
+ register ssize_t
i;
- CacheView
- *image_view,
- *reconstruct_view;
+ ssize_t
+ y;
status=MagickTrue;
image_view=AcquireCacheView(image);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
#pragma omp parallel for schedule(dynamic,4) shared(status)
#endif
- for (y=0; y < (long) image->rows; y++)
+ for (y=0; y < (ssize_t) image->rows; y++)
{
double
- channel_distortion[AllChannels+1];
+ channel_distortion[CompositeChannels+1];
register const IndexPacket
*restrict indexes,
*restrict p,
*restrict q;
- register long
+ register ssize_t
i,
x;
indexes=GetCacheViewVirtualIndexQueue(image_view);
reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
(void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
- for (x=0; x < (long) image->columns; x++)
+ for (x=0; x < (ssize_t) image->columns; x++)
{
MagickRealType
distance;
if ((channel & RedChannel) != 0)
{
- distance=QuantumScale*(p->red-(MagickRealType) q->red);
+ distance=QuantumScale*(GetRedPixelComponent(p)-(MagickRealType)
+ GetRedPixelComponent(q));
channel_distortion[RedChannel]+=distance*distance;
- channel_distortion[AllChannels]+=distance*distance;
+ channel_distortion[CompositeChannels]+=distance*distance;
}
if ((channel & GreenChannel) != 0)
{
- distance=QuantumScale*(p->green-(MagickRealType) q->green);
+ distance=QuantumScale*(GetGreenPixelComponent(p)-(MagickRealType)
+ GetGreenPixelComponent(q));
channel_distortion[GreenChannel]+=distance*distance;
- channel_distortion[AllChannels]+=distance*distance;
+ channel_distortion[CompositeChannels]+=distance*distance;
}
if ((channel & BlueChannel) != 0)
{
- distance=QuantumScale*(p->blue-(MagickRealType) q->blue);
+ distance=QuantumScale*(GetBluePixelComponent(p)-(MagickRealType)
+ GetBluePixelComponent(q));
channel_distortion[BlueChannel]+=distance*distance;
- channel_distortion[AllChannels]+=distance*distance;
+ channel_distortion[CompositeChannels]+=distance*distance;
}
if (((channel & OpacityChannel) != 0) &&
(image->matte != MagickFalse))
{
- distance=QuantumScale*(p->opacity-(MagickRealType) q->opacity);
+ distance=QuantumScale*(GetOpacityPixelComponent(p)-(MagickRealType)
+ GetOpacityPixelComponent(q));
channel_distortion[OpacityChannel]+=distance*distance;
- channel_distortion[AllChannels]+=distance*distance;
+ channel_distortion[CompositeChannels]+=distance*distance;
}
if (((channel & IndexChannel) != 0) &&
(image->colorspace == CMYKColorspace) &&
(reconstruct_image->colorspace == CMYKColorspace))
{
- distance=QuantumScale*(indexes[x]-(MagickRealType)
- reconstruct_indexes[x]);
+ distance=QuantumScale*(GetIndexPixelComponent(indexes+x)-
+ (MagickRealType) GetIndexPixelComponent(reconstruct_indexes+x));
channel_distortion[BlackChannel]+=distance*distance;
- channel_distortion[AllChannels]+=distance*distance;
+ channel_distortion[CompositeChannels]+=distance*distance;
}
p++;
q++;
#if defined(MAGICKCORE_OPENMP_SUPPORT)
#pragma omp critical (MagickCore_GetMeanSquaredError)
#endif
- for (i=0; i <= (long) AllChannels; i++)
+ for (i=0; i <= (ssize_t) CompositeChannels; i++)
distortion[i]+=channel_distortion[i];
}
reconstruct_view=DestroyCacheView(reconstruct_view);
image_view=DestroyCacheView(image_view);
- for (i=0; i <= (long) AllChannels; i++)
+ for (i=0; i <= (ssize_t) CompositeChannels; i++)
distortion[i]/=((double) image->columns*image->rows);
- distortion[AllChannels]/=(double) GetNumberChannels(image,channel);
+ distortion[CompositeChannels]/=(double) GetNumberChannels(image,channel);
return(status);
}
-static MagickBooleanType GetPeakAbsoluteError(const Image *image,
- const Image *reconstruct_image,const ChannelType channel,
+static MagickBooleanType GetNormalizedCrossCorrelationDistortion(
+ const Image *image,const Image *reconstruct_image,const ChannelType channel,
double *distortion,ExceptionInfo *exception)
{
- long
- y;
+#define SimilarityImageTag "Similarity/Image"
+
+ CacheView
+ *image_view,
+ *reconstruct_view;
+
+ ChannelStatistics
+ *image_statistics,
+ *reconstruct_statistics;
MagickBooleanType
status;
+ MagickOffsetType
+ progress;
+
+ MagickRealType
+ area;
+
+ register ssize_t
+ i;
+
+ ssize_t
+ y;
+
+ /*
+ Normalize to account for variation due to lighting and exposure condition.
+ */
+ image_statistics=GetImageChannelStatistics(image,exception);
+ reconstruct_statistics=GetImageChannelStatistics(reconstruct_image,exception);
+ status=MagickTrue;
+ progress=0;
+ for (i=0; i <= (ssize_t) CompositeChannels; i++)
+ distortion[i]=0.0;
+ area=1.0/((MagickRealType) image->columns*image->rows);
+ image_view=AcquireCacheView(image);
+ reconstruct_view=AcquireCacheView(reconstruct_image);
+ for (y=0; y < (ssize_t) image->rows; y++)
+ {
+ register const IndexPacket
+ *restrict indexes,
+ *restrict reconstruct_indexes;
+
+ register const PixelPacket
+ *restrict p,
+ *restrict q;
+
+ register ssize_t
+ x;
+
+ if (status == MagickFalse)
+ continue;
+ p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
+ q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
+ 1,exception);
+ if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
+ {
+ status=MagickFalse;
+ continue;
+ }
+ indexes=GetCacheViewVirtualIndexQueue(image_view);
+ reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
+ for (x=0; x < (ssize_t) image->columns; x++)
+ {
+ if ((channel & RedChannel) != 0)
+ distortion[RedChannel]+=area*QuantumScale*(GetRedPixelComponent(p)-
+ image_statistics[RedChannel].mean)*(GetRedPixelComponent(q)-
+ reconstruct_statistics[RedChannel].mean);
+ if ((channel & GreenChannel) != 0)
+ distortion[GreenChannel]+=area*QuantumScale*(GetGreenPixelComponent(p)-
+ image_statistics[GreenChannel].mean)*(GetGreenPixelComponent(q)-
+ reconstruct_statistics[GreenChannel].mean);
+ if ((channel & BlueChannel) != 0)
+ distortion[BlueChannel]+=area*QuantumScale*(GetBluePixelComponent(p)-
+ image_statistics[BlueChannel].mean)*(GetBluePixelComponent(q)-
+ reconstruct_statistics[BlueChannel].mean);
+ if (((channel & OpacityChannel) != 0) &&
+ (image->matte != MagickFalse))
+ distortion[OpacityChannel]+=area*QuantumScale*(
+ GetOpacityPixelComponent(p)-image_statistics[OpacityChannel].mean)*
+ (GetOpacityPixelComponent(q)-
+ reconstruct_statistics[OpacityChannel].mean);
+ if (((channel & IndexChannel) != 0) &&
+ (image->colorspace == CMYKColorspace) &&
+ (reconstruct_image->colorspace == CMYKColorspace))
+ distortion[BlackChannel]+=area*QuantumScale*(
+ GetIndexPixelComponent(indexes+x)-
+ image_statistics[OpacityChannel].mean)*(
+ GetIndexPixelComponent(reconstruct_indexes+x)-
+ reconstruct_statistics[OpacityChannel].mean);
+ p++;
+ q++;
+ }
+ if (image->progress_monitor != (MagickProgressMonitor) NULL)
+ {
+ MagickBooleanType
+ proceed;
+
+ proceed=SetImageProgress(image,SimilarityImageTag,progress++,
+ image->rows);
+ if (proceed == MagickFalse)
+ status=MagickFalse;
+ }
+ }
+ reconstruct_view=DestroyCacheView(reconstruct_view);
+ image_view=DestroyCacheView(image_view);
+ /*
+ Divide by the standard deviation.
+ */
+ for (i=0; i < (ssize_t) CompositeChannels; i++)
+ {
+ MagickRealType
+ gamma;
+
+ gamma=image_statistics[i].standard_deviation*
+ reconstruct_statistics[i].standard_deviation;
+ gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
+ distortion[i]=QuantumRange*gamma*distortion[i];
+ }
+ distortion[CompositeChannels]=0.0;
+ if ((channel & RedChannel) != 0)
+ distortion[CompositeChannels]+=distortion[RedChannel]*distortion[RedChannel];
+ if ((channel & GreenChannel) != 0)
+ distortion[CompositeChannels]+=distortion[GreenChannel]*distortion[GreenChannel];
+ if ((channel & BlueChannel) != 0)
+ distortion[CompositeChannels]+=distortion[BlueChannel]*distortion[BlueChannel];
+ if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
+ distortion[CompositeChannels]+=distortion[OpacityChannel]*
+ distortion[OpacityChannel];
+ if (((channel & IndexChannel) != 0) &&
+ (image->colorspace == CMYKColorspace))
+ distortion[CompositeChannels]+=distortion[BlackChannel]*distortion[BlackChannel];
+ distortion[CompositeChannels]=sqrt(distortion[CompositeChannels]/GetNumberChannels(image,
+ channel));
+ /*
+ Free resources.
+ */
+ reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
+ reconstruct_statistics);
+ image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
+ image_statistics);
+ return(status);
+}
+
+static MagickBooleanType GetPeakAbsoluteDistortion(const Image *image,
+ const Image *reconstruct_image,const ChannelType channel,
+ double *distortion,ExceptionInfo *exception)
+{
CacheView
*image_view,
*reconstruct_view;
+ MagickBooleanType
+ status;
+
+ ssize_t
+ y;
+
status=MagickTrue;
image_view=AcquireCacheView(image);
reconstruct_view=AcquireCacheView(reconstruct_image);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
#pragma omp parallel for schedule(dynamic,4) shared(status)
#endif
- for (y=0; y < (long) image->rows; y++)
+ for (y=0; y < (ssize_t) image->rows; y++)
{
double
- channel_distortion[AllChannels+1];
+ channel_distortion[CompositeChannels+1];
register const IndexPacket
*restrict indexes,
*restrict p,
*restrict q;
- register long
+ register ssize_t
i,
x;
indexes=GetCacheViewVirtualIndexQueue(image_view);
reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
(void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
- for (x=0; x < (long) image->columns; x++)
+ for (x=0; x < (ssize_t) image->columns; x++)
{
MagickRealType
distance;
if ((channel & RedChannel) != 0)
{
- distance=QuantumScale*fabs(p->red-(double) q->red);
+ distance=QuantumScale*fabs(GetRedPixelComponent(p)-(double)
+ GetRedPixelComponent(q));
if (distance > channel_distortion[RedChannel])
channel_distortion[RedChannel]=distance;
- if (distance > channel_distortion[AllChannels])
- channel_distortion[AllChannels]=distance;
+ if (distance > channel_distortion[CompositeChannels])
+ channel_distortion[CompositeChannels]=distance;
}
if ((channel & GreenChannel) != 0)
{
- distance=QuantumScale*fabs(p->green-(double) q->green);
+ distance=QuantumScale*fabs(GetGreenPixelComponent(p)-(double)
+ GetGreenPixelComponent(q));
if (distance > channel_distortion[GreenChannel])
channel_distortion[GreenChannel]=distance;
- if (distance > channel_distortion[AllChannels])
- channel_distortion[AllChannels]=distance;
+ if (distance > channel_distortion[CompositeChannels])
+ channel_distortion[CompositeChannels]=distance;
}
if ((channel & BlueChannel) != 0)
{
- distance=QuantumScale*fabs(p->blue-(double) q->blue);
+ distance=QuantumScale*fabs(GetBluePixelComponent(p)-(double)
+ GetBluePixelComponent(q));
if (distance > channel_distortion[BlueChannel])
channel_distortion[BlueChannel]=distance;
- if (distance > channel_distortion[AllChannels])
- channel_distortion[AllChannels]=distance;
+ if (distance > channel_distortion[CompositeChannels])
+ channel_distortion[CompositeChannels]=distance;
}
if (((channel & OpacityChannel) != 0) &&
(image->matte != MagickFalse))
{
- distance=QuantumScale*fabs(p->opacity-(double) q->opacity);
+ distance=QuantumScale*fabs(GetOpacityPixelComponent(p)-(double)
+ GetOpacityPixelComponent(q));
if (distance > channel_distortion[OpacityChannel])
channel_distortion[OpacityChannel]=distance;
- if (distance > channel_distortion[AllChannels])
- channel_distortion[AllChannels]=distance;
+ if (distance > channel_distortion[CompositeChannels])
+ channel_distortion[CompositeChannels]=distance;
}
if (((channel & IndexChannel) != 0) &&
(image->colorspace == CMYKColorspace) &&
(reconstruct_image->colorspace == CMYKColorspace))
{
- distance=QuantumScale*fabs(indexes[x]-(double)
- reconstruct_indexes[x]);
+ distance=QuantumScale*fabs(GetIndexPixelComponent(indexes+x)-(double)
+ GetIndexPixelComponent(reconstruct_indexes+x));
if (distance > channel_distortion[BlackChannel])
channel_distortion[BlackChannel]=distance;
- if (distance > channel_distortion[AllChannels])
- channel_distortion[AllChannels]=distance;
+ if (distance > channel_distortion[CompositeChannels])
+ channel_distortion[CompositeChannels]=distance;
}
p++;
q++;
#if defined(MAGICKCORE_OPENMP_SUPPORT)
#pragma omp critical (MagickCore_GetPeakAbsoluteError)
#endif
- for (i=0; i <= (long) AllChannels; i++)
+ for (i=0; i <= (ssize_t) CompositeChannels; i++)
if (channel_distortion[i] > distortion[i])
distortion[i]=channel_distortion[i];
}
MagickBooleanType
status;
- status=GetMeanSquaredError(image,reconstruct_image,channel,distortion,
+ status=GetMeanSquaredDistortion(image,reconstruct_image,channel,distortion,
exception);
if ((channel & RedChannel) != 0)
distortion[RedChannel]=20.0*log10((double) 1.0/sqrt(
(image->colorspace == CMYKColorspace))
distortion[BlackChannel]=20.0*log10((double) 1.0/sqrt(
distortion[BlackChannel]));
- distortion[AllChannels]=20.0*log10((double) 1.0/sqrt(
- distortion[AllChannels]));
+ distortion[CompositeChannels]=20.0*log10((double) 1.0/sqrt(
+ distortion[CompositeChannels]));
return(status);
}
-static MagickBooleanType GetRootMeanSquaredError(const Image *image,
+static MagickBooleanType GetRootMeanSquaredDistortion(const Image *image,
const Image *reconstruct_image,const ChannelType channel,
double *distortion,ExceptionInfo *exception)
{
MagickBooleanType
status;
- status=GetMeanSquaredError(image,reconstruct_image,channel,distortion,
+ status=GetMeanSquaredDistortion(image,reconstruct_image,channel,distortion,
exception);
if ((channel & RedChannel) != 0)
distortion[RedChannel]=sqrt(distortion[RedChannel]);
if (((channel & IndexChannel) != 0) &&
(image->colorspace == CMYKColorspace))
distortion[BlackChannel]=sqrt(distortion[BlackChannel]);
- distortion[AllChannels]=sqrt(distortion[AllChannels]);
+ distortion[CompositeChannels]=sqrt(distortion[CompositeChannels]);
return(status);
}
/*
Get image distortion.
*/
- length=AllChannels+1UL;
+ length=CompositeChannels+1UL;
channel_distortion=(double *) AcquireQuantumMemory(length,
sizeof(*channel_distortion));
if (channel_distortion == (double *) NULL)
{
case AbsoluteErrorMetric:
{
- status=GetAbsoluteError(image,reconstruct_image,channel,
+ status=GetAbsoluteDistortion(image,reconstruct_image,channel,
+ channel_distortion,exception);
+ break;
+ }
+ case FuzzErrorMetric:
+ {
+ status=GetFuzzDistortion(image,reconstruct_image,channel,
channel_distortion,exception);
break;
}
case MeanAbsoluteErrorMetric:
{
- status=GetMeanAbsoluteError(image,reconstruct_image,channel,
+ status=GetMeanAbsoluteDistortion(image,reconstruct_image,channel,
channel_distortion,exception);
break;
}
}
case MeanSquaredErrorMetric:
{
- status=GetMeanSquaredError(image,reconstruct_image,channel,
+ status=GetMeanSquaredDistortion(image,reconstruct_image,channel,
channel_distortion,exception);
break;
}
- case PeakAbsoluteErrorMetric:
+ case NormalizedCrossCorrelationErrorMetric:
default:
{
- status=GetPeakAbsoluteError(image,reconstruct_image,channel,
+ status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
+ channel,channel_distortion,exception);
+ break;
+ }
+ case PeakAbsoluteErrorMetric:
+ {
+ status=GetPeakAbsoluteDistortion(image,reconstruct_image,channel,
channel_distortion,exception);
break;
}
}
case RootMeanSquaredErrorMetric:
{
- status=GetRootMeanSquaredError(image,reconstruct_image,channel,
+ status=GetRootMeanSquaredDistortion(image,reconstruct_image,channel,
channel_distortion,exception);
break;
}
}
- *distortion=channel_distortion[AllChannels];
+ *distortion=channel_distortion[CompositeChannels];
channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
return(status);
}
/*
Get image distortion.
*/
- length=AllChannels+1UL;
+ length=CompositeChannels+1UL;
channel_distortion=(double *) AcquireQuantumMemory(length,
sizeof(*channel_distortion));
if (channel_distortion == (double *) NULL)
ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
(void) ResetMagickMemory(channel_distortion,0,length*
sizeof(*channel_distortion));
+ status=MagickTrue;
switch (metric)
{
case AbsoluteErrorMetric:
{
- status=GetAbsoluteError(image,reconstruct_image,AllChannels,
+ status=GetAbsoluteDistortion(image,reconstruct_image,CompositeChannels,
+ channel_distortion,exception);
+ break;
+ }
+ case FuzzErrorMetric:
+ {
+ status=GetFuzzDistortion(image,reconstruct_image,CompositeChannels,
channel_distortion,exception);
break;
}
case MeanAbsoluteErrorMetric:
{
- status=GetMeanAbsoluteError(image,reconstruct_image,AllChannels,
+ status=GetMeanAbsoluteDistortion(image,reconstruct_image,CompositeChannels,
channel_distortion,exception);
break;
}
case MeanErrorPerPixelMetric:
{
- status=GetMeanErrorPerPixel(image,reconstruct_image,AllChannels,
+ status=GetMeanErrorPerPixel(image,reconstruct_image,CompositeChannels,
channel_distortion,exception);
break;
}
case MeanSquaredErrorMetric:
{
- status=GetMeanSquaredError(image,reconstruct_image,AllChannels,
+ status=GetMeanSquaredDistortion(image,reconstruct_image,CompositeChannels,
channel_distortion,exception);
break;
}
- case PeakAbsoluteErrorMetric:
+ case NormalizedCrossCorrelationErrorMetric:
default:
{
- status=GetPeakAbsoluteError(image,reconstruct_image,AllChannels,
+ status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
+ CompositeChannels,channel_distortion,exception);
+ break;
+ }
+ case PeakAbsoluteErrorMetric:
+ {
+ status=GetPeakAbsoluteDistortion(image,reconstruct_image,CompositeChannels,
channel_distortion,exception);
break;
}
case PeakSignalToNoiseRatioMetric:
{
- status=GetPeakSignalToNoiseRatio(image,reconstruct_image,AllChannels,
+ status=GetPeakSignalToNoiseRatio(image,reconstruct_image,CompositeChannels,
channel_distortion,exception);
break;
}
case RootMeanSquaredErrorMetric:
{
- status=GetRootMeanSquaredError(image,reconstruct_image,AllChannels,
+ status=GetRootMeanSquaredDistortion(image,reconstruct_image,CompositeChannels,
channel_distortion,exception);
break;
}
}
+ if (status == MagickFalse)
+ {
+ channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
+ return((double *) NULL);
+ }
return(channel_distortion);
}
\f
MagickExport MagickBooleanType IsImagesEqual(Image *image,
const Image *reconstruct_image)
{
+ CacheView
+ *image_view,
+ *reconstruct_view;
+
ExceptionInfo
*exception;
- long
- y;
-
MagickBooleanType
status;
mean_error,
mean_error_per_pixel;
- CacheView
- *image_view,
- *reconstruct_view;
+ ssize_t
+ y;
assert(image != (Image *) NULL);
assert(image->signature == MagickSignature);
exception=(&image->exception);
image_view=AcquireCacheView(image);
reconstruct_view=AcquireCacheView(reconstruct_image);
- for (y=0; y < (long) image->rows; y++)
+ for (y=0; y < (ssize_t) image->rows; y++)
{
register const IndexPacket
*restrict indexes,
*restrict p,
*restrict q;
- register long
+ register ssize_t
x;
p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
break;
indexes=GetCacheViewVirtualIndexQueue(image_view);
reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
- for (x=0; x < (long) image->columns; x++)
+ for (x=0; x < (ssize_t) image->columns; x++)
{
MagickRealType
distance;
- distance=fabs(p->red-(double) q->red);
+ distance=fabs(GetRedPixelComponent(p)-(double)
+ GetRedPixelComponent(q));
mean_error_per_pixel+=distance;
mean_error+=distance*distance;
if (distance > maximum_error)
maximum_error=distance;
area++;
- distance=fabs(p->green-(double) q->green);
+ distance=fabs(GetGreenPixelComponent(p)-(double)
+ GetGreenPixelComponent(q));
mean_error_per_pixel+=distance;
mean_error+=distance*distance;
if (distance > maximum_error)
maximum_error=distance;
area++;
- distance=fabs(p->blue-(double) q->blue);
+ distance=fabs(GetBluePixelComponent(p)-(double)
+ GetBluePixelComponent(q));
mean_error_per_pixel+=distance;
mean_error+=distance*distance;
if (distance > maximum_error)
area++;
if (image->matte != MagickFalse)
{
- distance=fabs(p->opacity-(double) q->opacity);
+ distance=fabs(GetOpacityPixelComponent(p)-(double)
+ GetOpacityPixelComponent(q));
mean_error_per_pixel+=distance;
mean_error+=distance*distance;
if (distance > maximum_error)
if ((image->colorspace == CMYKColorspace) &&
(reconstruct_image->colorspace == CMYKColorspace))
{
- distance=fabs(indexes[x]-(double) reconstruct_indexes[x]);
+ distance=fabs(GetIndexPixelComponent(indexes+x)-(double)
+ GetIndexPixelComponent(reconstruct_indexes+x));
mean_error_per_pixel+=distance;
mean_error+=distance*distance;
if (distance > maximum_error)
%
*/
-static double GetSimilarityMetric(const Image *image,const Image *reference,
- const long x_offset,const long y_offset,ExceptionInfo *exception)
+static double GetNCCDistortion(const Image *image,
+ const Image *reconstruct_image,
+ const ChannelStatistics *reconstruct_statistics,ExceptionInfo *exception)
{
- double
- similarity;
-
- long
- y;
+#define SimilarityImageTag "Similarity/Image"
CacheView
*image_view,
- *reference_view;
+ *reconstruct_view;
+
+ ChannelStatistics
+ *image_statistics;
+
+ double
+ distortion;
MagickBooleanType
status;
+ MagickRealType
+ area,
+ gamma;
+
+ ssize_t
+ y;
+
+ unsigned long
+ number_channels;
+
/*
- Compute the similarity in pixels between two images.
+ Normalize to account for variation due to lighting and exposure condition.
*/
+ image_statistics=GetImageChannelStatistics(image,exception);
status=MagickTrue;
- similarity=0.0;
+ distortion=0.0;
+ area=1.0/((MagickRealType) image->columns*image->rows);
image_view=AcquireCacheView(image);
- reference_view=AcquireCacheView(reference);
-#if defined(MAGICKCORE_OPENMP_SUPPORT)
- #pragma omp parallel for schedule(dynamic,4) shared(status)
-#endif
- for (y=0; y < (long) reference->rows; y++)
+ reconstruct_view=AcquireCacheView(reconstruct_image);
+ for (y=0; y < (ssize_t) image->rows; y++)
{
register const IndexPacket
*restrict indexes,
- *restrict reference_indexes;
+ *restrict reconstruct_indexes;
register const PixelPacket
*restrict p,
*restrict q;
- register long
+ register ssize_t
x;
if (status == MagickFalse)
continue;
- p=GetCacheViewVirtualPixels(image_view,x_offset,y_offset+y,
- reference->columns,1,exception);
- q=GetCacheViewVirtualPixels(reference_view,0,y,reference->columns,1,
- exception);
+ p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
+ q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
+ 1,exception);
if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
{
status=MagickFalse;
continue;
}
indexes=GetCacheViewVirtualIndexQueue(image_view);
- reference_indexes=GetCacheViewVirtualIndexQueue(reference_view);
- for (x=0; x < (long) reference->columns; x++)
+ reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
+ for (x=0; x < (ssize_t) image->columns; x++)
{
- double
- thread_similarity;
-
- MagickRealType
- distance;
-
- thread_similarity=0.0;
- distance=QuantumScale*(p->red-(MagickRealType) q->red);
- thread_similarity+=distance*distance;
- distance=QuantumScale*(p->green-(MagickRealType) q->green);
- thread_similarity+=distance*distance;
- distance=QuantumScale*(p->blue-(MagickRealType) q->blue);
- thread_similarity+=distance*distance;
- if ((image->matte != MagickFalse) && (reference->matte != MagickFalse))
- {
- distance=QuantumScale*(p->opacity-(MagickRealType) q->opacity);
- thread_similarity+=distance*distance;
- }
+ distortion+=area*QuantumScale*(GetRedPixelComponent(p)-
+ image_statistics[RedChannel].mean)*(GetRedPixelComponent(q)-
+ reconstruct_statistics[RedChannel].mean);
+ distortion+=area*QuantumScale*(GetGreenPixelComponent(p)-
+ image_statistics[GreenChannel].mean)*(GetGreenPixelComponent(q)-
+ reconstruct_statistics[GreenChannel].mean);
+ distortion+=area*QuantumScale*(GetBluePixelComponent(p)-
+ image_statistics[BlueChannel].mean)*(q->blue-
+ reconstruct_statistics[BlueChannel].mean);
+ if (image->matte != MagickFalse)
+ distortion+=area*QuantumScale*(GetOpacityPixelComponent(p)-
+ image_statistics[OpacityChannel].mean)*(GetOpacityPixelComponent(q)-
+ reconstruct_statistics[OpacityChannel].mean);
if ((image->colorspace == CMYKColorspace) &&
- (reference->colorspace == CMYKColorspace))
- {
- distance=QuantumScale*(indexes[x]-(MagickRealType)
- reference_indexes[x]);
- thread_similarity+=distance*distance;
- }
-#if defined(MAGICKCORE_OPENMP_SUPPORT)
- #pragma omp critical (MagickCore_GetSimilarityMetric)
-#endif
- similarity+=thread_similarity;
+ (reconstruct_image->colorspace == CMYKColorspace))
+ distortion+=area*QuantumScale*(GetIndexPixelComponent(indexes+x)-
+ image_statistics[OpacityChannel].mean)*(GetIndexPixelComponent(
+ reconstruct_indexes+x)-reconstruct_statistics[OpacityChannel].mean);
p++;
q++;
}
}
- reference_view=DestroyCacheView(reference_view);
+ reconstruct_view=DestroyCacheView(reconstruct_view);
image_view=DestroyCacheView(image_view);
- if (status == MagickFalse)
+ /*
+ Divide by the standard deviation.
+ */
+ gamma=image_statistics[CompositeChannels].standard_deviation*
+ reconstruct_statistics[CompositeChannels].standard_deviation;
+ gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
+ distortion=QuantumRange*gamma*distortion;
+ number_channels=3;
+ if (image->matte != MagickFalse)
+ number_channels++;
+ if (image->colorspace == CMYKColorspace)
+ number_channels++;
+ distortion=sqrt(distortion/number_channels);
+ /*
+ Free resources.
+ */
+ image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
+ image_statistics);
+ return(1.0-distortion);
+}
+
+static double GetSimilarityMetric(const Image *image,const Image *reference,
+ const ChannelStatistics *reference_statistics,const ssize_t x_offset,
+ const ssize_t y_offset,ExceptionInfo *exception)
+{
+ double
+ distortion;
+
+ Image
+ *similarity_image;
+
+ RectangleInfo
+ geometry;
+
+ SetGeometry(reference,&geometry);
+ geometry.x=x_offset;
+ geometry.y=y_offset;
+ similarity_image=CropImage(image,&geometry,exception);
+ if (similarity_image == (Image *) NULL)
return(0.0);
- similarity/=((double) reference->columns*reference->rows);
- similarity/=(double) GetNumberChannels(reference,AllChannels);
- return(sqrt(similarity));
+ distortion=GetNCCDistortion(reference,similarity_image,reference_statistics,
+ exception);
+ similarity_image=DestroyImage(similarity_image);
+ return(distortion);
}
MagickExport Image *SimilarityImage(Image *image,const Image *reference,
{
#define SimilarityImageTag "Similarity/Image"
- long
- progress,
- y;
+ CacheView
+ *similarity_view;
+
+ ChannelStatistics
+ *reference_statistics;
Image
*similarity_image;
MagickBooleanType
status;
- CacheView
- *similarity_view;
+ MagickOffsetType
+ progress;
+
+ ssize_t
+ y;
assert(image != (const Image *) NULL);
assert(image->signature == MagickSignature);
assert(offset != (RectangleInfo *) NULL);
SetGeometry(reference,offset);
*similarity_metric=1.0;
- if ((reference->columns > image->columns) ||
- (reference->rows > image->rows))
+ if ((reference->columns > image->columns) || (reference->rows > image->rows))
ThrowImageException(ImageError,"ImageSizeDiffers");
similarity_image=CloneImage(image,image->columns-reference->columns+1,
image->rows-reference->rows+1,MagickTrue,exception);
*/
status=MagickTrue;
progress=0;
+ reference_statistics=GetImageChannelStatistics(reference,exception);
similarity_view=AcquireCacheView(similarity_image);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
#pragma omp parallel for schedule(dynamic,4) shared(progress,status)
#endif
- for (y=0; y < (long) (image->rows-reference->rows+1); y++)
+ for (y=0; y < (ssize_t) (image->rows-reference->rows+1); y++)
{
double
similarity;
- register long
+ register ssize_t
x;
register PixelPacket
if (status == MagickFalse)
continue;
- q=GetCacheViewAuthenticPixels(similarity_view,0,y,
- similarity_image->columns,1,exception);
+ q=GetCacheViewAuthenticPixels(similarity_view,0,y,similarity_image->columns,
+ 1,exception);
if (q == (const PixelPacket *) NULL)
{
status=MagickFalse;
continue;
}
- for (x=0; x < (long) (image->columns-reference->columns+1); x++)
+ for (x=0; x < (ssize_t) (image->columns-reference->columns+1); x++)
{
- similarity=GetSimilarityMetric(image,reference,x,y,exception);
+ similarity=GetSimilarityMetric(image,reference,reference_statistics,x,y,
+ exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
#pragma omp critical (MagickCore_SimilarityImage)
#endif
offset->x=x;
offset->y=y;
}
- q->red=RoundToQuantum(QuantumRange-QuantumRange*similarity);
- q->green=q->red;
- q->blue=q->red;
+ SetRedPixelComponent(q,ClampToQuantum(QuantumRange-QuantumRange*
+ similarity));
+ SetGreenPixelComponent(q,GetRedPixelComponent(q));
+ SetBluePixelComponent(q,GetRedPixelComponent(q));
q++;
}
if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
}
}
similarity_view=DestroyCacheView(similarity_view);
+ reference_statistics=(ChannelStatistics *) RelinquishMagickMemory(
+ reference_statistics);
return(similarity_image);
}