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-2013 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 % http://www.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/cache-view.h"
46 #include "MagickCore/client.h"
47 #include "MagickCore/color.h"
48 #include "MagickCore/color-private.h"
49 #include "MagickCore/colorspace.h"
50 #include "MagickCore/colorspace-private.h"
51 #include "MagickCore/compare.h"
52 #include "MagickCore/composite-private.h"
53 #include "MagickCore/constitute.h"
54 #include "MagickCore/exception-private.h"
55 #include "MagickCore/geometry.h"
56 #include "MagickCore/image-private.h"
57 #include "MagickCore/list.h"
58 #include "MagickCore/log.h"
59 #include "MagickCore/memory_.h"
60 #include "MagickCore/monitor.h"
61 #include "MagickCore/monitor-private.h"
62 #include "MagickCore/option.h"
63 #include "MagickCore/pixel-accessor.h"
64 #include "MagickCore/resource_.h"
65 #include "MagickCore/string_.h"
66 #include "MagickCore/statistic.h"
67 #include "MagickCore/thread-private.h"
68 #include "MagickCore/transform.h"
69 #include "MagickCore/utility.h"
70 #include "MagickCore/version.h"
73 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
77 % C o m p a r e I m a g e %
81 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
83 % CompareImages() compares one or more pixel channels of an image to a
84 % reconstructed image and returns the difference image.
86 % The format of the CompareImages method is:
88 % Image *CompareImages(const Image *image,const Image *reconstruct_image,
89 % const MetricType metric,double *distortion,ExceptionInfo *exception)
91 % A description of each parameter follows:
95 % o reconstruct_image: the reconstruct image.
97 % o metric: the metric.
99 % o distortion: the computed distortion between the images.
101 % o exception: return any errors or warnings in this structure.
104 MagickExport Image *CompareImages(Image *image,const Image *reconstruct_image,
105 const MetricType metric,double *distortion,ExceptionInfo *exception)
129 assert(image != (Image *) NULL);
130 assert(image->signature == MagickSignature);
131 if (image->debug != MagickFalse)
132 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
133 assert(reconstruct_image != (const Image *) NULL);
134 assert(reconstruct_image->signature == MagickSignature);
135 assert(distortion != (double *) NULL);
137 if (image->debug != MagickFalse)
138 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
139 if ((reconstruct_image->columns != image->columns) ||
140 (reconstruct_image->rows != image->rows))
141 ThrowImageException(ImageError,"ImageSizeDiffers");
142 status=GetImageDistortion(image,reconstruct_image,metric,distortion,
144 if (status == MagickFalse)
145 return((Image *) NULL);
146 difference_image=CloneImage(image,0,0,MagickTrue,exception);
147 if (difference_image == (Image *) NULL)
148 return((Image *) NULL);
149 (void) SetImageAlphaChannel(difference_image,OpaqueAlphaChannel,exception);
150 highlight_image=CloneImage(image,image->columns,image->rows,MagickTrue,
152 if (highlight_image == (Image *) NULL)
154 difference_image=DestroyImage(difference_image);
155 return((Image *) NULL);
157 status=SetImageStorageClass(highlight_image,DirectClass,exception);
158 if (status == MagickFalse)
160 difference_image=DestroyImage(difference_image);
161 highlight_image=DestroyImage(highlight_image);
162 return((Image *) NULL);
164 (void) SetImageAlphaChannel(highlight_image,OpaqueAlphaChannel,exception);
165 (void) QueryColorCompliance("#f1001ecc",AllCompliance,&highlight,exception);
166 artifact=GetImageArtifact(image,"highlight-color");
167 if (artifact != (const char *) NULL)
168 (void) QueryColorCompliance(artifact,AllCompliance,&highlight,exception);
169 (void) QueryColorCompliance("#ffffffcc",AllCompliance,&lowlight,exception);
170 artifact=GetImageArtifact(image,"lowlight-color");
171 if (artifact != (const char *) NULL)
172 (void) QueryColorCompliance(artifact,AllCompliance,&lowlight,exception);
174 Generate difference image.
177 image_view=AcquireVirtualCacheView(image,exception);
178 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
179 highlight_view=AcquireAuthenticCacheView(highlight_image,exception);
180 #if defined(MAGICKCORE_OPENMP_SUPPORT)
181 #pragma omp parallel for schedule(static,4) shared(status) \
182 magick_threads(image,highlight_image,image->rows,1)
184 for (y=0; y < (ssize_t) image->rows; y++)
189 register const Quantum
199 if (status == MagickFalse)
201 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
202 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
204 r=QueueCacheViewAuthenticPixels(highlight_view,0,y,highlight_image->columns,
206 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL) ||
207 (r == (Quantum *) NULL))
212 for (x=0; x < (ssize_t) image->columns; x++)
224 if (GetPixelReadMask(image,p) == 0)
226 SetPixelInfoPixel(highlight_image,&lowlight,r);
227 p+=GetPixelChannels(image);
228 q+=GetPixelChannels(reconstruct_image);
229 r+=GetPixelChannels(highlight_image);
232 difference=MagickFalse;
233 Sa=QuantumScale*GetPixelAlpha(image,p);
234 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
235 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
240 PixelChannel channel=GetPixelChannelChannel(image,i);
241 PixelTrait traits=GetPixelChannelTraits(image,channel);
242 PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
244 if ((traits == UndefinedPixelTrait) ||
245 (reconstruct_traits == UndefinedPixelTrait) ||
246 ((reconstruct_traits & UpdatePixelTrait) == 0))
248 distance=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
249 if (fabs((double) distance) >= MagickEpsilon)
250 difference=MagickTrue;
252 if (difference == MagickFalse)
253 SetPixelInfoPixel(highlight_image,&lowlight,r);
255 SetPixelInfoPixel(highlight_image,&highlight,r);
256 p+=GetPixelChannels(image);
257 q+=GetPixelChannels(reconstruct_image);
258 r+=GetPixelChannels(highlight_image);
260 sync=SyncCacheViewAuthenticPixels(highlight_view,exception);
261 if (sync == MagickFalse)
264 highlight_view=DestroyCacheView(highlight_view);
265 reconstruct_view=DestroyCacheView(reconstruct_view);
266 image_view=DestroyCacheView(image_view);
267 (void) CompositeImage(difference_image,highlight_image,image->compose,
268 MagickTrue,0,0,exception);
269 highlight_image=DestroyImage(highlight_image);
270 if (status == MagickFalse)
271 difference_image=DestroyImage(difference_image);
272 return(difference_image);
276 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
280 % G e t I m a g e D i s t o r t i o n %
284 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
286 % GetImageDistortion() compares one or more pixel channels of an image to a
287 % reconstructed image and returns the specified distortion metric.
289 % The format of the GetImageDistortion method is:
291 % MagickBooleanType GetImageDistortion(const Image *image,
292 % const Image *reconstruct_image,const MetricType metric,
293 % double *distortion,ExceptionInfo *exception)
295 % A description of each parameter follows:
297 % o image: the image.
299 % o reconstruct_image: the reconstruct image.
301 % o metric: the metric.
303 % o distortion: the computed distortion between the images.
305 % o exception: return any errors or warnings in this structure.
309 static inline double MagickMax(const double x,const double y)
316 static MagickBooleanType GetAbsoluteDistortion(const Image *image,
317 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
333 Compute the absolute difference in pixels between two images.
336 if (image->fuzz == 0.0)
337 fuzz=MagickMax(reconstruct_image->fuzz,MagickSQ1_2)*
338 MagickMax(reconstruct_image->fuzz,MagickSQ1_2);
340 if (reconstruct_image->fuzz == 0.0)
341 fuzz=MagickMax(image->fuzz,MagickSQ1_2)*
342 MagickMax(image->fuzz,MagickSQ1_2);
344 fuzz=MagickMax(image->fuzz,MagickSQ1_2)*
345 MagickMax(reconstruct_image->fuzz,MagickSQ1_2);
346 image_view=AcquireVirtualCacheView(image,exception);
347 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
348 #if defined(MAGICKCORE_OPENMP_SUPPORT)
349 #pragma omp parallel for schedule(static,4) shared(status) \
350 magick_threads(image,image,image->rows,1)
352 for (y=0; y < (ssize_t) image->rows; y++)
355 channel_distortion[MaxPixelChannels+1];
357 register const Quantum
365 if (status == MagickFalse)
367 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
368 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
370 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
375 (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
376 for (x=0; x < (ssize_t) image->columns; x++)
388 if (GetPixelReadMask(image,p) == 0)
390 p+=GetPixelChannels(image);
391 q+=GetPixelChannels(reconstruct_image);
394 difference=MagickFalse;
395 Sa=QuantumScale*GetPixelAlpha(image,p);
396 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
397 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
402 PixelChannel channel=GetPixelChannelChannel(image,i);
403 PixelTrait traits=GetPixelChannelTraits(image,channel);
404 PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
406 if ((traits == UndefinedPixelTrait) ||
407 (reconstruct_traits == UndefinedPixelTrait) ||
408 ((reconstruct_traits & UpdatePixelTrait) == 0))
410 distance=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
411 if ((distance*distance) > fuzz)
413 difference=MagickTrue;
417 if (difference != MagickFalse)
419 channel_distortion[i]++;
420 channel_distortion[CompositePixelChannel]++;
422 p+=GetPixelChannels(image);
423 q+=GetPixelChannels(reconstruct_image);
425 #if defined(MAGICKCORE_OPENMP_SUPPORT)
426 #pragma omp critical (MagickCore_GetAbsoluteError)
428 for (i=0; i <= MaxPixelChannels; i++)
429 distortion[i]+=channel_distortion[i];
431 reconstruct_view=DestroyCacheView(reconstruct_view);
432 image_view=DestroyCacheView(image_view);
436 static size_t GetImageChannels(const Image *image)
445 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
447 PixelChannel channel=GetPixelChannelChannel(image,i);
448 PixelTrait traits=GetPixelChannelTraits(image,channel);
449 if ((traits & UpdatePixelTrait) != 0)
455 static MagickBooleanType GetFuzzDistortion(const Image *image,
456 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
472 image_view=AcquireVirtualCacheView(image,exception);
473 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
474 #if defined(MAGICKCORE_OPENMP_SUPPORT)
475 #pragma omp parallel for schedule(static,4) shared(status) \
476 magick_threads(image,image,image->rows,1)
478 for (y=0; y < (ssize_t) image->rows; y++)
481 channel_distortion[MaxPixelChannels+1];
483 register const Quantum
491 if (status == MagickFalse)
493 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
494 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
496 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
501 (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
502 for (x=0; x < (ssize_t) image->columns; x++)
511 if (GetPixelReadMask(image,p) == 0)
513 p+=GetPixelChannels(image);
514 q+=GetPixelChannels(reconstruct_image);
517 Sa=QuantumScale*GetPixelAlpha(image,p);
518 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
519 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
524 PixelChannel channel=GetPixelChannelChannel(image,i);
525 PixelTrait traits=GetPixelChannelTraits(image,channel);
526 PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
528 if ((traits == UndefinedPixelTrait) ||
529 (reconstruct_traits == UndefinedPixelTrait) ||
530 ((reconstruct_traits & UpdatePixelTrait) == 0))
532 distance=QuantumScale*(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
534 channel_distortion[i]+=distance*distance;
535 channel_distortion[CompositePixelChannel]+=distance*distance;
537 p+=GetPixelChannels(image);
538 q+=GetPixelChannels(reconstruct_image);
540 #if defined(MAGICKCORE_OPENMP_SUPPORT)
541 #pragma omp critical (MagickCore_GetFuzzDistortion)
543 for (i=0; i <= MaxPixelChannels; i++)
544 distortion[i]+=channel_distortion[i];
546 reconstruct_view=DestroyCacheView(reconstruct_view);
547 image_view=DestroyCacheView(image_view);
548 for (i=0; i <= MaxPixelChannels; i++)
549 distortion[i]/=((double) image->columns*image->rows);
550 distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
551 distortion[CompositePixelChannel]=sqrt(distortion[CompositePixelChannel]);
555 static MagickBooleanType GetMeanAbsoluteDistortion(const Image *image,
556 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
572 image_view=AcquireVirtualCacheView(image,exception);
573 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
574 #if defined(MAGICKCORE_OPENMP_SUPPORT)
575 #pragma omp parallel for schedule(static,4) shared(status) \
576 magick_threads(image,image,image->rows,1)
578 for (y=0; y < (ssize_t) image->rows; y++)
581 channel_distortion[MaxPixelChannels+1];
583 register const Quantum
591 if (status == MagickFalse)
593 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
594 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
596 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
601 (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
602 for (x=0; x < (ssize_t) image->columns; x++)
611 if (GetPixelReadMask(image,p) == 0)
613 p+=GetPixelChannels(image);
614 q+=GetPixelChannels(reconstruct_image);
617 Sa=QuantumScale*GetPixelAlpha(image,p);
618 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
619 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
624 PixelChannel channel=GetPixelChannelChannel(image,i);
625 PixelTrait traits=GetPixelChannelTraits(image,channel);
626 PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
628 if ((traits == UndefinedPixelTrait) ||
629 (reconstruct_traits == UndefinedPixelTrait) ||
630 ((reconstruct_traits & UpdatePixelTrait) == 0))
632 distance=QuantumScale*fabs(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
634 channel_distortion[i]+=distance;
635 channel_distortion[CompositePixelChannel]+=distance;
637 p+=GetPixelChannels(image);
638 q+=GetPixelChannels(reconstruct_image);
640 #if defined(MAGICKCORE_OPENMP_SUPPORT)
641 #pragma omp critical (MagickCore_GetMeanAbsoluteError)
643 for (i=0; i <= MaxPixelChannels; i++)
644 distortion[i]+=channel_distortion[i];
646 reconstruct_view=DestroyCacheView(reconstruct_view);
647 image_view=DestroyCacheView(image_view);
648 for (i=0; i <= MaxPixelChannels; i++)
649 distortion[i]/=((double) image->columns*image->rows);
650 distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
654 static MagickBooleanType GetMeanErrorPerPixel(Image *image,
655 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
676 image_view=AcquireVirtualCacheView(image,exception);
677 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
678 for (y=0; y < (ssize_t) image->rows; y++)
680 register const Quantum
687 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
688 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
690 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
695 for (x=0; x < (ssize_t) image->columns; x++)
704 if (GetPixelReadMask(image,p) == 0)
706 p+=GetPixelChannels(image);
707 q+=GetPixelChannels(reconstruct_image);
710 Sa=QuantumScale*GetPixelAlpha(image,p);
711 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
712 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
717 PixelChannel channel=GetPixelChannelChannel(image,i);
718 PixelTrait traits=GetPixelChannelTraits(image,channel);
719 PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
721 if ((traits == UndefinedPixelTrait) ||
722 (reconstruct_traits == UndefinedPixelTrait) ||
723 ((reconstruct_traits & UpdatePixelTrait) == 0))
725 distance=fabs((double) (Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
727 distortion[i]+=distance;
728 distortion[CompositePixelChannel]+=distance;
729 mean_error+=distance*distance;
730 if (distance > maximum_error)
731 maximum_error=distance;
734 p+=GetPixelChannels(image);
735 q+=GetPixelChannels(reconstruct_image);
738 reconstruct_view=DestroyCacheView(reconstruct_view);
739 image_view=DestroyCacheView(image_view);
740 image->error.mean_error_per_pixel=distortion[CompositePixelChannel]/area;
741 image->error.normalized_mean_error=QuantumScale*QuantumScale*mean_error/area;
742 image->error.normalized_maximum_error=QuantumScale*maximum_error;
746 static MagickBooleanType GetMeanSquaredDistortion(const Image *image,
747 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
763 image_view=AcquireVirtualCacheView(image,exception);
764 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
765 #if defined(MAGICKCORE_OPENMP_SUPPORT)
766 #pragma omp parallel for schedule(static,4) shared(status) \
767 magick_threads(image,image,image->rows,1)
769 for (y=0; y < (ssize_t) image->rows; y++)
772 channel_distortion[MaxPixelChannels+1];
774 register const Quantum
782 if (status == MagickFalse)
784 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
785 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
787 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
792 (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
793 for (x=0; x < (ssize_t) image->columns; x++)
802 if (GetPixelReadMask(image,p) == 0)
804 p+=GetPixelChannels(image);
805 q+=GetPixelChannels(reconstruct_image);
808 Sa=QuantumScale*GetPixelAlpha(image,p);
809 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
810 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
815 PixelChannel channel=GetPixelChannelChannel(image,i);
816 PixelTrait traits=GetPixelChannelTraits(image,channel);
817 PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
819 if ((traits == UndefinedPixelTrait) ||
820 (reconstruct_traits == UndefinedPixelTrait) ||
821 ((reconstruct_traits & UpdatePixelTrait) == 0))
823 distance=QuantumScale*(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
825 channel_distortion[i]+=distance*distance;
826 channel_distortion[CompositePixelChannel]+=distance*distance;
828 p+=GetPixelChannels(image);
829 q+=GetPixelChannels(reconstruct_image);
831 #if defined(MAGICKCORE_OPENMP_SUPPORT)
832 #pragma omp critical (MagickCore_GetMeanSquaredError)
834 for (i=0; i <= MaxPixelChannels; i++)
835 distortion[i]+=channel_distortion[i];
837 reconstruct_view=DestroyCacheView(reconstruct_view);
838 image_view=DestroyCacheView(image_view);
839 for (i=0; i <= MaxPixelChannels; i++)
840 distortion[i]/=((double) image->columns*image->rows);
841 distortion[CompositePixelChannel]/=GetImageChannels(image);
845 static MagickBooleanType GetNormalizedCrossCorrelationDistortion(
846 const Image *image,const Image *reconstruct_image,double *distortion,
847 ExceptionInfo *exception)
849 #define SimilarityImageTag "Similarity/Image"
857 *reconstruct_statistics;
875 Normalize to account for variation due to lighting and exposure condition.
877 image_statistics=GetImageStatistics(image,exception);
878 reconstruct_statistics=GetImageStatistics(reconstruct_image,exception);
881 for (i=0; i <= MaxPixelChannels; i++)
883 area=1.0/((double) image->columns*image->rows-1);
884 image_view=AcquireVirtualCacheView(image,exception);
885 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
886 for (y=0; y < (ssize_t) image->rows; y++)
888 register const Quantum
895 if (status == MagickFalse)
897 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
898 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
900 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
905 for (x=0; x < (ssize_t) image->columns; x++)
914 if (GetPixelReadMask(image,p) == 0)
916 p+=GetPixelChannels(image);
917 q+=GetPixelChannels(reconstruct_image);
920 Sa=QuantumScale*GetPixelAlpha(image,p);
921 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
922 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
924 PixelChannel channel=GetPixelChannelChannel(image,i);
925 PixelTrait traits=GetPixelChannelTraits(image,channel);
926 PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
928 if ((traits == UndefinedPixelTrait) ||
929 (reconstruct_traits == UndefinedPixelTrait) ||
930 ((reconstruct_traits & UpdatePixelTrait) == 0))
932 distortion[i]+=area*QuantumScale*(Sa*p[i]-image_statistics[i].mean)*
933 (Da*GetPixelChannel(reconstruct_image,channel,q)-
934 reconstruct_statistics[channel].mean);
936 p+=GetPixelChannels(image);
937 q+=GetPixelChannels(reconstruct_image);
939 if (image->progress_monitor != (MagickProgressMonitor) NULL)
944 proceed=SetImageProgress(image,SimilarityImageTag,progress++,
946 if (proceed == MagickFalse)
950 reconstruct_view=DestroyCacheView(reconstruct_view);
951 image_view=DestroyCacheView(image_view);
953 Divide by the standard deviation.
955 distortion[CompositePixelChannel]=0.0;
956 for (i=0; i < MaxPixelChannels; i++)
961 PixelChannel channel=GetPixelChannelChannel(image,i);
962 gamma=image_statistics[i].standard_deviation*
963 reconstruct_statistics[channel].standard_deviation;
964 gamma=PerceptibleReciprocal(gamma);
965 distortion[i]=QuantumRange*gamma*distortion[i];
966 distortion[CompositePixelChannel]+=distortion[i]*distortion[i];
968 distortion[CompositePixelChannel]=sqrt(distortion[CompositePixelChannel]/
969 GetImageChannels(image));
973 reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
974 reconstruct_statistics);
975 image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
980 static MagickBooleanType GetPeakAbsoluteDistortion(const Image *image,
981 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
994 image_view=AcquireVirtualCacheView(image,exception);
995 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
996 #if defined(MAGICKCORE_OPENMP_SUPPORT)
997 #pragma omp parallel for schedule(static,4) shared(status) \
998 magick_threads(image,image,image->rows,1)
1000 for (y=0; y < (ssize_t) image->rows; y++)
1003 channel_distortion[MaxPixelChannels+1];
1005 register const Quantum
1013 if (status == MagickFalse)
1015 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1016 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
1018 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1023 (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
1024 for (x=0; x < (ssize_t) image->columns; x++)
1033 if (GetPixelReadMask(image,p) == 0)
1035 p+=GetPixelChannels(image);
1036 q+=GetPixelChannels(reconstruct_image);
1039 Sa=QuantumScale*GetPixelAlpha(image,p);
1040 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
1041 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1046 PixelChannel channel=GetPixelChannelChannel(image,i);
1047 PixelTrait traits=GetPixelChannelTraits(image,channel);
1048 PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
1050 if ((traits == UndefinedPixelTrait) ||
1051 (reconstruct_traits == UndefinedPixelTrait) ||
1052 ((reconstruct_traits & UpdatePixelTrait) == 0))
1054 distance=QuantumScale*fabs(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
1056 if (distance > channel_distortion[i])
1057 channel_distortion[i]=distance;
1058 if (distance > channel_distortion[CompositePixelChannel])
1059 channel_distortion[CompositePixelChannel]=distance;
1061 p+=GetPixelChannels(image);
1062 q+=GetPixelChannels(reconstruct_image);
1064 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1065 #pragma omp critical (MagickCore_GetPeakAbsoluteError)
1067 for (i=0; i <= MaxPixelChannels; i++)
1068 if (channel_distortion[i] > distortion[i])
1069 distortion[i]=channel_distortion[i];
1071 reconstruct_view=DestroyCacheView(reconstruct_view);
1072 image_view=DestroyCacheView(image_view);
1076 static MagickBooleanType GetPeakSignalToNoiseRatio(const Image *image,
1077 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1085 status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1086 for (i=0; i <= MaxPixelChannels; i++)
1087 distortion[i]=20.0*log10((double) 1.0/sqrt(distortion[i]));
1091 static MagickBooleanType GetRootMeanSquaredDistortion(const Image *image,
1092 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1100 status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1101 for (i=0; i <= MaxPixelChannels; i++)
1102 distortion[i]=sqrt(distortion[i]);
1106 MagickExport MagickBooleanType GetImageDistortion(Image *image,
1107 const Image *reconstruct_image,const MetricType metric,double *distortion,
1108 ExceptionInfo *exception)
1111 *channel_distortion;
1119 assert(image != (Image *) NULL);
1120 assert(image->signature == MagickSignature);
1121 if (image->debug != MagickFalse)
1122 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1123 assert(reconstruct_image != (const Image *) NULL);
1124 assert(reconstruct_image->signature == MagickSignature);
1125 assert(distortion != (double *) NULL);
1127 if (image->debug != MagickFalse)
1128 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1129 if ((reconstruct_image->columns != image->columns) ||
1130 (reconstruct_image->rows != image->rows))
1131 ThrowBinaryException(ImageError,"ImageSizeDiffers",image->filename);
1133 Get image distortion.
1135 length=MaxPixelChannels+1;
1136 channel_distortion=(double *) AcquireQuantumMemory(length,
1137 sizeof(*channel_distortion));
1138 if (channel_distortion == (double *) NULL)
1139 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1140 (void) ResetMagickMemory(channel_distortion,0,length*
1141 sizeof(*channel_distortion));
1144 case AbsoluteErrorMetric:
1146 status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1150 case FuzzErrorMetric:
1152 status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1156 case MeanAbsoluteErrorMetric:
1158 status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1159 channel_distortion,exception);
1162 case MeanErrorPerPixelMetric:
1164 status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1168 case MeanSquaredErrorMetric:
1170 status=GetMeanSquaredDistortion(image,reconstruct_image,
1171 channel_distortion,exception);
1174 case NormalizedCrossCorrelationErrorMetric:
1177 status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1178 channel_distortion,exception);
1181 case PeakAbsoluteErrorMetric:
1183 status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1184 channel_distortion,exception);
1187 case PeakSignalToNoiseRatioMetric:
1189 status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1190 channel_distortion,exception);
1193 case RootMeanSquaredErrorMetric:
1195 status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1196 channel_distortion,exception);
1200 *distortion=channel_distortion[CompositePixelChannel];
1201 channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1206 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1210 % G e t I m a g e D i s t o r t i o n s %
1214 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1216 % GetImageDistortions() compares the pixel channels of an image to a
1217 % reconstructed image and returns the specified distortion metric for each
1220 % The format of the GetImageDistortions method is:
1222 % double *GetImageDistortions(const Image *image,
1223 % const Image *reconstruct_image,const MetricType metric,
1224 % ExceptionInfo *exception)
1226 % A description of each parameter follows:
1228 % o image: the image.
1230 % o reconstruct_image: the reconstruct image.
1232 % o metric: the metric.
1234 % o exception: return any errors or warnings in this structure.
1237 MagickExport double *GetImageDistortions(Image *image,
1238 const Image *reconstruct_image,const MetricType metric,
1239 ExceptionInfo *exception)
1242 *channel_distortion;
1250 assert(image != (Image *) NULL);
1251 assert(image->signature == MagickSignature);
1252 if (image->debug != MagickFalse)
1253 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1254 assert(reconstruct_image != (const Image *) NULL);
1255 assert(reconstruct_image->signature == MagickSignature);
1256 if (image->debug != MagickFalse)
1257 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1258 if ((reconstruct_image->columns != image->columns) ||
1259 (reconstruct_image->rows != image->rows))
1261 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1262 "ImageSizeDiffers","`%s'",image->filename);
1263 return((double *) NULL);
1266 Get image distortion.
1268 length=MaxPixelChannels+1UL;
1269 channel_distortion=(double *) AcquireQuantumMemory(length,
1270 sizeof(*channel_distortion));
1271 if (channel_distortion == (double *) NULL)
1272 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1273 (void) ResetMagickMemory(channel_distortion,0,length*
1274 sizeof(*channel_distortion));
1278 case AbsoluteErrorMetric:
1280 status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1284 case FuzzErrorMetric:
1286 status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1290 case MeanAbsoluteErrorMetric:
1292 status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1293 channel_distortion,exception);
1296 case MeanErrorPerPixelMetric:
1298 status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1302 case MeanSquaredErrorMetric:
1304 status=GetMeanSquaredDistortion(image,reconstruct_image,
1305 channel_distortion,exception);
1308 case NormalizedCrossCorrelationErrorMetric:
1311 status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1312 channel_distortion,exception);
1315 case PeakAbsoluteErrorMetric:
1317 status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1318 channel_distortion,exception);
1321 case PeakSignalToNoiseRatioMetric:
1323 status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1324 channel_distortion,exception);
1327 case RootMeanSquaredErrorMetric:
1329 status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1330 channel_distortion,exception);
1334 if (status == MagickFalse)
1336 channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1337 return((double *) NULL);
1339 return(channel_distortion);
1343 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1347 % I s I m a g e s E q u a l %
1351 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1353 % IsImagesEqual() measures the difference between colors at each pixel
1354 % location of two images. A value other than 0 means the colors match
1355 % exactly. Otherwise an error measure is computed by summing over all
1356 % pixels in an image the distance squared in RGB space between each image
1357 % pixel and its corresponding pixel in the reconstruct image. The error
1358 % measure is assigned to these image members:
1360 % o mean_error_per_pixel: The mean error for any single pixel in
1363 % o normalized_mean_error: The normalized mean quantization error for
1364 % any single pixel in the image. This distance measure is normalized to
1365 % a range between 0 and 1. It is independent of the range of red, green,
1366 % and blue values in the image.
1368 % o normalized_maximum_error: The normalized maximum quantization
1369 % error for any single pixel in the image. This distance measure is
1370 % normalized to a range between 0 and 1. It is independent of the range
1371 % of red, green, and blue values in your image.
1373 % A small normalized mean square error, accessed as
1374 % image->normalized_mean_error, suggests the images are very similar in
1375 % spatial layout and color.
1377 % The format of the IsImagesEqual method is:
1379 % MagickBooleanType IsImagesEqual(Image *image,
1380 % const Image *reconstruct_image,ExceptionInfo *exception)
1382 % A description of each parameter follows.
1384 % o image: the image.
1386 % o reconstruct_image: the reconstruct image.
1388 % o exception: return any errors or warnings in this structure.
1391 MagickExport MagickBooleanType IsImagesEqual(Image *image,
1392 const Image *reconstruct_image,ExceptionInfo *exception)
1405 mean_error_per_pixel;
1410 assert(image != (Image *) NULL);
1411 assert(image->signature == MagickSignature);
1412 assert(reconstruct_image != (const Image *) NULL);
1413 assert(reconstruct_image->signature == MagickSignature);
1414 if ((reconstruct_image->columns != image->columns) ||
1415 (reconstruct_image->rows != image->rows))
1416 ThrowBinaryException(ImageError,"ImageSizeDiffers",image->filename);
1419 mean_error_per_pixel=0.0;
1421 image_view=AcquireVirtualCacheView(image,exception);
1422 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1423 for (y=0; y < (ssize_t) image->rows; y++)
1425 register const Quantum
1432 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1433 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
1435 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
1437 for (x=0; x < (ssize_t) image->columns; x++)
1442 if (GetPixelReadMask(image,p) == 0)
1444 p+=GetPixelChannels(image);
1445 q+=GetPixelChannels(reconstruct_image);
1448 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1453 PixelChannel channel=GetPixelChannelChannel(image,i);
1454 PixelTrait traits=GetPixelChannelTraits(image,channel);
1455 PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
1457 if ((traits == UndefinedPixelTrait) ||
1458 (reconstruct_traits == UndefinedPixelTrait) ||
1459 ((reconstruct_traits & UpdatePixelTrait) == 0))
1461 distance=fabs(p[i]-(double) GetPixelChannel(reconstruct_image,
1463 mean_error_per_pixel+=distance;
1464 mean_error+=distance*distance;
1465 if (distance > maximum_error)
1466 maximum_error=distance;
1469 p+=GetPixelChannels(image);
1470 q+=GetPixelChannels(reconstruct_image);
1473 reconstruct_view=DestroyCacheView(reconstruct_view);
1474 image_view=DestroyCacheView(image_view);
1475 image->error.mean_error_per_pixel=(double) (mean_error_per_pixel/area);
1476 image->error.normalized_mean_error=(double) (QuantumScale*QuantumScale*
1478 image->error.normalized_maximum_error=(double) (QuantumScale*maximum_error);
1479 status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
1484 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1488 % S i m i l a r i t y I m a g e %
1492 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1494 % SimilarityImage() compares the reference image of the image and returns the
1495 % best match offset. In addition, it returns a similarity image such that an
1496 % exact match location is completely white and if none of the pixels match,
1497 % black, otherwise some gray level in-between.
1499 % The format of the SimilarityImageImage method is:
1501 % Image *SimilarityImage(const Image *image,const Image *reference,
1502 % const MetricType metric,const double similarity_threshold,
1503 % RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
1505 % A description of each parameter follows:
1507 % o image: the image.
1509 % o reference: find an area of the image that closely resembles this image.
1511 % o metric: the metric.
1513 % o similarity_threshold: minimum distortion for (sub)image match.
1515 % o offset: the best match offset of the reference image within the image.
1517 % o similarity: the computed similarity between the images.
1519 % o exception: return any errors or warnings in this structure.
1523 static double GetSimilarityMetric(const Image *image,const Image *reference,
1524 const MetricType metric,const ssize_t x_offset,const ssize_t y_offset,
1525 ExceptionInfo *exception)
1539 SetGeometry(reference,&geometry);
1540 geometry.x=x_offset;
1541 geometry.y=y_offset;
1542 similarity_image=CropImage(image,&geometry,exception);
1543 if (similarity_image == (Image *) NULL)
1546 status=GetImageDistortion(similarity_image,reference,metric,&distortion,
1548 similarity_image=DestroyImage(similarity_image);
1549 if (status == MagickFalse)
1554 MagickExport Image *SimilarityImage(Image *image,const Image *reference,
1555 const MetricType metric,const double similarity_threshold,
1556 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
1558 #define SimilarityImageTag "Similarity/Image"
1575 assert(image != (const Image *) NULL);
1576 assert(image->signature == MagickSignature);
1577 if (image->debug != MagickFalse)
1578 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1579 assert(exception != (ExceptionInfo *) NULL);
1580 assert(exception->signature == MagickSignature);
1581 assert(offset != (RectangleInfo *) NULL);
1582 SetGeometry(reference,offset);
1583 *similarity_metric=1.0;
1584 if ((reference->columns > image->columns) || (reference->rows > image->rows))
1585 ThrowImageException(ImageError,"ImageSizeDiffers");
1586 similarity_image=CloneImage(image,image->columns-reference->columns+1,
1587 image->rows-reference->rows+1,MagickTrue,exception);
1588 if (similarity_image == (Image *) NULL)
1589 return((Image *) NULL);
1590 status=SetImageStorageClass(similarity_image,DirectClass,exception);
1591 if (status == MagickFalse)
1593 similarity_image=DestroyImage(similarity_image);
1594 return((Image *) NULL);
1596 (void) SetImageAlphaChannel(similarity_image,DeactivateAlphaChannel,
1599 Measure similarity of reference image against image.
1603 similarity_view=AcquireAuthenticCacheView(similarity_image,exception);
1604 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1605 #pragma omp parallel for schedule(static,4) \
1606 shared(progress,status,similarity_metric) \
1607 magick_threads(image,image,image->rows,1)
1609 for (y=0; y < (ssize_t) (image->rows-reference->rows+1); y++)
1620 if (status == MagickFalse)
1622 if (*similarity_metric <= similarity_threshold)
1624 q=GetCacheViewAuthenticPixels(similarity_view,0,y,similarity_image->columns,
1626 if (q == (Quantum *) NULL)
1631 for (x=0; x < (ssize_t) (image->columns-reference->columns+1); x++)
1636 if (*similarity_metric <= similarity_threshold)
1638 similarity=GetSimilarityMetric(image,reference,metric,x,y,exception);
1639 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1640 #pragma omp critical (MagickCore_SimilarityImage)
1642 if (similarity < *similarity_metric)
1646 *similarity_metric=similarity;
1648 if (GetPixelReadMask(similarity_image,q) == 0)
1650 q+=GetPixelChannels(similarity_image);
1653 for (i=0; i < (ssize_t) GetPixelChannels(similarity_image); i++)
1655 PixelChannel channel=GetPixelChannelChannel(image,i);
1656 PixelTrait traits=GetPixelChannelTraits(image,channel);
1657 PixelTrait similarity_traits=GetPixelChannelTraits(similarity_image,
1659 if ((traits == UndefinedPixelTrait) ||
1660 (similarity_traits == UndefinedPixelTrait) ||
1661 ((similarity_traits & UpdatePixelTrait) == 0))
1663 SetPixelChannel(similarity_image,channel,ClampToQuantum(QuantumRange-
1664 QuantumRange*similarity),q);
1666 q+=GetPixelChannels(similarity_image);
1668 if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
1670 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1675 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1676 #pragma omp critical (MagickCore_SimilarityImage)
1678 proceed=SetImageProgress(image,SimilarityImageTag,progress++,
1680 if (proceed == MagickFalse)
1684 similarity_view=DestroyCacheView(similarity_view);
1685 if (status == MagickFalse)
1686 similarity_image=DestroyImage(similarity_image);
1687 return(similarity_image);