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-2011 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 "magick/studio.h"
44 #include "magick/artifact.h"
45 #include "magick/cache-view.h"
46 #include "magick/client.h"
47 #include "magick/color.h"
48 #include "magick/color-private.h"
49 #include "magick/colorspace.h"
50 #include "magick/colorspace-private.h"
51 #include "magick/compare.h"
52 #include "magick/composite-private.h"
53 #include "magick/constitute.h"
54 #include "magick/exception-private.h"
55 #include "magick/geometry.h"
56 #include "magick/image-private.h"
57 #include "magick/list.h"
58 #include "magick/log.h"
59 #include "magick/memory_.h"
60 #include "magick/monitor.h"
61 #include "magick/monitor-private.h"
62 #include "magick/option.h"
63 #include "magick/pixel-private.h"
64 #include "magick/resource_.h"
65 #include "magick/string_.h"
66 #include "magick/statistic.h"
67 #include "magick/transform.h"
68 #include "magick/utility.h"
69 #include "magick/version.h"
72 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
76 % C o m p a r e I m a g e C h a n n e l s %
80 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
82 % CompareImageChannels() compares one or more image channels of an image
83 % to a reconstructed image and returns the difference image.
85 % The format of the CompareImageChannels method is:
87 % Image *CompareImageChannels(const Image *image,
88 % const Image *reconstruct_image,const ChannelType channel,
89 % const MetricType metric,double *distortion,ExceptionInfo *exception)
91 % A description of each parameter follows:
95 % o reconstruct_image: the reconstruct image.
97 % o channel: the channel.
99 % o metric: the metric.
101 % o distortion: the computed distortion between the images.
103 % o exception: return any errors or warnings in this structure.
107 MagickExport Image *CompareImages(Image *image,const Image *reconstruct_image,
108 const MetricType metric,double *distortion,ExceptionInfo *exception)
113 highlight_image=CompareImageChannels(image,reconstruct_image,AllChannels,
114 metric,distortion,exception);
115 return(highlight_image);
118 MagickExport Image *CompareImageChannels(Image *image,
119 const Image *reconstruct_image,const ChannelType channel,
120 const MetricType metric,double *distortion,ExceptionInfo *exception)
145 assert(image != (Image *) NULL);
146 assert(image->signature == MagickSignature);
147 if (image->debug != MagickFalse)
148 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
149 assert(reconstruct_image != (const Image *) NULL);
150 assert(reconstruct_image->signature == MagickSignature);
151 assert(distortion != (double *) NULL);
153 if (image->debug != MagickFalse)
154 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
155 if ((reconstruct_image->columns != image->columns) ||
156 (reconstruct_image->rows != image->rows))
157 ThrowImageException(ImageError,"ImageSizeDiffers");
158 status=GetImageChannelDistortion(image,reconstruct_image,channel,metric,
159 distortion,exception);
160 if (status == MagickFalse)
161 return((Image *) NULL);
162 difference_image=CloneImage(image,0,0,MagickTrue,exception);
163 if (difference_image == (Image *) NULL)
164 return((Image *) NULL);
165 (void) SetImageAlphaChannel(difference_image,OpaqueAlphaChannel);
166 highlight_image=CloneImage(image,image->columns,image->rows,MagickTrue,
168 if (highlight_image == (Image *) NULL)
170 difference_image=DestroyImage(difference_image);
171 return((Image *) NULL);
173 if (SetImageStorageClass(highlight_image,DirectClass) == MagickFalse)
175 InheritException(exception,&highlight_image->exception);
176 difference_image=DestroyImage(difference_image);
177 highlight_image=DestroyImage(highlight_image);
178 return((Image *) NULL);
180 (void) SetImageAlphaChannel(highlight_image,OpaqueAlphaChannel);
181 (void) QueryMagickColor("#f1001ecc",&highlight,exception);
182 artifact=GetImageArtifact(image,"highlight-color");
183 if (artifact != (const char *) NULL)
184 (void) QueryMagickColor(artifact,&highlight,exception);
185 (void) QueryMagickColor("#ffffffcc",&lowlight,exception);
186 artifact=GetImageArtifact(image,"lowlight-color");
187 if (artifact != (const char *) NULL)
188 (void) QueryMagickColor(artifact,&lowlight,exception);
189 if (highlight_image->colorspace == CMYKColorspace)
191 ConvertRGBToCMYK(&highlight);
192 ConvertRGBToCMYK(&lowlight);
195 Generate difference image.
198 GetMagickPixelPacket(image,&zero);
199 image_view=AcquireCacheView(image);
200 reconstruct_view=AcquireCacheView(reconstruct_image);
201 highlight_view=AcquireCacheView(highlight_image);
202 #if defined(MAGICKCORE_OPENMP_SUPPORT)
203 #pragma omp parallel for schedule(dynamic,4) shared(status)
205 for (y=0; y < (ssize_t) image->rows; y++)
214 register const IndexPacket
216 *restrict reconstruct_indexes;
218 register const PixelPacket
223 *restrict highlight_indexes;
231 if (status == MagickFalse)
233 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
234 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
236 r=QueueCacheViewAuthenticPixels(highlight_view,0,y,highlight_image->columns,
238 if ((p == (const PixelPacket *) NULL) ||
239 (q == (const PixelPacket *) NULL) || (r == (PixelPacket *) NULL))
244 indexes=GetCacheViewVirtualIndexQueue(image_view);
245 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
246 highlight_indexes=GetCacheViewAuthenticIndexQueue(highlight_view);
248 reconstruct_pixel=zero;
249 for (x=0; x < (ssize_t) image->columns; x++)
254 SetMagickPixelPacket(image,p,indexes+x,&pixel);
255 SetMagickPixelPacket(reconstruct_image,q,reconstruct_indexes+x,
257 difference=MagickFalse;
258 if (channel == AllChannels)
260 if (IsMagickColorSimilar(&pixel,&reconstruct_pixel) == MagickFalse)
261 difference=MagickTrue;
265 if (((channel & RedChannel) != 0) && (p->red != q->red))
266 difference=MagickTrue;
267 if (((channel & GreenChannel) != 0) && (p->green != q->green))
268 difference=MagickTrue;
269 if (((channel & BlueChannel) != 0) && (p->blue != q->blue))
270 difference=MagickTrue;
271 if (((channel & OpacityChannel) != 0) &&
272 (image->matte != MagickFalse) && (p->opacity != q->opacity))
273 difference=MagickTrue;
274 if ((((channel & IndexChannel) != 0) &&
275 (image->colorspace == CMYKColorspace) &&
276 (reconstruct_image->colorspace == CMYKColorspace)) &&
277 (indexes[x] != reconstruct_indexes[x]))
278 difference=MagickTrue;
280 if (difference != MagickFalse)
281 SetPixelPacket(highlight_image,&highlight,r,highlight_indexes+x);
283 SetPixelPacket(highlight_image,&lowlight,r,highlight_indexes+x);
288 sync=SyncCacheViewAuthenticPixels(highlight_view,exception);
289 if (sync == MagickFalse)
292 highlight_view=DestroyCacheView(highlight_view);
293 reconstruct_view=DestroyCacheView(reconstruct_view);
294 image_view=DestroyCacheView(image_view);
295 (void) CompositeImage(difference_image,image->compose,highlight_image,0,0);
296 highlight_image=DestroyImage(highlight_image);
297 if (status == MagickFalse)
298 difference_image=DestroyImage(difference_image);
299 return(difference_image);
303 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
307 % G e t I m a g e C h a n n e l D i s t o r t i o n %
311 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
313 % GetImageChannelDistortion() compares one or more image channels of an image
314 % to a reconstructed image and returns the specified distortion metric.
316 % The format of the CompareImageChannels method is:
318 % MagickBooleanType GetImageChannelDistortion(const Image *image,
319 % const Image *reconstruct_image,const ChannelType channel,
320 % const MetricType metric,double *distortion,ExceptionInfo *exception)
322 % A description of each parameter follows:
324 % o image: the image.
326 % o reconstruct_image: the reconstruct image.
328 % o channel: the channel.
330 % o metric: the metric.
332 % o distortion: the computed distortion between the images.
334 % o exception: return any errors or warnings in this structure.
338 MagickExport MagickBooleanType GetImageDistortion(Image *image,
339 const Image *reconstruct_image,const MetricType metric,double *distortion,
340 ExceptionInfo *exception)
345 status=GetImageChannelDistortion(image,reconstruct_image,AllChannels,
346 metric,distortion,exception);
350 static MagickBooleanType GetAbsoluteDistortion(const Image *image,
351 const Image *reconstruct_image,const ChannelType channel,double *distortion,
352 ExceptionInfo *exception)
368 Compute the absolute difference in pixels between two images.
371 GetMagickPixelPacket(image,&zero);
372 image_view=AcquireCacheView(image);
373 reconstruct_view=AcquireCacheView(reconstruct_image);
374 #if defined(MAGICKCORE_OPENMP_SUPPORT)
375 #pragma omp parallel for schedule(dynamic,4) shared(status)
377 for (y=0; y < (ssize_t) image->rows; y++)
380 channel_distortion[AllChannels+1];
386 register const IndexPacket
388 *restrict reconstruct_indexes;
390 register const PixelPacket
398 if (status == MagickFalse)
400 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
401 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
403 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
408 indexes=GetCacheViewVirtualIndexQueue(image_view);
409 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
411 reconstruct_pixel=pixel;
412 (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
413 for (x=0; x < (ssize_t) image->columns; x++)
415 SetMagickPixelPacket(image,p,indexes+x,&pixel);
416 SetMagickPixelPacket(reconstruct_image,q,reconstruct_indexes+x,
418 if (IsMagickColorSimilar(&pixel,&reconstruct_pixel) == MagickFalse)
420 if ((channel & RedChannel) != 0)
421 channel_distortion[RedChannel]++;
422 if ((channel & GreenChannel) != 0)
423 channel_distortion[GreenChannel]++;
424 if ((channel & BlueChannel) != 0)
425 channel_distortion[BlueChannel]++;
426 if (((channel & OpacityChannel) != 0) &&
427 (image->matte != MagickFalse))
428 channel_distortion[OpacityChannel]++;
429 if (((channel & IndexChannel) != 0) &&
430 (image->colorspace == CMYKColorspace))
431 channel_distortion[BlackChannel]++;
432 channel_distortion[AllChannels]++;
437 #if defined(MAGICKCORE_OPENMP_SUPPORT)
438 #pragma omp critical (MagickCore_GetAbsoluteError)
440 for (i=0; i <= (ssize_t) AllChannels; i++)
441 distortion[i]+=channel_distortion[i];
443 reconstruct_view=DestroyCacheView(reconstruct_view);
444 image_view=DestroyCacheView(image_view);
448 static size_t GetNumberChannels(const Image *image,
449 const ChannelType channel)
455 if ((channel & RedChannel) != 0)
457 if ((channel & GreenChannel) != 0)
459 if ((channel & BlueChannel) != 0)
461 if (((channel & OpacityChannel) != 0) &&
462 (image->matte != MagickFalse))
464 if (((channel & IndexChannel) != 0) &&
465 (image->colorspace == CMYKColorspace))
470 static MagickBooleanType GetFuzzDistortion(const Image *image,
471 const Image *reconstruct_image,const ChannelType channel,
472 double *distortion,ExceptionInfo *exception)
488 image_view=AcquireCacheView(image);
489 reconstruct_view=AcquireCacheView(reconstruct_image);
490 #if defined(MAGICKCORE_OPENMP_SUPPORT)
491 #pragma omp parallel for schedule(dynamic,4) shared(status)
493 for (y=0; y < (ssize_t) image->rows; y++)
496 channel_distortion[AllChannels+1];
498 register const IndexPacket
500 *restrict reconstruct_indexes;
502 register const PixelPacket
510 if (status == MagickFalse)
512 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
513 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
515 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
520 indexes=GetCacheViewVirtualIndexQueue(image_view);
521 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
522 (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
523 for (x=0; x < (ssize_t) image->columns; x++)
528 if ((channel & RedChannel) != 0)
530 distance=QuantumScale*(p->red-(MagickRealType) q->red);
531 channel_distortion[RedChannel]+=distance*distance;
532 channel_distortion[AllChannels]+=distance*distance;
534 if ((channel & GreenChannel) != 0)
536 distance=QuantumScale*(p->green-(MagickRealType) q->green);
537 channel_distortion[GreenChannel]+=distance*distance;
538 channel_distortion[AllChannels]+=distance*distance;
540 if ((channel & BlueChannel) != 0)
542 distance=QuantumScale*(p->blue-(MagickRealType) q->blue);
543 channel_distortion[BlueChannel]+=distance*distance;
544 channel_distortion[AllChannels]+=distance*distance;
546 if (((channel & OpacityChannel) != 0) && ((image->matte != MagickFalse) ||
547 (reconstruct_image->matte != MagickFalse)))
549 distance=QuantumScale*((image->matte != MagickFalse ? p->opacity :
550 OpaqueOpacity)-(reconstruct_image->matte != MagickFalse ?
551 q->opacity : OpaqueOpacity));
552 channel_distortion[OpacityChannel]+=distance*distance;
553 channel_distortion[AllChannels]+=distance*distance;
555 if (((channel & IndexChannel) != 0) &&
556 (image->colorspace == CMYKColorspace) &&
557 (reconstruct_image->colorspace == CMYKColorspace))
559 distance=QuantumScale*(indexes[x]-(MagickRealType)
560 reconstruct_indexes[x]);
561 channel_distortion[BlackChannel]+=distance*distance;
562 channel_distortion[AllChannels]+=distance*distance;
567 #if defined(MAGICKCORE_OPENMP_SUPPORT)
568 #pragma omp critical (MagickCore_GetMeanSquaredError)
570 for (i=0; i <= (ssize_t) AllChannels; i++)
571 distortion[i]+=channel_distortion[i];
573 reconstruct_view=DestroyCacheView(reconstruct_view);
574 image_view=DestroyCacheView(image_view);
575 for (i=0; i <= (ssize_t) AllChannels; i++)
576 distortion[i]/=((double) image->columns*image->rows);
577 if (((channel & OpacityChannel) != 0) && ((image->matte != MagickFalse) ||
578 (reconstruct_image->matte != MagickFalse)))
579 distortion[AllChannels]/=(double) (GetNumberChannels(image,channel)-1);
581 distortion[AllChannels]/=(double) GetNumberChannels(image,channel);
582 distortion[AllChannels]=sqrt(distortion[AllChannels]);
586 static MagickBooleanType GetMeanAbsoluteDistortion(const Image *image,
587 const Image *reconstruct_image,const ChannelType channel,
588 double *distortion,ExceptionInfo *exception)
604 image_view=AcquireCacheView(image);
605 reconstruct_view=AcquireCacheView(reconstruct_image);
606 #if defined(MAGICKCORE_OPENMP_SUPPORT)
607 #pragma omp parallel for schedule(dynamic,4) shared(status)
609 for (y=0; y < (ssize_t) image->rows; y++)
612 channel_distortion[AllChannels+1];
614 register const IndexPacket
616 *restrict reconstruct_indexes;
618 register const PixelPacket
626 if (status == MagickFalse)
628 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
629 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,
630 reconstruct_image->columns,1,exception);
631 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
636 indexes=GetCacheViewVirtualIndexQueue(image_view);
637 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
638 (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
639 for (x=0; x < (ssize_t) image->columns; x++)
644 if ((channel & RedChannel) != 0)
646 distance=QuantumScale*fabs(p->red-(double) q->red);
647 channel_distortion[RedChannel]+=distance;
648 channel_distortion[AllChannels]+=distance;
650 if ((channel & GreenChannel) != 0)
652 distance=QuantumScale*fabs(p->green-(double) q->green);
653 channel_distortion[GreenChannel]+=distance;
654 channel_distortion[AllChannels]+=distance;
656 if ((channel & BlueChannel) != 0)
658 distance=QuantumScale*fabs(p->blue-(double) q->blue);
659 channel_distortion[BlueChannel]+=distance;
660 channel_distortion[AllChannels]+=distance;
662 if (((channel & OpacityChannel) != 0) &&
663 (image->matte != MagickFalse))
665 distance=QuantumScale*fabs(p->opacity-(double) q->opacity);
666 channel_distortion[OpacityChannel]+=distance;
667 channel_distortion[AllChannels]+=distance;
669 if (((channel & IndexChannel) != 0) &&
670 (image->colorspace == CMYKColorspace))
672 distance=QuantumScale*fabs(indexes[x]-(double)
673 reconstruct_indexes[x]);
674 channel_distortion[BlackChannel]+=distance;
675 channel_distortion[AllChannels]+=distance;
680 #if defined(MAGICKCORE_OPENMP_SUPPORT)
681 #pragma omp critical (MagickCore_GetMeanAbsoluteError)
683 for (i=0; i <= (ssize_t) AllChannels; i++)
684 distortion[i]+=channel_distortion[i];
686 reconstruct_view=DestroyCacheView(reconstruct_view);
687 image_view=DestroyCacheView(image_view);
688 for (i=0; i <= (ssize_t) AllChannels; i++)
689 distortion[i]/=((double) image->columns*image->rows);
690 distortion[AllChannels]/=(double) GetNumberChannels(image,channel);
694 static MagickBooleanType GetMeanErrorPerPixel(Image *image,
695 const Image *reconstruct_image,const ChannelType channel,double *distortion,
696 ExceptionInfo *exception)
721 image_view=AcquireCacheView(image);
722 reconstruct_view=AcquireCacheView(reconstruct_image);
723 for (y=0; y < (ssize_t) image->rows; y++)
725 register const IndexPacket
727 *restrict reconstruct_indexes;
729 register const PixelPacket
736 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
737 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
739 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
744 indexes=GetCacheViewVirtualIndexQueue(image_view);
745 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
746 for (x=0; x < (ssize_t) image->columns; x++)
751 if ((channel & OpacityChannel) != 0)
753 if (image->matte != MagickFalse)
754 alpha=(MagickRealType) (QuantumScale*(GetAlphaPixelComponent(p)));
755 if (reconstruct_image->matte != MagickFalse)
756 beta=(MagickRealType) (QuantumScale*GetAlphaPixelComponent(q));
758 if ((channel & RedChannel) != 0)
760 distance=fabs(alpha*p->red-beta*q->red);
761 distortion[RedChannel]+=distance;
762 distortion[AllChannels]+=distance;
763 mean_error+=distance*distance;
764 if (distance > maximum_error)
765 maximum_error=distance;
768 if ((channel & GreenChannel) != 0)
770 distance=fabs(alpha*p->green-beta*q->green);
771 distortion[GreenChannel]+=distance;
772 distortion[AllChannels]+=distance;
773 mean_error+=distance*distance;
774 if (distance > maximum_error)
775 maximum_error=distance;
778 if ((channel & BlueChannel) != 0)
780 distance=fabs(alpha*p->blue-beta*q->blue);
781 distortion[BlueChannel]+=distance;
782 distortion[AllChannels]+=distance;
783 mean_error+=distance*distance;
784 if (distance > maximum_error)
785 maximum_error=distance;
788 if (((channel & OpacityChannel) != 0) &&
789 (image->matte != MagickFalse))
791 distance=fabs((double) p->opacity-q->opacity);
792 distortion[OpacityChannel]+=distance;
793 distortion[AllChannels]+=distance;
794 mean_error+=distance*distance;
795 if (distance > maximum_error)
796 maximum_error=distance;
799 if (((channel & IndexChannel) != 0) &&
800 (image->colorspace == CMYKColorspace) &&
801 (reconstruct_image->colorspace == CMYKColorspace))
803 distance=fabs(alpha*indexes[x]-beta*reconstruct_indexes[x]);
804 distortion[BlackChannel]+=distance;
805 distortion[AllChannels]+=distance;
806 mean_error+=distance*distance;
807 if (distance > maximum_error)
808 maximum_error=distance;
815 reconstruct_view=DestroyCacheView(reconstruct_view);
816 image_view=DestroyCacheView(image_view);
817 image->error.mean_error_per_pixel=distortion[AllChannels]/area;
818 image->error.normalized_mean_error=QuantumScale*QuantumScale*mean_error/area;
819 image->error.normalized_maximum_error=QuantumScale*maximum_error;
823 static MagickBooleanType GetMeanSquaredDistortion(const Image *image,
824 const Image *reconstruct_image,const ChannelType channel,
825 double *distortion,ExceptionInfo *exception)
841 image_view=AcquireCacheView(image);
842 reconstruct_view=AcquireCacheView(reconstruct_image);
843 #if defined(MAGICKCORE_OPENMP_SUPPORT)
844 #pragma omp parallel for schedule(dynamic,4) shared(status)
846 for (y=0; y < (ssize_t) image->rows; y++)
849 channel_distortion[AllChannels+1];
851 register const IndexPacket
853 *restrict reconstruct_indexes;
855 register const PixelPacket
863 if (status == MagickFalse)
865 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
866 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,
867 reconstruct_image->columns,1,exception);
868 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
873 indexes=GetCacheViewVirtualIndexQueue(image_view);
874 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
875 (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
876 for (x=0; x < (ssize_t) image->columns; x++)
881 if ((channel & RedChannel) != 0)
883 distance=QuantumScale*(p->red-(MagickRealType) q->red);
884 channel_distortion[RedChannel]+=distance*distance;
885 channel_distortion[AllChannels]+=distance*distance;
887 if ((channel & GreenChannel) != 0)
889 distance=QuantumScale*(p->green-(MagickRealType) q->green);
890 channel_distortion[GreenChannel]+=distance*distance;
891 channel_distortion[AllChannels]+=distance*distance;
893 if ((channel & BlueChannel) != 0)
895 distance=QuantumScale*(p->blue-(MagickRealType) q->blue);
896 channel_distortion[BlueChannel]+=distance*distance;
897 channel_distortion[AllChannels]+=distance*distance;
899 if (((channel & OpacityChannel) != 0) &&
900 (image->matte != MagickFalse))
902 distance=QuantumScale*(p->opacity-(MagickRealType) q->opacity);
903 channel_distortion[OpacityChannel]+=distance*distance;
904 channel_distortion[AllChannels]+=distance*distance;
906 if (((channel & IndexChannel) != 0) &&
907 (image->colorspace == CMYKColorspace) &&
908 (reconstruct_image->colorspace == CMYKColorspace))
910 distance=QuantumScale*(indexes[x]-(MagickRealType)
911 reconstruct_indexes[x]);
912 channel_distortion[BlackChannel]+=distance*distance;
913 channel_distortion[AllChannels]+=distance*distance;
918 #if defined(MAGICKCORE_OPENMP_SUPPORT)
919 #pragma omp critical (MagickCore_GetMeanSquaredError)
921 for (i=0; i <= (ssize_t) AllChannels; i++)
922 distortion[i]+=channel_distortion[i];
924 reconstruct_view=DestroyCacheView(reconstruct_view);
925 image_view=DestroyCacheView(image_view);
926 for (i=0; i <= (ssize_t) AllChannels; i++)
927 distortion[i]/=((double) image->columns*image->rows);
928 distortion[AllChannels]/=(double) GetNumberChannels(image,channel);
932 static MagickBooleanType GetNormalizedCrossCorrelationDistortion(
933 const Image *image,const Image *reconstruct_image,const ChannelType channel,
934 double *distortion,ExceptionInfo *exception)
936 #define SimilarityImageTag "Similarity/Image"
944 *reconstruct_statistics;
962 Normalize to account for variation due to lighting and exposure condition.
964 image_statistics=GetImageChannelStatistics(image,exception);
965 reconstruct_statistics=GetImageChannelStatistics(reconstruct_image,exception);
968 for (i=0; i <= (ssize_t) AllChannels; i++)
970 area=1.0/((MagickRealType) image->columns*image->rows);
971 image_view=AcquireCacheView(image);
972 reconstruct_view=AcquireCacheView(reconstruct_image);
973 for (y=0; y < (ssize_t) image->rows; y++)
975 register const IndexPacket
977 *restrict reconstruct_indexes;
979 register const PixelPacket
986 if (status == MagickFalse)
988 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
989 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
991 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
996 indexes=GetCacheViewVirtualIndexQueue(image_view);
997 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
998 for (x=0; x < (ssize_t) image->columns; x++)
1000 if ((channel & RedChannel) != 0)
1001 distortion[RedChannel]+=area*QuantumScale*(p->red-
1002 image_statistics[RedChannel].mean)*(q->red-
1003 reconstruct_statistics[RedChannel].mean);
1004 if ((channel & GreenChannel) != 0)
1005 distortion[GreenChannel]+=area*QuantumScale*(p->green-
1006 image_statistics[GreenChannel].mean)*(q->green-
1007 reconstruct_statistics[GreenChannel].mean);
1008 if ((channel & BlueChannel) != 0)
1009 distortion[BlueChannel]+=area*QuantumScale*(p->blue-
1010 image_statistics[BlueChannel].mean)*(q->blue-
1011 reconstruct_statistics[BlueChannel].mean);
1012 if (((channel & OpacityChannel) != 0) &&
1013 (image->matte != MagickFalse))
1014 distortion[OpacityChannel]+=area*QuantumScale*(p->opacity-
1015 image_statistics[OpacityChannel].mean)*(q->opacity-
1016 reconstruct_statistics[OpacityChannel].mean);
1017 if (((channel & IndexChannel) != 0) &&
1018 (image->colorspace == CMYKColorspace) &&
1019 (reconstruct_image->colorspace == CMYKColorspace))
1020 distortion[BlackChannel]+=area*QuantumScale*(indexes[x]-
1021 image_statistics[OpacityChannel].mean)*(reconstruct_indexes[x]-
1022 reconstruct_statistics[OpacityChannel].mean);
1026 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1031 proceed=SetImageProgress(image,SimilarityImageTag,progress++,
1033 if (proceed == MagickFalse)
1037 reconstruct_view=DestroyCacheView(reconstruct_view);
1038 image_view=DestroyCacheView(image_view);
1040 Divide by the standard deviation.
1042 for (i=0; i < (ssize_t) AllChannels; i++)
1047 gamma=image_statistics[i].standard_deviation*
1048 reconstruct_statistics[i].standard_deviation;
1049 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1050 distortion[i]=QuantumRange*gamma*distortion[i];
1052 distortion[AllChannels]=0.0;
1053 if ((channel & RedChannel) != 0)
1054 distortion[AllChannels]+=distortion[RedChannel]*distortion[RedChannel];
1055 if ((channel & GreenChannel) != 0)
1056 distortion[AllChannels]+=distortion[GreenChannel]*distortion[GreenChannel];
1057 if ((channel & BlueChannel) != 0)
1058 distortion[AllChannels]+=distortion[BlueChannel]*distortion[BlueChannel];
1059 if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
1060 distortion[AllChannels]+=distortion[OpacityChannel]*
1061 distortion[OpacityChannel];
1062 if (((channel & IndexChannel) != 0) &&
1063 (image->colorspace == CMYKColorspace))
1064 distortion[AllChannels]+=distortion[BlackChannel]*distortion[BlackChannel];
1065 distortion[AllChannels]=sqrt(distortion[AllChannels]/GetNumberChannels(image,
1070 reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1071 reconstruct_statistics);
1072 image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1077 static MagickBooleanType GetPeakAbsoluteDistortion(const Image *image,
1078 const Image *reconstruct_image,const ChannelType channel,
1079 double *distortion,ExceptionInfo *exception)
1092 image_view=AcquireCacheView(image);
1093 reconstruct_view=AcquireCacheView(reconstruct_image);
1094 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1095 #pragma omp parallel for schedule(dynamic,4) shared(status)
1097 for (y=0; y < (ssize_t) image->rows; y++)
1100 channel_distortion[AllChannels+1];
1102 register const IndexPacket
1104 *restrict reconstruct_indexes;
1106 register const PixelPacket
1114 if (status == MagickFalse)
1116 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1117 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,
1118 reconstruct_image->columns,1,exception);
1119 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1124 indexes=GetCacheViewVirtualIndexQueue(image_view);
1125 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
1126 (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
1127 for (x=0; x < (ssize_t) image->columns; x++)
1132 if ((channel & RedChannel) != 0)
1134 distance=QuantumScale*fabs(p->red-(double) q->red);
1135 if (distance > channel_distortion[RedChannel])
1136 channel_distortion[RedChannel]=distance;
1137 if (distance > channel_distortion[AllChannels])
1138 channel_distortion[AllChannels]=distance;
1140 if ((channel & GreenChannel) != 0)
1142 distance=QuantumScale*fabs(p->green-(double) q->green);
1143 if (distance > channel_distortion[GreenChannel])
1144 channel_distortion[GreenChannel]=distance;
1145 if (distance > channel_distortion[AllChannels])
1146 channel_distortion[AllChannels]=distance;
1148 if ((channel & BlueChannel) != 0)
1150 distance=QuantumScale*fabs(p->blue-(double) q->blue);
1151 if (distance > channel_distortion[BlueChannel])
1152 channel_distortion[BlueChannel]=distance;
1153 if (distance > channel_distortion[AllChannels])
1154 channel_distortion[AllChannels]=distance;
1156 if (((channel & OpacityChannel) != 0) &&
1157 (image->matte != MagickFalse))
1159 distance=QuantumScale*fabs(p->opacity-(double) q->opacity);
1160 if (distance > channel_distortion[OpacityChannel])
1161 channel_distortion[OpacityChannel]=distance;
1162 if (distance > channel_distortion[AllChannels])
1163 channel_distortion[AllChannels]=distance;
1165 if (((channel & IndexChannel) != 0) &&
1166 (image->colorspace == CMYKColorspace) &&
1167 (reconstruct_image->colorspace == CMYKColorspace))
1169 distance=QuantumScale*fabs(indexes[x]-(double)
1170 reconstruct_indexes[x]);
1171 if (distance > channel_distortion[BlackChannel])
1172 channel_distortion[BlackChannel]=distance;
1173 if (distance > channel_distortion[AllChannels])
1174 channel_distortion[AllChannels]=distance;
1179 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1180 #pragma omp critical (MagickCore_GetPeakAbsoluteError)
1182 for (i=0; i <= (ssize_t) AllChannels; i++)
1183 if (channel_distortion[i] > distortion[i])
1184 distortion[i]=channel_distortion[i];
1186 reconstruct_view=DestroyCacheView(reconstruct_view);
1187 image_view=DestroyCacheView(image_view);
1191 static MagickBooleanType GetPeakSignalToNoiseRatio(const Image *image,
1192 const Image *reconstruct_image,const ChannelType channel,
1193 double *distortion,ExceptionInfo *exception)
1198 status=GetMeanSquaredDistortion(image,reconstruct_image,channel,distortion,
1200 if ((channel & RedChannel) != 0)
1201 distortion[RedChannel]=20.0*log10((double) 1.0/sqrt(
1202 distortion[RedChannel]));
1203 if ((channel & GreenChannel) != 0)
1204 distortion[GreenChannel]=20.0*log10((double) 1.0/sqrt(
1205 distortion[GreenChannel]));
1206 if ((channel & BlueChannel) != 0)
1207 distortion[BlueChannel]=20.0*log10((double) 1.0/sqrt(
1208 distortion[BlueChannel]));
1209 if (((channel & OpacityChannel) != 0) &&
1210 (image->matte != MagickFalse))
1211 distortion[OpacityChannel]=20.0*log10((double) 1.0/sqrt(
1212 distortion[OpacityChannel]));
1213 if (((channel & IndexChannel) != 0) &&
1214 (image->colorspace == CMYKColorspace))
1215 distortion[BlackChannel]=20.0*log10((double) 1.0/sqrt(
1216 distortion[BlackChannel]));
1217 distortion[AllChannels]=20.0*log10((double) 1.0/sqrt(
1218 distortion[AllChannels]));
1222 static MagickBooleanType GetRootMeanSquaredDistortion(const Image *image,
1223 const Image *reconstruct_image,const ChannelType channel,
1224 double *distortion,ExceptionInfo *exception)
1229 status=GetMeanSquaredDistortion(image,reconstruct_image,channel,distortion,
1231 if ((channel & RedChannel) != 0)
1232 distortion[RedChannel]=sqrt(distortion[RedChannel]);
1233 if ((channel & GreenChannel) != 0)
1234 distortion[GreenChannel]=sqrt(distortion[GreenChannel]);
1235 if ((channel & BlueChannel) != 0)
1236 distortion[BlueChannel]=sqrt(distortion[BlueChannel]);
1237 if (((channel & OpacityChannel) != 0) &&
1238 (image->matte != MagickFalse))
1239 distortion[OpacityChannel]=sqrt(distortion[OpacityChannel]);
1240 if (((channel & IndexChannel) != 0) &&
1241 (image->colorspace == CMYKColorspace))
1242 distortion[BlackChannel]=sqrt(distortion[BlackChannel]);
1243 distortion[AllChannels]=sqrt(distortion[AllChannels]);
1247 MagickExport MagickBooleanType GetImageChannelDistortion(Image *image,
1248 const Image *reconstruct_image,const ChannelType channel,
1249 const MetricType metric,double *distortion,ExceptionInfo *exception)
1252 *channel_distortion;
1260 assert(image != (Image *) NULL);
1261 assert(image->signature == MagickSignature);
1262 if (image->debug != MagickFalse)
1263 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1264 assert(reconstruct_image != (const Image *) NULL);
1265 assert(reconstruct_image->signature == MagickSignature);
1266 assert(distortion != (double *) NULL);
1268 if (image->debug != MagickFalse)
1269 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1270 if ((reconstruct_image->columns != image->columns) ||
1271 (reconstruct_image->rows != image->rows))
1272 ThrowBinaryException(ImageError,"ImageSizeDiffers",image->filename);
1274 Get image distortion.
1276 length=AllChannels+1UL;
1277 channel_distortion=(double *) AcquireQuantumMemory(length,
1278 sizeof(*channel_distortion));
1279 if (channel_distortion == (double *) NULL)
1280 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1281 (void) ResetMagickMemory(channel_distortion,0,length*
1282 sizeof(*channel_distortion));
1285 case AbsoluteErrorMetric:
1287 status=GetAbsoluteDistortion(image,reconstruct_image,channel,
1288 channel_distortion,exception);
1291 case FuzzErrorMetric:
1293 status=GetFuzzDistortion(image,reconstruct_image,channel,
1294 channel_distortion,exception);
1297 case MeanAbsoluteErrorMetric:
1299 status=GetMeanAbsoluteDistortion(image,reconstruct_image,channel,
1300 channel_distortion,exception);
1303 case MeanErrorPerPixelMetric:
1305 status=GetMeanErrorPerPixel(image,reconstruct_image,channel,
1306 channel_distortion,exception);
1309 case MeanSquaredErrorMetric:
1311 status=GetMeanSquaredDistortion(image,reconstruct_image,channel,
1312 channel_distortion,exception);
1315 case NormalizedCrossCorrelationErrorMetric:
1318 status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1319 channel,channel_distortion,exception);
1322 case PeakAbsoluteErrorMetric:
1324 status=GetPeakAbsoluteDistortion(image,reconstruct_image,channel,
1325 channel_distortion,exception);
1328 case PeakSignalToNoiseRatioMetric:
1330 status=GetPeakSignalToNoiseRatio(image,reconstruct_image,channel,
1331 channel_distortion,exception);
1334 case RootMeanSquaredErrorMetric:
1336 status=GetRootMeanSquaredDistortion(image,reconstruct_image,channel,
1337 channel_distortion,exception);
1341 *distortion=channel_distortion[AllChannels];
1342 channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1347 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1351 % G e t I m a g e C h a n n e l D i s t o r t i o n s %
1355 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1357 % GetImageChannelDistrortion() compares the image channels of an image to a
1358 % reconstructed image and returns the specified distortion metric for each
1361 % The format of the CompareImageChannels method is:
1363 % double *GetImageChannelDistortions(const Image *image,
1364 % const Image *reconstruct_image,const MetricType metric,
1365 % ExceptionInfo *exception)
1367 % A description of each parameter follows:
1369 % o image: the image.
1371 % o reconstruct_image: the reconstruct image.
1373 % o metric: the metric.
1375 % o exception: return any errors or warnings in this structure.
1378 MagickExport double *GetImageChannelDistortions(Image *image,
1379 const Image *reconstruct_image,const MetricType metric,
1380 ExceptionInfo *exception)
1383 *channel_distortion;
1391 assert(image != (Image *) NULL);
1392 assert(image->signature == MagickSignature);
1393 if (image->debug != MagickFalse)
1394 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1395 assert(reconstruct_image != (const Image *) NULL);
1396 assert(reconstruct_image->signature == MagickSignature);
1397 if (image->debug != MagickFalse)
1398 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1399 if ((reconstruct_image->columns != image->columns) ||
1400 (reconstruct_image->rows != image->rows))
1402 (void) ThrowMagickException(&image->exception,GetMagickModule(),
1403 ImageError,"ImageSizeDiffers","`%s'",image->filename);
1404 return((double *) NULL);
1407 Get image distortion.
1409 length=AllChannels+1UL;
1410 channel_distortion=(double *) AcquireQuantumMemory(length,
1411 sizeof(*channel_distortion));
1412 if (channel_distortion == (double *) NULL)
1413 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1414 (void) ResetMagickMemory(channel_distortion,0,length*
1415 sizeof(*channel_distortion));
1419 case AbsoluteErrorMetric:
1421 status=GetAbsoluteDistortion(image,reconstruct_image,AllChannels,
1422 channel_distortion,exception);
1425 case FuzzErrorMetric:
1427 status=GetFuzzDistortion(image,reconstruct_image,AllChannels,
1428 channel_distortion,exception);
1431 case MeanAbsoluteErrorMetric:
1433 status=GetMeanAbsoluteDistortion(image,reconstruct_image,AllChannels,
1434 channel_distortion,exception);
1437 case MeanErrorPerPixelMetric:
1439 status=GetMeanErrorPerPixel(image,reconstruct_image,AllChannels,
1440 channel_distortion,exception);
1443 case MeanSquaredErrorMetric:
1445 status=GetMeanSquaredDistortion(image,reconstruct_image,AllChannels,
1446 channel_distortion,exception);
1449 case NormalizedCrossCorrelationErrorMetric:
1452 status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1453 AllChannels,channel_distortion,exception);
1456 case PeakAbsoluteErrorMetric:
1458 status=GetPeakAbsoluteDistortion(image,reconstruct_image,AllChannels,
1459 channel_distortion,exception);
1462 case PeakSignalToNoiseRatioMetric:
1464 status=GetPeakSignalToNoiseRatio(image,reconstruct_image,AllChannels,
1465 channel_distortion,exception);
1468 case RootMeanSquaredErrorMetric:
1470 status=GetRootMeanSquaredDistortion(image,reconstruct_image,AllChannels,
1471 channel_distortion,exception);
1475 if (status == MagickFalse)
1477 channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1478 return((double *) NULL);
1480 return(channel_distortion);
1484 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1488 % I s I m a g e s E q u a l %
1492 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1494 % IsImagesEqual() measures the difference between colors at each pixel
1495 % location of two images. A value other than 0 means the colors match
1496 % exactly. Otherwise an error measure is computed by summing over all
1497 % pixels in an image the distance squared in RGB space between each image
1498 % pixel and its corresponding pixel in the reconstruct image. The error
1499 % measure is assigned to these image members:
1501 % o mean_error_per_pixel: The mean error for any single pixel in
1504 % o normalized_mean_error: The normalized mean quantization error for
1505 % any single pixel in the image. This distance measure is normalized to
1506 % a range between 0 and 1. It is independent of the range of red, green,
1507 % and blue values in the image.
1509 % o normalized_maximum_error: The normalized maximum quantization
1510 % error for any single pixel in the image. This distance measure is
1511 % normalized to a range between 0 and 1. It is independent of the range
1512 % of red, green, and blue values in your image.
1514 % A small normalized mean square error, accessed as
1515 % image->normalized_mean_error, suggests the images are very similar in
1516 % spatial layout and color.
1518 % The format of the IsImagesEqual method is:
1520 % MagickBooleanType IsImagesEqual(Image *image,
1521 % const Image *reconstruct_image)
1523 % A description of each parameter follows.
1525 % o image: the image.
1527 % o reconstruct_image: the reconstruct image.
1530 MagickExport MagickBooleanType IsImagesEqual(Image *image,
1531 const Image *reconstruct_image)
1547 mean_error_per_pixel;
1552 assert(image != (Image *) NULL);
1553 assert(image->signature == MagickSignature);
1554 assert(reconstruct_image != (const Image *) NULL);
1555 assert(reconstruct_image->signature == MagickSignature);
1556 if ((reconstruct_image->columns != image->columns) ||
1557 (reconstruct_image->rows != image->rows))
1558 ThrowBinaryException(ImageError,"ImageSizeDiffers",image->filename);
1561 mean_error_per_pixel=0.0;
1563 exception=(&image->exception);
1564 image_view=AcquireCacheView(image);
1565 reconstruct_view=AcquireCacheView(reconstruct_image);
1566 for (y=0; y < (ssize_t) image->rows; y++)
1568 register const IndexPacket
1570 *restrict reconstruct_indexes;
1572 register const PixelPacket
1579 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1580 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
1582 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1584 indexes=GetCacheViewVirtualIndexQueue(image_view);
1585 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
1586 for (x=0; x < (ssize_t) image->columns; x++)
1591 distance=fabs(p->red-(double) q->red);
1592 mean_error_per_pixel+=distance;
1593 mean_error+=distance*distance;
1594 if (distance > maximum_error)
1595 maximum_error=distance;
1597 distance=fabs(p->green-(double) q->green);
1598 mean_error_per_pixel+=distance;
1599 mean_error+=distance*distance;
1600 if (distance > maximum_error)
1601 maximum_error=distance;
1603 distance=fabs(p->blue-(double) q->blue);
1604 mean_error_per_pixel+=distance;
1605 mean_error+=distance*distance;
1606 if (distance > maximum_error)
1607 maximum_error=distance;
1609 if (image->matte != MagickFalse)
1611 distance=fabs(p->opacity-(double) q->opacity);
1612 mean_error_per_pixel+=distance;
1613 mean_error+=distance*distance;
1614 if (distance > maximum_error)
1615 maximum_error=distance;
1618 if ((image->colorspace == CMYKColorspace) &&
1619 (reconstruct_image->colorspace == CMYKColorspace))
1621 distance=fabs(indexes[x]-(double) reconstruct_indexes[x]);
1622 mean_error_per_pixel+=distance;
1623 mean_error+=distance*distance;
1624 if (distance > maximum_error)
1625 maximum_error=distance;
1632 reconstruct_view=DestroyCacheView(reconstruct_view);
1633 image_view=DestroyCacheView(image_view);
1634 image->error.mean_error_per_pixel=(double) (mean_error_per_pixel/area);
1635 image->error.normalized_mean_error=(double) (QuantumScale*QuantumScale*
1637 image->error.normalized_maximum_error=(double) (QuantumScale*maximum_error);
1638 status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
1643 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1647 % S i m i l a r i t y I m a g e %
1651 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1653 % SimilarityImage() compares the reference image of the image and returns the
1654 % best match offset. In addition, it returns a similarity image such that an
1655 % exact match location is completely white and if none of the pixels match,
1656 % black, otherwise some gray level in-between.
1658 % The format of the SimilarityImageImage method is:
1660 % Image *SimilarityImage(const Image *image,const Image *reference,
1661 % RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
1663 % A description of each parameter follows:
1665 % o image: the image.
1667 % o reference: find an area of the image that closely resembles this image.
1669 % o the best match offset of the reference image within the image.
1671 % o similarity: the computed similarity between the images.
1673 % o exception: return any errors or warnings in this structure.
1677 static double GetNCCDistortion(const Image *image,
1678 const Image *reconstruct_image,
1679 const ChannelStatistics *reconstruct_statistics,ExceptionInfo *exception)
1681 #define SimilarityImageTag "Similarity/Image"
1707 Normalize to account for variation due to lighting and exposure condition.
1709 image_statistics=GetImageChannelStatistics(image,exception);
1712 area=1.0/((MagickRealType) image->columns*image->rows);
1713 image_view=AcquireCacheView(image);
1714 reconstruct_view=AcquireCacheView(reconstruct_image);
1715 for (y=0; y < (ssize_t) image->rows; y++)
1717 register const IndexPacket
1719 *restrict reconstruct_indexes;
1721 register const PixelPacket
1728 if (status == MagickFalse)
1730 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1731 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
1733 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1738 indexes=GetCacheViewVirtualIndexQueue(image_view);
1739 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
1740 for (x=0; x < (ssize_t) image->columns; x++)
1742 distortion+=area*QuantumScale*(p->red-
1743 image_statistics[RedChannel].mean)*(q->red-
1744 reconstruct_statistics[RedChannel].mean);
1745 distortion+=area*QuantumScale*(p->green-
1746 image_statistics[GreenChannel].mean)*(q->green-
1747 reconstruct_statistics[GreenChannel].mean);
1748 distortion+=area*QuantumScale*(p->blue-
1749 image_statistics[BlueChannel].mean)*(q->blue-
1750 reconstruct_statistics[BlueChannel].mean);
1751 if (image->matte != MagickFalse)
1752 distortion+=area*QuantumScale*(p->opacity-
1753 image_statistics[OpacityChannel].mean)*(q->opacity-
1754 reconstruct_statistics[OpacityChannel].mean);
1755 if ((image->colorspace == CMYKColorspace) &&
1756 (reconstruct_image->colorspace == CMYKColorspace))
1757 distortion+=area*QuantumScale*(indexes[x]-
1758 image_statistics[OpacityChannel].mean)*(reconstruct_indexes[x]-
1759 reconstruct_statistics[OpacityChannel].mean);
1764 reconstruct_view=DestroyCacheView(reconstruct_view);
1765 image_view=DestroyCacheView(image_view);
1767 Divide by the standard deviation.
1769 gamma=image_statistics[AllChannels].standard_deviation*
1770 reconstruct_statistics[AllChannels].standard_deviation;
1771 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1772 distortion=QuantumRange*gamma*distortion;
1774 if (image->matte != MagickFalse)
1776 if (image->colorspace == CMYKColorspace)
1778 distortion=sqrt(distortion/number_channels);
1782 image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1784 return(1.0-distortion);
1787 static double GetSimilarityMetric(const Image *image,const Image *reference,
1788 const ChannelStatistics *reference_statistics,const ssize_t x_offset,
1789 const ssize_t y_offset,ExceptionInfo *exception)
1800 SetGeometry(reference,&geometry);
1801 geometry.x=x_offset;
1802 geometry.y=y_offset;
1803 similarity_image=CropImage(image,&geometry,exception);
1804 if (similarity_image == (Image *) NULL)
1806 distortion=GetNCCDistortion(reference,similarity_image,reference_statistics,
1808 similarity_image=DestroyImage(similarity_image);
1812 MagickExport Image *SimilarityImage(Image *image,const Image *reference,
1813 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
1815 #define SimilarityImageTag "Similarity/Image"
1821 *reference_statistics;
1835 assert(image != (const Image *) NULL);
1836 assert(image->signature == MagickSignature);
1837 if (image->debug != MagickFalse)
1838 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1839 assert(exception != (ExceptionInfo *) NULL);
1840 assert(exception->signature == MagickSignature);
1841 assert(offset != (RectangleInfo *) NULL);
1842 SetGeometry(reference,offset);
1843 *similarity_metric=1.0;
1844 if ((reference->columns > image->columns) || (reference->rows > image->rows))
1845 ThrowImageException(ImageError,"ImageSizeDiffers");
1846 similarity_image=CloneImage(image,image->columns-reference->columns+1,
1847 image->rows-reference->rows+1,MagickTrue,exception);
1848 if (similarity_image == (Image *) NULL)
1849 return((Image *) NULL);
1850 if (SetImageStorageClass(similarity_image,DirectClass) == MagickFalse)
1852 InheritException(exception,&similarity_image->exception);
1853 similarity_image=DestroyImage(similarity_image);
1854 return((Image *) NULL);
1857 Measure similarity of reference image against image.
1861 reference_statistics=GetImageChannelStatistics(reference,exception);
1862 similarity_view=AcquireCacheView(similarity_image);
1863 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1864 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1866 for (y=0; y < (ssize_t) (image->rows-reference->rows+1); y++)
1874 register PixelPacket
1877 if (status == MagickFalse)
1879 q=GetCacheViewAuthenticPixels(similarity_view,0,y,similarity_image->columns,
1881 if (q == (const PixelPacket *) NULL)
1886 for (x=0; x < (ssize_t) (image->columns-reference->columns+1); x++)
1888 similarity=GetSimilarityMetric(image,reference,reference_statistics,x,y,
1890 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1891 #pragma omp critical (MagickCore_SimilarityImage)
1893 if (similarity < *similarity_metric)
1895 *similarity_metric=similarity;
1899 q->red=ClampToQuantum(QuantumRange-QuantumRange*similarity);
1904 if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
1906 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1911 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1912 #pragma omp critical (MagickCore_SimilarityImage)
1914 proceed=SetImageProgress(image,SimilarityImageTag,progress++,
1916 if (proceed == MagickFalse)
1920 similarity_view=DestroyCacheView(similarity_view);
1921 reference_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1922 reference_statistics);
1923 return(similarity_image);