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,
514 reconstruct_image->columns,1,exception);
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++)
530 if ((channel & OpacityChannel) != 0)
532 if (image->matte != MagickFalse)
533 scale*=QuantumScale*GetAlphaPixelComponent(p);
534 if (reconstruct_image->matte != MagickFalse)
535 scale*=QuantumScale*GetAlphaPixelComponent(q);
537 if (((channel & IndexChannel) != 0) &&
538 (image->colorspace == CMYKColorspace) &&
539 (reconstruct_image->colorspace == CMYKColorspace))
541 scale*=(MagickRealType) (QuantumScale*(QuantumRange-indexes[x]));
542 scale*=(MagickRealType) (QuantumScale*(QuantumRange-
543 reconstruct_indexes[x]));
545 if ((channel & RedChannel) != 0)
547 distance=QuantumScale*(p->red-(MagickRealType) q->red);
548 channel_distortion[RedChannel]+=scale*distance*distance;
549 channel_distortion[AllChannels]+=scale*distance*distance;
551 if ((channel & GreenChannel) != 0)
553 distance=QuantumScale*(p->green-(MagickRealType) q->green);
554 channel_distortion[GreenChannel]+=scale*distance*distance;
555 channel_distortion[AllChannels]+=scale*distance*distance;
557 if ((channel & BlueChannel) != 0)
559 distance=QuantumScale*(p->blue-(MagickRealType) q->blue);
560 channel_distortion[BlueChannel]+=scale*distance*distance;
561 channel_distortion[AllChannels]+=scale*distance*distance;
563 if (((channel & OpacityChannel) != 0) && ((image->matte != MagickFalse) || (reconstruct_image->matte != MagickFalse)))
565 distance=QuantumScale*((image->matte != MagickFalse ? p->opacity :
566 OpaqueOpacity)-(reconstruct_image->matte != MagickFalse ?
567 q->opacity : OpaqueOpacity));
568 channel_distortion[OpacityChannel]+=QuantumScale*distance*distance;
569 channel_distortion[AllChannels]+=QuantumScale*distance*distance;
571 if (((channel & IndexChannel) != 0) &&
572 (image->colorspace == CMYKColorspace) &&
573 (reconstruct_image->colorspace == CMYKColorspace))
575 distance=QuantumScale*(indexes[x]-(MagickRealType)
576 reconstruct_indexes[x]);
577 channel_distortion[BlackChannel]+=scale*distance*distance;
578 channel_distortion[AllChannels]+=scale*distance*distance;
583 #if defined(MAGICKCORE_OPENMP_SUPPORT)
584 #pragma omp critical (MagickCore_GetMeanSquaredError)
586 for (i=0; i <= (ssize_t) AllChannels; i++)
587 distortion[i]+=channel_distortion[i];
589 reconstruct_view=DestroyCacheView(reconstruct_view);
590 image_view=DestroyCacheView(image_view);
591 for (i=0; i <= (ssize_t) AllChannels; i++)
592 distortion[i]/=((double) image->columns*image->rows);
593 distortion[AllChannels]/=(double) GetNumberChannels(image,channel);
597 static MagickBooleanType GetMeanAbsoluteDistortion(const Image *image,
598 const Image *reconstruct_image,const ChannelType channel,
599 double *distortion,ExceptionInfo *exception)
615 image_view=AcquireCacheView(image);
616 reconstruct_view=AcquireCacheView(reconstruct_image);
617 #if defined(MAGICKCORE_OPENMP_SUPPORT)
618 #pragma omp parallel for schedule(dynamic,4) shared(status)
620 for (y=0; y < (ssize_t) image->rows; y++)
623 channel_distortion[AllChannels+1];
625 register const IndexPacket
627 *restrict reconstruct_indexes;
629 register const PixelPacket
637 if (status == MagickFalse)
639 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
640 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,
641 reconstruct_image->columns,1,exception);
642 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
647 indexes=GetCacheViewVirtualIndexQueue(image_view);
648 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
649 (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
650 for (x=0; x < (ssize_t) image->columns; x++)
655 if ((channel & RedChannel) != 0)
657 distance=QuantumScale*fabs(p->red-(double) q->red);
658 channel_distortion[RedChannel]+=distance;
659 channel_distortion[AllChannels]+=distance;
661 if ((channel & GreenChannel) != 0)
663 distance=QuantumScale*fabs(p->green-(double) q->green);
664 channel_distortion[GreenChannel]+=distance;
665 channel_distortion[AllChannels]+=distance;
667 if ((channel & BlueChannel) != 0)
669 distance=QuantumScale*fabs(p->blue-(double) q->blue);
670 channel_distortion[BlueChannel]+=distance;
671 channel_distortion[AllChannels]+=distance;
673 if (((channel & OpacityChannel) != 0) &&
674 (image->matte != MagickFalse))
676 distance=QuantumScale*fabs(p->opacity-(double) q->opacity);
677 channel_distortion[OpacityChannel]+=distance;
678 channel_distortion[AllChannels]+=distance;
680 if (((channel & IndexChannel) != 0) &&
681 (image->colorspace == CMYKColorspace))
683 distance=QuantumScale*fabs(indexes[x]-(double)
684 reconstruct_indexes[x]);
685 channel_distortion[BlackChannel]+=distance;
686 channel_distortion[AllChannels]+=distance;
691 #if defined(MAGICKCORE_OPENMP_SUPPORT)
692 #pragma omp critical (MagickCore_GetMeanAbsoluteError)
694 for (i=0; i <= (ssize_t) AllChannels; i++)
695 distortion[i]+=channel_distortion[i];
697 reconstruct_view=DestroyCacheView(reconstruct_view);
698 image_view=DestroyCacheView(image_view);
699 for (i=0; i <= (ssize_t) AllChannels; i++)
700 distortion[i]/=((double) image->columns*image->rows);
701 distortion[AllChannels]/=(double) GetNumberChannels(image,channel);
705 static MagickBooleanType GetMeanErrorPerPixel(Image *image,
706 const Image *reconstruct_image,const ChannelType channel,double *distortion,
707 ExceptionInfo *exception)
732 image_view=AcquireCacheView(image);
733 reconstruct_view=AcquireCacheView(reconstruct_image);
734 for (y=0; y < (ssize_t) image->rows; y++)
736 register const IndexPacket
738 *restrict reconstruct_indexes;
740 register const PixelPacket
747 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
748 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
750 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
755 indexes=GetCacheViewVirtualIndexQueue(image_view);
756 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
757 for (x=0; x < (ssize_t) image->columns; x++)
762 if ((channel & OpacityChannel) != 0)
764 if (image->matte != MagickFalse)
765 alpha=(MagickRealType) (QuantumScale*(GetAlphaPixelComponent(p)));
766 if (reconstruct_image->matte != MagickFalse)
767 beta=(MagickRealType) (QuantumScale*GetAlphaPixelComponent(q));
769 if ((channel & RedChannel) != 0)
771 distance=fabs(alpha*p->red-beta*q->red);
772 distortion[RedChannel]+=distance;
773 distortion[AllChannels]+=distance;
774 mean_error+=distance*distance;
775 if (distance > maximum_error)
776 maximum_error=distance;
779 if ((channel & GreenChannel) != 0)
781 distance=fabs(alpha*p->green-beta*q->green);
782 distortion[GreenChannel]+=distance;
783 distortion[AllChannels]+=distance;
784 mean_error+=distance*distance;
785 if (distance > maximum_error)
786 maximum_error=distance;
789 if ((channel & BlueChannel) != 0)
791 distance=fabs(alpha*p->blue-beta*q->blue);
792 distortion[BlueChannel]+=distance;
793 distortion[AllChannels]+=distance;
794 mean_error+=distance*distance;
795 if (distance > maximum_error)
796 maximum_error=distance;
799 if (((channel & OpacityChannel) != 0) &&
800 (image->matte != MagickFalse))
802 distance=fabs((double) p->opacity-q->opacity);
803 distortion[OpacityChannel]+=distance;
804 distortion[AllChannels]+=distance;
805 mean_error+=distance*distance;
806 if (distance > maximum_error)
807 maximum_error=distance;
810 if (((channel & IndexChannel) != 0) &&
811 (image->colorspace == CMYKColorspace) &&
812 (reconstruct_image->colorspace == CMYKColorspace))
814 distance=fabs(alpha*indexes[x]-beta*reconstruct_indexes[x]);
815 distortion[BlackChannel]+=distance;
816 distortion[AllChannels]+=distance;
817 mean_error+=distance*distance;
818 if (distance > maximum_error)
819 maximum_error=distance;
826 reconstruct_view=DestroyCacheView(reconstruct_view);
827 image_view=DestroyCacheView(image_view);
828 image->error.mean_error_per_pixel=distortion[AllChannels]/area;
829 image->error.normalized_mean_error=QuantumScale*QuantumScale*mean_error/area;
830 image->error.normalized_maximum_error=QuantumScale*maximum_error;
834 static MagickBooleanType GetMeanSquaredDistortion(const Image *image,
835 const Image *reconstruct_image,const ChannelType channel,
836 double *distortion,ExceptionInfo *exception)
852 image_view=AcquireCacheView(image);
853 reconstruct_view=AcquireCacheView(reconstruct_image);
854 #if defined(MAGICKCORE_OPENMP_SUPPORT)
855 #pragma omp parallel for schedule(dynamic,4) shared(status)
857 for (y=0; y < (ssize_t) image->rows; y++)
860 channel_distortion[AllChannels+1];
862 register const IndexPacket
864 *restrict reconstruct_indexes;
866 register const PixelPacket
874 if (status == MagickFalse)
876 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
877 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,
878 reconstruct_image->columns,1,exception);
879 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
884 indexes=GetCacheViewVirtualIndexQueue(image_view);
885 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
886 (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
887 for (x=0; x < (ssize_t) image->columns; x++)
892 if ((channel & RedChannel) != 0)
894 distance=QuantumScale*(p->red-(MagickRealType) q->red);
895 channel_distortion[RedChannel]+=distance*distance;
896 channel_distortion[AllChannels]+=distance*distance;
898 if ((channel & GreenChannel) != 0)
900 distance=QuantumScale*(p->green-(MagickRealType) q->green);
901 channel_distortion[GreenChannel]+=distance*distance;
902 channel_distortion[AllChannels]+=distance*distance;
904 if ((channel & BlueChannel) != 0)
906 distance=QuantumScale*(p->blue-(MagickRealType) q->blue);
907 channel_distortion[BlueChannel]+=distance*distance;
908 channel_distortion[AllChannels]+=distance*distance;
910 if (((channel & OpacityChannel) != 0) &&
911 (image->matte != MagickFalse))
913 distance=QuantumScale*(p->opacity-(MagickRealType) q->opacity);
914 channel_distortion[OpacityChannel]+=distance*distance;
915 channel_distortion[AllChannels]+=distance*distance;
917 if (((channel & IndexChannel) != 0) &&
918 (image->colorspace == CMYKColorspace) &&
919 (reconstruct_image->colorspace == CMYKColorspace))
921 distance=QuantumScale*(indexes[x]-(MagickRealType)
922 reconstruct_indexes[x]);
923 channel_distortion[BlackChannel]+=distance*distance;
924 channel_distortion[AllChannels]+=distance*distance;
929 #if defined(MAGICKCORE_OPENMP_SUPPORT)
930 #pragma omp critical (MagickCore_GetMeanSquaredError)
932 for (i=0; i <= (ssize_t) AllChannels; i++)
933 distortion[i]+=channel_distortion[i];
935 reconstruct_view=DestroyCacheView(reconstruct_view);
936 image_view=DestroyCacheView(image_view);
937 for (i=0; i <= (ssize_t) AllChannels; i++)
938 distortion[i]/=((double) image->columns*image->rows);
939 distortion[AllChannels]/=(double) GetNumberChannels(image,channel);
943 static MagickBooleanType GetNormalizedCrossCorrelationDistortion(
944 const Image *image,const Image *reconstruct_image,const ChannelType channel,
945 double *distortion,ExceptionInfo *exception)
947 #define SimilarityImageTag "Similarity/Image"
955 *reconstruct_statistics;
973 Normalize to account for variation due to lighting and exposure condition.
975 image_statistics=GetImageChannelStatistics(image,exception);
976 reconstruct_statistics=GetImageChannelStatistics(reconstruct_image,exception);
979 for (i=0; i <= (ssize_t) AllChannels; i++)
981 area=1.0/((MagickRealType) image->columns*image->rows);
982 image_view=AcquireCacheView(image);
983 reconstruct_view=AcquireCacheView(reconstruct_image);
984 for (y=0; y < (ssize_t) image->rows; y++)
986 register const IndexPacket
988 *restrict reconstruct_indexes;
990 register const PixelPacket
997 if (status == MagickFalse)
999 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1000 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
1002 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1007 indexes=GetCacheViewVirtualIndexQueue(image_view);
1008 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
1009 for (x=0; x < (ssize_t) image->columns; x++)
1011 if ((channel & RedChannel) != 0)
1012 distortion[RedChannel]+=area*QuantumScale*(p->red-
1013 image_statistics[RedChannel].mean)*(q->red-
1014 reconstruct_statistics[RedChannel].mean);
1015 if ((channel & GreenChannel) != 0)
1016 distortion[GreenChannel]+=area*QuantumScale*(p->green-
1017 image_statistics[GreenChannel].mean)*(q->green-
1018 reconstruct_statistics[GreenChannel].mean);
1019 if ((channel & BlueChannel) != 0)
1020 distortion[BlueChannel]+=area*QuantumScale*(p->blue-
1021 image_statistics[BlueChannel].mean)*(q->blue-
1022 reconstruct_statistics[BlueChannel].mean);
1023 if (((channel & OpacityChannel) != 0) &&
1024 (image->matte != MagickFalse))
1025 distortion[OpacityChannel]+=area*QuantumScale*(p->opacity-
1026 image_statistics[OpacityChannel].mean)*(q->opacity-
1027 reconstruct_statistics[OpacityChannel].mean);
1028 if (((channel & IndexChannel) != 0) &&
1029 (image->colorspace == CMYKColorspace) &&
1030 (reconstruct_image->colorspace == CMYKColorspace))
1031 distortion[BlackChannel]+=area*QuantumScale*(indexes[x]-
1032 image_statistics[OpacityChannel].mean)*(reconstruct_indexes[x]-
1033 reconstruct_statistics[OpacityChannel].mean);
1037 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1042 proceed=SetImageProgress(image,SimilarityImageTag,progress++,
1044 if (proceed == MagickFalse)
1048 reconstruct_view=DestroyCacheView(reconstruct_view);
1049 image_view=DestroyCacheView(image_view);
1051 Divide by the standard deviation.
1053 for (i=0; i < (ssize_t) AllChannels; i++)
1058 gamma=image_statistics[i].standard_deviation*
1059 reconstruct_statistics[i].standard_deviation;
1060 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1061 distortion[i]=QuantumRange*gamma*distortion[i];
1063 distortion[AllChannels]=0.0;
1064 if ((channel & RedChannel) != 0)
1065 distortion[AllChannels]+=distortion[RedChannel]*distortion[RedChannel];
1066 if ((channel & GreenChannel) != 0)
1067 distortion[AllChannels]+=distortion[GreenChannel]*distortion[GreenChannel];
1068 if ((channel & BlueChannel) != 0)
1069 distortion[AllChannels]+=distortion[BlueChannel]*distortion[BlueChannel];
1070 if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
1071 distortion[AllChannels]+=distortion[OpacityChannel]*
1072 distortion[OpacityChannel];
1073 if (((channel & IndexChannel) != 0) &&
1074 (image->colorspace == CMYKColorspace))
1075 distortion[AllChannels]+=distortion[BlackChannel]*distortion[BlackChannel];
1076 distortion[AllChannels]=sqrt(distortion[AllChannels]/GetNumberChannels(image,
1081 reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1082 reconstruct_statistics);
1083 image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1088 static MagickBooleanType GetPeakAbsoluteDistortion(const Image *image,
1089 const Image *reconstruct_image,const ChannelType channel,
1090 double *distortion,ExceptionInfo *exception)
1103 image_view=AcquireCacheView(image);
1104 reconstruct_view=AcquireCacheView(reconstruct_image);
1105 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1106 #pragma omp parallel for schedule(dynamic,4) shared(status)
1108 for (y=0; y < (ssize_t) image->rows; y++)
1111 channel_distortion[AllChannels+1];
1113 register const IndexPacket
1115 *restrict reconstruct_indexes;
1117 register const PixelPacket
1125 if (status == MagickFalse)
1127 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1128 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,
1129 reconstruct_image->columns,1,exception);
1130 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1135 indexes=GetCacheViewVirtualIndexQueue(image_view);
1136 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
1137 (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
1138 for (x=0; x < (ssize_t) image->columns; x++)
1143 if ((channel & RedChannel) != 0)
1145 distance=QuantumScale*fabs(p->red-(double) q->red);
1146 if (distance > channel_distortion[RedChannel])
1147 channel_distortion[RedChannel]=distance;
1148 if (distance > channel_distortion[AllChannels])
1149 channel_distortion[AllChannels]=distance;
1151 if ((channel & GreenChannel) != 0)
1153 distance=QuantumScale*fabs(p->green-(double) q->green);
1154 if (distance > channel_distortion[GreenChannel])
1155 channel_distortion[GreenChannel]=distance;
1156 if (distance > channel_distortion[AllChannels])
1157 channel_distortion[AllChannels]=distance;
1159 if ((channel & BlueChannel) != 0)
1161 distance=QuantumScale*fabs(p->blue-(double) q->blue);
1162 if (distance > channel_distortion[BlueChannel])
1163 channel_distortion[BlueChannel]=distance;
1164 if (distance > channel_distortion[AllChannels])
1165 channel_distortion[AllChannels]=distance;
1167 if (((channel & OpacityChannel) != 0) &&
1168 (image->matte != MagickFalse))
1170 distance=QuantumScale*fabs(p->opacity-(double) q->opacity);
1171 if (distance > channel_distortion[OpacityChannel])
1172 channel_distortion[OpacityChannel]=distance;
1173 if (distance > channel_distortion[AllChannels])
1174 channel_distortion[AllChannels]=distance;
1176 if (((channel & IndexChannel) != 0) &&
1177 (image->colorspace == CMYKColorspace) &&
1178 (reconstruct_image->colorspace == CMYKColorspace))
1180 distance=QuantumScale*fabs(indexes[x]-(double)
1181 reconstruct_indexes[x]);
1182 if (distance > channel_distortion[BlackChannel])
1183 channel_distortion[BlackChannel]=distance;
1184 if (distance > channel_distortion[AllChannels])
1185 channel_distortion[AllChannels]=distance;
1190 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1191 #pragma omp critical (MagickCore_GetPeakAbsoluteError)
1193 for (i=0; i <= (ssize_t) AllChannels; i++)
1194 if (channel_distortion[i] > distortion[i])
1195 distortion[i]=channel_distortion[i];
1197 reconstruct_view=DestroyCacheView(reconstruct_view);
1198 image_view=DestroyCacheView(image_view);
1202 static MagickBooleanType GetPeakSignalToNoiseRatio(const Image *image,
1203 const Image *reconstruct_image,const ChannelType channel,
1204 double *distortion,ExceptionInfo *exception)
1209 status=GetMeanSquaredDistortion(image,reconstruct_image,channel,distortion,
1211 if ((channel & RedChannel) != 0)
1212 distortion[RedChannel]=20.0*log10((double) 1.0/sqrt(
1213 distortion[RedChannel]));
1214 if ((channel & GreenChannel) != 0)
1215 distortion[GreenChannel]=20.0*log10((double) 1.0/sqrt(
1216 distortion[GreenChannel]));
1217 if ((channel & BlueChannel) != 0)
1218 distortion[BlueChannel]=20.0*log10((double) 1.0/sqrt(
1219 distortion[BlueChannel]));
1220 if (((channel & OpacityChannel) != 0) &&
1221 (image->matte != MagickFalse))
1222 distortion[OpacityChannel]=20.0*log10((double) 1.0/sqrt(
1223 distortion[OpacityChannel]));
1224 if (((channel & IndexChannel) != 0) &&
1225 (image->colorspace == CMYKColorspace))
1226 distortion[BlackChannel]=20.0*log10((double) 1.0/sqrt(
1227 distortion[BlackChannel]));
1228 distortion[AllChannels]=20.0*log10((double) 1.0/sqrt(
1229 distortion[AllChannels]));
1233 static MagickBooleanType GetRootMeanSquaredDistortion(const Image *image,
1234 const Image *reconstruct_image,const ChannelType channel,
1235 double *distortion,ExceptionInfo *exception)
1240 status=GetMeanSquaredDistortion(image,reconstruct_image,channel,distortion,
1242 if ((channel & RedChannel) != 0)
1243 distortion[RedChannel]=sqrt(distortion[RedChannel]);
1244 if ((channel & GreenChannel) != 0)
1245 distortion[GreenChannel]=sqrt(distortion[GreenChannel]);
1246 if ((channel & BlueChannel) != 0)
1247 distortion[BlueChannel]=sqrt(distortion[BlueChannel]);
1248 if (((channel & OpacityChannel) != 0) &&
1249 (image->matte != MagickFalse))
1250 distortion[OpacityChannel]=sqrt(distortion[OpacityChannel]);
1251 if (((channel & IndexChannel) != 0) &&
1252 (image->colorspace == CMYKColorspace))
1253 distortion[BlackChannel]=sqrt(distortion[BlackChannel]);
1254 distortion[AllChannels]=sqrt(distortion[AllChannels]);
1258 MagickExport MagickBooleanType GetImageChannelDistortion(Image *image,
1259 const Image *reconstruct_image,const ChannelType channel,
1260 const MetricType metric,double *distortion,ExceptionInfo *exception)
1263 *channel_distortion;
1271 assert(image != (Image *) NULL);
1272 assert(image->signature == MagickSignature);
1273 if (image->debug != MagickFalse)
1274 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1275 assert(reconstruct_image != (const Image *) NULL);
1276 assert(reconstruct_image->signature == MagickSignature);
1277 assert(distortion != (double *) NULL);
1279 if (image->debug != MagickFalse)
1280 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1281 if ((reconstruct_image->columns != image->columns) ||
1282 (reconstruct_image->rows != image->rows))
1283 ThrowBinaryException(ImageError,"ImageSizeDiffers",image->filename);
1285 Get image distortion.
1287 length=AllChannels+1UL;
1288 channel_distortion=(double *) AcquireQuantumMemory(length,
1289 sizeof(*channel_distortion));
1290 if (channel_distortion == (double *) NULL)
1291 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1292 (void) ResetMagickMemory(channel_distortion,0,length*
1293 sizeof(*channel_distortion));
1296 case AbsoluteErrorMetric:
1298 status=GetAbsoluteDistortion(image,reconstruct_image,channel,
1299 channel_distortion,exception);
1302 case FuzzErrorMetric:
1304 status=GetFuzzDistortion(image,reconstruct_image,channel,
1305 channel_distortion,exception);
1308 case MeanAbsoluteErrorMetric:
1310 status=GetMeanAbsoluteDistortion(image,reconstruct_image,channel,
1311 channel_distortion,exception);
1314 case MeanErrorPerPixelMetric:
1316 status=GetMeanErrorPerPixel(image,reconstruct_image,channel,
1317 channel_distortion,exception);
1320 case MeanSquaredErrorMetric:
1322 status=GetMeanSquaredDistortion(image,reconstruct_image,channel,
1323 channel_distortion,exception);
1326 case NormalizedCrossCorrelationErrorMetric:
1329 status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1330 channel,channel_distortion,exception);
1333 case PeakAbsoluteErrorMetric:
1335 status=GetPeakAbsoluteDistortion(image,reconstruct_image,channel,
1336 channel_distortion,exception);
1339 case PeakSignalToNoiseRatioMetric:
1341 status=GetPeakSignalToNoiseRatio(image,reconstruct_image,channel,
1342 channel_distortion,exception);
1345 case RootMeanSquaredErrorMetric:
1347 status=GetRootMeanSquaredDistortion(image,reconstruct_image,channel,
1348 channel_distortion,exception);
1352 *distortion=channel_distortion[AllChannels];
1353 channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1358 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1362 % 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 %
1366 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1368 % GetImageChannelDistrortion() compares the image channels of an image to a
1369 % reconstructed image and returns the specified distortion metric for each
1372 % The format of the CompareImageChannels method is:
1374 % double *GetImageChannelDistortions(const Image *image,
1375 % const Image *reconstruct_image,const MetricType metric,
1376 % ExceptionInfo *exception)
1378 % A description of each parameter follows:
1380 % o image: the image.
1382 % o reconstruct_image: the reconstruct image.
1384 % o metric: the metric.
1386 % o exception: return any errors or warnings in this structure.
1389 MagickExport double *GetImageChannelDistortions(Image *image,
1390 const Image *reconstruct_image,const MetricType metric,
1391 ExceptionInfo *exception)
1394 *channel_distortion;
1402 assert(image != (Image *) NULL);
1403 assert(image->signature == MagickSignature);
1404 if (image->debug != MagickFalse)
1405 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1406 assert(reconstruct_image != (const Image *) NULL);
1407 assert(reconstruct_image->signature == MagickSignature);
1408 if (image->debug != MagickFalse)
1409 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1410 if ((reconstruct_image->columns != image->columns) ||
1411 (reconstruct_image->rows != image->rows))
1413 (void) ThrowMagickException(&image->exception,GetMagickModule(),
1414 ImageError,"ImageSizeDiffers","`%s'",image->filename);
1415 return((double *) NULL);
1418 Get image distortion.
1420 length=AllChannels+1UL;
1421 channel_distortion=(double *) AcquireQuantumMemory(length,
1422 sizeof(*channel_distortion));
1423 if (channel_distortion == (double *) NULL)
1424 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1425 (void) ResetMagickMemory(channel_distortion,0,length*
1426 sizeof(*channel_distortion));
1429 case AbsoluteErrorMetric:
1431 status=GetAbsoluteDistortion(image,reconstruct_image,AllChannels,
1432 channel_distortion,exception);
1435 case FuzzErrorMetric:
1437 status=GetFuzzDistortion(image,reconstruct_image,AllChannels,
1438 channel_distortion,exception);
1441 case MeanAbsoluteErrorMetric:
1443 status=GetMeanAbsoluteDistortion(image,reconstruct_image,AllChannels,
1444 channel_distortion,exception);
1447 case MeanErrorPerPixelMetric:
1449 status=GetMeanErrorPerPixel(image,reconstruct_image,AllChannels,
1450 channel_distortion,exception);
1453 case MeanSquaredErrorMetric:
1455 status=GetMeanSquaredDistortion(image,reconstruct_image,AllChannels,
1456 channel_distortion,exception);
1459 case NormalizedCrossCorrelationErrorMetric:
1462 status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1463 AllChannels,channel_distortion,exception);
1466 case PeakAbsoluteErrorMetric:
1468 status=GetPeakAbsoluteDistortion(image,reconstruct_image,AllChannels,
1469 channel_distortion,exception);
1472 case PeakSignalToNoiseRatioMetric:
1474 status=GetPeakSignalToNoiseRatio(image,reconstruct_image,AllChannels,
1475 channel_distortion,exception);
1478 case RootMeanSquaredErrorMetric:
1480 status=GetRootMeanSquaredDistortion(image,reconstruct_image,AllChannels,
1481 channel_distortion,exception);
1485 return(channel_distortion);
1489 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1493 % I s I m a g e s E q u a l %
1497 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1499 % IsImagesEqual() measures the difference between colors at each pixel
1500 % location of two images. A value other than 0 means the colors match
1501 % exactly. Otherwise an error measure is computed by summing over all
1502 % pixels in an image the distance squared in RGB space between each image
1503 % pixel and its corresponding pixel in the reconstruct image. The error
1504 % measure is assigned to these image members:
1506 % o mean_error_per_pixel: The mean error for any single pixel in
1509 % o normalized_mean_error: The normalized mean quantization error for
1510 % any single pixel in the image. This distance measure is normalized to
1511 % a range between 0 and 1. It is independent of the range of red, green,
1512 % and blue values in the image.
1514 % o normalized_maximum_error: The normalized maximum quantization
1515 % error for any single pixel in the image. This distance measure is
1516 % normalized to a range between 0 and 1. It is independent of the range
1517 % of red, green, and blue values in your image.
1519 % A small normalized mean square error, accessed as
1520 % image->normalized_mean_error, suggests the images are very similar in
1521 % spatial layout and color.
1523 % The format of the IsImagesEqual method is:
1525 % MagickBooleanType IsImagesEqual(Image *image,
1526 % const Image *reconstruct_image)
1528 % A description of each parameter follows.
1530 % o image: the image.
1532 % o reconstruct_image: the reconstruct image.
1535 MagickExport MagickBooleanType IsImagesEqual(Image *image,
1536 const Image *reconstruct_image)
1555 mean_error_per_pixel;
1557 assert(image != (Image *) NULL);
1558 assert(image->signature == MagickSignature);
1559 assert(reconstruct_image != (const Image *) NULL);
1560 assert(reconstruct_image->signature == MagickSignature);
1561 if ((reconstruct_image->columns != image->columns) ||
1562 (reconstruct_image->rows != image->rows))
1563 ThrowBinaryException(ImageError,"ImageSizeDiffers",image->filename);
1566 mean_error_per_pixel=0.0;
1568 exception=(&image->exception);
1569 image_view=AcquireCacheView(image);
1570 reconstruct_view=AcquireCacheView(reconstruct_image);
1571 for (y=0; y < (ssize_t) image->rows; y++)
1573 register const IndexPacket
1575 *restrict reconstruct_indexes;
1577 register const PixelPacket
1584 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1585 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
1587 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1589 indexes=GetCacheViewVirtualIndexQueue(image_view);
1590 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
1591 for (x=0; x < (ssize_t) image->columns; x++)
1596 distance=fabs(p->red-(double) q->red);
1597 mean_error_per_pixel+=distance;
1598 mean_error+=distance*distance;
1599 if (distance > maximum_error)
1600 maximum_error=distance;
1602 distance=fabs(p->green-(double) q->green);
1603 mean_error_per_pixel+=distance;
1604 mean_error+=distance*distance;
1605 if (distance > maximum_error)
1606 maximum_error=distance;
1608 distance=fabs(p->blue-(double) q->blue);
1609 mean_error_per_pixel+=distance;
1610 mean_error+=distance*distance;
1611 if (distance > maximum_error)
1612 maximum_error=distance;
1614 if (image->matte != MagickFalse)
1616 distance=fabs(p->opacity-(double) q->opacity);
1617 mean_error_per_pixel+=distance;
1618 mean_error+=distance*distance;
1619 if (distance > maximum_error)
1620 maximum_error=distance;
1623 if ((image->colorspace == CMYKColorspace) &&
1624 (reconstruct_image->colorspace == CMYKColorspace))
1626 distance=fabs(indexes[x]-(double) reconstruct_indexes[x]);
1627 mean_error_per_pixel+=distance;
1628 mean_error+=distance*distance;
1629 if (distance > maximum_error)
1630 maximum_error=distance;
1637 reconstruct_view=DestroyCacheView(reconstruct_view);
1638 image_view=DestroyCacheView(image_view);
1639 image->error.mean_error_per_pixel=(double) (mean_error_per_pixel/area);
1640 image->error.normalized_mean_error=(double) (QuantumScale*QuantumScale*
1642 image->error.normalized_maximum_error=(double) (QuantumScale*maximum_error);
1643 status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
1648 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1652 % S i m i l a r i t y I m a g e %
1656 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1658 % SimilarityImage() compares the reference image of the image and returns the
1659 % best match offset. In addition, it returns a similarity image such that an
1660 % exact match location is completely white and if none of the pixels match,
1661 % black, otherwise some gray level in-between.
1663 % The format of the SimilarityImageImage method is:
1665 % Image *SimilarityImage(const Image *image,const Image *reference,
1666 % RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
1668 % A description of each parameter follows:
1670 % o image: the image.
1672 % o reference: find an area of the image that closely resembles this image.
1674 % o the best match offset of the reference image within the image.
1676 % o similarity: the computed similarity between the images.
1678 % o exception: return any errors or warnings in this structure.
1682 static double GetNCCDistortion(const Image *image,
1683 const Image *reconstruct_image,
1684 const ChannelStatistics *reconstruct_statistics,ExceptionInfo *exception)
1686 #define SimilarityImageTag "Similarity/Image"
1712 Normalize to account for variation due to lighting and exposure condition.
1714 image_statistics=GetImageChannelStatistics(image,exception);
1717 area=1.0/((MagickRealType) image->columns*image->rows);
1718 image_view=AcquireCacheView(image);
1719 reconstruct_view=AcquireCacheView(reconstruct_image);
1720 for (y=0; y < (ssize_t) image->rows; y++)
1722 register const IndexPacket
1724 *restrict reconstruct_indexes;
1726 register const PixelPacket
1733 if (status == MagickFalse)
1735 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1736 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
1738 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1743 indexes=GetCacheViewVirtualIndexQueue(image_view);
1744 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
1745 for (x=0; x < (ssize_t) image->columns; x++)
1747 distortion+=area*QuantumScale*(p->red-
1748 image_statistics[RedChannel].mean)*(q->red-
1749 reconstruct_statistics[RedChannel].mean);
1750 distortion+=area*QuantumScale*(p->green-
1751 image_statistics[GreenChannel].mean)*(q->green-
1752 reconstruct_statistics[GreenChannel].mean);
1753 distortion+=area*QuantumScale*(p->blue-
1754 image_statistics[BlueChannel].mean)*(q->blue-
1755 reconstruct_statistics[BlueChannel].mean);
1756 if (image->matte != MagickFalse)
1757 distortion+=area*QuantumScale*(p->opacity-
1758 image_statistics[OpacityChannel].mean)*(q->opacity-
1759 reconstruct_statistics[OpacityChannel].mean);
1760 if ((image->colorspace == CMYKColorspace) &&
1761 (reconstruct_image->colorspace == CMYKColorspace))
1762 distortion+=area*QuantumScale*(indexes[x]-
1763 image_statistics[OpacityChannel].mean)*(reconstruct_indexes[x]-
1764 reconstruct_statistics[OpacityChannel].mean);
1769 reconstruct_view=DestroyCacheView(reconstruct_view);
1770 image_view=DestroyCacheView(image_view);
1772 Divide by the standard deviation.
1774 gamma=image_statistics[AllChannels].standard_deviation*
1775 reconstruct_statistics[AllChannels].standard_deviation;
1776 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1777 distortion=QuantumRange*gamma*distortion;
1779 if (image->matte != MagickFalse)
1781 if (image->colorspace == CMYKColorspace)
1783 distortion=sqrt(distortion/number_channels);
1787 image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1789 return(1.0-distortion);
1792 static double GetSimilarityMetric(const Image *image,const Image *reference,
1793 const ChannelStatistics *reference_statistics,const ssize_t x_offset,
1794 const ssize_t y_offset,ExceptionInfo *exception)
1805 SetGeometry(reference,&geometry);
1806 geometry.x=x_offset;
1807 geometry.y=y_offset;
1808 similarity_image=CropImage(image,&geometry,exception);
1809 if (similarity_image == (Image *) NULL)
1811 distortion=GetNCCDistortion(reference,similarity_image,reference_statistics,
1813 similarity_image=DestroyImage(similarity_image);
1817 MagickExport Image *SimilarityImage(Image *image,const Image *reference,
1818 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
1820 #define SimilarityImageTag "Similarity/Image"
1826 *reference_statistics;
1840 assert(image != (const Image *) NULL);
1841 assert(image->signature == MagickSignature);
1842 if (image->debug != MagickFalse)
1843 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1844 assert(exception != (ExceptionInfo *) NULL);
1845 assert(exception->signature == MagickSignature);
1846 assert(offset != (RectangleInfo *) NULL);
1847 SetGeometry(reference,offset);
1848 *similarity_metric=1.0;
1849 if ((reference->columns > image->columns) || (reference->rows > image->rows))
1850 ThrowImageException(ImageError,"ImageSizeDiffers");
1851 similarity_image=CloneImage(image,image->columns-reference->columns+1,
1852 image->rows-reference->rows+1,MagickTrue,exception);
1853 if (similarity_image == (Image *) NULL)
1854 return((Image *) NULL);
1855 if (SetImageStorageClass(similarity_image,DirectClass) == MagickFalse)
1857 InheritException(exception,&similarity_image->exception);
1858 similarity_image=DestroyImage(similarity_image);
1859 return((Image *) NULL);
1862 Measure similarity of reference image against image.
1866 reference_statistics=GetImageChannelStatistics(reference,exception);
1867 similarity_view=AcquireCacheView(similarity_image);
1868 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1869 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1871 for (y=0; y < (ssize_t) (image->rows-reference->rows+1); y++)
1879 register PixelPacket
1882 if (status == MagickFalse)
1884 q=GetCacheViewAuthenticPixels(similarity_view,0,y,similarity_image->columns,
1886 if (q == (const PixelPacket *) NULL)
1891 for (x=0; x < (ssize_t) (image->columns-reference->columns+1); x++)
1893 similarity=GetSimilarityMetric(image,reference,reference_statistics,x,y,
1895 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1896 #pragma omp critical (MagickCore_SimilarityImage)
1898 if (similarity < *similarity_metric)
1900 *similarity_metric=similarity;
1904 q->red=ClampToQuantum(QuantumRange-QuantumRange*similarity);
1909 if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
1911 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1916 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1917 #pragma omp critical (MagickCore_SimilarityImage)
1919 proceed=SetImageProgress(image,SimilarityImageTag,progress++,
1921 if (proceed == MagickFalse)
1925 similarity_view=DestroyCacheView(similarity_view);
1926 reference_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1927 reference_statistics);
1928 return(similarity_image);