43#include "MagickCore/studio.h"
44#include "MagickCore/artifact.h"
45#include "MagickCore/attribute.h"
46#include "MagickCore/cache-view.h"
47#include "MagickCore/channel.h"
48#include "MagickCore/client.h"
49#include "MagickCore/color.h"
50#include "MagickCore/color-private.h"
51#include "MagickCore/colorspace.h"
52#include "MagickCore/colorspace-private.h"
53#include "MagickCore/compare.h"
54#include "MagickCore/composite-private.h"
55#include "MagickCore/constitute.h"
56#include "MagickCore/distort.h"
57#include "MagickCore/exception-private.h"
58#include "MagickCore/enhance.h"
59#include "MagickCore/fourier.h"
60#include "MagickCore/geometry.h"
61#include "MagickCore/image-private.h"
62#include "MagickCore/list.h"
63#include "MagickCore/log.h"
64#include "MagickCore/memory_.h"
65#include "MagickCore/monitor.h"
66#include "MagickCore/monitor-private.h"
67#include "MagickCore/option.h"
68#include "MagickCore/pixel-accessor.h"
69#include "MagickCore/property.h"
70#include "MagickCore/registry.h"
71#include "MagickCore/resource_.h"
72#include "MagickCore/string_.h"
73#include "MagickCore/statistic.h"
74#include "MagickCore/statistic-private.h"
75#include "MagickCore/string-private.h"
76#include "MagickCore/thread-private.h"
77#include "MagickCore/threshold.h"
78#include "MagickCore/transform.h"
79#include "MagickCore/utility.h"
80#include "MagickCore/version.h"
114MagickExport
Image *CompareImages(
Image *image,
const Image *reconstruct_image,
115 const MetricType metric,
double *distortion,
ExceptionInfo *exception)
151 assert(image != (
Image *) NULL);
152 assert(image->signature == MagickCoreSignature);
153 assert(reconstruct_image != (
const Image *) NULL);
154 assert(reconstruct_image->signature == MagickCoreSignature);
155 assert(distortion != (
double *) NULL);
156 if (IsEventLogging() != MagickFalse)
157 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
159 status=GetImageDistortion(image,reconstruct_image,metric,distortion,
161 if (status == MagickFalse)
162 return((
Image *) NULL);
163 columns=MagickMax(image->columns,reconstruct_image->columns);
164 rows=MagickMax(image->rows,reconstruct_image->rows);
165 SetGeometry(image,&geometry);
166 geometry.width=columns;
167 geometry.height=rows;
168 clone_image=CloneImage(image,0,0,MagickTrue,exception);
169 if (clone_image == (
Image *) NULL)
170 return((
Image *) NULL);
171 (void) SetImageMask(clone_image,ReadPixelMask,(
Image *) NULL,exception);
172 difference_image=ExtentImage(clone_image,&geometry,exception);
173 clone_image=DestroyImage(clone_image);
174 if (difference_image == (
Image *) NULL)
175 return((
Image *) NULL);
176 (void) ResetImagePage(difference_image,
"0x0+0+0");
177 (void) SetImageAlphaChannel(difference_image,OpaqueAlphaChannel,exception);
178 highlight_image=CloneImage(image,columns,rows,MagickTrue,exception);
179 if (highlight_image == (
Image *) NULL)
181 difference_image=DestroyImage(difference_image);
182 return((
Image *) NULL);
184 status=SetImageStorageClass(highlight_image,DirectClass,exception);
185 if (status == MagickFalse)
187 difference_image=DestroyImage(difference_image);
188 highlight_image=DestroyImage(highlight_image);
189 return((
Image *) NULL);
191 (void) SetImageMask(highlight_image,ReadPixelMask,(
Image *) NULL,exception);
192 (void) SetImageAlphaChannel(highlight_image,OpaqueAlphaChannel,exception);
193 (void) QueryColorCompliance(
"#f1001ecc",AllCompliance,&highlight,exception);
194 artifact=GetImageArtifact(image,
"compare:highlight-color");
195 if (artifact != (
const char *) NULL)
196 (void) QueryColorCompliance(artifact,AllCompliance,&highlight,exception);
197 (void) QueryColorCompliance(
"#ffffffcc",AllCompliance,&lowlight,exception);
198 artifact=GetImageArtifact(image,
"compare:lowlight-color");
199 if (artifact != (
const char *) NULL)
200 (void) QueryColorCompliance(artifact,AllCompliance,&lowlight,exception);
201 (void) QueryColorCompliance(
"#888888cc",AllCompliance,&masklight,exception);
202 artifact=GetImageArtifact(image,
"compare:masklight-color");
203 if (artifact != (
const char *) NULL)
204 (void) QueryColorCompliance(artifact,AllCompliance,&masklight,exception);
209 fuzz=GetFuzzyColorDistance(image,reconstruct_image);
210 image_view=AcquireVirtualCacheView(image,exception);
211 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
212 highlight_view=AcquireAuthenticCacheView(highlight_image,exception);
213#if defined(MAGICKCORE_OPENMP_SUPPORT)
214 #pragma omp parallel for schedule(static) shared(status) \
215 magick_number_threads(image,highlight_image,rows,1)
217 for (y=0; y < (ssize_t) rows; y++)
232 if (status == MagickFalse)
234 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
235 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
236 r=QueueCacheViewAuthenticPixels(highlight_view,0,y,columns,1,exception);
237 if ((p == (
const Quantum *) NULL) || (q == (
const Quantum *) NULL) ||
238 (r == (Quantum *) NULL))
243 for (x=0; x < (ssize_t) columns; x++)
255 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
256 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
258 SetPixelViaPixelInfo(highlight_image,&masklight,r);
259 p+=(ptrdiff_t) GetPixelChannels(image);
260 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
261 r+=(ptrdiff_t) GetPixelChannels(highlight_image);
264 difference=MagickFalse;
265 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
266 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
267 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
273 PixelChannel channel = GetPixelChannelChannel(image,i);
274 PixelTrait traits = GetPixelChannelTraits(image,channel);
275 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
277 if ((traits == UndefinedPixelTrait) ||
278 (reconstruct_traits == UndefinedPixelTrait) ||
279 ((reconstruct_traits & UpdatePixelTrait) == 0))
281 if (channel == AlphaPixelChannel)
282 pixel=(double) p[i]-(
double) GetPixelChannel(reconstruct_image,
285 pixel=Sa*(double) p[i]-Da*(
double)
286 GetPixelChannel(reconstruct_image,channel,q);
287 distance=pixel*pixel;
288 if (distance >= fuzz)
290 difference=MagickTrue;
294 if (difference == MagickFalse)
295 SetPixelViaPixelInfo(highlight_image,&lowlight,r);
297 SetPixelViaPixelInfo(highlight_image,&highlight,r);
298 p+=(ptrdiff_t) GetPixelChannels(image);
299 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
300 r+=(ptrdiff_t) GetPixelChannels(highlight_image);
302 sync=SyncCacheViewAuthenticPixels(highlight_view,exception);
303 if (sync == MagickFalse)
306 highlight_view=DestroyCacheView(highlight_view);
307 reconstruct_view=DestroyCacheView(reconstruct_view);
308 image_view=DestroyCacheView(image_view);
309 (void) CompositeImage(difference_image,highlight_image,image->compose,
310 MagickTrue,0,0,exception);
311 highlight_image=DestroyImage(highlight_image);
312 if (status == MagickFalse)
313 difference_image=DestroyImage(difference_image);
314 return(difference_image);
351static MagickBooleanType GetAbsoluteDistortion(
const Image *image,
375 fuzz=GetFuzzyColorDistance(image,reconstruct_image);
376 rows=MagickMax(image->rows,reconstruct_image->rows);
377 columns=MagickMax(image->columns,reconstruct_image->columns);
378 image_view=AcquireVirtualCacheView(image,exception);
379 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
380#if defined(MAGICKCORE_OPENMP_SUPPORT)
381 #pragma omp parallel for schedule(static) shared(status) \
382 magick_number_threads(image,image,rows,1)
384 for (y=0; y < (ssize_t) rows; y++)
387 channel_distortion[MaxPixelChannels+1];
397 if (status == MagickFalse)
399 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
400 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
401 if ((p == (
const Quantum *) NULL) || (q == (
const Quantum *) NULL))
406 (void) memset(channel_distortion,0,
sizeof(channel_distortion));
407 for (x=0; x < (ssize_t) columns; x++)
419 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
420 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
422 p+=(ptrdiff_t) GetPixelChannels(image);
423 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
426 difference=MagickFalse;
427 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
428 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
429 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
435 PixelChannel channel = GetPixelChannelChannel(image,i);
436 PixelTrait traits = GetPixelChannelTraits(image,channel);
437 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
439 if ((traits == UndefinedPixelTrait) ||
440 (reconstruct_traits == UndefinedPixelTrait) ||
441 ((reconstruct_traits & UpdatePixelTrait) == 0))
443 if (channel == AlphaPixelChannel)
444 pixel=(double) p[i]-(
double) GetPixelChannel(reconstruct_image,
447 pixel=Sa*(double) p[i]-Da*(
double) GetPixelChannel(reconstruct_image,
449 distance=pixel*pixel;
450 if (distance >= fuzz)
452 channel_distortion[i]++;
453 difference=MagickTrue;
456 if (difference != MagickFalse)
457 channel_distortion[CompositePixelChannel]++;
458 p+=(ptrdiff_t) GetPixelChannels(image);
459 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
461#if defined(MAGICKCORE_OPENMP_SUPPORT)
462 #pragma omp critical (MagickCore_GetAbsoluteDistortion)
464 for (j=0; j <= MaxPixelChannels; j++)
465 distortion[j]+=channel_distortion[j];
467 reconstruct_view=DestroyCacheView(reconstruct_view);
468 image_view=DestroyCacheView(image_view);
472static MagickBooleanType GetFuzzDistortion(
const Image *image,
496 rows=MagickMax(image->rows,reconstruct_image->rows);
497 columns=MagickMax(image->columns,reconstruct_image->columns);
499 image_view=AcquireVirtualCacheView(image,exception);
500 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
501#if defined(MAGICKCORE_OPENMP_SUPPORT)
502 #pragma omp parallel for schedule(static) shared(status) \
503 magick_number_threads(image,image,rows,1)
505 for (y=0; y < (ssize_t) rows; y++)
508 channel_distortion[MaxPixelChannels+1];
520 if (status == MagickFalse)
522 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
523 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
524 if ((p == (
const Quantum *) NULL) || (q == (Quantum *) NULL))
529 (void) memset(channel_distortion,0,
sizeof(channel_distortion));
530 for (x=0; x < (ssize_t) columns; x++)
539 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
540 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
542 p+=(ptrdiff_t) GetPixelChannels(image);
543 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
546 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
547 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
548 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
553 PixelChannel channel = GetPixelChannelChannel(image,i);
554 PixelTrait traits = GetPixelChannelTraits(image,channel);
555 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
557 if ((traits == UndefinedPixelTrait) ||
558 (reconstruct_traits == UndefinedPixelTrait) ||
559 ((reconstruct_traits & UpdatePixelTrait) == 0))
561 if (channel == AlphaPixelChannel)
562 distance=QuantumScale*((double) p[i]-(
double) GetPixelChannel(
563 reconstruct_image,channel,q));
565 distance=QuantumScale*(Sa*(double) p[i]-Da*(
double) GetPixelChannel(
566 reconstruct_image,channel,q));
567 channel_distortion[i]+=distance*distance;
568 channel_distortion[CompositePixelChannel]+=distance*distance;
571 p+=(ptrdiff_t) GetPixelChannels(image);
572 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
574#if defined(MAGICKCORE_OPENMP_SUPPORT)
575 #pragma omp critical (MagickCore_GetFuzzDistortion)
579 for (j=0; j <= MaxPixelChannels; j++)
580 distortion[j]+=channel_distortion[j];
583 reconstruct_view=DestroyCacheView(reconstruct_view);
584 image_view=DestroyCacheView(image_view);
585 area=PerceptibleReciprocal(area);
586 for (j=0; j <= MaxPixelChannels; j++)
588 distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
589 distortion[CompositePixelChannel]=sqrt(distortion[CompositePixelChannel]);
593static MagickBooleanType GetMeanAbsoluteDistortion(
const Image *image,
617 rows=MagickMax(image->rows,reconstruct_image->rows);
618 columns=MagickMax(image->columns,reconstruct_image->columns);
620 image_view=AcquireVirtualCacheView(image,exception);
621 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
622#if defined(MAGICKCORE_OPENMP_SUPPORT)
623 #pragma omp parallel for schedule(static) shared(status) \
624 magick_number_threads(image,image,rows,1)
626 for (y=0; y < (ssize_t) rows; y++)
629 channel_distortion[MaxPixelChannels+1];
641 if (status == MagickFalse)
643 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
644 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
645 if ((p == (
const Quantum *) NULL) || (q == (
const Quantum *) NULL))
650 (void) memset(channel_distortion,0,
sizeof(channel_distortion));
651 for (x=0; x < (ssize_t) columns; x++)
660 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
661 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
663 p+=(ptrdiff_t) GetPixelChannels(image);
664 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
667 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
668 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
669 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
674 PixelChannel channel = GetPixelChannelChannel(image,i);
675 PixelTrait traits = GetPixelChannelTraits(image,channel);
676 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
678 if ((traits == UndefinedPixelTrait) ||
679 (reconstruct_traits == UndefinedPixelTrait) ||
680 ((reconstruct_traits & UpdatePixelTrait) == 0))
682 if (channel == AlphaPixelChannel)
683 distance=QuantumScale*fabs((
double) p[i]-(
double)
684 GetPixelChannel(reconstruct_image,channel,q));
686 distance=QuantumScale*fabs(Sa*(
double) p[i]-Da*(
double)
687 GetPixelChannel(reconstruct_image,channel,q));
688 channel_distortion[i]+=distance;
689 channel_distortion[CompositePixelChannel]+=distance;
692 p+=(ptrdiff_t) GetPixelChannels(image);
693 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
695#if defined(MAGICKCORE_OPENMP_SUPPORT)
696 #pragma omp critical (MagickCore_GetMeanAbsoluteError)
700 for (j=0; j <= MaxPixelChannels; j++)
701 distortion[j]+=channel_distortion[j];
704 reconstruct_view=DestroyCacheView(reconstruct_view);
705 image_view=DestroyCacheView(image_view);
706 area=PerceptibleReciprocal(area);
707 for (j=0; j <= MaxPixelChannels; j++)
709 distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
713static MagickBooleanType GetMeanErrorPerPixel(
Image *image,
739 rows=MagickMax(image->rows,reconstruct_image->rows);
740 columns=MagickMax(image->columns,reconstruct_image->columns);
741 image_view=AcquireVirtualCacheView(image,exception);
742 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
743 for (y=0; y < (ssize_t) rows; y++)
752 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
753 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
754 if ((p == (
const Quantum *) NULL) || (q == (
const Quantum *) NULL))
759 for (x=0; x < (ssize_t) columns; x++)
768 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
769 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
771 p+=(ptrdiff_t) GetPixelChannels(image);
772 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
775 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
776 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
777 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
782 PixelChannel channel = GetPixelChannelChannel(image,i);
783 PixelTrait traits = GetPixelChannelTraits(image,channel);
784 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
786 if ((traits == UndefinedPixelTrait) ||
787 (reconstruct_traits == UndefinedPixelTrait) ||
788 ((reconstruct_traits & UpdatePixelTrait) == 0))
790 if (channel == AlphaPixelChannel)
791 distance=fabs((
double) p[i]-(
double)
792 GetPixelChannel(reconstruct_image,channel,q));
794 distance=fabs(Sa*(
double) p[i]-Da*(
double)
795 GetPixelChannel(reconstruct_image,channel,q));
796 distortion[i]+=distance;
797 distortion[CompositePixelChannel]+=distance;
798 mean_error+=distance*distance;
799 if (distance > maximum_error)
800 maximum_error=distance;
803 p+=(ptrdiff_t) GetPixelChannels(image);
804 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
807 reconstruct_view=DestroyCacheView(reconstruct_view);
808 image_view=DestroyCacheView(image_view);
809 area=PerceptibleReciprocal(area);
810 image->error.mean_error_per_pixel=area*distortion[CompositePixelChannel];
811 image->error.normalized_mean_error=area*QuantumScale*QuantumScale*mean_error;
812 image->error.normalized_maximum_error=QuantumScale*maximum_error;
816static MagickBooleanType GetMeanSquaredDistortion(
const Image *image,
838 rows=MagickMax(image->rows,reconstruct_image->rows);
839 columns=MagickMax(image->columns,reconstruct_image->columns);
841 image_view=AcquireVirtualCacheView(image,exception);
842 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
843#if defined(MAGICKCORE_OPENMP_SUPPORT)
844 #pragma omp parallel for schedule(static) shared(status) \
845 magick_number_threads(image,image,rows,1)
847 for (y=0; y < (ssize_t) rows; y++)
854 channel_distortion[MaxPixelChannels+1];
862 if (status == MagickFalse)
864 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
865 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
866 if ((p == (
const Quantum *) NULL) || (q == (
const Quantum *) NULL))
871 (void) memset(channel_distortion,0,
sizeof(channel_distortion));
872 for (x=0; x < (ssize_t) columns; x++)
881 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
882 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
884 p+=(ptrdiff_t) GetPixelChannels(image);
885 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
888 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
889 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
890 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
895 PixelChannel channel = GetPixelChannelChannel(image,i);
896 PixelTrait traits = GetPixelChannelTraits(image,channel);
897 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
899 if ((traits == UndefinedPixelTrait) ||
900 (reconstruct_traits == UndefinedPixelTrait) ||
901 ((reconstruct_traits & UpdatePixelTrait) == 0))
903 if (channel == AlphaPixelChannel)
904 distance=QuantumScale*((double) p[i]-(
double) GetPixelChannel(
905 reconstruct_image,channel,q));
907 distance=QuantumScale*(Sa*(double) p[i]-Da*(
double) GetPixelChannel(
908 reconstruct_image,channel,q));
909 channel_distortion[i]+=distance*distance;
910 channel_distortion[CompositePixelChannel]+=distance*distance;
913 p+=(ptrdiff_t) GetPixelChannels(image);
914 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
916#if defined(MAGICKCORE_OPENMP_SUPPORT)
917 #pragma omp critical (MagickCore_GetMeanSquaredError)
921 for (j=0; j <= MaxPixelChannels; j++)
922 distortion[j]+=channel_distortion[j];
925 reconstruct_view=DestroyCacheView(reconstruct_view);
926 image_view=DestroyCacheView(image_view);
927 area=PerceptibleReciprocal(area);
928 for (j=0; j <= MaxPixelChannels; j++)
930 distortion[CompositePixelChannel]/=GetImageChannels(image);
934static MagickBooleanType GetNormalizedCrossCorrelationDistortion(
935 const Image *image,
const Image *reconstruct_image,
double *distortion,
938#define SimilarityImageTag "Similarity/Image"
946 *reconstruct_statistics;
949 alpha_variance[MaxPixelChannels+1],
950 beta_variance[MaxPixelChannels+1];
969 image_statistics=GetImageStatistics(image,exception);
970 reconstruct_statistics=GetImageStatistics(reconstruct_image,exception);
979 reconstruct_statistics);
982 (void) memset(distortion,0,(MaxPixelChannels+1)*
sizeof(*distortion));
983 (void) memset(alpha_variance,0,(MaxPixelChannels+1)*
sizeof(*alpha_variance));
984 (void) memset(beta_variance,0,(MaxPixelChannels+1)*
sizeof(*beta_variance));
987 columns=MagickMax(image->columns,reconstruct_image->columns);
988 rows=MagickMax(image->rows,reconstruct_image->rows);
989 image_view=AcquireVirtualCacheView(image,exception);
990 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
991 for (y=0; y < (ssize_t) rows; y++)
1000 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1001 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1002 if ((p == (
const Quantum *) NULL) || (q == (
const Quantum *) NULL))
1007 for (x=0; x < (ssize_t) columns; x++)
1013 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1014 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1016 p+=(ptrdiff_t) GetPixelChannels(image);
1017 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1020 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
1021 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
1022 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1028 PixelChannel channel = GetPixelChannelChannel(image,i);
1029 PixelTrait traits = GetPixelChannelTraits(image,channel);
1030 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1032 if ((traits == UndefinedPixelTrait) ||
1033 (reconstruct_traits == UndefinedPixelTrait) ||
1034 ((reconstruct_traits & UpdatePixelTrait) == 0))
1036 if (channel == AlphaPixelChannel)
1038 alpha=QuantumScale*((double) p[i]-image_statistics[channel].mean);
1039 beta=QuantumScale*((double) GetPixelChannel(reconstruct_image,
1040 channel,q)-reconstruct_statistics[channel].mean);
1044 alpha=QuantumScale*(Sa*(double) p[i]-
1045 image_statistics[channel].mean);
1046 beta=QuantumScale*(Da*(double) GetPixelChannel(reconstruct_image,
1047 channel,q)-reconstruct_statistics[channel].mean);
1049 distortion[i]+=alpha*beta;
1050 alpha_variance[i]+=alpha*alpha;
1051 beta_variance[i]+=beta*beta;
1053 p+=(ptrdiff_t) GetPixelChannels(image);
1054 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1056 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1061#if defined(MAGICKCORE_OPENMP_SUPPORT)
1065 proceed=SetImageProgress(image,SimilarityImageTag,progress,rows);
1066 if (proceed == MagickFalse)
1073 reconstruct_view=DestroyCacheView(reconstruct_view);
1074 image_view=DestroyCacheView(image_view);
1078 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1080 distortion[i]/=sqrt(alpha_variance[i]*beta_variance[i]);
1081 if (fabs(distortion[i]) > MagickEpsilon)
1082 distortion[CompositePixelChannel]+=distortion[i];
1084 distortion[CompositePixelChannel]=distortion[CompositePixelChannel]/
1085 GetImageChannels(image);
1090 reconstruct_statistics);
1096static MagickBooleanType GetPeakAbsoluteDistortion(
const Image *image,
1114 rows=MagickMax(image->rows,reconstruct_image->rows);
1115 columns=MagickMax(image->columns,reconstruct_image->columns);
1116 image_view=AcquireVirtualCacheView(image,exception);
1117 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1118#if defined(MAGICKCORE_OPENMP_SUPPORT)
1119 #pragma omp parallel for schedule(static) shared(status) \
1120 magick_number_threads(image,image,rows,1)
1122 for (y=0; y < (ssize_t) rows; y++)
1125 channel_distortion[MaxPixelChannels+1];
1135 if (status == MagickFalse)
1137 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1138 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1139 if ((p == (
const Quantum *) NULL) || (q == (
const Quantum *) NULL))
1144 (void) memset(channel_distortion,0,
sizeof(channel_distortion));
1145 for (x=0; x < (ssize_t) columns; x++)
1154 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1155 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1157 p+=(ptrdiff_t) GetPixelChannels(image);
1158 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1161 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
1162 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
1163 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1168 PixelChannel channel = GetPixelChannelChannel(image,i);
1169 PixelTrait traits = GetPixelChannelTraits(image,channel);
1170 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1172 if ((traits == UndefinedPixelTrait) ||
1173 (reconstruct_traits == UndefinedPixelTrait) ||
1174 ((reconstruct_traits & UpdatePixelTrait) == 0))
1176 if (channel == AlphaPixelChannel)
1177 distance=QuantumScale*fabs((
double) p[i]-(
double)
1178 GetPixelChannel(reconstruct_image,channel,q));
1180 distance=QuantumScale*fabs(Sa*(
double) p[i]-Da*(
double)
1181 GetPixelChannel(reconstruct_image,channel,q));
1182 if (distance > channel_distortion[i])
1183 channel_distortion[i]=distance;
1184 if (distance > channel_distortion[CompositePixelChannel])
1185 channel_distortion[CompositePixelChannel]=distance;
1187 p+=(ptrdiff_t) GetPixelChannels(image);
1188 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1190#if defined(MAGICKCORE_OPENMP_SUPPORT)
1191 #pragma omp critical (MagickCore_GetPeakAbsoluteError)
1193 for (j=0; j <= MaxPixelChannels; j++)
1194 if (channel_distortion[j] > distortion[j])
1195 distortion[j]=channel_distortion[j];
1197 reconstruct_view=DestroyCacheView(reconstruct_view);
1198 image_view=DestroyCacheView(image_view);
1202static MagickBooleanType GetPeakSignalToNoiseRatio(
const Image *image,
1211 status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1212 for (i=0; i <= MaxPixelChannels; i++)
1213 if (fabs(distortion[i]) >= MagickEpsilon)
1214 distortion[i]=(10.0*MagickLog10(distortion[i]))/48.1647;
1218static MagickBooleanType GetPerceptualHashDistortion(
const Image *image,
1234 channel_phash=GetImagePerceptualHash(image,exception);
1236 return(MagickFalse);
1237 reconstruct_phash=GetImagePerceptualHash(reconstruct_image,exception);
1242 return(MagickFalse);
1244#if defined(MAGICKCORE_OPENMP_SUPPORT)
1245 #pragma omp parallel for schedule(static)
1247 for (channel=0; channel < MaxPixelChannels; channel++)
1256 for (i=0; i < MaximumNumberOfImageMoments; i++)
1265 for (j=0; j < (ssize_t) channel_phash[0].number_colorspaces; j++)
1267 alpha=channel_phash[channel].phash[j][i];
1268 beta=reconstruct_phash[channel].phash[j][i];
1269 difference+=(beta-alpha)*(beta-alpha);
1272 distortion[channel]+=difference;
1273#if defined(MAGICKCORE_OPENMP_SUPPORT)
1274 #pragma omp critical (MagickCore_GetPerceptualHashDistortion)
1276 distortion[CompositePixelChannel]+=difference;
1278 artifact=GetImageArtifact(image,
"phash:normalize");
1279 if ((artifact != (
const char *) NULL) &&
1280 (IsStringTrue(artifact) != MagickFalse))
1285 for (j=0; j <= MaxPixelChannels; j++)
1286 distortion[j]=sqrt(distortion[j]/channel_phash[0].number_channels);
1297static MagickBooleanType GetRootMeanSquaredDistortion(
const Image *image,
1306 status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1307 for (i=0; i <= MaxPixelChannels; i++)
1308 distortion[i]=sqrt(distortion[i]);
1312static MagickBooleanType GetStructuralSimilarityDistortion(
const Image *image,
1315#define SSIMRadius 5.0
1316#define SSIMSigma 1.5
1317#define SSIMBlocksize 8
1327 geometry[MagickPathExtent];
1358 artifact=GetImageArtifact(image,
"compare:ssim-radius");
1359 if (artifact != (
const char *) NULL)
1360 radius=StringToDouble(artifact,(
char **) NULL);
1362 artifact=GetImageArtifact(image,
"compare:ssim-sigma");
1363 if (artifact != (
const char *) NULL)
1364 sigma=StringToDouble(artifact,(
char **) NULL);
1365 (void) FormatLocaleString(geometry,MagickPathExtent,
"gaussian:%.20gx%.20g",
1367 kernel_info=AcquireKernelInfo(geometry,exception);
1369 ThrowBinaryException(ResourceLimitError,
"MemoryAllocationFailed",
1371 c1=pow(SSIMK1*SSIML,2.0);
1372 artifact=GetImageArtifact(image,
"compare:ssim-k1");
1373 if (artifact != (
const char *) NULL)
1374 c1=pow(StringToDouble(artifact,(
char **) NULL)*SSIML,2.0);
1375 c2=pow(SSIMK2*SSIML,2.0);
1376 artifact=GetImageArtifact(image,
"compare:ssim-k2");
1377 if (artifact != (
const char *) NULL)
1378 c2=pow(StringToDouble(artifact,(
char **) NULL)*SSIML,2.0);
1381 rows=MagickMax(image->rows,reconstruct_image->rows);
1382 columns=MagickMax(image->columns,reconstruct_image->columns);
1383 image_view=AcquireVirtualCacheView(image,exception);
1384 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1385#if defined(MAGICKCORE_OPENMP_SUPPORT)
1386 #pragma omp parallel for schedule(static,1) shared(area,distortion,status) \
1387 magick_number_threads(image,reconstruct_image,rows,1)
1389 for (y=0; y < (ssize_t) rows; y++)
1396 channel_distortion[MaxPixelChannels+1];
1405 if (status == MagickFalse)
1407 p=GetCacheViewVirtualPixels(image_view,-((ssize_t) kernel_info->width/2L),y-
1408 ((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
1409 kernel_info->height,exception);
1410 q=GetCacheViewVirtualPixels(reconstruct_view,-((ssize_t) kernel_info->width/
1411 2L),y-((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
1412 kernel_info->height,exception);
1413 if ((p == (
const Quantum *) NULL) || (q == (
const Quantum *) NULL))
1418 (void) memset(channel_distortion,0,
sizeof(channel_distortion));
1419 for (x=0; x < (ssize_t) columns; x++)
1422 *magick_restrict reconstruct,
1423 *magick_restrict test;
1426 x_pixel_mu[MaxPixelChannels+1],
1427 x_pixel_sigma_squared[MaxPixelChannels+1],
1428 xy_sigma[MaxPixelChannels+1],
1429 y_pixel_mu[MaxPixelChannels+1],
1430 y_pixel_sigma_squared[MaxPixelChannels+1];
1438 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1439 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1441 p+=(ptrdiff_t) GetPixelChannels(image);
1442 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1445 (void) memset(x_pixel_mu,0,
sizeof(x_pixel_mu));
1446 (void) memset(x_pixel_sigma_squared,0,
sizeof(x_pixel_sigma_squared));
1447 (void) memset(xy_sigma,0,
sizeof(xy_sigma));
1448 (void) memset(y_pixel_mu,0,
sizeof(y_pixel_mu));
1449 (void) memset(y_pixel_sigma_squared,0,
sizeof(y_pixel_sigma_squared));
1450 k=kernel_info->values;
1453 for (v=0; v < (ssize_t) kernel_info->height; v++)
1458 for (u=0; u < (ssize_t) kernel_info->width; u++)
1460 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1466 PixelChannel channel = GetPixelChannelChannel(image,i);
1467 PixelTrait traits = GetPixelChannelTraits(image,channel);
1468 PixelTrait reconstruct_traits = GetPixelChannelTraits(
1469 reconstruct_image,channel);
1470 if ((traits == UndefinedPixelTrait) ||
1471 (reconstruct_traits == UndefinedPixelTrait) ||
1472 ((reconstruct_traits & UpdatePixelTrait) == 0))
1474 x_pixel=QuantumScale*(double) reconstruct[i];
1475 x_pixel_mu[i]+=(*k)*x_pixel;
1476 x_pixel_sigma_squared[i]+=(*k)*x_pixel*x_pixel;
1477 y_pixel=QuantumScale*(double)
1478 GetPixelChannel(reconstruct_image,channel,test);
1479 y_pixel_mu[i]+=(*k)*y_pixel;
1480 y_pixel_sigma_squared[i]+=(*k)*y_pixel*y_pixel;
1481 xy_sigma[i]+=(*k)*x_pixel*y_pixel;
1484 reconstruct+=(ptrdiff_t) GetPixelChannels(image);
1485 test+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1487 reconstruct+=(ptrdiff_t) GetPixelChannels(image)*columns;
1488 test+=(ptrdiff_t) GetPixelChannels(reconstruct_image)*columns;
1490 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1495 x_pixel_sigmas_squared,
1499 y_pixel_sigmas_squared;
1501 PixelChannel channel = GetPixelChannelChannel(image,i);
1502 PixelTrait traits = GetPixelChannelTraits(image,channel);
1503 PixelTrait reconstruct_traits = GetPixelChannelTraits(
1504 reconstruct_image,channel);
1505 if ((traits == UndefinedPixelTrait) ||
1506 (reconstruct_traits == UndefinedPixelTrait) ||
1507 ((reconstruct_traits & UpdatePixelTrait) == 0))
1509 x_pixel_mu_squared=x_pixel_mu[i]*x_pixel_mu[i];
1510 y_pixel_mu_squared=y_pixel_mu[i]*y_pixel_mu[i];
1511 xy_mu=x_pixel_mu[i]*y_pixel_mu[i];
1512 xy_sigmas=xy_sigma[i]-xy_mu;
1513 x_pixel_sigmas_squared=x_pixel_sigma_squared[i]-x_pixel_mu_squared;
1514 y_pixel_sigmas_squared=y_pixel_sigma_squared[i]-y_pixel_mu_squared;
1515 ssim=((2.0*xy_mu+c1)*(2.0*xy_sigmas+c2))/
1516 ((x_pixel_mu_squared+y_pixel_mu_squared+c1)*
1517 (x_pixel_sigmas_squared+y_pixel_sigmas_squared+c2));
1518 channel_distortion[i]+=ssim;
1519 channel_distortion[CompositePixelChannel]+=ssim;
1521 p+=(ptrdiff_t) GetPixelChannels(image);
1522 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1525#if defined(MAGICKCORE_OPENMP_SUPPORT)
1526 #pragma omp critical (MagickCore_GetStructuralSimilarityDistortion)
1530 for (i=0; i <= MaxPixelChannels; i++)
1531 distortion[i]+=channel_distortion[i];
1534 image_view=DestroyCacheView(image_view);
1535 reconstruct_view=DestroyCacheView(reconstruct_view);
1536 for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
1538 PixelChannel channel = GetPixelChannelChannel(image,j);
1539 PixelTrait traits = GetPixelChannelTraits(image,channel);
1540 if ((traits == UndefinedPixelTrait) || ((traits & UpdatePixelTrait) == 0))
1542 distortion[j]*=PerceptibleReciprocal(area);
1544 distortion[CompositePixelChannel]*=PerceptibleReciprocal(area);
1545 distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
1546 kernel_info=DestroyKernelInfo(kernel_info);
1550static MagickBooleanType GetStructuralDisimilarityDistortion(
const Image *image,
1559 status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1560 distortion,exception);
1561 for (i=0; i <= MaxPixelChannels; i++)
1562 distortion[i]=(1.0-distortion[i])/2.0;
1566MagickExport MagickBooleanType GetImageDistortion(
Image *image,
1567 const Image *reconstruct_image,
const MetricType metric,
double *distortion,
1571 *channel_distortion;
1579 assert(image != (
Image *) NULL);
1580 assert(image->signature == MagickCoreSignature);
1581 assert(reconstruct_image != (
const Image *) NULL);
1582 assert(reconstruct_image->signature == MagickCoreSignature);
1583 assert(distortion != (
double *) NULL);
1584 if (IsEventLogging() != MagickFalse)
1585 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
1590 length=MaxPixelChannels+1UL;
1591 channel_distortion=(
double *) AcquireQuantumMemory(length,
1592 sizeof(*channel_distortion));
1593 if (channel_distortion == (
double *) NULL)
1594 ThrowFatalException(ResourceLimitFatalError,
"MemoryAllocationFailed");
1595 (void) memset(channel_distortion,0,length*
1596 sizeof(*channel_distortion));
1599 case AbsoluteErrorMetric:
1601 status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1605 case FuzzErrorMetric:
1607 status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1611 case MeanAbsoluteErrorMetric:
1613 status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1614 channel_distortion,exception);
1617 case MeanErrorPerPixelErrorMetric:
1619 status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1623 case MeanSquaredErrorMetric:
1625 status=GetMeanSquaredDistortion(image,reconstruct_image,
1626 channel_distortion,exception);
1629 case NormalizedCrossCorrelationErrorMetric:
1632 status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1633 channel_distortion,exception);
1636 case PeakAbsoluteErrorMetric:
1638 status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1639 channel_distortion,exception);
1642 case PeakSignalToNoiseRatioErrorMetric:
1644 status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1645 channel_distortion,exception);
1648 case PerceptualHashErrorMetric:
1650 status=GetPerceptualHashDistortion(image,reconstruct_image,
1651 channel_distortion,exception);
1654 case RootMeanSquaredErrorMetric:
1656 status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1657 channel_distortion,exception);
1660 case StructuralSimilarityErrorMetric:
1662 status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1663 channel_distortion,exception);
1666 case StructuralDissimilarityErrorMetric:
1668 status=GetStructuralDisimilarityDistortion(image,reconstruct_image,
1669 channel_distortion,exception);
1673 *distortion=channel_distortion[CompositePixelChannel];
1674 channel_distortion=(
double *) RelinquishMagickMemory(channel_distortion);
1675 (void) FormatImageProperty(image,
"distortion",
"%.*g",GetMagickPrecision(),
1712MagickExport
double *GetImageDistortions(
Image *image,
1713 const Image *reconstruct_image,
const MetricType metric,
1717 *channel_distortion;
1725 assert(image != (
Image *) NULL);
1726 assert(image->signature == MagickCoreSignature);
1727 assert(reconstruct_image != (
const Image *) NULL);
1728 assert(reconstruct_image->signature == MagickCoreSignature);
1729 if (IsEventLogging() != MagickFalse)
1730 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
1734 length=MaxPixelChannels+1UL;
1735 channel_distortion=(
double *) AcquireQuantumMemory(length,
1736 sizeof(*channel_distortion));
1737 if (channel_distortion == (
double *) NULL)
1738 ThrowFatalException(ResourceLimitFatalError,
"MemoryAllocationFailed");
1739 (void) memset(channel_distortion,0,length*
1740 sizeof(*channel_distortion));
1744 case AbsoluteErrorMetric:
1746 status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1750 case FuzzErrorMetric:
1752 status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1756 case MeanAbsoluteErrorMetric:
1758 status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1759 channel_distortion,exception);
1762 case MeanErrorPerPixelErrorMetric:
1764 status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1768 case MeanSquaredErrorMetric:
1770 status=GetMeanSquaredDistortion(image,reconstruct_image,
1771 channel_distortion,exception);
1774 case NormalizedCrossCorrelationErrorMetric:
1777 status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1778 channel_distortion,exception);
1781 case PeakAbsoluteErrorMetric:
1783 status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1784 channel_distortion,exception);
1787 case PeakSignalToNoiseRatioErrorMetric:
1789 status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1790 channel_distortion,exception);
1793 case PerceptualHashErrorMetric:
1795 status=GetPerceptualHashDistortion(image,reconstruct_image,
1796 channel_distortion,exception);
1799 case RootMeanSquaredErrorMetric:
1801 status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1802 channel_distortion,exception);
1805 case StructuralSimilarityErrorMetric:
1807 status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1808 channel_distortion,exception);
1811 case StructuralDissimilarityErrorMetric:
1813 status=GetStructuralDisimilarityDistortion(image,reconstruct_image,
1814 channel_distortion,exception);
1818 if (status == MagickFalse)
1820 channel_distortion=(
double *) RelinquishMagickMemory(channel_distortion);
1821 return((
double *) NULL);
1823 return(channel_distortion);
1854MagickExport MagickBooleanType IsImagesEqual(
const Image *image,
1868 assert(image != (
Image *) NULL);
1869 assert(image->signature == MagickCoreSignature);
1870 assert(reconstruct_image != (
const Image *) NULL);
1871 assert(reconstruct_image->signature == MagickCoreSignature);
1872 rows=MagickMax(image->rows,reconstruct_image->rows);
1873 columns=MagickMax(image->columns,reconstruct_image->columns);
1874 image_view=AcquireVirtualCacheView(image,exception);
1875 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1876 for (y=0; y < (ssize_t) rows; y++)
1885 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1886 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1887 if ((p == (
const Quantum *) NULL) || (q == (
const Quantum *) NULL))
1889 for (x=0; x < (ssize_t) columns; x++)
1894 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1899 PixelChannel channel = GetPixelChannelChannel(image,i);
1900 PixelTrait traits = GetPixelChannelTraits(image,channel);
1901 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1903 if ((traits == UndefinedPixelTrait) ||
1904 (reconstruct_traits == UndefinedPixelTrait) ||
1905 ((reconstruct_traits & UpdatePixelTrait) == 0))
1907 distance=fabs((
double) p[i]-(
double) GetPixelChannel(reconstruct_image,
1909 if (distance >= MagickEpsilon)
1912 if (i < (ssize_t) GetPixelChannels(image))
1914 p+=(ptrdiff_t) GetPixelChannels(image);
1915 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1917 if (x < (ssize_t) columns)
1920 reconstruct_view=DestroyCacheView(reconstruct_view);
1921 image_view=DestroyCacheView(image_view);
1922 return(y < (ssize_t) rows ? MagickFalse : MagickTrue);
1974MagickExport MagickBooleanType SetImageColorMetric(
Image *image,
1985 mean_error_per_pixel;
1997 assert(image != (
Image *) NULL);
1998 assert(image->signature == MagickCoreSignature);
1999 assert(reconstruct_image != (
const Image *) NULL);
2000 assert(reconstruct_image->signature == MagickCoreSignature);
2003 mean_error_per_pixel=0.0;
2005 rows=MagickMax(image->rows,reconstruct_image->rows);
2006 columns=MagickMax(image->columns,reconstruct_image->columns);
2007 image_view=AcquireVirtualCacheView(image,exception);
2008 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
2009 for (y=0; y < (ssize_t) rows; y++)
2018 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
2019 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
2020 if ((p == (
const Quantum *) NULL) || (q == (Quantum *) NULL))
2022 for (x=0; x < (ssize_t) columns; x++)
2027 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2032 PixelChannel channel = GetPixelChannelChannel(image,i);
2033 PixelTrait traits = GetPixelChannelTraits(image,channel);
2034 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2036 if ((traits == UndefinedPixelTrait) ||
2037 (reconstruct_traits == UndefinedPixelTrait) ||
2038 ((reconstruct_traits & UpdatePixelTrait) == 0))
2040 distance=fabs((
double) p[i]-(
double) GetPixelChannel(reconstruct_image,
2042 if (distance >= MagickEpsilon)
2044 mean_error_per_pixel+=distance;
2045 mean_error+=distance*distance;
2046 if (distance > maximum_error)
2047 maximum_error=distance;
2051 p+=(ptrdiff_t) GetPixelChannels(image);
2052 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
2055 reconstruct_view=DestroyCacheView(reconstruct_view);
2056 image_view=DestroyCacheView(image_view);
2057 image->error.mean_error_per_pixel=(double) (mean_error_per_pixel/area);
2058 image->error.normalized_mean_error=(double) (QuantumScale*QuantumScale*
2060 image->error.normalized_maximum_error=(double) (QuantumScale*maximum_error);
2061 status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
2108#if defined(MAGICKCORE_HDRI_SUPPORT) && defined(MAGICKCORE_FFTW_DELEGATE)
2109static Image *SIMCrossCorrelationImage(
const Image *alpha_image,
2113 *alpha_fft = (
Image *) NULL,
2114 *beta_fft = (
Image *) NULL,
2115 *complex_conjugate = (
Image *) NULL,
2116 *complex_multiplication = (
Image *) NULL,
2117 *cross_correlation = (
Image *) NULL,
2118 *temp_image = (
Image *) NULL;
2123 temp_image=CloneImage(beta_image,0,0,MagickTrue,exception);
2124 if (temp_image == (
Image *) NULL)
2125 return((
Image *) NULL);
2126 (void) SetImageArtifact(temp_image,
"fourier:normalize",
"inverse");
2127 beta_fft=ForwardFourierTransformImage(temp_image,MagickFalse,exception);
2128 temp_image=DestroyImageList(temp_image);
2129 if (beta_fft == (
Image *) NULL)
2130 return((
Image *) NULL);
2134 complex_conjugate=ComplexImages(beta_fft,ConjugateComplexOperator,exception);
2135 beta_fft=DestroyImageList(beta_fft);
2136 if (complex_conjugate == (
Image *) NULL)
2137 return((
Image *) NULL);
2141 temp_image=CloneImage(alpha_image,0,0,MagickTrue,exception);
2142 if (temp_image == (
Image *) NULL)
2144 complex_conjugate=DestroyImageList(complex_conjugate);
2145 return((
Image *) NULL);
2147 (void) SetImageArtifact(temp_image,
"fourier:normalize",
"inverse");
2148 alpha_fft=ForwardFourierTransformImage(temp_image,MagickFalse,exception);
2149 temp_image=DestroyImageList(temp_image);
2150 if (alpha_fft == (
Image *) NULL)
2152 complex_conjugate=DestroyImageList(complex_conjugate);
2153 return((
Image *) NULL);
2158 DisableCompositeClampUnlessSpecified(complex_conjugate);
2159 DisableCompositeClampUnlessSpecified(complex_conjugate->next);
2160 complex_conjugate->next->next=alpha_fft;
2161 complex_multiplication=ComplexImages(complex_conjugate,
2162 MultiplyComplexOperator,exception);
2163 complex_conjugate=DestroyImageList(complex_conjugate);
2164 if (complex_multiplication == (
Image *) NULL)
2165 return((
Image *) NULL);
2169 cross_correlation=InverseFourierTransformImage(complex_multiplication,
2170 complex_multiplication->next,MagickFalse,exception);
2171 complex_multiplication=DestroyImageList(complex_multiplication);
2172 return(cross_correlation);
2175static Image *SIMDerivativeImage(
const Image *image,
const char *kernel,
2184 kernel_info=AcquireKernelInfo(kernel,exception);
2186 return((
Image *) NULL);
2187 derivative_image=MorphologyImage(image,ConvolveMorphology,1,kernel_info,
2189 kernel_info=DestroyKernelInfo(kernel_info);
2190 return(derivative_image);
2193static Image *SIMDivideImage(
const Image *numerator_image,
2212 divide_image=CloneImage(numerator_image,0,0,MagickTrue,exception);
2213 if (divide_image == (
Image *) NULL)
2214 return(divide_image);
2216 numerator_view=AcquireAuthenticCacheView(divide_image,exception);
2217 denominator_view=AcquireVirtualCacheView(denominator_image,exception);
2218#if defined(MAGICKCORE_OPENMP_SUPPORT)
2219 #pragma omp parallel for schedule(static) shared(status) \
2220 magick_number_threads(denominator_image,divide_image,divide_image->rows,1)
2222 for (y=0; y < (ssize_t) divide_image->rows; y++)
2233 if (status == MagickFalse)
2235 p=GetCacheViewVirtualPixels(denominator_view,0,y,
2236 denominator_image->columns,1,exception);
2237 q=GetCacheViewAuthenticPixels(numerator_view,0,y,divide_image->columns,1,
2239 if ((p == (
const Quantum *) NULL) || (q == (Quantum *) NULL))
2244 for (x=0; x < (ssize_t) divide_image->columns; x++)
2249 for (i=0; i < (ssize_t) GetPixelChannels(divide_image); i++)
2251 PixelChannel channel = GetPixelChannelChannel(divide_image,i);
2252 PixelTrait traits = GetPixelChannelTraits(divide_image,channel);
2253 if ((traits & UpdatePixelTrait) == 0)
2255 if (fabs(p[i]) >= MagickEpsilon)
2256 q[i]=(Quantum) ((
double) q[i]*PerceptibleReciprocal(QuantumScale*
2259 p+=(ptrdiff_t) GetPixelChannels(denominator_image);
2260 q+=(ptrdiff_t) GetPixelChannels(divide_image);
2262 if (SyncCacheViewAuthenticPixels(numerator_view,exception) == MagickFalse)
2265 denominator_view=DestroyCacheView(denominator_view);
2266 numerator_view=DestroyCacheView(numerator_view);
2267 if (status == MagickFalse)
2268 divide_image=DestroyImage(divide_image);
2269 return(divide_image);
2272static Image *SIMDivideByMagnitude(
Image *image,
Image *magnitude_image,
2282 divide_image=SIMDivideImage(image,magnitude_image,exception);
2283 if (divide_image == (
Image *) NULL)
2284 return((
Image *) NULL);
2285 GetPixelInfoRGBA(0,0,0,0,÷_image->background_color);
2286 SetGeometry(source_image,&geometry);
2287 geometry.width=MagickMax(source_image->columns,divide_image->columns);
2288 geometry.height=MagickMax(source_image->rows,divide_image->rows);
2289 result_image=ExtentImage(divide_image,&geometry,exception);
2290 divide_image=DestroyImage(divide_image);
2291 return(result_image);
2309 image_view=AcquireAuthenticCacheView(image,exception);
2310#if defined(MAGICKCORE_OPENMP_SUPPORT)
2311 #pragma omp parallel for schedule(static) shared(status) \
2312 magick_number_threads(image,image,image->rows,1)
2314 for (y=0; y < (ssize_t) image->rows; y++)
2322 if (status == MagickFalse)
2324 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
2325 if (q == (Quantum *) NULL)
2330 for (x=0; x < (ssize_t) image->columns; x++)
2335 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2337 PixelChannel channel = GetPixelChannelChannel(image,i);
2338 PixelTrait traits = GetPixelChannelTraits(image,channel);
2339 if ((traits & UpdatePixelTrait) == 0)
2341 q[i]=(Quantum) (QuantumRange*10.0*MagickLog10((
double) q[i]));
2343 q+=(ptrdiff_t) GetPixelChannels(image);
2345 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2348 image_view=DestroyCacheView(image_view);
2369 square_image=CloneImage(image,0,0,MagickTrue,exception);
2370 if (square_image == (
Image *) NULL)
2371 return(square_image);
2373 image_view=AcquireAuthenticCacheView(square_image,exception);
2374#if defined(MAGICKCORE_OPENMP_SUPPORT)
2375 #pragma omp parallel for schedule(static) shared(status) \
2376 magick_number_threads(square_image,square_image,square_image->rows,1)
2378 for (y=0; y < (ssize_t) square_image->rows; y++)
2386 if (status == MagickFalse)
2388 q=GetCacheViewAuthenticPixels(image_view,0,y,square_image->columns,1,
2390 if (q == (Quantum *) NULL)
2395 for (x=0; x < (ssize_t) square_image->columns; x++)
2400 for (i=0; i < (ssize_t) GetPixelChannels(square_image); i++)
2402 PixelChannel channel = GetPixelChannelChannel(square_image,i);
2403 PixelTrait traits = GetPixelChannelTraits(square_image,channel);
2404 if ((traits & UpdatePixelTrait) == 0)
2406 q[i]=(Quantum) ((
double) q[i]*QuantumScale*(double) q[i]);
2408 q+=(ptrdiff_t) GetPixelChannels(square_image);
2410 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2413 image_view=DestroyCacheView(image_view);
2414 if (status == MagickFalse)
2415 square_image=DestroyImage(square_image);
2416 return(square_image);
2430 (void) SetImageArtifact(alpha_image,
"compose:clamp",
"False");
2431 xsq_image=SIMSquareImage(alpha_image,exception);
2432 if (xsq_image == (
Image *) NULL)
2433 return((
Image *) NULL);
2434 (void) SetImageArtifact(beta_image,
"compose:clamp",
"False");
2435 ysq_image=SIMSquareImage(beta_image,exception);
2436 if (ysq_image == (
Image *) NULL)
2438 xsq_image=DestroyImage(xsq_image);
2439 return((
Image *) NULL);
2441 status=CompositeImage(xsq_image,ysq_image,PlusCompositeOp,MagickTrue,0,0,
2443 magnitude_image=xsq_image;
2444 ysq_image=DestroyImage(ysq_image);
2445 if (status == MagickFalse)
2447 magnitude_image=DestroyImage(magnitude_image);
2448 return((
Image *) NULL);
2450 status=EvaluateImage(magnitude_image,PowEvaluateOperator,0.5,exception);
2451 if (status == MagickFalse)
2453 magnitude_image=DestroyImage(magnitude_image);
2454 return (
Image *) NULL;
2456 return(magnitude_image);
2459static MagickBooleanType SIMMaximaImage(
const Image *image,
double *maxima,
2475 *maxima=MagickMinimumValue;
2478 image_view=AcquireVirtualCacheView(image,exception);
2479#if defined(MAGICKCORE_OPENMP_SUPPORT)
2480 #pragma omp parallel for schedule(static) shared(status) \
2481 magick_number_threads(image,image,image->rows,1)
2483 for (y=0; y < (ssize_t) image->rows; y++)
2495 if (status == MagickFalse)
2497 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2498 if (p == (
const Quantum *) NULL)
2503 row_maxima=(double) p[0];
2505 for (x=0; x < (ssize_t) image->columns; x++)
2510 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2512 PixelChannel channel = GetPixelChannelChannel(image,i);
2513 PixelTrait traits = GetPixelChannelTraits(image,channel);
2514 if ((traits & UpdatePixelTrait) == 0)
2516 if ((
double) p[i] > row_maxima)
2518 row_maxima=(double) p[i];
2522 p+=(ptrdiff_t) GetPixelChannels(image);
2524#if defined(MAGICKCORE_OPENMP_SUPPORT)
2525 #pragma omp critical (MagickCore_SIMMaximaImage)
2527 if (row_maxima > *maxima)
2533 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2536 image_view=DestroyCacheView(image_view);
2540static MagickBooleanType SIMMinimaImage(
const Image *image,
double *minima,
2556 *minima=MagickMaximumValue;
2559 image_view=AcquireVirtualCacheView(image,exception);
2560#if defined(MAGICKCORE_OPENMP_SUPPORT)
2561 #pragma omp parallel for schedule(static) shared(status) \
2562 magick_number_threads(image,image,image->rows,1)
2564 for (y=0; y < (ssize_t) image->rows; y++)
2576 if (status == MagickFalse)
2578 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2579 if (p == (
const Quantum *) NULL)
2584 row_minima=(double) p[0];
2586 for (x=0; x < (ssize_t) image->columns; x++)
2591 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2593 PixelChannel channel = GetPixelChannelChannel(image,i);
2594 PixelTrait traits = GetPixelChannelTraits(image,channel);
2595 if ((traits & UpdatePixelTrait) == 0)
2597 if ((
double) p[i] < row_minima)
2599 row_minima=(double) p[i];
2603 p+=(ptrdiff_t) GetPixelChannels(image);
2605#if defined(MAGICKCORE_OPENMP_SUPPORT)
2606 #pragma omp critical (MagickCore_SIMMinimaImage)
2608 if (row_minima < *minima)
2614 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2617 image_view=DestroyCacheView(image_view);
2621static MagickBooleanType SIMMultiplyImage(
Image *image,
const double factor,
2637 image_view=AcquireAuthenticCacheView(image,exception);
2638#if defined(MAGICKCORE_OPENMP_SUPPORT)
2639 #pragma omp parallel for schedule(static) shared(status) \
2640 magick_number_threads(image,image,image->rows,1)
2642 for (y=0; y < (ssize_t) image->rows; y++)
2650 if (status == MagickFalse)
2652 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
2653 if (q == (Quantum *) NULL)
2658 for (x=0; x < (ssize_t) image->columns; x++)
2663 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2665 PixelChannel channel = GetPixelChannelChannel(image,i);
2666 PixelTrait traits = GetPixelChannelTraits(image,channel);
2667 if ((traits & UpdatePixelTrait) == 0)
2670 q[i]=(Quantum) ((
double) q[i]*QuantumScale*
2671 channel_statistics[channel].standard_deviation);
2672 q[i]=(Quantum) ((
double) q[i]*factor);
2674 q+=(ptrdiff_t) GetPixelChannels(image);
2676 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2679 image_view=DestroyCacheView(image_view);
2683static Image *SIMPhaseCorrelationImage(
const Image *alpha_image,
2687 *alpha_fft = (
Image *) NULL,
2688 *beta_fft = (
Image *) NULL,
2689 *complex_multiplication = (
Image *) NULL,
2690 *cross_correlation = (
Image *) NULL;
2695 beta_fft=CloneImage(beta_image,0,0,MagickTrue,exception);
2696 if (beta_fft == NULL)
2697 return((
Image *) NULL);
2698 (void) SetImageArtifact(beta_fft,
"fourier:normalize",
"inverse");
2699 beta_fft=ForwardFourierTransformImage(beta_fft,MagickFalse,exception);
2700 if (beta_fft == NULL)
2701 return((
Image *) NULL);
2705 alpha_fft=CloneImage(alpha_image,0,0,MagickTrue,exception);
2706 if (alpha_fft == (
Image *) NULL)
2708 beta_fft=DestroyImageList(beta_fft);
2709 return((
Image *) NULL);
2711 (void) SetImageArtifact(alpha_fft,
"fourier:normalize",
"inverse");
2712 alpha_fft=ForwardFourierTransformImage(alpha_fft,MagickFalse,exception);
2713 if (alpha_fft == (
Image *) NULL)
2715 beta_fft=DestroyImageList(beta_fft);
2716 return((
Image *) NULL);
2721 beta_fft=ComplexImages(beta_fft,ConjugateComplexOperator,exception);
2722 if (beta_fft == (
Image *) NULL)
2724 alpha_fft=DestroyImageList(alpha_fft);
2725 return((
Image *) NULL);
2730 beta_fft->next->next=alpha_fft;
2731 DisableCompositeClampUnlessSpecified(beta_fft);
2732 DisableCompositeClampUnlessSpecified(beta_fft->next);
2733 complex_multiplication=ComplexImages(beta_fft,MultiplyComplexOperator,
2735 beta_fft=DestroyImageList(beta_fft);
2736 if (complex_multiplication == (
Image *) NULL)
2737 return((
Image *) NULL);
2741 CompositeLayers(complex_multiplication,DivideSrcCompositeOp,(
Image *)
2742 magnitude_image,0,0,exception);
2746 (void) SetImageArtifact(complex_multiplication,
"fourier:normalize",
"inverse");
2747 cross_correlation=InverseFourierTransformImage(complex_multiplication,
2748 complex_multiplication->next,MagickFalse,exception);
2749 complex_multiplication=DestroyImageList(complex_multiplication);
2750 return(cross_correlation);
2753static MagickBooleanType SIMSetImageMean(
const Image *image,
2769 image_view=AcquireAuthenticCacheView(image,exception);
2770#if defined(MAGICKCORE_OPENMP_SUPPORT)
2771 #pragma omp parallel for schedule(static) shared(status) \
2772 magick_number_threads(image,image,image->rows,1)
2774 for (y=0; y < (ssize_t) image->rows; y++)
2782 if (status == MagickFalse)
2784 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
2785 if (q == (Quantum *) NULL)
2790 for (x=0; x < (ssize_t) image->columns; x++)
2795 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2797 PixelChannel channel = GetPixelChannelChannel(image,i);
2798 PixelTrait traits = GetPixelChannelTraits(image,channel);
2799 if ((traits & UpdatePixelTrait) == 0)
2801 q[i]=(Quantum) channel_statistics[channel].mean;
2803 q+=(ptrdiff_t) GetPixelChannels(image);
2805 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2808 image_view=DestroyCacheView(image_view);
2812static Image *SIMSubtractImageMean(
const Image *alpha_image,
2832 subtract_image=CloneImage(beta_image,alpha_image->columns,alpha_image->rows,
2833 MagickTrue,exception);
2834 if (subtract_image == (
Image *) NULL)
2835 return(subtract_image);
2837 image_view=AcquireAuthenticCacheView(subtract_image,exception);
2838 beta_view=AcquireVirtualCacheView(beta_image,exception);
2839#if defined(MAGICKCORE_OPENMP_SUPPORT)
2840 #pragma omp parallel for schedule(static) shared(status) \
2841 magick_number_threads(beta_image,subtract_image,subtract_image->rows,1)
2843 for (y=0; y < (ssize_t) subtract_image->rows; y++)
2854 if (status == MagickFalse)
2856 p=GetCacheViewVirtualPixels(beta_view,0,y,beta_image->columns,1,exception);
2857 q=GetCacheViewAuthenticPixels(image_view,0,y,subtract_image->columns,1,
2859 if ((p == (
const Quantum *) NULL) || (q == (Quantum *) NULL))
2864 for (x=0; x < (ssize_t) subtract_image->columns; x++)
2869 for (i=0; i < (ssize_t) GetPixelChannels(subtract_image); i++)
2871 PixelChannel channel = GetPixelChannelChannel(subtract_image,i);
2872 PixelTrait traits = GetPixelChannelTraits(subtract_image,channel);
2873 if ((traits & UpdatePixelTrait) == 0)
2875 if ((x >= (ssize_t) beta_image->columns) ||
2876 (y >= (ssize_t) beta_image->rows))
2879 q[i]=(Quantum) ((
double) p[i]-channel_statistics[channel].mean);
2881 p+=(ptrdiff_t) GetPixelChannels(beta_image);
2882 q+=(ptrdiff_t) GetPixelChannels(subtract_image);
2884 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2887 beta_view=DestroyCacheView(beta_view);
2888 image_view=DestroyCacheView(image_view);
2889 if (status == MagickFalse)
2890 subtract_image=DestroyImage(subtract_image);
2891 return(subtract_image);
2894static Image *SIMUnityImage(
const Image *alpha_image,
const Image *beta_image,
2912 unity_image=CloneImage(alpha_image,alpha_image->columns,alpha_image->rows,
2913 MagickTrue,exception);
2914 if (unity_image == (
Image *) NULL)
2915 return(unity_image);
2916 if (SetImageStorageClass(unity_image,DirectClass,exception) == MagickFalse)
2917 return(DestroyImage(unity_image));
2919 image_view=AcquireAuthenticCacheView(unity_image,exception);
2920#if defined(MAGICKCORE_OPENMP_SUPPORT)
2921 #pragma omp parallel for schedule(static) shared(status) \
2922 magick_number_threads(unity_image,unity_image,unity_image->rows,1)
2924 for (y=0; y < (ssize_t) unity_image->rows; y++)
2932 if (status == MagickFalse)
2934 q=GetCacheViewAuthenticPixels(image_view,0,y,unity_image->columns,1,
2936 if (q == (Quantum *) NULL)
2941 for (x=0; x < (ssize_t) unity_image->columns; x++)
2946 for (i=0; i < (ssize_t) GetPixelChannels(unity_image); i++)
2948 PixelChannel channel = GetPixelChannelChannel(unity_image,i);
2949 PixelTrait traits = GetPixelChannelTraits(unity_image,channel);
2950 if ((traits & UpdatePixelTrait) == 0)
2952 if ((x >= (ssize_t) beta_image->columns) ||
2953 (y >= (ssize_t) beta_image->rows))
2958 q+=(ptrdiff_t) GetPixelChannels(unity_image);
2960 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2963 image_view=DestroyCacheView(image_view);
2964 if (status == MagickFalse)
2965 unity_image=DestroyImage(unity_image);
2966 return(unity_image);
2969static Image *SIMVarianceImage(
Image *alpha_image,
const Image *beta_image,
2988 variance_image=CloneImage(alpha_image,0,0,MagickTrue,exception);
2989 if (variance_image == (
Image *) NULL)
2990 return(variance_image);
2992 image_view=AcquireAuthenticCacheView(variance_image,exception);
2993 beta_view=AcquireVirtualCacheView(beta_image,exception);
2994#if defined(MAGICKCORE_OPENMP_SUPPORT)
2995 #pragma omp parallel for schedule(static) shared(status) \
2996 magick_number_threads(beta_image,variance_image,variance_image->rows,1)
2998 for (y=0; y < (ssize_t) variance_image->rows; y++)
3009 if (status == MagickFalse)
3011 p=GetCacheViewVirtualPixels(beta_view,0,y,beta_image->columns,1,
3013 q=GetCacheViewAuthenticPixels(image_view,0,y,variance_image->columns,1,
3015 if ((p == (
const Quantum *) NULL) || (q == (Quantum *) NULL))
3020 for (x=0; x < (ssize_t) variance_image->columns; x++)
3025 for (i=0; i < (ssize_t) GetPixelChannels(variance_image); i++)
3027 PixelChannel channel = GetPixelChannelChannel(variance_image,i);
3028 PixelTrait traits = GetPixelChannelTraits(variance_image,channel);
3029 if ((traits & UpdatePixelTrait) == 0)
3031 q[i]=(Quantum) ((
double) ClampToQuantum((
double) QuantumRange*
3032 sqrt(fabs(QuantumScale*((
double) q[i]-(
double) p[i]))))/
3033 sqrt((
double) QuantumRange));
3035 p+=(ptrdiff_t) GetPixelChannels(beta_image);
3036 q+=(ptrdiff_t) GetPixelChannels(variance_image);
3038 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
3041 beta_view=DestroyCacheView(beta_view);
3042 image_view=DestroyCacheView(image_view);
3043 if (status == MagickFalse)
3044 variance_image=DestroyImage(variance_image);
3045 return(variance_image);
3048static Image *DPCSimilarityImage(
const Image *image,
const Image *reconstruct,
3051#define ThrowDPCSimilarityException() \
3053 if (dot_product_image != (Image *) NULL) \
3054 dot_product_image=DestroyImage(dot_product_image); \
3055 if (magnitude_image != (Image *) NULL) \
3056 magnitude_image=DestroyImage(magnitude_image); \
3057 if (reconstruct_image != (Image *) NULL) \
3058 reconstruct_image=DestroyImage(reconstruct_image); \
3059 if (rx_image != (Image *) NULL) \
3060 rx_image=DestroyImage(rx_image); \
3061 if (ry_image != (Image *) NULL) \
3062 ry_image=DestroyImage(ry_image); \
3063 if (test_image != (Image *) NULL) \
3064 test_image=DestroyImage(test_image); \
3065 if (threshold_image != (Image *) NULL) \
3066 threshold_image=DestroyImage(threshold_image); \
3067 if (trx_image != (Image *) NULL) \
3068 trx_image=DestroyImage(trx_image); \
3069 if (try_image != (Image *) NULL) \
3070 try_image=DestroyImage(try_image); \
3071 if (tx_image != (Image *) NULL) \
3072 tx_image=DestroyImage(tx_image); \
3073 if (ty_image != (Image *) NULL) \
3074 ty_image=DestroyImage(ty_image); \
3075 return((Image *) NULL); \
3082 standard_deviation = 0.0;
3085 *dot_product_image = (
Image *) NULL,
3086 *magnitude_image = (
Image *) NULL,
3087 *reconstruct_image = (
Image *) NULL,
3088 *rx_image = (
Image *) NULL,
3089 *ry_image = (
Image *) NULL,
3090 *trx_image = (
Image *) NULL,
3091 *temp_image = (
Image *) NULL,
3092 *test_image = (
Image *) NULL,
3093 *threshold_image = (
Image *) NULL,
3094 *try_image = (
Image *) NULL,
3095 *tx_image = (
Image *) NULL,
3096 *ty_image = (
Image *) NULL;
3107 test_image=CloneImage(image,0,0,MagickTrue,exception);
3108 if (test_image == (
Image *) NULL)
3109 return((
Image *) NULL);
3110 GetPixelInfoRGBA(0,0,0,0,&test_image->background_color);
3111 (void) ResetImagePage(test_image,
"0x0+0+0");
3112 status=SetImageExtent(test_image,2*(
size_t) ceil((
double) image->columns/2.0),
3113 2*(
size_t) ceil((
double) image->rows/2.0),exception);
3114 if (status == MagickFalse)
3115 ThrowDPCSimilarityException();
3116 (void) SetImageAlphaChannel(test_image,OffAlphaChannel,exception);
3120 reconstruct_image=CloneImage(reconstruct,0,0,MagickTrue,exception);
3121 if (reconstruct_image == (
Image *) NULL)
3122 ThrowDPCSimilarityException();
3123 GetPixelInfoRGBA(0,0,0,0,&reconstruct_image->background_color);
3124 (void) ResetImagePage(reconstruct_image,
"0x0+0+0");
3125 (void) SetImageAlphaChannel(reconstruct_image,OffAlphaChannel,exception);
3129 (void) SetImageVirtualPixelMethod(reconstruct_image,EdgeVirtualPixelMethod,
3131 rx_image=SIMDerivativeImage(reconstruct_image,
"Sobel",exception);
3132 if (rx_image == (
Image *) NULL)
3133 ThrowDPCSimilarityException();
3134 ry_image=SIMDerivativeImage(reconstruct_image,
"Sobel:90",exception);
3135 reconstruct_image=DestroyImage(reconstruct_image);
3136 if (ry_image == (
Image *) NULL)
3137 ThrowDPCSimilarityException();
3141 magnitude_image=SIMMagnitudeImage(rx_image,ry_image,exception);
3142 if (magnitude_image == (
Image *) NULL)
3143 ThrowDPCSimilarityException();
3147 threshold_image=CloneImage(magnitude_image,0,0,MagickTrue,exception);
3148 if (threshold_image == (
Image *) NULL)
3149 ThrowDPCSimilarityException();
3150 status=BilevelImage(threshold_image,0.0,exception);
3151 if (status == MagickFalse)
3152 ThrowDPCSimilarityException();
3153 status=GetImageMean(threshold_image,&mean,&standard_deviation,exception);
3154 threshold_image=DestroyImage(threshold_image);
3155 if (status == MagickFalse)
3156 ThrowDPCSimilarityException();
3157 edge_factor=1.0/(QuantumScale*mean)/reconstruct->columns/reconstruct->rows;
3161 temp_image=SIMDivideByMagnitude(rx_image,magnitude_image,image,exception);
3162 rx_image=DestroyImage(rx_image);
3163 if (temp_image == (
Image *) NULL)
3164 ThrowDPCSimilarityException();
3165 rx_image=temp_image;
3166 try_image=SIMDivideByMagnitude(ry_image,magnitude_image,image,exception);
3167 magnitude_image=DestroyImage(magnitude_image);
3168 ry_image=DestroyImage(ry_image);
3169 if (try_image == (
Image *) NULL)
3170 ThrowDPCSimilarityException();
3175 (void) SetImageVirtualPixelMethod(test_image,EdgeVirtualPixelMethod,
3177 tx_image=SIMDerivativeImage(test_image,
"Sobel",exception);
3178 if (tx_image == (
Image *) NULL)
3179 ThrowDPCSimilarityException();
3180 ty_image=SIMDerivativeImage(test_image,
"Sobel:90",exception);
3181 test_image=DestroyImage(test_image);
3182 if (ty_image == (
Image *) NULL)
3183 ThrowDPCSimilarityException();
3187 magnitude_image=SIMMagnitudeImage(tx_image,ty_image,exception);
3188 if (magnitude_image == (
Image *) NULL)
3189 ThrowDPCSimilarityException();
3193 temp_image=SIMDivideByMagnitude(tx_image,magnitude_image,image,exception);
3194 tx_image=DestroyImage(tx_image);
3195 if (temp_image == (
Image *) NULL)
3196 ThrowDPCSimilarityException();
3197 tx_image=temp_image;
3198 try_image=SIMDivideByMagnitude(ty_image,magnitude_image,image,exception);
3199 ty_image=DestroyImage(ty_image);
3200 magnitude_image=DestroyImage(magnitude_image);
3201 if (try_image == (
Image *) NULL)
3202 ThrowDPCSimilarityException();
3207 trx_image=SIMCrossCorrelationImage(tx_image,rx_image,exception);
3208 rx_image=DestroyImage(rx_image);
3209 tx_image=DestroyImage(tx_image);
3210 if (trx_image == (
Image *) NULL)
3211 ThrowDPCSimilarityException();
3212 try_image=SIMCrossCorrelationImage(ty_image,ry_image,exception);
3213 ry_image=DestroyImage(ry_image);
3214 ty_image=DestroyImage(ty_image);
3215 if (try_image == (
Image *) NULL)
3216 ThrowDPCSimilarityException();
3220 (void) SetImageArtifact(try_image,
"compose:clamp",
"False");
3221 status=CompositeImage(trx_image,try_image,PlusCompositeOp,MagickTrue,0,0,
3223 try_image=DestroyImage(try_image);
3224 if (status == MagickFalse)
3225 ThrowDPCSimilarityException();
3226 status=SIMMultiplyImage(trx_image,edge_factor,
3228 if (status == MagickFalse)
3229 ThrowDPCSimilarityException();
3233 SetGeometry(image,&geometry);
3234 geometry.width=image->columns;
3235 geometry.height=image->rows;
3236 (void) ResetImagePage(trx_image,
"0x0+0+0");
3238 dot_product_image=CropImage(trx_image,&geometry,exception);
3239 trx_image=DestroyImage(trx_image);
3240 if (dot_product_image == (
Image *) NULL)
3241 ThrowDPCSimilarityException();
3242 (void) ResetImagePage(dot_product_image,
"0x0+0+0");
3246 status=GrayscaleImage(dot_product_image,AveragePixelIntensityMethod,
3248 if (status == MagickFalse)
3249 ThrowDPCSimilarityException();
3250 dot_product_image->depth=MAGICKCORE_QUANTUM_DEPTH;
3251 status=SIMMaximaImage(dot_product_image,&maxima,offset,exception);
3252 if (status == MagickFalse)
3253 ThrowDPCSimilarityException();
3254 *similarity_metric=1.0-QuantumScale*maxima;
3255 return(dot_product_image);
3258static Image *MSESimilarityImage(
const Image *image,
const Image *reconstruct,
3261#define ThrowMSESimilarityException() \
3263 if (alpha_image != (Image *) NULL) \
3264 alpha_image=DestroyImage(alpha_image); \
3265 if (beta_image != (Image *) NULL) \
3266 beta_image=DestroyImage(beta_image); \
3267 if (channel_statistics != (ChannelStatistics *) NULL) \
3268 channel_statistics=(ChannelStatistics *) \
3269 RelinquishMagickMemory(channel_statistics); \
3270 if (mean_image != (Image *) NULL) \
3271 mean_image=DestroyImage(mean_image); \
3272 if (mse_image != (Image *) NULL) \
3273 mse_image=DestroyImage(mse_image); \
3274 if (reconstruct_image != (Image *) NULL) \
3275 reconstruct_image=DestroyImage(reconstruct_image); \
3276 if (sum_image != (Image *) NULL) \
3277 sum_image=DestroyImage(sum_image); \
3278 if (alpha_image != (Image *) NULL) \
3279 alpha_image=DestroyImage(alpha_image); \
3280 return((Image *) NULL); \
3290 *alpha_image = (
Image *) NULL,
3291 *beta_image = (
Image *) NULL,
3292 *mean_image = (
Image *) NULL,
3293 *mse_image = (
Image *) NULL,
3294 *reconstruct_image = (
Image *) NULL,
3295 *sum_image = (
Image *) NULL,
3296 *test_image = (
Image *) NULL;
3307 test_image=SIMSquareImage(image,exception);
3308 if (test_image == (
Image *) NULL)
3309 ThrowMSESimilarityException();
3310 reconstruct_image=SIMUnityImage(image,reconstruct,exception);
3311 if (reconstruct_image == (
Image *) NULL)
3312 ThrowMSESimilarityException();
3316 alpha_image=SIMCrossCorrelationImage(test_image,reconstruct_image,exception);
3317 test_image=DestroyImage(test_image);
3318 if (alpha_image == (
Image *) NULL)
3319 ThrowMSESimilarityException();
3320 status=SIMMultiplyImage(alpha_image,1.0/reconstruct->columns/(
double)
3322 if (status == MagickFalse)
3323 ThrowMSESimilarityException();
3327 (void) CompositeImage(reconstruct_image,reconstruct,CopyCompositeOp,
3328 MagickTrue,0,0,exception);
3329 beta_image=SIMCrossCorrelationImage(image,reconstruct_image,exception);
3330 reconstruct_image=DestroyImage(reconstruct_image);
3331 if (beta_image == (
Image *) NULL)
3332 ThrowMSESimilarityException();
3333 status=SIMMultiplyImage(beta_image,-2.0/reconstruct->columns/(
double)
3335 if (status == MagickFalse)
3336 ThrowMSESimilarityException();
3340 sum_image=SIMSquareImage(reconstruct,exception);
3341 if (sum_image == (
Image *) NULL)
3342 ThrowMSESimilarityException();
3343 channel_statistics=GetImageStatistics(sum_image,exception);
3345 ThrowMSESimilarityException();
3346 status=SetImageExtent(sum_image,image->columns,image->rows,exception);
3347 if (status == MagickFalse)
3348 ThrowMSESimilarityException();
3349 status=SetImageStorageClass(sum_image,DirectClass,exception);
3350 if (status == MagickFalse)
3351 ThrowMSESimilarityException();
3352 status=SIMSetImageMean(sum_image,channel_statistics,exception);
3354 RelinquishMagickMemory(channel_statistics);
3355 if (status == MagickFalse)
3356 ThrowMSESimilarityException();
3360 AppendImageToList(&sum_image,alpha_image);
3361 AppendImageToList(&sum_image,beta_image);
3362 mean_image=EvaluateImages(sum_image,SumEvaluateOperator,exception);
3363 if (mean_image == (
Image *) NULL)
3364 ThrowMSESimilarityException();
3365 sum_image=DestroyImage(sum_image);
3366 status=GrayscaleImage(mean_image,AveragePixelIntensityMethod,exception);
3367 if (status == MagickFalse)
3368 ThrowMSESimilarityException();
3369 mean_image->depth=MAGICKCORE_QUANTUM_DEPTH;
3373 SetGeometry(image,&geometry);
3374 geometry.width=image->columns;
3375 geometry.height=image->rows;
3376 (void) ResetImagePage(mean_image,
"0x0+0+0");
3377 mse_image=CropImage(mean_image,&geometry,exception);
3378 mean_image=DestroyImage(mean_image);
3379 if (mse_image == (
Image *) NULL)
3380 ThrowMSESimilarityException();
3384 (void) ResetImagePage(mse_image,
"0x0+0+0");
3385 (void) ClampImage(mse_image,exception);
3386 status=SIMMinimaImage(mse_image,&minima,offset,exception);
3387 if (status == MagickFalse)
3388 ThrowMSESimilarityException();
3389 status=NegateImage(mse_image,MagickFalse,exception);
3390 if (status == MagickFalse)
3391 ThrowMSESimilarityException();
3392 *similarity_metric=QuantumScale*minima;
3393 alpha_image=DestroyImage(alpha_image);
3394 beta_image=DestroyImage(beta_image);
3398static Image *NCCSimilarityImage(
const Image *image,
const Image *reconstruct,
3401#define ThrowNCCSimilarityException() \
3403 if (alpha_image != (Image *) NULL) \
3404 alpha_image=DestroyImage(alpha_image); \
3405 if (beta_image != (Image *) NULL) \
3406 beta_image=DestroyImage(beta_image); \
3407 if (channel_statistics != (ChannelStatistics *) NULL) \
3408 channel_statistics=(ChannelStatistics *) \
3409 RelinquishMagickMemory(channel_statistics); \
3410 if (correlation_image != (Image *) NULL) \
3411 correlation_image=DestroyImage(correlation_image); \
3412 if (divide_image != (Image *) NULL) \
3413 divide_image=DestroyImage(divide_image); \
3414 if (ncc_image != (Image *) NULL) \
3415 ncc_image=DestroyImage(ncc_image); \
3416 if (normalize_image != (Image *) NULL) \
3417 normalize_image=DestroyImage(normalize_image); \
3418 if (reconstruct_image != (Image *) NULL) \
3419 reconstruct_image=DestroyImage(reconstruct_image); \
3420 if (test_image != (Image *) NULL) \
3421 test_image=DestroyImage(test_image); \
3422 if (variance_image != (Image *) NULL) \
3423 variance_image=DestroyImage(variance_image); \
3424 return((Image *) NULL); \
3434 *alpha_image = (
Image *) NULL,
3435 *beta_image = (
Image *) NULL,
3436 *correlation_image = (
Image *) NULL,
3437 *divide_image = (
Image *) NULL,
3438 *ncc_image = (
Image *) NULL,
3439 *normalize_image = (
Image *) NULL,
3440 *reconstruct_image = (
Image *) NULL,
3441 *test_image = (
Image *) NULL,
3442 *variance_image = (
Image *) NULL;
3453 test_image=SIMSquareImage(image,exception);
3454 if (test_image == (
Image *) NULL)
3455 ThrowNCCSimilarityException();
3456 reconstruct_image=SIMUnityImage(image,reconstruct,exception);
3457 if (reconstruct_image == (
Image *) NULL)
3458 ThrowNCCSimilarityException();
3462 alpha_image=SIMCrossCorrelationImage(test_image,reconstruct_image,exception);
3463 test_image=DestroyImage(test_image);
3464 if (alpha_image == (
Image *) NULL)
3465 ThrowNCCSimilarityException();
3466 status=SIMMultiplyImage(alpha_image,(
double) QuantumRange*
3469 if (status == MagickFalse)
3470 ThrowNCCSimilarityException();
3474 beta_image=SIMCrossCorrelationImage(image,reconstruct_image,exception);
3475 reconstruct_image=DestroyImage(reconstruct_image);
3476 if (beta_image == (
Image *) NULL)
3477 ThrowNCCSimilarityException();
3478 test_image=SIMSquareImage(beta_image,exception);
3479 beta_image=DestroyImage(beta_image);
3480 if (test_image == (
Image *) NULL)
3481 ThrowNCCSimilarityException();
3482 status=SIMMultiplyImage(test_image,(
double) QuantumRange,
3484 if (status == MagickFalse)
3485 ThrowNCCSimilarityException();
3489 variance_image=SIMVarianceImage(alpha_image,test_image,exception);
3490 test_image=DestroyImage(test_image);
3491 alpha_image=DestroyImage(alpha_image);
3492 if (variance_image == (
Image *) NULL)
3493 ThrowNCCSimilarityException();
3497 channel_statistics=GetImageStatistics(reconstruct,exception);
3499 ThrowNCCSimilarityException();
3500 status=SIMMultiplyImage(variance_image,1.0,channel_statistics,exception);
3501 if (status == MagickFalse)
3502 ThrowNCCSimilarityException();
3503 normalize_image=SIMSubtractImageMean(image,reconstruct,channel_statistics,
3506 RelinquishMagickMemory(channel_statistics);
3507 if (normalize_image == (
Image *) NULL)
3508 ThrowNCCSimilarityException();
3509 correlation_image=SIMCrossCorrelationImage(image,normalize_image,exception);
3510 normalize_image=DestroyImage(normalize_image);
3511 if (correlation_image == (
Image *) NULL)
3512 ThrowNCCSimilarityException();
3516 divide_image=SIMDivideImage(correlation_image,variance_image,exception);
3517 correlation_image=DestroyImage(correlation_image);
3518 variance_image=DestroyImage(variance_image);
3519 if (divide_image == (
Image *) NULL)
3520 ThrowNCCSimilarityException();
3524 SetGeometry(image,&geometry);
3525 geometry.width=image->columns;
3526 geometry.height=image->rows;
3527 (void) ResetImagePage(divide_image,
"0x0+0+0");
3528 ncc_image=CropImage(divide_image,&geometry,exception);
3529 divide_image=DestroyImage(divide_image);
3530 if (ncc_image == (
Image *) NULL)
3531 ThrowNCCSimilarityException();
3535 (void) ResetImagePage(ncc_image,
"0x0+0+0");
3536 status=GrayscaleImage(ncc_image,AveragePixelIntensityMethod,exception);
3537 if (status == MagickFalse)
3538 ThrowNCCSimilarityException();
3539 ncc_image->depth=MAGICKCORE_QUANTUM_DEPTH;
3540 status=SIMMaximaImage(ncc_image,&maxima,offset,exception);
3541 if (status == MagickFalse)
3542 ThrowNCCSimilarityException();
3543 *similarity_metric=1.0-QuantumScale*maxima;
3547static Image *PhaseSimilarityImage(
const Image *image,
const Image *reconstruct,
3550#define ThrowPhaseSimilarityException() \
3552 if (correlation_image != (Image *) NULL) \
3553 correlation_image=DestroyImage(correlation_image); \
3554 if (fft_images != (Image *) NULL) \
3555 fft_images=DestroyImageList(fft_images); \
3556 if (gamma_image != (Image *) NULL) \
3557 gamma_image=DestroyImage(gamma_image); \
3558 if (magnitude_image != (Image *) NULL) \
3559 magnitude_image=DestroyImage(magnitude_image); \
3560 if (phase_image != (Image *) NULL) \
3561 phase_image=DestroyImage(phase_image); \
3562 if (reconstruct_image != (Image *) NULL) \
3563 reconstruct_image=DestroyImage(reconstruct_image); \
3564 if (reconstruct_magnitude != (Image *) NULL) \
3565 reconstruct_magnitude=DestroyImage(reconstruct_magnitude); \
3566 if (test_image != (Image *) NULL) \
3567 test_image=DestroyImage(test_image); \
3568 if (test_magnitude != (Image *) NULL) \
3569 test_magnitude=DestroyImage(test_magnitude); \
3570 return((Image *) NULL); \
3577 *correlation_image = (
Image *) NULL,
3578 *fft_images = (
Image *) NULL,
3579 *gamma_image = (
Image *) NULL,
3580 *magnitude_image = (
Image *) NULL,
3581 *phase_image = (
Image *) NULL,
3582 *reconstruct_image = (
Image *) NULL,
3583 *reconstruct_magnitude = (
Image *) NULL,
3584 *test_image = (
Image *) NULL,
3585 *test_magnitude = (
Image *) NULL;
3596 test_image=CloneImage(image,0,0,MagickTrue,exception);
3597 if (test_image == (
Image *) NULL)
3598 ThrowPhaseSimilarityException();
3599 GetPixelInfoRGBA(0,0,0,0,&test_image->background_color);
3600 (void) ResetImagePage(test_image,
"0x0+0+0");
3601 status=SetImageExtent(test_image,2*(
size_t) ceil((
double) image->columns/2.0),
3602 2*(
size_t) ceil((
double) image->rows/2.0),exception);
3603 if (status == MagickFalse)
3604 ThrowPhaseSimilarityException();
3605 (void) SetImageAlphaChannel(test_image,OffAlphaChannel,exception);
3609 reconstruct_image=CloneImage(reconstruct,0,0,MagickTrue,exception);
3610 if (reconstruct_image == (
Image *) NULL)
3611 ThrowPhaseSimilarityException();
3612 GetPixelInfoRGBA(0,0,0,0,&reconstruct_image->background_color);
3613 (void) ResetImagePage(reconstruct_image,
"0x0+0+0");
3614 status=SetImageExtent(reconstruct_image,2*(
size_t) ceil((
double)
3615 image->columns/2.0),2*(
size_t) ceil((
double) image->rows/2.0),exception);
3616 if (status == MagickFalse)
3617 ThrowPhaseSimilarityException();
3618 (void) SetImageAlphaChannel(reconstruct_image,OffAlphaChannel,exception);
3619 (void) SetImageArtifact(test_image,
"fourier:normalize",
"inverse");
3620 fft_images=ForwardFourierTransformImage(test_image,MagickTrue,exception);
3621 if (fft_images == (
Image *) NULL)
3622 ThrowPhaseSimilarityException();
3623 test_magnitude=CloneImage(fft_images,0,0,MagickTrue,exception);
3624 fft_images=DestroyImageList(fft_images);
3625 if (test_magnitude == (
Image *) NULL)
3626 ThrowPhaseSimilarityException();
3627 (void) SetImageArtifact(reconstruct_image,
"fourier:normalize",
"inverse");
3628 fft_images=ForwardFourierTransformImage(reconstruct_image,MagickTrue,
3630 if (fft_images == (
Image *) NULL)
3631 ThrowPhaseSimilarityException();
3632 reconstruct_magnitude=CloneImage(fft_images,0,0,MagickTrue,exception);
3633 fft_images=DestroyImageList(fft_images);
3634 if (reconstruct_magnitude == (
Image *) NULL)
3635 ThrowPhaseSimilarityException();
3636 magnitude_image=CloneImage(reconstruct_magnitude,0,0,MagickTrue,exception);
3637 if (magnitude_image == (
Image *) NULL)
3638 ThrowPhaseSimilarityException();
3639 DisableCompositeClampUnlessSpecified(magnitude_image);
3640 (void) CompositeImage(magnitude_image,test_magnitude,MultiplyCompositeOp,
3641 MagickTrue,0,0,exception);
3645 correlation_image=SIMPhaseCorrelationImage(test_image,reconstruct_image,
3646 magnitude_image,exception);
3647 test_image=DestroyImage(test_image);
3648 reconstruct_image=DestroyImage(reconstruct_image);
3649 test_magnitude=DestroyImage(test_magnitude);
3650 reconstruct_magnitude=DestroyImage(reconstruct_magnitude);
3651 if (correlation_image == (
Image *) NULL)
3652 ThrowPhaseSimilarityException();
3656 gamma_image=CloneImage(correlation_image,0,0,MagickTrue,exception);
3657 correlation_image=DestroyImage(correlation_image);
3658 if (gamma_image == (
Image *) NULL)
3659 ThrowPhaseSimilarityException();
3663 SetGeometry(image,&geometry);
3664 geometry.width=image->columns;
3665 geometry.height=image->rows;
3666 (void) ResetImagePage(gamma_image,
"0x0+0+0");
3667 phase_image=CropImage(gamma_image,&geometry,exception);
3668 gamma_image=DestroyImage(gamma_image);
3669 if (phase_image == (
Image *) NULL)
3670 ThrowPhaseSimilarityException();
3671 (void) ResetImagePage(phase_image,
"0x0+0+0");
3675 status=GrayscaleImage(phase_image,AveragePixelIntensityMethod,exception);
3676 if (status == MagickFalse)
3677 ThrowPhaseSimilarityException();
3678 phase_image->depth=MAGICKCORE_QUANTUM_DEPTH;
3679 status=SIMMaximaImage(phase_image,&maxima,offset,exception);
3680 if (status == MagickFalse)
3681 ThrowPhaseSimilarityException();
3682 *similarity_metric=QuantumScale*maxima;
3683 magnitude_image=DestroyImage(magnitude_image);
3684 return(phase_image);
3687static Image *PSNRSimilarityImage(
const Image *image,
const Image *reconstruct,
3690#define ThrowPSNRSimilarityException() \
3692 if (alpha_image != (Image *) NULL) \
3693 alpha_image=DestroyImage(alpha_image); \
3694 if (beta_image != (Image *) NULL) \
3695 beta_image=DestroyImage(beta_image); \
3696 if (channel_statistics != (ChannelStatistics *) NULL) \
3697 channel_statistics=(ChannelStatistics *) \
3698 RelinquishMagickMemory(channel_statistics); \
3699 if (mean_image != (Image *) NULL) \
3700 mean_image=DestroyImage(mean_image); \
3701 if (psnr_image != (Image *) NULL) \
3702 psnr_image=DestroyImage(psnr_image); \
3703 if (reconstruct_image != (Image *) NULL) \
3704 reconstruct_image=DestroyImage(reconstruct_image); \
3705 if (sum_image != (Image *) NULL) \
3706 sum_image=DestroyImage(sum_image); \
3707 if (test_image != (Image *) NULL) \
3708 test_image=DestroyImage(test_image); \
3709 return((Image *) NULL); \
3719 *alpha_image = (
Image *) NULL,
3720 *beta_image = (
Image *) NULL,
3721 *mean_image = (
Image *) NULL,
3722 *psnr_image = (
Image *) NULL,
3723 *reconstruct_image = (
Image *) NULL,
3724 *sum_image = (
Image *) NULL,
3725 *test_image = (
Image *) NULL;
3736 test_image=SIMSquareImage(image,exception);
3737 if (test_image == (
Image *) NULL)
3738 ThrowPSNRSimilarityException();
3739 reconstruct_image=SIMUnityImage(image,reconstruct,exception);
3740 if (reconstruct_image == (
Image *) NULL)
3741 ThrowPSNRSimilarityException();
3745 alpha_image=SIMCrossCorrelationImage(test_image,reconstruct_image,exception);
3746 test_image=DestroyImage(test_image);
3747 if (alpha_image == (
Image *) NULL)
3748 ThrowPSNRSimilarityException();
3749 status=SIMMultiplyImage(alpha_image,1.0/reconstruct->columns/(
double)
3751 if (status == MagickFalse)
3752 ThrowPSNRSimilarityException();
3756 (void) CompositeImage(reconstruct_image,reconstruct,CopyCompositeOp,
3757 MagickTrue,0,0,exception);
3758 beta_image=SIMCrossCorrelationImage(image,reconstruct_image,exception);
3759 reconstruct_image=DestroyImage(reconstruct_image);
3760 if (beta_image == (
Image *) NULL)
3761 ThrowPSNRSimilarityException();
3762 status=SIMMultiplyImage(beta_image,-2.0/reconstruct->columns/(
double)
3764 if (status == MagickFalse)
3765 ThrowPSNRSimilarityException();
3769 sum_image=SIMSquareImage(reconstruct,exception);
3770 if (sum_image == (
Image *) NULL)
3771 ThrowPSNRSimilarityException();
3772 channel_statistics=GetImageStatistics(sum_image,exception);
3774 ThrowPSNRSimilarityException();
3775 status=SetImageExtent(sum_image,image->columns,image->rows,exception);
3776 if (status == MagickFalse)
3777 ThrowPSNRSimilarityException();
3778 status=SetImageStorageClass(sum_image,DirectClass,exception);
3779 if (status == MagickFalse)
3780 ThrowPSNRSimilarityException();
3781 status=SIMSetImageMean(sum_image,channel_statistics,exception);
3783 RelinquishMagickMemory(channel_statistics);
3784 if (status == MagickFalse)
3785 ThrowPSNRSimilarityException();
3789 AppendImageToList(&sum_image,alpha_image);
3790 AppendImageToList(&sum_image,beta_image);
3791 mean_image=EvaluateImages(sum_image,SumEvaluateOperator,exception);
3792 if (mean_image == (
Image *) NULL)
3793 ThrowPSNRSimilarityException();
3794 sum_image=DestroyImage(sum_image);
3795 status=SIMLogImage(mean_image,exception);
3796 if (status == MagickFalse)
3797 ThrowPSNRSimilarityException();
3798 status=GrayscaleImage(mean_image,AveragePixelIntensityMethod,exception);
3799 if (status == MagickFalse)
3800 ThrowPSNRSimilarityException();
3801 mean_image->depth=MAGICKCORE_QUANTUM_DEPTH;
3802 status=SIMMultiplyImage(mean_image,1.0/48.1647,
3804 if (status == MagickFalse)
3805 ThrowPSNRSimilarityException();
3809 SetGeometry(image,&geometry);
3810 geometry.width=image->columns;
3811 geometry.height=image->rows;
3812 (void) ResetImagePage(mean_image,
"0x0+0+0");
3813 psnr_image=CropImage(mean_image,&geometry,exception);
3814 mean_image=DestroyImage(mean_image);
3815 if (psnr_image == (
Image *) NULL)
3816 ThrowPSNRSimilarityException();
3820 (void) ResetImagePage(psnr_image,
"0x0+0+0");
3821 (void) EvaluateImage(psnr_image,MaxEvaluateOperator,0.0,exception);
3822 status=SIMMinimaImage(psnr_image,&minima,offset,exception);
3823 if (status == MagickFalse)
3824 ThrowPSNRSimilarityException();
3825 status=NegateImage(psnr_image,MagickFalse,exception);
3826 if (status == MagickFalse)
3827 ThrowPSNRSimilarityException();
3828 *similarity_metric=(-1.0+QuantumScale*minima);
3829 alpha_image=DestroyImage(alpha_image);
3830 beta_image=DestroyImage(beta_image);
3834static Image *RMSESimilarityImage(
const Image *image,
const Image *reconstruct,
3837#define ThrowRMSESimilarityException() \
3839 if (alpha_image != (Image *) NULL) \
3840 alpha_image=DestroyImage(alpha_image); \
3841 if (beta_image != (Image *) NULL) \
3842 beta_image=DestroyImage(beta_image); \
3843 if (channel_statistics != (ChannelStatistics *) NULL) \
3844 channel_statistics=(ChannelStatistics *) \
3845 RelinquishMagickMemory(channel_statistics); \
3846 if (mean_image != (Image *) NULL) \
3847 mean_image=DestroyImage(mean_image); \
3848 if (rmse_image != (Image *) NULL) \
3849 rmse_image=DestroyImage(rmse_image); \
3850 if (reconstruct_image != (Image *) NULL) \
3851 reconstruct_image=DestroyImage(reconstruct_image); \
3852 if (sum_image != (Image *) NULL) \
3853 sum_image=DestroyImage(sum_image); \
3854 if (alpha_image != (Image *) NULL) \
3855 alpha_image=DestroyImage(alpha_image); \
3856 return((Image *) NULL); \
3866 *alpha_image = (
Image *) NULL,
3867 *beta_image = (
Image *) NULL,
3868 *mean_image = (
Image *) NULL,
3869 *reconstruct_image = (
Image *) NULL,
3870 *rmse_image = (
Image *) NULL,
3871 *sum_image = (
Image *) NULL,
3872 *test_image = (
Image *) NULL;
3883 test_image=SIMSquareImage(image,exception);
3884 if (test_image == (
Image *) NULL)
3885 ThrowRMSESimilarityException();
3886 reconstruct_image=SIMUnityImage(image,reconstruct,exception);
3887 if (reconstruct_image == (
Image *) NULL)
3888 ThrowRMSESimilarityException();
3892 alpha_image=SIMCrossCorrelationImage(test_image,reconstruct_image,exception);
3893 test_image=DestroyImage(test_image);
3894 if (alpha_image == (
Image *) NULL)
3895 ThrowRMSESimilarityException();
3896 status=SIMMultiplyImage(alpha_image,1.0/reconstruct->columns/(
double)
3898 if (status == MagickFalse)
3899 ThrowRMSESimilarityException();
3903 (void) CompositeImage(reconstruct_image,reconstruct,CopyCompositeOp,
3904 MagickTrue,0,0,exception);
3905 beta_image=SIMCrossCorrelationImage(image,reconstruct_image,exception);
3906 reconstruct_image=DestroyImage(reconstruct_image);
3907 if (beta_image == (
Image *) NULL)
3908 ThrowRMSESimilarityException();
3909 status=SIMMultiplyImage(beta_image,-2.0/reconstruct->columns/(
double)
3911 if (status == MagickFalse)
3912 ThrowRMSESimilarityException();
3916 sum_image=SIMSquareImage(reconstruct,exception);
3917 if (sum_image == (
Image *) NULL)
3918 ThrowRMSESimilarityException();
3919 channel_statistics=GetImageStatistics(sum_image,exception);
3921 ThrowRMSESimilarityException();
3922 status=SetImageExtent(sum_image,image->columns,image->rows,exception);
3923 if (status == MagickFalse)
3924 ThrowRMSESimilarityException();
3925 status=SetImageStorageClass(sum_image,DirectClass,exception);
3926 if (status == MagickFalse)
3927 ThrowRMSESimilarityException();
3928 status=SIMSetImageMean(sum_image,channel_statistics,exception);
3930 RelinquishMagickMemory(channel_statistics);
3931 if (status == MagickFalse)
3932 ThrowRMSESimilarityException();
3936 AppendImageToList(&sum_image,alpha_image);
3937 AppendImageToList(&sum_image,beta_image);
3938 mean_image=EvaluateImages(sum_image,SumEvaluateOperator,exception);
3939 if (mean_image == (
Image *) NULL)
3940 ThrowRMSESimilarityException();
3941 status=EvaluateImage(mean_image,PowEvaluateOperator,0.5,exception);
3942 if (mean_image == (
Image *) NULL)
3943 ThrowRMSESimilarityException();
3944 sum_image=DestroyImage(sum_image);
3945 status=GrayscaleImage(mean_image,AveragePixelIntensityMethod,exception);
3946 if (status == MagickFalse)
3947 ThrowRMSESimilarityException();
3948 mean_image->depth=MAGICKCORE_QUANTUM_DEPTH;
3952 SetGeometry(image,&geometry);
3953 geometry.width=image->columns;
3954 geometry.height=image->rows;
3955 (void) ResetImagePage(mean_image,
"0x0+0+0");
3956 rmse_image=CropImage(mean_image,&geometry,exception);
3957 mean_image=DestroyImage(mean_image);
3958 if (rmse_image == (
Image *) NULL)
3959 ThrowRMSESimilarityException();
3963 (void) ResetImagePage(rmse_image,
"0x0+0+0");
3964 status=SIMMinimaImage(rmse_image,&minima,offset,exception);
3965 if (status == MagickFalse)
3966 ThrowRMSESimilarityException();
3967 status=NegateImage(rmse_image,MagickFalse,exception);
3968 if (status == MagickFalse)
3969 ThrowRMSESimilarityException();
3970 *similarity_metric=QuantumScale*minima;
3971 alpha_image=DestroyImage(alpha_image);
3972 beta_image=DestroyImage(beta_image);
3978static double GetSimilarityMetric(
const Image *image,
const Image *reconstruct,
3979 const MetricType metric,
const ssize_t x_offset,
const ssize_t y_offset,
3986 *sans_exception = AcquireExceptionInfo();
3997 SetGeometry(reconstruct,&geometry);
3998 geometry.x=x_offset;
3999 geometry.y=y_offset;
4000 similarity_image=CropImage(image,&geometry,sans_exception);
4001 sans_exception=DestroyExceptionInfo(sans_exception);
4002 if (similarity_image == (
Image *) NULL)
4005 status=GetImageDistortion(similarity_image,reconstruct,metric,&distortion,
4007 similarity_image=DestroyImage(similarity_image);
4008 if (status == MagickFalse)
4013MagickExport
Image *SimilarityImage(
const Image *image,
const Image *reconstruct,
4014 const MetricType metric,
const double similarity_threshold,
4017#define SimilarityImageTag "Similarity/Image"
4023 *similarity_image = (
Image *) NULL;
4037 assert(image != (
const Image *) NULL);
4038 assert(image->signature == MagickCoreSignature);
4040 assert(exception->signature == MagickCoreSignature);
4042 if (IsEventLogging() != MagickFalse)
4043 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
4044 SetGeometry(reconstruct,offset);
4045 *similarity_metric=MagickMaximumValue;
4046#if defined(MAGICKCORE_HDRI_SUPPORT) && defined(MAGICKCORE_FFTW_DELEGATE)
4048 const char *artifact = GetImageArtifact(image,
"compare:frequency-domain");
4049 if (artifact == (
const char *) NULL)
4050 artifact=GetImageArtifact(image,
"compare:accelerate-ncc");
4051 if ((image->channels & ReadMaskChannel) == 0)
4054 case DotProductCorrelationErrorMetric:
4056 similarity_image=DPCSimilarityImage(image,reconstruct,offset,
4057 similarity_metric,exception);
4058 return(similarity_image);
4060 case MeanSquaredErrorMetric:
4062 if ((artifact != (
const char *) NULL) &&
4063 (IsStringTrue(artifact) == MagickFalse))
4065 similarity_image=MSESimilarityImage(image,reconstruct,offset,
4066 similarity_metric,exception);
4067 return(similarity_image);
4069 case NormalizedCrossCorrelationErrorMetric:
4071 if ((artifact != (
const char *) NULL) &&
4072 (IsStringTrue(artifact) == MagickFalse))
4074 similarity_image=NCCSimilarityImage(image,reconstruct,offset,
4075 similarity_metric,exception);
4076 return(similarity_image);
4078 case PeakSignalToNoiseRatioErrorMetric:
4080 if ((artifact != (
const char *) NULL) &&
4081 (IsStringTrue(artifact) == MagickFalse))
4083 similarity_image=PSNRSimilarityImage(image,reconstruct,offset,
4084 similarity_metric,exception);
4085 return(similarity_image);
4087 case PhaseCorrelationErrorMetric:
4089 similarity_image=PhaseSimilarityImage(image,reconstruct,offset,
4090 similarity_metric,exception);
4091 return(similarity_image);
4093 case RootMeanSquaredErrorMetric:
4095 if ((artifact != (
const char *) NULL) &&
4096 (IsStringTrue(artifact) == MagickFalse))
4098 similarity_image=RMSESimilarityImage(image,reconstruct,offset,
4099 similarity_metric,exception);
4100 return(similarity_image);
4106 if ((metric == DotProductCorrelationErrorMetric) ||
4107 (metric == PhaseCorrelationErrorMetric))
4109 (void) ThrowMagickException(exception,GetMagickModule(),
4110 MissingDelegateError,
"DelegateLibrarySupportNotBuiltIn",
4111 "'%s' (HDRI, FFT)",image->filename);
4112 return((
Image *) NULL);
4115 if ((image->columns < reconstruct->columns) ||
4116 (image->rows < reconstruct->rows))
4118 (void) ThrowMagickException(exception,GetMagickModule(),OptionWarning,
4119 "GeometryDoesNotContainImage",
"`%s'",image->filename);
4120 return((
Image *) NULL);
4122 similarity_image=CloneImage(image,image->columns,image->rows,MagickTrue,
4124 if (similarity_image == (
Image *) NULL)
4125 return((
Image *) NULL);
4126 similarity_image->depth=MAGICKCORE_QUANTUM_DEPTH;
4127 similarity_image->alpha_trait=UndefinedPixelTrait;
4128 similarity_image->type=GrayscaleType;
4129 status=SetImageStorageClass(similarity_image,DirectClass,exception);
4130 if (status == MagickFalse)
4131 return(DestroyImage(similarity_image));
4137 similarity_view=AcquireAuthenticCacheView(similarity_image,exception);
4138 rows=similarity_image->rows;
4139#if defined(MAGICKCORE_OPENMP_SUPPORT)
4140 #pragma omp parallel for schedule(static,1) \
4141 shared(progress,similarity_metric,status) \
4142 magick_number_threads(similarity_image,similarity_image,rows << 3,1)
4144 for (y=0; y < (ssize_t) rows; y++)
4155 if (status == MagickFalse)
4157#if defined(MAGICKCORE_OPENMP_SUPPORT)
4158 #pragma omp flush(similarity_metric)
4160 if (*similarity_metric <= similarity_threshold)
4162 q=QueueCacheViewAuthenticPixels(similarity_view,0,y,
4163 similarity_image->columns,1,exception);
4164 if (q == (Quantum *) NULL)
4169 for (x=0; x < (ssize_t) similarity_image->columns; x++)
4174#if defined(MAGICKCORE_OPENMP_SUPPORT)
4175 #pragma omp flush(similarity_metric)
4177 if (*similarity_metric <= similarity_threshold)
4179 similarity=GetSimilarityMetric(image,reconstruct,metric,x,y,exception);
4180 if ((metric == DotProductCorrelationErrorMetric) ||
4181 (metric == PhaseCorrelationErrorMetric) ||
4182 (metric == NormalizedCrossCorrelationErrorMetric) ||
4183 (metric == StructuralSimilarityErrorMetric) ||
4184 (metric == UndefinedErrorMetric))
4185 similarity=1.0-similarity;
4186 if (metric == PerceptualHashErrorMetric)
4187 similarity=MagickMin(0.01*similarity,1.0);
4188#if defined(MAGICKCORE_OPENMP_SUPPORT)
4189 #pragma omp critical (MagickCore_SimilarityImage)
4191 if (similarity < *similarity_metric)
4195 *similarity_metric=similarity;
4197 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
4199 PixelChannel channel = GetPixelChannelChannel(image,i);
4200 PixelTrait traits = GetPixelChannelTraits(image,channel);
4201 PixelTrait similarity_traits = GetPixelChannelTraits(similarity_image,
4203 if ((traits == UndefinedPixelTrait) ||
4204 (similarity_traits == UndefinedPixelTrait) ||
4205 ((similarity_traits & UpdatePixelTrait) == 0))
4209 case FuzzErrorMetric:
4210 case MeanAbsoluteErrorMetric:
4211 case MeanSquaredErrorMetric:
4212 case NormalizedCrossCorrelationErrorMetric:
4213 case PerceptualHashErrorMetric:
4214 case RootMeanSquaredErrorMetric:
4215 case StructuralSimilarityErrorMetric:
4216 case StructuralDissimilarityErrorMetric:
4218 SetPixelChannel(similarity_image,channel,ClampToQuantum((
double)
4219 QuantumRange-QuantumRange*similarity),q);
4224 SetPixelChannel(similarity_image,channel,ClampToQuantum((
double)
4225 QuantumRange*similarity),q);
4230 q+=(ptrdiff_t) GetPixelChannels(similarity_image);
4232 if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
4234 if (image->progress_monitor != (MagickProgressMonitor) NULL)
4239#if defined(MAGICKCORE_OPENMP_SUPPORT)
4243 proceed=SetImageProgress(image,SimilarityImageTag,progress,image->rows);
4244 if (proceed == MagickFalse)
4248 similarity_view=DestroyCacheView(similarity_view);
4249 if (status == MagickFalse)
4250 similarity_image=DestroyImage(similarity_image);
4251 return(similarity_image);