2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
6 % CCCC OOO M M PPPP AAA RRRR EEEEE %
7 % C O O MM MM P P A A R R E %
8 % C O O M M M PPPP AAAAA RRRR EEE %
9 % C O O M M P A A R R E %
10 % CCCC OOO M M P A A R R EEEEE %
13 % MagickCore Image Comparison Methods %
20 % Copyright 1999-2019 ImageMagick Studio LLC, a non-profit organization %
21 % dedicated to making software imaging solutions freely available. %
23 % You may not use this file except in compliance with the License. You may %
24 % obtain a copy of the License at %
26 % https://imagemagick.org/script/license.php %
28 % Unless required by applicable law or agreed to in writing, software %
29 % distributed under the License is distributed on an "AS IS" BASIS, %
30 % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. %
31 % See the License for the specific language governing permissions and %
32 % limitations under the License. %
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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/exception-private.h"
57 #include "MagickCore/geometry.h"
58 #include "MagickCore/image-private.h"
59 #include "MagickCore/list.h"
60 #include "MagickCore/log.h"
61 #include "MagickCore/memory_.h"
62 #include "MagickCore/monitor.h"
63 #include "MagickCore/monitor-private.h"
64 #include "MagickCore/option.h"
65 #include "MagickCore/pixel-accessor.h"
66 #include "MagickCore/property.h"
67 #include "MagickCore/resource_.h"
68 #include "MagickCore/string_.h"
69 #include "MagickCore/statistic.h"
70 #include "MagickCore/string-private.h"
71 #include "MagickCore/thread-private.h"
72 #include "MagickCore/transform.h"
73 #include "MagickCore/utility.h"
74 #include "MagickCore/version.h"
77 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
81 % C o m p a r e I m a g e %
85 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
87 % CompareImages() compares one or more pixel channels of an image to a
88 % reconstructed image and returns the difference image.
90 % The format of the CompareImages method is:
92 % Image *CompareImages(const Image *image,const Image *reconstruct_image,
93 % const MetricType metric,double *distortion,ExceptionInfo *exception)
95 % A description of each parameter follows:
99 % o reconstruct_image: the reconstruct image.
101 % o metric: the metric.
103 % o distortion: the computed distortion between the images.
105 % o exception: return any errors or warnings in this structure.
109 static size_t GetImageChannels(const Image *image)
118 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
120 PixelChannel channel = GetPixelChannelChannel(image,i);
121 PixelTrait traits = GetPixelChannelTraits(image,channel);
122 if ((traits & UpdatePixelTrait) != 0)
125 return(channels == 0 ? (size_t) 1 : channels);
128 MagickExport Image *CompareImages(Image *image,const Image *reconstruct_image,
129 const MetricType metric,double *distortion,ExceptionInfo *exception)
165 assert(image != (Image *) NULL);
166 assert(image->signature == MagickCoreSignature);
167 if (image->debug != MagickFalse)
168 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
169 assert(reconstruct_image != (const Image *) NULL);
170 assert(reconstruct_image->signature == MagickCoreSignature);
171 assert(distortion != (double *) NULL);
173 if (image->debug != MagickFalse)
174 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
175 status=GetImageDistortion(image,reconstruct_image,metric,distortion,
177 if (status == MagickFalse)
178 return((Image *) NULL);
179 columns=MagickMax(image->columns,reconstruct_image->columns);
180 rows=MagickMax(image->rows,reconstruct_image->rows);
181 SetGeometry(image,&geometry);
182 geometry.width=columns;
183 geometry.height=rows;
184 clone_image=CloneImage(image,0,0,MagickTrue,exception);
185 if (clone_image == (Image *) NULL)
186 return((Image *) NULL);
187 (void) SetImageMask(clone_image,ReadPixelMask,(Image *) NULL,exception);
188 difference_image=ExtentImage(clone_image,&geometry,exception);
189 clone_image=DestroyImage(clone_image);
190 if (difference_image == (Image *) NULL)
191 return((Image *) NULL);
192 (void) SetImageAlphaChannel(difference_image,OpaqueAlphaChannel,exception);
193 highlight_image=CloneImage(image,columns,rows,MagickTrue,exception);
194 if (highlight_image == (Image *) NULL)
196 difference_image=DestroyImage(difference_image);
197 return((Image *) NULL);
199 status=SetImageStorageClass(highlight_image,DirectClass,exception);
200 if (status == MagickFalse)
202 difference_image=DestroyImage(difference_image);
203 highlight_image=DestroyImage(highlight_image);
204 return((Image *) NULL);
206 (void) SetImageMask(highlight_image,ReadPixelMask,(Image *) NULL,exception);
207 (void) SetImageAlphaChannel(highlight_image,OpaqueAlphaChannel,exception);
208 (void) QueryColorCompliance("#f1001ecc",AllCompliance,&highlight,exception);
209 artifact=GetImageArtifact(image,"compare:highlight-color");
210 if (artifact != (const char *) NULL)
211 (void) QueryColorCompliance(artifact,AllCompliance,&highlight,exception);
212 (void) QueryColorCompliance("#ffffffcc",AllCompliance,&lowlight,exception);
213 artifact=GetImageArtifact(image,"compare:lowlight-color");
214 if (artifact != (const char *) NULL)
215 (void) QueryColorCompliance(artifact,AllCompliance,&lowlight,exception);
216 (void) QueryColorCompliance("#888888cc",AllCompliance,&masklight,exception);
217 artifact=GetImageArtifact(image,"compare:masklight-color");
218 if (artifact != (const char *) NULL)
219 (void) QueryColorCompliance(artifact,AllCompliance,&masklight,exception);
221 Generate difference image.
224 fuzz=GetFuzzyColorDistance(image,reconstruct_image);
225 image_view=AcquireVirtualCacheView(image,exception);
226 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
227 highlight_view=AcquireAuthenticCacheView(highlight_image,exception);
228 #if defined(MAGICKCORE_OPENMP_SUPPORT)
229 #pragma omp parallel for schedule(static) shared(status) \
230 magick_number_threads(image,highlight_image,rows,1)
232 for (y=0; y < (ssize_t) rows; y++)
237 register const Quantum
247 if (status == MagickFalse)
249 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
250 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
251 r=QueueCacheViewAuthenticPixels(highlight_view,0,y,columns,1,exception);
252 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL) ||
253 (r == (Quantum *) NULL))
258 for (x=0; x < (ssize_t) columns; x++)
270 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
271 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
273 SetPixelViaPixelInfo(highlight_image,&masklight,r);
274 p+=GetPixelChannels(image);
275 q+=GetPixelChannels(reconstruct_image);
276 r+=GetPixelChannels(highlight_image);
279 difference=MagickFalse;
280 Sa=QuantumScale*GetPixelAlpha(image,p);
281 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
282 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
287 PixelChannel channel = GetPixelChannelChannel(image,i);
288 PixelTrait traits = GetPixelChannelTraits(image,channel);
289 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
291 if ((traits == UndefinedPixelTrait) ||
292 (reconstruct_traits == UndefinedPixelTrait) ||
293 ((reconstruct_traits & UpdatePixelTrait) == 0))
295 if (channel == AlphaPixelChannel)
296 distance=(double) p[i]-GetPixelChannel(reconstruct_image,channel,q);
298 distance=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
299 if ((distance*distance) > fuzz)
301 difference=MagickTrue;
305 if (difference == MagickFalse)
306 SetPixelViaPixelInfo(highlight_image,&lowlight,r);
308 SetPixelViaPixelInfo(highlight_image,&highlight,r);
309 p+=GetPixelChannels(image);
310 q+=GetPixelChannels(reconstruct_image);
311 r+=GetPixelChannels(highlight_image);
313 sync=SyncCacheViewAuthenticPixels(highlight_view,exception);
314 if (sync == MagickFalse)
317 highlight_view=DestroyCacheView(highlight_view);
318 reconstruct_view=DestroyCacheView(reconstruct_view);
319 image_view=DestroyCacheView(image_view);
320 (void) CompositeImage(difference_image,highlight_image,image->compose,
321 MagickTrue,0,0,exception);
322 highlight_image=DestroyImage(highlight_image);
323 if (status == MagickFalse)
324 difference_image=DestroyImage(difference_image);
325 return(difference_image);
329 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
333 % G e t I m a g e D i s t o r t i o n %
337 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
339 % GetImageDistortion() compares one or more pixel channels of an image to a
340 % reconstructed image and returns the specified distortion metric.
342 % The format of the GetImageDistortion method is:
344 % MagickBooleanType GetImageDistortion(const Image *image,
345 % const Image *reconstruct_image,const MetricType metric,
346 % double *distortion,ExceptionInfo *exception)
348 % A description of each parameter follows:
350 % o image: the image.
352 % o reconstruct_image: the reconstruct image.
354 % o metric: the metric.
356 % o distortion: the computed distortion between the images.
358 % o exception: return any errors or warnings in this structure.
362 static MagickBooleanType GetAbsoluteDistortion(const Image *image,
363 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
383 Compute the absolute difference in pixels between two images.
386 fuzz=(double) MagickMin(GetPixelChannels(image),
387 GetPixelChannels(reconstruct_image))*
388 GetFuzzyColorDistance(image,reconstruct_image);
389 rows=MagickMax(image->rows,reconstruct_image->rows);
390 columns=MagickMax(image->columns,reconstruct_image->columns);
391 image_view=AcquireVirtualCacheView(image,exception);
392 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
393 #if defined(MAGICKCORE_OPENMP_SUPPORT)
394 #pragma omp parallel for schedule(static) shared(status) \
395 magick_number_threads(image,image,rows,1)
397 for (y=0; y < (ssize_t) rows; y++)
400 channel_distortion[MaxPixelChannels+1];
402 register const Quantum
410 if (status == MagickFalse)
412 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
413 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
414 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
419 (void) memset(channel_distortion,0,sizeof(channel_distortion));
420 for (x=0; x < (ssize_t) columns; x++)
433 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
434 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
436 p+=GetPixelChannels(image);
437 q+=GetPixelChannels(reconstruct_image);
440 difference=MagickFalse;
442 Sa=QuantumScale*GetPixelAlpha(image,p);
443 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
444 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
449 PixelChannel channel = GetPixelChannelChannel(image,i);
450 PixelTrait traits = GetPixelChannelTraits(image,channel);
451 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
453 if ((traits == UndefinedPixelTrait) ||
454 (reconstruct_traits == UndefinedPixelTrait) ||
455 ((reconstruct_traits & UpdatePixelTrait) == 0))
457 if (channel == AlphaPixelChannel)
458 pixel=(double) p[i]-GetPixelChannel(reconstruct_image,channel,q);
460 pixel=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
461 distance+=pixel*pixel;
464 channel_distortion[i]++;
465 difference=MagickTrue;
468 if (difference != MagickFalse)
469 channel_distortion[CompositePixelChannel]++;
470 p+=GetPixelChannels(image);
471 q+=GetPixelChannels(reconstruct_image);
473 #if defined(MAGICKCORE_OPENMP_SUPPORT)
474 #pragma omp critical (MagickCore_GetAbsoluteDistortion)
476 for (j=0; j <= MaxPixelChannels; j++)
477 distortion[j]+=channel_distortion[j];
479 reconstruct_view=DestroyCacheView(reconstruct_view);
480 image_view=DestroyCacheView(image_view);
484 static MagickBooleanType GetFuzzDistortion(const Image *image,
485 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
508 rows=MagickMax(image->rows,reconstruct_image->rows);
509 columns=MagickMax(image->columns,reconstruct_image->columns);
511 image_view=AcquireVirtualCacheView(image,exception);
512 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
513 #if defined(MAGICKCORE_OPENMP_SUPPORT)
514 #pragma omp parallel for schedule(static) shared(status) \
515 magick_number_threads(image,image,rows,1) reduction(+:area)
517 for (y=0; y < (ssize_t) rows; y++)
520 channel_distortion[MaxPixelChannels+1];
522 register const Quantum
529 if (status == MagickFalse)
531 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
532 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
533 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
538 (void) memset(channel_distortion,0,sizeof(channel_distortion));
539 for (x=0; x < (ssize_t) columns; x++)
548 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
549 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
551 p+=GetPixelChannels(image);
552 q+=GetPixelChannels(reconstruct_image);
555 Sa=QuantumScale*GetPixelAlpha(image,p);
556 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
557 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
562 PixelChannel channel = GetPixelChannelChannel(image,i);
563 PixelTrait traits = GetPixelChannelTraits(image,channel);
564 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
566 if ((traits == UndefinedPixelTrait) ||
567 (reconstruct_traits == UndefinedPixelTrait) ||
568 ((reconstruct_traits & UpdatePixelTrait) == 0))
570 if (channel == AlphaPixelChannel)
571 distance=QuantumScale*(p[i]-GetPixelChannel(reconstruct_image,
574 distance=QuantumScale*(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
576 channel_distortion[i]+=distance*distance;
577 channel_distortion[CompositePixelChannel]+=distance*distance;
580 p+=GetPixelChannels(image);
581 q+=GetPixelChannels(reconstruct_image);
583 #if defined(MAGICKCORE_OPENMP_SUPPORT)
584 #pragma omp critical (MagickCore_GetFuzzDistortion)
586 for (j=0; j <= MaxPixelChannels; j++)
587 distortion[j]+=channel_distortion[j];
589 reconstruct_view=DestroyCacheView(reconstruct_view);
590 image_view=DestroyCacheView(image_view);
591 area=PerceptibleReciprocal(area);
592 for (j=0; j <= MaxPixelChannels; j++)
594 distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
595 distortion[CompositePixelChannel]=sqrt(distortion[CompositePixelChannel]);
599 static MagickBooleanType GetMeanAbsoluteDistortion(const Image *image,
600 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
623 rows=MagickMax(image->rows,reconstruct_image->rows);
624 columns=MagickMax(image->columns,reconstruct_image->columns);
626 image_view=AcquireVirtualCacheView(image,exception);
627 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
628 #if defined(MAGICKCORE_OPENMP_SUPPORT)
629 #pragma omp parallel for schedule(static) shared(status) \
630 magick_number_threads(image,image,rows,1) reduction(+:area)
632 for (y=0; y < (ssize_t) rows; y++)
635 channel_distortion[MaxPixelChannels+1];
637 register const Quantum
644 if (status == MagickFalse)
646 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
647 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
648 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
653 (void) memset(channel_distortion,0,sizeof(channel_distortion));
654 for (x=0; x < (ssize_t) columns; x++)
663 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
664 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
666 p+=GetPixelChannels(image);
667 q+=GetPixelChannels(reconstruct_image);
670 Sa=QuantumScale*GetPixelAlpha(image,p);
671 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
672 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
677 PixelChannel channel = GetPixelChannelChannel(image,i);
678 PixelTrait traits = GetPixelChannelTraits(image,channel);
679 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
681 if ((traits == UndefinedPixelTrait) ||
682 (reconstruct_traits == UndefinedPixelTrait) ||
683 ((reconstruct_traits & UpdatePixelTrait) == 0))
685 if (channel == AlphaPixelChannel)
686 distance=QuantumScale*fabs((double) p[i]-
687 GetPixelChannel(reconstruct_image,channel,q));
689 distance=QuantumScale*fabs(Sa*p[i]-Da*
690 GetPixelChannel(reconstruct_image,channel,q));
691 channel_distortion[i]+=distance;
692 channel_distortion[CompositePixelChannel]+=distance;
695 p+=GetPixelChannels(image);
696 q+=GetPixelChannels(reconstruct_image);
698 #if defined(MAGICKCORE_OPENMP_SUPPORT)
699 #pragma omp critical (MagickCore_GetMeanAbsoluteError)
701 for (j=0; j <= MaxPixelChannels; j++)
702 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);
713 static MagickBooleanType GetMeanErrorPerPixel(Image *image,
714 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
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++)
745 register const Quantum
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+=GetPixelChannels(image);
772 q+=GetPixelChannels(reconstruct_image);
775 Sa=QuantumScale*GetPixelAlpha(image,p);
776 Da=QuantumScale*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]-
792 GetPixelChannel(reconstruct_image,channel,q));
794 distance=fabs(Sa*p[i]-Da*
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+=GetPixelChannels(image);
804 q+=GetPixelChannels(reconstruct_image);
807 reconstruct_view=DestroyCacheView(reconstruct_view);
808 image_view=DestroyCacheView(image_view);
809 image->error.mean_error_per_pixel=distortion[CompositePixelChannel]/area;
810 image->error.normalized_mean_error=QuantumScale*QuantumScale*mean_error/area;
811 image->error.normalized_maximum_error=QuantumScale*maximum_error;
815 static MagickBooleanType GetMeanSquaredDistortion(const Image *image,
816 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
839 rows=MagickMax(image->rows,reconstruct_image->rows);
840 columns=MagickMax(image->columns,reconstruct_image->columns);
842 image_view=AcquireVirtualCacheView(image,exception);
843 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
844 #if defined(MAGICKCORE_OPENMP_SUPPORT)
845 #pragma omp parallel for schedule(static) shared(status) \
846 magick_number_threads(image,image,rows,1) reduction(+:area)
848 for (y=0; y < (ssize_t) rows; y++)
851 channel_distortion[MaxPixelChannels+1];
853 register const Quantum
860 if (status == MagickFalse)
862 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
863 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
864 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
869 (void) memset(channel_distortion,0,sizeof(channel_distortion));
870 for (x=0; x < (ssize_t) columns; x++)
879 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
880 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
882 p+=GetPixelChannels(image);
883 q+=GetPixelChannels(reconstruct_image);
886 Sa=QuantumScale*GetPixelAlpha(image,p);
887 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
888 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
893 PixelChannel channel = GetPixelChannelChannel(image,i);
894 PixelTrait traits = GetPixelChannelTraits(image,channel);
895 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
897 if ((traits == UndefinedPixelTrait) ||
898 (reconstruct_traits == UndefinedPixelTrait) ||
899 ((reconstruct_traits & UpdatePixelTrait) == 0))
901 if (channel == AlphaPixelChannel)
902 distance=QuantumScale*(p[i]-GetPixelChannel(reconstruct_image,
905 distance=QuantumScale*(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
907 channel_distortion[i]+=distance*distance;
908 channel_distortion[CompositePixelChannel]+=distance*distance;
911 p+=GetPixelChannels(image);
912 q+=GetPixelChannels(reconstruct_image);
914 #if defined(MAGICKCORE_OPENMP_SUPPORT)
915 #pragma omp critical (MagickCore_GetMeanSquaredError)
917 for (j=0; j <= MaxPixelChannels; j++)
918 distortion[j]+=channel_distortion[j];
920 reconstruct_view=DestroyCacheView(reconstruct_view);
921 image_view=DestroyCacheView(image_view);
922 area=PerceptibleReciprocal(area);
923 for (j=0; j <= MaxPixelChannels; j++)
925 distortion[CompositePixelChannel]/=GetImageChannels(image);
929 static MagickBooleanType GetNormalizedCrossCorrelationDistortion(
930 const Image *image,const Image *reconstruct_image,double *distortion,
931 ExceptionInfo *exception)
933 #define SimilarityImageTag "Similarity/Image"
941 *reconstruct_statistics;
963 Normalize to account for variation due to lighting and exposure condition.
965 image_statistics=GetImageStatistics(image,exception);
966 reconstruct_statistics=GetImageStatistics(reconstruct_image,exception);
967 if ((image_statistics == (ChannelStatistics *) NULL) ||
968 (reconstruct_statistics == (ChannelStatistics *) NULL))
970 if (image_statistics != (ChannelStatistics *) NULL)
971 image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
973 if (reconstruct_statistics != (ChannelStatistics *) NULL)
974 reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
975 reconstruct_statistics);
980 for (i=0; i <= MaxPixelChannels; i++)
982 rows=MagickMax(image->rows,reconstruct_image->rows);
983 columns=MagickMax(image->columns,reconstruct_image->columns);
985 image_view=AcquireVirtualCacheView(image,exception);
986 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
987 for (y=0; y < (ssize_t) rows; y++)
989 register const Quantum
996 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
997 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
998 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1003 for (x=0; x < (ssize_t) columns; x++)
1005 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1006 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1008 p+=GetPixelChannels(image);
1009 q+=GetPixelChannels(reconstruct_image);
1013 p+=GetPixelChannels(image);
1014 q+=GetPixelChannels(reconstruct_image);
1017 area=PerceptibleReciprocal(area);
1018 for (y=0; y < (ssize_t) rows; y++)
1020 register const Quantum
1027 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1028 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1029 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1034 for (x=0; x < (ssize_t) columns; x++)
1040 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1041 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1043 p+=GetPixelChannels(image);
1044 q+=GetPixelChannels(reconstruct_image);
1047 Sa=QuantumScale*GetPixelAlpha(image,p);
1048 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
1049 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1051 PixelChannel channel = GetPixelChannelChannel(image,i);
1052 PixelTrait traits = GetPixelChannelTraits(image,channel);
1053 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1055 if ((traits == UndefinedPixelTrait) ||
1056 (reconstruct_traits == UndefinedPixelTrait) ||
1057 ((reconstruct_traits & UpdatePixelTrait) == 0))
1059 if (channel == AlphaPixelChannel)
1061 distortion[i]+=area*QuantumScale*(p[i]-
1062 image_statistics[channel].mean)*(GetPixelChannel(
1063 reconstruct_image,channel,q)-
1064 reconstruct_statistics[channel].mean);
1068 distortion[i]+=area*QuantumScale*(Sa*p[i]-
1069 image_statistics[channel].mean)*(Da*GetPixelChannel(
1070 reconstruct_image,channel,q)-
1071 reconstruct_statistics[channel].mean);
1074 p+=GetPixelChannels(image);
1075 q+=GetPixelChannels(reconstruct_image);
1077 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1082 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1086 proceed=SetImageProgress(image,SimilarityImageTag,progress,rows);
1087 if (proceed == MagickFalse)
1094 reconstruct_view=DestroyCacheView(reconstruct_view);
1095 image_view=DestroyCacheView(image_view);
1097 Divide by the standard deviation.
1099 distortion[CompositePixelChannel]=0.0;
1100 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1105 PixelChannel channel = GetPixelChannelChannel(image,i);
1106 gamma=image_statistics[channel].standard_deviation*
1107 reconstruct_statistics[channel].standard_deviation;
1108 gamma=PerceptibleReciprocal(gamma);
1109 distortion[i]=QuantumRange*gamma*distortion[i];
1110 distortion[CompositePixelChannel]+=distortion[i]*distortion[i];
1112 distortion[CompositePixelChannel]=sqrt(distortion[CompositePixelChannel]/
1113 GetImageChannels(image));
1117 reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1118 reconstruct_statistics);
1119 image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1124 static MagickBooleanType GetPeakAbsoluteDistortion(const Image *image,
1125 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1142 rows=MagickMax(image->rows,reconstruct_image->rows);
1143 columns=MagickMax(image->columns,reconstruct_image->columns);
1144 image_view=AcquireVirtualCacheView(image,exception);
1145 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1146 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1147 #pragma omp parallel for schedule(static) shared(status) \
1148 magick_number_threads(image,image,rows,1)
1150 for (y=0; y < (ssize_t) rows; y++)
1153 channel_distortion[MaxPixelChannels+1];
1155 register const Quantum
1163 if (status == MagickFalse)
1165 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1166 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1167 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1172 (void) memset(channel_distortion,0,sizeof(channel_distortion));
1173 for (x=0; x < (ssize_t) columns; x++)
1182 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1183 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1185 p+=GetPixelChannels(image);
1186 q+=GetPixelChannels(reconstruct_image);
1189 Sa=QuantumScale*GetPixelAlpha(image,p);
1190 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
1191 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1196 PixelChannel channel = GetPixelChannelChannel(image,i);
1197 PixelTrait traits = GetPixelChannelTraits(image,channel);
1198 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1200 if ((traits == UndefinedPixelTrait) ||
1201 (reconstruct_traits == UndefinedPixelTrait) ||
1202 ((reconstruct_traits & UpdatePixelTrait) == 0))
1204 if (channel == AlphaPixelChannel)
1205 distance=QuantumScale*fabs((double) p[i]-
1206 GetPixelChannel(reconstruct_image,channel,q));
1208 distance=QuantumScale*fabs(Sa*p[i]-Da*
1209 GetPixelChannel(reconstruct_image,channel,q));
1210 if (distance > channel_distortion[i])
1211 channel_distortion[i]=distance;
1212 if (distance > channel_distortion[CompositePixelChannel])
1213 channel_distortion[CompositePixelChannel]=distance;
1215 p+=GetPixelChannels(image);
1216 q+=GetPixelChannels(reconstruct_image);
1218 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1219 #pragma omp critical (MagickCore_GetPeakAbsoluteError)
1221 for (j=0; j <= MaxPixelChannels; j++)
1222 if (channel_distortion[j] > distortion[j])
1223 distortion[j]=channel_distortion[j];
1225 reconstruct_view=DestroyCacheView(reconstruct_view);
1226 image_view=DestroyCacheView(image_view);
1230 static inline double MagickLog10(const double x)
1232 #define Log10Epsilon (1.0e-11)
1234 if (fabs(x) < Log10Epsilon)
1235 return(log10(Log10Epsilon));
1236 return(log10(fabs(x)));
1239 static MagickBooleanType GetPeakSignalToNoiseRatio(const Image *image,
1240 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1248 status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1249 for (i=0; i <= MaxPixelChannels; i++)
1250 if (fabs(distortion[i]) < MagickEpsilon)
1251 distortion[i]=INFINITY;
1253 distortion[i]=10.0*MagickLog10(1.0)-10.0*MagickLog10(distortion[i]);
1257 static MagickBooleanType GetPerceptualHashDistortion(const Image *image,
1258 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1260 ChannelPerceptualHash
1274 Compute perceptual hash in the sRGB colorspace.
1276 channel_phash=GetImagePerceptualHash(image,exception);
1277 if (channel_phash == (ChannelPerceptualHash *) NULL)
1278 return(MagickFalse);
1279 reconstruct_phash=GetImagePerceptualHash(reconstruct_image,exception);
1280 if (reconstruct_phash == (ChannelPerceptualHash *) NULL)
1282 channel_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1284 return(MagickFalse);
1286 artifact=GetImageArtifact(image,"phash:normalize");
1287 normalize=(artifact == (const char *) NULL) ||
1288 (IsStringTrue(artifact) == MagickFalse) ? MagickFalse : MagickTrue;
1289 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1290 #pragma omp parallel for schedule(static)
1292 for (channel=0; channel < MaxPixelChannels; channel++)
1301 for (i=0; i < MaximumNumberOfImageMoments; i++)
1310 for (j=0; j < (ssize_t) channel_phash[0].number_colorspaces; j++)
1312 alpha=channel_phash[channel].phash[j][i];
1313 beta=reconstruct_phash[channel].phash[j][i];
1314 if (normalize == MagickFalse)
1315 difference+=(beta-alpha)*(beta-alpha);
1317 difference=sqrt((beta-alpha)*(beta-alpha)/
1318 channel_phash[0].number_channels);
1321 distortion[channel]+=difference;
1322 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1323 #pragma omp critical (MagickCore_GetPerceptualHashDistortion)
1325 distortion[CompositePixelChannel]+=difference;
1330 reconstruct_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1332 channel_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(channel_phash);
1336 static MagickBooleanType GetRootMeanSquaredDistortion(const Image *image,
1337 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1345 status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1346 for (i=0; i <= MaxPixelChannels; i++)
1347 distortion[i]=sqrt(distortion[i]);
1351 static MagickBooleanType GetStructuralSimilarityDistortion(const Image *image,
1352 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1354 #define SSIMRadius 5.0
1355 #define SSIMSigma 1.5
1356 #define SSIMBlocksize 8
1366 geometry[MagickPathExtent];
1394 Compute structural similarity index @
1395 https://en.wikipedia.org/wiki/Structural_similarity.
1398 artifact=GetImageArtifact(image,"compare:ssim-radius");
1399 if (artifact != (const char *) NULL)
1400 radius=StringToDouble(artifact,(char **) NULL);
1402 artifact=GetImageArtifact(image,"compare:ssim-sigma");
1403 if (artifact != (const char *) NULL)
1404 sigma=StringToDouble(artifact,(char **) NULL);
1405 (void) FormatLocaleString(geometry,MagickPathExtent,"gaussian:%.20gx%.20g",
1407 kernel_info=AcquireKernelInfo(geometry,exception);
1408 if (kernel_info == (KernelInfo *) NULL)
1409 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
1411 c1=pow(SSIMK1*SSIML,2.0);
1412 artifact=GetImageArtifact(image,"compare:ssim-k1");
1413 if (artifact != (const char *) NULL)
1414 c1=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
1415 c2=pow(SSIMK2*SSIML,2.0);
1416 artifact=GetImageArtifact(image,"compare:ssim-k2");
1417 if (artifact != (const char *) NULL)
1418 c2=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
1420 rows=MagickMax(image->rows,reconstruct_image->rows);
1421 columns=MagickMax(image->columns,reconstruct_image->columns);
1422 image_view=AcquireVirtualCacheView(image,exception);
1423 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1424 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1425 #pragma omp parallel for schedule(static) shared(status) \
1426 magick_number_threads(image,reconstruct_image,rows,1)
1428 for (y=0; y < (ssize_t) rows; y++)
1431 channel_distortion[MaxPixelChannels+1];
1433 register const Quantum
1441 if (status == MagickFalse)
1443 p=GetCacheViewVirtualPixels(image_view,-((ssize_t) kernel_info->width/2L),y-
1444 ((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
1445 kernel_info->height,exception);
1446 q=GetCacheViewVirtualPixels(reconstruct_view,-((ssize_t) kernel_info->width/
1447 2L),y-((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
1448 kernel_info->height,exception);
1449 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1454 (void) memset(channel_distortion,0,sizeof(channel_distortion));
1455 for (x=0; x < (ssize_t) columns; x++)
1458 x_pixel_mu[MaxPixelChannels+1],
1459 x_pixel_sigma_squared[MaxPixelChannels+1],
1460 xy_sigma[MaxPixelChannels+1],
1461 y_pixel_mu[MaxPixelChannels+1],
1462 y_pixel_sigma_squared[MaxPixelChannels+1];
1464 register const Quantum
1465 *magick_restrict reference,
1466 *magick_restrict target;
1468 register MagickRealType
1474 (void) memset(x_pixel_mu,0,sizeof(x_pixel_mu));
1475 (void) memset(x_pixel_sigma_squared,0,sizeof(x_pixel_sigma_squared));
1476 (void) memset(xy_sigma,0,sizeof(xy_sigma));
1477 (void) memset(x_pixel_sigma_squared,0,sizeof(y_pixel_sigma_squared));
1478 (void) memset(y_pixel_mu,0,sizeof(y_pixel_mu));
1479 (void) memset(y_pixel_sigma_squared,0,sizeof(y_pixel_sigma_squared));
1480 k=kernel_info->values;
1483 for (v=0; v < (ssize_t) kernel_info->height; v++)
1488 for (u=0; u < (ssize_t) kernel_info->width; u++)
1490 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1496 PixelChannel channel = GetPixelChannelChannel(image,i);
1497 PixelTrait traits = GetPixelChannelTraits(image,channel);
1498 PixelTrait reconstruct_traits = GetPixelChannelTraits(
1499 reconstruct_image,channel);
1500 if ((traits == UndefinedPixelTrait) ||
1501 (reconstruct_traits == UndefinedPixelTrait) ||
1502 ((reconstruct_traits & UpdatePixelTrait) == 0))
1504 x_pixel=QuantumScale*reference[i];
1505 x_pixel_mu[i]+=(*k)*x_pixel;
1506 x_pixel_sigma_squared[i]+=(*k)*x_pixel*x_pixel;
1507 y_pixel=QuantumScale*
1508 GetPixelChannel(reconstruct_image,channel,target);
1509 y_pixel_mu[i]+=(*k)*y_pixel;
1510 y_pixel_sigma_squared[i]+=(*k)*y_pixel*y_pixel;
1511 xy_sigma[i]+=(*k)*x_pixel*y_pixel;
1514 reference+=GetPixelChannels(image);
1515 target+=GetPixelChannels(reconstruct_image);
1517 reference+=GetPixelChannels(image)*columns;
1518 target+=GetPixelChannels(reconstruct_image)*columns;
1520 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1525 x_pixel_sigmas_squared,
1529 y_pixel_sigmas_squared;
1531 PixelChannel channel = GetPixelChannelChannel(image,i);
1532 PixelTrait traits = GetPixelChannelTraits(image,channel);
1533 PixelTrait reconstruct_traits = GetPixelChannelTraits(
1534 reconstruct_image,channel);
1535 if ((traits == UndefinedPixelTrait) ||
1536 (reconstruct_traits == UndefinedPixelTrait) ||
1537 ((reconstruct_traits & UpdatePixelTrait) == 0))
1539 x_pixel_mu_squared=x_pixel_mu[i]*x_pixel_mu[i];
1540 y_pixel_mu_squared=y_pixel_mu[i]*y_pixel_mu[i];
1541 xy_mu=x_pixel_mu[i]*y_pixel_mu[i];
1542 xy_sigmas=xy_sigma[i]-xy_mu;
1543 x_pixel_sigmas_squared=x_pixel_sigma_squared[i]-x_pixel_mu_squared;
1544 y_pixel_sigmas_squared=y_pixel_sigma_squared[i]-y_pixel_mu_squared;
1545 ssim=((2.0*xy_mu+c1)*(2.0*xy_sigmas+c2))/
1546 ((x_pixel_mu_squared+y_pixel_mu_squared+c1)*
1547 (x_pixel_sigmas_squared+y_pixel_sigmas_squared+c2));
1548 channel_distortion[i]+=ssim;
1549 channel_distortion[CompositePixelChannel]+=ssim;
1551 p+=GetPixelChannels(image);
1552 q+=GetPixelChannels(reconstruct_image);
1554 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1555 #pragma omp critical (MagickCore_GetStructuralSimilarityDistortion)
1557 for (i=0; i <= MaxPixelChannels; i++)
1558 distortion[i]+=channel_distortion[i];
1560 image_view=DestroyCacheView(image_view);
1561 reconstruct_view=DestroyCacheView(reconstruct_view);
1562 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1564 PixelChannel channel = GetPixelChannelChannel(image,i);
1565 PixelTrait traits = GetPixelChannelTraits(image,channel);
1566 if ((traits == UndefinedPixelTrait) || ((traits & UpdatePixelTrait) == 0))
1568 distortion[i]/=((double) columns*rows);
1570 distortion[CompositePixelChannel]/=((double) columns*rows);
1571 distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
1572 kernel_info=DestroyKernelInfo(kernel_info);
1576 static MagickBooleanType GetStructuralDisimilarityDistortion(const Image *image,
1577 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1585 status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1586 distortion,exception);
1587 for (i=0; i <= MaxPixelChannels; i++)
1588 distortion[i]=(1.0-(distortion[i]))/2.0;
1592 MagickExport MagickBooleanType GetImageDistortion(Image *image,
1593 const Image *reconstruct_image,const MetricType metric,double *distortion,
1594 ExceptionInfo *exception)
1597 *channel_distortion;
1605 assert(image != (Image *) NULL);
1606 assert(image->signature == MagickCoreSignature);
1607 if (image->debug != MagickFalse)
1608 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1609 assert(reconstruct_image != (const Image *) NULL);
1610 assert(reconstruct_image->signature == MagickCoreSignature);
1611 assert(distortion != (double *) NULL);
1613 if (image->debug != MagickFalse)
1614 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1616 Get image distortion.
1618 length=MaxPixelChannels+1;
1619 channel_distortion=(double *) AcquireQuantumMemory(length,
1620 sizeof(*channel_distortion));
1621 if (channel_distortion == (double *) NULL)
1622 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1623 (void) memset(channel_distortion,0,length*
1624 sizeof(*channel_distortion));
1627 case AbsoluteErrorMetric:
1629 status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1633 case FuzzErrorMetric:
1635 status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1639 case MeanAbsoluteErrorMetric:
1641 status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1642 channel_distortion,exception);
1645 case MeanErrorPerPixelErrorMetric:
1647 status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1651 case MeanSquaredErrorMetric:
1653 status=GetMeanSquaredDistortion(image,reconstruct_image,
1654 channel_distortion,exception);
1657 case NormalizedCrossCorrelationErrorMetric:
1660 status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1661 channel_distortion,exception);
1664 case PeakAbsoluteErrorMetric:
1666 status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1667 channel_distortion,exception);
1670 case PeakSignalToNoiseRatioErrorMetric:
1672 status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1673 channel_distortion,exception);
1676 case PerceptualHashErrorMetric:
1678 status=GetPerceptualHashDistortion(image,reconstruct_image,
1679 channel_distortion,exception);
1682 case RootMeanSquaredErrorMetric:
1684 status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1685 channel_distortion,exception);
1688 case StructuralSimilarityErrorMetric:
1690 status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1691 channel_distortion,exception);
1694 case StructuralDissimilarityErrorMetric:
1696 status=GetStructuralDisimilarityDistortion(image,reconstruct_image,
1697 channel_distortion,exception);
1701 *distortion=channel_distortion[CompositePixelChannel];
1702 channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1703 (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
1709 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1713 % G e t I m a g e D i s t o r t i o n s %
1717 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1719 % GetImageDistortions() compares the pixel channels of an image to a
1720 % reconstructed image and returns the specified distortion metric for each
1723 % The format of the GetImageDistortions method is:
1725 % double *GetImageDistortions(const Image *image,
1726 % const Image *reconstruct_image,const MetricType metric,
1727 % ExceptionInfo *exception)
1729 % A description of each parameter follows:
1731 % o image: the image.
1733 % o reconstruct_image: the reconstruct image.
1735 % o metric: the metric.
1737 % o exception: return any errors or warnings in this structure.
1740 MagickExport double *GetImageDistortions(Image *image,
1741 const Image *reconstruct_image,const MetricType metric,
1742 ExceptionInfo *exception)
1745 *channel_distortion;
1753 assert(image != (Image *) NULL);
1754 assert(image->signature == MagickCoreSignature);
1755 if (image->debug != MagickFalse)
1756 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1757 assert(reconstruct_image != (const Image *) NULL);
1758 assert(reconstruct_image->signature == MagickCoreSignature);
1759 if (image->debug != MagickFalse)
1760 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1762 Get image distortion.
1764 length=MaxPixelChannels+1UL;
1765 channel_distortion=(double *) AcquireQuantumMemory(length,
1766 sizeof(*channel_distortion));
1767 if (channel_distortion == (double *) NULL)
1768 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1769 (void) memset(channel_distortion,0,length*
1770 sizeof(*channel_distortion));
1774 case AbsoluteErrorMetric:
1776 status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1780 case FuzzErrorMetric:
1782 status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1786 case MeanAbsoluteErrorMetric:
1788 status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1789 channel_distortion,exception);
1792 case MeanErrorPerPixelErrorMetric:
1794 status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1798 case MeanSquaredErrorMetric:
1800 status=GetMeanSquaredDistortion(image,reconstruct_image,
1801 channel_distortion,exception);
1804 case NormalizedCrossCorrelationErrorMetric:
1807 status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1808 channel_distortion,exception);
1811 case PeakAbsoluteErrorMetric:
1813 status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1814 channel_distortion,exception);
1817 case PeakSignalToNoiseRatioErrorMetric:
1819 status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1820 channel_distortion,exception);
1823 case PerceptualHashErrorMetric:
1825 status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1826 channel_distortion,exception);
1829 case RootMeanSquaredErrorMetric:
1831 status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1832 channel_distortion,exception);
1835 case StructuralSimilarityErrorMetric:
1837 status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1838 channel_distortion,exception);
1841 case StructuralDissimilarityErrorMetric:
1843 status=GetStructuralDisimilarityDistortion(image,reconstruct_image,
1844 channel_distortion,exception);
1848 if (status == MagickFalse)
1850 channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1851 return((double *) NULL);
1853 return(channel_distortion);
1857 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1861 % I s I m a g e s E q u a l %
1865 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1867 % IsImagesEqual() compare the pixels of two images and returns immediately
1868 % if any pixel is not identical.
1870 % The format of the IsImagesEqual method is:
1872 % MagickBooleanType IsImagesEqual(const Image *image,
1873 % const Image *reconstruct_image,ExceptionInfo *exception)
1875 % A description of each parameter follows.
1877 % o image: the image.
1879 % o reconstruct_image: the reconstruct image.
1881 % o exception: return any errors or warnings in this structure.
1884 MagickExport MagickBooleanType IsImagesEqual(const Image *image,
1885 const Image *reconstruct_image,ExceptionInfo *exception)
1898 assert(image != (Image *) NULL);
1899 assert(image->signature == MagickCoreSignature);
1900 assert(reconstruct_image != (const Image *) NULL);
1901 assert(reconstruct_image->signature == MagickCoreSignature);
1902 rows=MagickMax(image->rows,reconstruct_image->rows);
1903 columns=MagickMax(image->columns,reconstruct_image->columns);
1904 image_view=AcquireVirtualCacheView(image,exception);
1905 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1906 for (y=0; y < (ssize_t) rows; y++)
1908 register const Quantum
1915 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1916 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1917 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
1919 for (x=0; x < (ssize_t) columns; x++)
1924 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1929 PixelChannel channel = GetPixelChannelChannel(image,i);
1930 PixelTrait traits = GetPixelChannelTraits(image,channel);
1931 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1933 if ((traits == UndefinedPixelTrait) ||
1934 (reconstruct_traits == UndefinedPixelTrait) ||
1935 ((reconstruct_traits & UpdatePixelTrait) == 0))
1937 distance=fabs(p[i]-(double) GetPixelChannel(reconstruct_image,
1939 if (distance >= MagickEpsilon)
1942 if (i < (ssize_t) GetPixelChannels(image))
1944 p+=GetPixelChannels(image);
1945 q+=GetPixelChannels(reconstruct_image);
1947 if (x < (ssize_t) columns)
1950 reconstruct_view=DestroyCacheView(reconstruct_view);
1951 image_view=DestroyCacheView(image_view);
1952 return(y < (ssize_t) rows ? MagickFalse : MagickTrue);
1956 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1960 % S e t I m a g e C o l o r M e t r i c %
1964 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1966 % SetImageColorMetric() measures the difference between colors at each pixel
1967 % location of two images. A value other than 0 means the colors match
1968 % exactly. Otherwise an error measure is computed by summing over all
1969 % pixels in an image the distance squared in RGB space between each image
1970 % pixel and its corresponding pixel in the reconstruct image. The error
1971 % measure is assigned to these image members:
1973 % o mean_error_per_pixel: The mean error for any single pixel in
1976 % o normalized_mean_error: The normalized mean quantization error for
1977 % any single pixel in the image. This distance measure is normalized to
1978 % a range between 0 and 1. It is independent of the range of red, green,
1979 % and blue values in the image.
1981 % o normalized_maximum_error: The normalized maximum quantization
1982 % error for any single pixel in the image. This distance measure is
1983 % normalized to a range between 0 and 1. It is independent of the range
1984 % of red, green, and blue values in your image.
1986 % A small normalized mean square error, accessed as
1987 % image->normalized_mean_error, suggests the images are very similar in
1988 % spatial layout and color.
1990 % The format of the SetImageColorMetric method is:
1992 % MagickBooleanType SetImageColorMetric(Image *image,
1993 % const Image *reconstruct_image,ExceptionInfo *exception)
1995 % A description of each parameter follows.
1997 % o image: the image.
1999 % o reconstruct_image: the reconstruct image.
2001 % o exception: return any errors or warnings in this structure.
2004 MagickExport MagickBooleanType SetImageColorMetric(Image *image,
2005 const Image *reconstruct_image,ExceptionInfo *exception)
2015 mean_error_per_pixel;
2027 assert(image != (Image *) NULL);
2028 assert(image->signature == MagickCoreSignature);
2029 assert(reconstruct_image != (const Image *) NULL);
2030 assert(reconstruct_image->signature == MagickCoreSignature);
2033 mean_error_per_pixel=0.0;
2035 rows=MagickMax(image->rows,reconstruct_image->rows);
2036 columns=MagickMax(image->columns,reconstruct_image->columns);
2037 image_view=AcquireVirtualCacheView(image,exception);
2038 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
2039 for (y=0; y < (ssize_t) rows; y++)
2041 register const Quantum
2048 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
2049 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
2050 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2052 for (x=0; x < (ssize_t) columns; x++)
2057 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2062 PixelChannel channel = GetPixelChannelChannel(image,i);
2063 PixelTrait traits = GetPixelChannelTraits(image,channel);
2064 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2066 if ((traits == UndefinedPixelTrait) ||
2067 (reconstruct_traits == UndefinedPixelTrait) ||
2068 ((reconstruct_traits & UpdatePixelTrait) == 0))
2070 distance=fabs(p[i]-(double) GetPixelChannel(reconstruct_image,
2072 if (distance >= MagickEpsilon)
2074 mean_error_per_pixel+=distance;
2075 mean_error+=distance*distance;
2076 if (distance > maximum_error)
2077 maximum_error=distance;
2081 p+=GetPixelChannels(image);
2082 q+=GetPixelChannels(reconstruct_image);
2085 reconstruct_view=DestroyCacheView(reconstruct_view);
2086 image_view=DestroyCacheView(image_view);
2087 image->error.mean_error_per_pixel=(double) (mean_error_per_pixel/area);
2088 image->error.normalized_mean_error=(double) (QuantumScale*QuantumScale*
2090 image->error.normalized_maximum_error=(double) (QuantumScale*maximum_error);
2091 status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
2096 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2100 % S i m i l a r i t y I m a g e %
2104 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2106 % SimilarityImage() compares the reference image of the image and returns the
2107 % best match offset. In addition, it returns a similarity image such that an
2108 % exact match location is completely white and if none of the pixels match,
2109 % black, otherwise some gray level in-between.
2111 % The format of the SimilarityImageImage method is:
2113 % Image *SimilarityImage(const Image *image,const Image *reference,
2114 % const MetricType metric,const double similarity_threshold,
2115 % RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
2117 % A description of each parameter follows:
2119 % o image: the image.
2121 % o reference: find an area of the image that closely resembles this image.
2123 % o metric: the metric.
2125 % o similarity_threshold: minimum distortion for (sub)image match.
2127 % o offset: the best match offset of the reference image within the image.
2129 % o similarity: the computed similarity between the images.
2131 % o exception: return any errors or warnings in this structure.
2135 static double GetSimilarityMetric(const Image *image,const Image *reference,
2136 const MetricType metric,const ssize_t x_offset,const ssize_t y_offset,
2137 ExceptionInfo *exception)
2151 SetGeometry(reference,&geometry);
2152 geometry.x=x_offset;
2153 geometry.y=y_offset;
2154 similarity_image=CropImage(image,&geometry,exception);
2155 if (similarity_image == (Image *) NULL)
2158 status=GetImageDistortion(similarity_image,reference,metric,&distortion,
2160 similarity_image=DestroyImage(similarity_image);
2161 if (status == MagickFalse)
2166 MagickExport Image *SimilarityImage(const Image *image,const Image *reference,
2167 const MetricType metric,const double similarity_threshold,
2168 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
2170 #define SimilarityImageTag "Similarity/Image"
2187 assert(image != (const Image *) NULL);
2188 assert(image->signature == MagickCoreSignature);
2189 if (image->debug != MagickFalse)
2190 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2191 assert(exception != (ExceptionInfo *) NULL);
2192 assert(exception->signature == MagickCoreSignature);
2193 assert(offset != (RectangleInfo *) NULL);
2194 SetGeometry(reference,offset);
2195 *similarity_metric=MagickMaximumValue;
2196 similarity_image=CloneImage(image,image->columns-reference->columns+1,
2197 image->rows-reference->rows+1,MagickTrue,exception);
2198 if (similarity_image == (Image *) NULL)
2199 return((Image *) NULL);
2200 status=SetImageStorageClass(similarity_image,DirectClass,exception);
2201 if (status == MagickFalse)
2203 similarity_image=DestroyImage(similarity_image);
2204 return((Image *) NULL);
2206 (void) SetImageAlphaChannel(similarity_image,DeactivateAlphaChannel,
2209 Measure similarity of reference image against image.
2213 similarity_view=AcquireAuthenticCacheView(similarity_image,exception);
2214 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2215 #pragma omp parallel for schedule(static) \
2216 shared(progress,status,similarity_metric) \
2217 magick_number_threads(image,image,image->rows-reference->rows+1,1)
2219 for (y=0; y < (ssize_t) (image->rows-reference->rows+1); y++)
2230 if (status == MagickFalse)
2232 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2233 #pragma omp flush(similarity_metric)
2235 if (*similarity_metric <= similarity_threshold)
2237 q=GetCacheViewAuthenticPixels(similarity_view,0,y,similarity_image->columns,
2239 if (q == (Quantum *) NULL)
2244 for (x=0; x < (ssize_t) (image->columns-reference->columns+1); x++)
2249 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2250 #pragma omp flush(similarity_metric)
2252 if (*similarity_metric <= similarity_threshold)
2254 similarity=GetSimilarityMetric(image,reference,metric,x,y,exception);
2255 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2256 #pragma omp critical (MagickCore_SimilarityImage)
2258 if ((metric == NormalizedCrossCorrelationErrorMetric) ||
2259 (metric == UndefinedErrorMetric))
2260 similarity=1.0-similarity;
2261 if (similarity < *similarity_metric)
2265 *similarity_metric=similarity;
2267 if (metric == PerceptualHashErrorMetric)
2268 similarity=MagickMin(0.01*similarity,1.0);
2269 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2271 PixelChannel channel = GetPixelChannelChannel(image,i);
2272 PixelTrait traits = GetPixelChannelTraits(image,channel);
2273 PixelTrait similarity_traits=GetPixelChannelTraits(similarity_image,
2275 if ((traits == UndefinedPixelTrait) ||
2276 (similarity_traits == UndefinedPixelTrait) ||
2277 ((similarity_traits & UpdatePixelTrait) == 0))
2279 SetPixelChannel(similarity_image,channel,ClampToQuantum(QuantumRange-
2280 QuantumRange*similarity),q);
2282 q+=GetPixelChannels(similarity_image);
2284 if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
2286 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2291 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2295 proceed=SetImageProgress(image,SimilarityImageTag,progress,image->rows);
2296 if (proceed == MagickFalse)
2300 similarity_view=DestroyCacheView(similarity_view);
2301 if (status == MagickFalse)
2302 similarity_image=DestroyImage(similarity_image);
2303 return(similarity_image);