]> granicus.if.org Git - imagemagick/blob - MagickCore/compare.c
(no commit message)
[imagemagick] / MagickCore / compare.c
1 /*
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %                                                                             %
4 %                                                                             %
5 %                                                                             %
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               %
11 %                                                                             %
12 %                                                                             %
13 %                      MagickCore Image Comparison Methods                    %
14 %                                                                             %
15 %                              Software Design                                %
16 %                                  Cristy                                     %
17 %                               December 2003                                 %
18 %                                                                             %
19 %                                                                             %
20 %  Copyright 1999-2014 ImageMagick Studio LLC, a non-profit organization      %
21 %  dedicated to making software imaging solutions freely available.           %
22 %                                                                             %
23 %  You may not use this file except in compliance with the License.  You may  %
24 %  obtain a copy of the License at                                            %
25 %                                                                             %
26 %    http://www.imagemagick.org/script/license.php                            %
27 %                                                                             %
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.                                             %
33 %                                                                             %
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35 %
36 %
37 %
38 */
39 \f
40 /*
41   Include declarations.
42 */
43 #include "MagickCore/studio.h"
44 #include "MagickCore/artifact.h"
45 #include "MagickCore/attribute.h"
46 #include "MagickCore/cache-view.h"
47 #include "MagickCore/channel.h"
48 #include "MagickCore/client.h"
49 #include "MagickCore/color.h"
50 #include "MagickCore/color-private.h"
51 #include "MagickCore/colorspace.h"
52 #include "MagickCore/colorspace-private.h"
53 #include "MagickCore/compare.h"
54 #include "MagickCore/composite-private.h"
55 #include "MagickCore/constitute.h"
56 #include "MagickCore/exception-private.h"
57 #include "MagickCore/geometry.h"
58 #include "MagickCore/image-private.h"
59 #include "MagickCore/list.h"
60 #include "MagickCore/log.h"
61 #include "MagickCore/memory_.h"
62 #include "MagickCore/monitor.h"
63 #include "MagickCore/monitor-private.h"
64 #include "MagickCore/option.h"
65 #include "MagickCore/pixel-accessor.h"
66 #include "MagickCore/property.h"
67 #include "MagickCore/resource_.h"
68 #include "MagickCore/string_.h"
69 #include "MagickCore/statistic.h"
70 #include "MagickCore/thread-private.h"
71 #include "MagickCore/transform.h"
72 #include "MagickCore/utility.h"
73 #include "MagickCore/version.h"
74 \f
75 /*
76 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
77 %                                                                             %
78 %                                                                             %
79 %                                                                             %
80 %   C o m p a r e I m a g e                                                   %
81 %                                                                             %
82 %                                                                             %
83 %                                                                             %
84 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
85 %
86 %  CompareImages() compares one or more pixel channels of an image to a
87 %  reconstructed image and returns the difference image.
88 %
89 %  The format of the CompareImages method is:
90 %
91 %      Image *CompareImages(const Image *image,const Image *reconstruct_image,
92 %        const MetricType metric,double *distortion,ExceptionInfo *exception)
93 %
94 %  A description of each parameter follows:
95 %
96 %    o image: the image.
97 %
98 %    o reconstruct_image: the reconstruct image.
99 %
100 %    o metric: the metric.
101 %
102 %    o distortion: the computed distortion between the images.
103 %
104 %    o exception: return any errors or warnings in this structure.
105 %
106 */
107 MagickExport Image *CompareImages(Image *image,const Image *reconstruct_image,
108   const MetricType metric,double *distortion,ExceptionInfo *exception)
109 {
110   CacheView
111     *highlight_view,
112     *image_view,
113     *reconstruct_view;
114
115   const char
116     *artifact;
117
118   Image
119     *difference_image,
120     *highlight_image;
121
122   MagickBooleanType
123     status;
124
125   PixelInfo
126     highlight,
127     lowlight;
128
129   ssize_t
130     y;
131
132   assert(image != (Image *) NULL);
133   assert(image->signature == MagickSignature);
134   if (image->debug != MagickFalse)
135     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
136   assert(reconstruct_image != (const Image *) NULL);
137   assert(reconstruct_image->signature == MagickSignature);
138   assert(distortion != (double *) NULL);
139   *distortion=0.0;
140   if (image->debug != MagickFalse)
141     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
142   if (metric != PerceptualHashErrorMetric)
143     if ((reconstruct_image->columns != image->columns) ||
144         (reconstruct_image->rows != image->rows))
145       ThrowImageException(ImageError,"ImageSizeDiffers");
146   status=GetImageDistortion(image,reconstruct_image,metric,distortion,
147     exception);
148   if (status == MagickFalse)
149     return((Image *) NULL);
150   difference_image=CloneImage(image,0,0,MagickTrue,exception);
151   if (difference_image == (Image *) NULL)
152     return((Image *) NULL);
153   (void) SetImageAlphaChannel(difference_image,OpaqueAlphaChannel,exception);
154   highlight_image=CloneImage(image,image->columns,image->rows,MagickTrue,
155     exception);
156   if (highlight_image == (Image *) NULL)
157     {
158       difference_image=DestroyImage(difference_image);
159       return((Image *) NULL);
160     }
161   status=SetImageStorageClass(highlight_image,DirectClass,exception);
162   if (status == MagickFalse)
163     {
164       difference_image=DestroyImage(difference_image);
165       highlight_image=DestroyImage(highlight_image);
166       return((Image *) NULL);
167     }
168   (void) SetImageAlphaChannel(highlight_image,OpaqueAlphaChannel,exception);
169   (void) QueryColorCompliance("#f1001ecc",AllCompliance,&highlight,exception);
170   artifact=GetImageArtifact(image,"highlight-color");
171   if (artifact != (const char *) NULL)
172     (void) QueryColorCompliance(artifact,AllCompliance,&highlight,exception);
173   (void) QueryColorCompliance("#ffffffcc",AllCompliance,&lowlight,exception);
174   artifact=GetImageArtifact(image,"lowlight-color");
175   if (artifact != (const char *) NULL)
176     (void) QueryColorCompliance(artifact,AllCompliance,&lowlight,exception);
177   /*
178     Generate difference image.
179   */
180   status=MagickTrue;
181   image_view=AcquireVirtualCacheView(image,exception);
182   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
183   highlight_view=AcquireAuthenticCacheView(highlight_image,exception);
184 #if defined(MAGICKCORE_OPENMP_SUPPORT)
185   #pragma omp parallel for schedule(static,4) shared(status) \
186     magick_threads(image,highlight_image,image->rows,1)
187 #endif
188   for (y=0; y < (ssize_t) image->rows; y++)
189   {
190     MagickBooleanType
191       sync;
192
193     register const Quantum
194       *restrict p,
195       *restrict q;
196
197     register Quantum
198       *restrict r;
199
200     register ssize_t
201       x;
202
203     if (status == MagickFalse)
204       continue;
205     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
206     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,image->columns,1,
207       exception);
208     r=QueueCacheViewAuthenticPixels(highlight_view,0,y,highlight_image->columns,
209       1,exception);
210     if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL) ||
211         (r == (Quantum *) NULL))
212       {
213         status=MagickFalse;
214         continue;
215       }
216     for (x=0; x < (ssize_t) image->columns; x++)
217     {
218       double
219         Da,
220         Sa;
221
222       MagickStatusType
223         difference;
224
225       register ssize_t
226         i;
227
228       if (GetPixelReadMask(image,p) == 0)
229         {
230           SetPixelInfoPixel(highlight_image,&lowlight,r);
231           p+=GetPixelChannels(image);
232           q+=GetPixelChannels(reconstruct_image);
233           r+=GetPixelChannels(highlight_image);
234           continue;
235         }
236       difference=MagickFalse;
237       Sa=QuantumScale*GetPixelAlpha(image,p);
238       Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
239       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
240       {
241         double
242           distance;
243
244         PixelChannel channel=GetPixelChannelChannel(image,i);
245         PixelTrait traits=GetPixelChannelTraits(image,channel);
246         PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
247           channel);
248         if ((traits == UndefinedPixelTrait) ||
249             (reconstruct_traits == UndefinedPixelTrait) ||
250             ((reconstruct_traits & UpdatePixelTrait) == 0))
251           continue;
252         distance=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
253         if (fabs((double) distance) >= MagickEpsilon)
254           difference=MagickTrue;
255       }
256       if (difference == MagickFalse)
257         SetPixelInfoPixel(highlight_image,&lowlight,r);
258       else
259         SetPixelInfoPixel(highlight_image,&highlight,r);
260       p+=GetPixelChannels(image);
261       q+=GetPixelChannels(reconstruct_image);
262       r+=GetPixelChannels(highlight_image);
263     }
264     sync=SyncCacheViewAuthenticPixels(highlight_view,exception);
265     if (sync == MagickFalse)
266       status=MagickFalse;
267   }
268   highlight_view=DestroyCacheView(highlight_view);
269   reconstruct_view=DestroyCacheView(reconstruct_view);
270   image_view=DestroyCacheView(image_view);
271   (void) CompositeImage(difference_image,highlight_image,image->compose,
272     MagickTrue,0,0,exception);
273   highlight_image=DestroyImage(highlight_image);
274   if (status == MagickFalse)
275     difference_image=DestroyImage(difference_image);
276   return(difference_image);
277 }
278 \f
279 /*
280 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
281 %                                                                             %
282 %                                                                             %
283 %                                                                             %
284 %   G e t I m a g e D i s t o r t i o n                                       %
285 %                                                                             %
286 %                                                                             %
287 %                                                                             %
288 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
289 %
290 %  GetImageDistortion() compares one or more pixel channels of an image to a
291 %  reconstructed image and returns the specified distortion metric.
292 %
293 %  The format of the GetImageDistortion method is:
294 %
295 %      MagickBooleanType GetImageDistortion(const Image *image,
296 %        const Image *reconstruct_image,const MetricType metric,
297 %        double *distortion,ExceptionInfo *exception)
298 %
299 %  A description of each parameter follows:
300 %
301 %    o image: the image.
302 %
303 %    o reconstruct_image: the reconstruct image.
304 %
305 %    o metric: the metric.
306 %
307 %    o distortion: the computed distortion between the images.
308 %
309 %    o exception: return any errors or warnings in this structure.
310 %
311 */
312
313 static inline double MagickMax(const double x,const double y)
314 {
315   if (x > y)
316     return(x);
317   return(y);
318 }
319
320 static MagickBooleanType GetAbsoluteDistortion(const Image *image,
321   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
322 {
323   CacheView
324     *image_view,
325     *reconstruct_view;
326
327   double
328     fuzz;
329
330   MagickBooleanType
331     status;
332
333   ssize_t
334     y;
335
336   /*
337     Compute the absolute difference in pixels between two images.
338   */
339   status=MagickTrue;
340   if (image->fuzz == 0.0)
341     fuzz=MagickMax(reconstruct_image->fuzz,MagickSQ1_2)*
342       MagickMax(reconstruct_image->fuzz,MagickSQ1_2);
343   else
344     if (reconstruct_image->fuzz == 0.0)
345       fuzz=MagickMax(image->fuzz,MagickSQ1_2)*
346         MagickMax(image->fuzz,MagickSQ1_2);
347     else
348       fuzz=MagickMax(image->fuzz,MagickSQ1_2)*
349         MagickMax(reconstruct_image->fuzz,MagickSQ1_2);
350   image_view=AcquireVirtualCacheView(image,exception);
351   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
352 #if defined(MAGICKCORE_OPENMP_SUPPORT)
353   #pragma omp parallel for schedule(static,4) shared(status) \
354     magick_threads(image,image,image->rows,1)
355 #endif
356   for (y=0; y < (ssize_t) image->rows; y++)
357   {
358     double
359       channel_distortion[MaxPixelChannels+1];
360
361     register const Quantum
362       *restrict p,
363       *restrict q;
364
365     register ssize_t
366       i,
367       x;
368
369     if (status == MagickFalse)
370       continue;
371     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
372     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
373       1,exception);
374     if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
375       {
376         status=MagickFalse;
377         continue;
378       }
379     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
380     for (x=0; x < (ssize_t) image->columns; x++)
381     {
382       double
383         Da,
384         Sa;
385
386       MagickBooleanType
387         difference;
388
389       register ssize_t
390         i;
391
392       if (GetPixelReadMask(image,p) == 0)
393         {
394           p+=GetPixelChannels(image);
395           q+=GetPixelChannels(reconstruct_image);
396           continue;
397         }
398       difference=MagickFalse;
399       Sa=QuantumScale*GetPixelAlpha(image,p);
400       Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
401       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
402       {
403         double
404           distance;
405
406         PixelChannel channel=GetPixelChannelChannel(image,i);
407         PixelTrait traits=GetPixelChannelTraits(image,channel);
408         PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
409           channel);
410         if ((traits == UndefinedPixelTrait) ||
411             (reconstruct_traits == UndefinedPixelTrait) ||
412             ((reconstruct_traits & UpdatePixelTrait) == 0))
413           continue;
414         distance=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
415         if ((distance*distance) > fuzz)
416           {
417             difference=MagickTrue;
418             break;
419           }
420       }
421       if (difference != MagickFalse)
422         {
423           channel_distortion[i]++;
424           channel_distortion[CompositePixelChannel]++;
425         }
426       p+=GetPixelChannels(image);
427       q+=GetPixelChannels(reconstruct_image);
428     }
429 #if defined(MAGICKCORE_OPENMP_SUPPORT)
430     #pragma omp critical (MagickCore_GetAbsoluteError)
431 #endif
432     for (i=0; i <= MaxPixelChannels; i++)
433       distortion[i]+=channel_distortion[i];
434   }
435   reconstruct_view=DestroyCacheView(reconstruct_view);
436   image_view=DestroyCacheView(image_view);
437   return(status);
438 }
439
440 static size_t GetImageChannels(const Image *image)
441 {
442   register ssize_t
443     i;
444
445   size_t
446     channels;
447
448   channels=0;
449   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
450   {
451     PixelChannel channel=GetPixelChannelChannel(image,i);
452     PixelTrait traits=GetPixelChannelTraits(image,channel);
453     if ((traits & UpdatePixelTrait) != 0)
454       channels++;
455   }
456   return(channels);
457 }
458
459 static MagickBooleanType GetFuzzDistortion(const Image *image,
460   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
461 {
462   CacheView
463     *image_view,
464     *reconstruct_view;
465
466   MagickBooleanType
467     status;
468
469   register ssize_t
470     i;
471
472   ssize_t
473     y;
474
475   status=MagickTrue;
476   image_view=AcquireVirtualCacheView(image,exception);
477   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
478 #if defined(MAGICKCORE_OPENMP_SUPPORT)
479   #pragma omp parallel for schedule(static,4) shared(status) \
480     magick_threads(image,image,image->rows,1)
481 #endif
482   for (y=0; y < (ssize_t) image->rows; y++)
483   {
484     double
485       channel_distortion[MaxPixelChannels+1];
486
487     register const Quantum
488       *restrict p,
489       *restrict q;
490
491     register ssize_t
492       i,
493       x;
494
495     if (status == MagickFalse)
496       continue;
497     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
498     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
499       1,exception);
500     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
501       {
502         status=MagickFalse;
503         continue;
504       }
505     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
506     for (x=0; x < (ssize_t) image->columns; x++)
507     {
508       double
509         Da,
510         Sa;
511
512       register ssize_t
513         i;
514
515       if (GetPixelReadMask(image,p) == 0)
516         {
517           p+=GetPixelChannels(image);
518           q+=GetPixelChannels(reconstruct_image);
519           continue;
520         }
521       Sa=QuantumScale*GetPixelAlpha(image,p);
522       Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
523       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
524       {
525         double
526           distance;
527
528         PixelChannel channel=GetPixelChannelChannel(image,i);
529         PixelTrait traits=GetPixelChannelTraits(image,channel);
530         PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
531           channel);
532         if ((traits == UndefinedPixelTrait) ||
533             (reconstruct_traits == UndefinedPixelTrait) ||
534             ((reconstruct_traits & UpdatePixelTrait) == 0))
535           continue;
536         distance=QuantumScale*(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
537           channel,q));
538         channel_distortion[i]+=distance*distance;
539         channel_distortion[CompositePixelChannel]+=distance*distance;
540       }
541       p+=GetPixelChannels(image);
542       q+=GetPixelChannels(reconstruct_image);
543     }
544 #if defined(MAGICKCORE_OPENMP_SUPPORT)
545     #pragma omp critical (MagickCore_GetFuzzDistortion)
546 #endif
547     for (i=0; i <= MaxPixelChannels; i++)
548       distortion[i]+=channel_distortion[i];
549   }
550   reconstruct_view=DestroyCacheView(reconstruct_view);
551   image_view=DestroyCacheView(image_view);
552   for (i=0; i <= MaxPixelChannels; i++)
553     distortion[i]/=((double) image->columns*image->rows);
554   distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
555   distortion[CompositePixelChannel]=sqrt(distortion[CompositePixelChannel]);
556   return(status);
557 }
558
559 static MagickBooleanType GetMeanAbsoluteDistortion(const Image *image,
560   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
561 {
562   CacheView
563     *image_view,
564     *reconstruct_view;
565
566   MagickBooleanType
567     status;
568
569   register ssize_t
570     i;
571
572   ssize_t
573     y;
574
575   status=MagickTrue;
576   image_view=AcquireVirtualCacheView(image,exception);
577   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
578 #if defined(MAGICKCORE_OPENMP_SUPPORT)
579   #pragma omp parallel for schedule(static,4) shared(status) \
580     magick_threads(image,image,image->rows,1)
581 #endif
582   for (y=0; y < (ssize_t) image->rows; y++)
583   {
584     double
585       channel_distortion[MaxPixelChannels+1];
586
587     register const Quantum
588       *restrict p,
589       *restrict q;
590
591     register ssize_t
592       i,
593       x;
594
595     if (status == MagickFalse)
596       continue;
597     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
598     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
599       1,exception);
600     if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
601       {
602         status=MagickFalse;
603         continue;
604       }
605     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
606     for (x=0; x < (ssize_t) image->columns; x++)
607     {
608       double
609         Da,
610         Sa;
611
612       register ssize_t
613         i;
614
615       if (GetPixelReadMask(image,p) == 0)
616         {
617           p+=GetPixelChannels(image);
618           q+=GetPixelChannels(reconstruct_image);
619           continue;
620         }
621       Sa=QuantumScale*GetPixelAlpha(image,p);
622       Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
623       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
624       {
625         double
626           distance;
627
628         PixelChannel channel=GetPixelChannelChannel(image,i);
629         PixelTrait traits=GetPixelChannelTraits(image,channel);
630         PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
631           channel);
632         if ((traits == UndefinedPixelTrait) ||
633             (reconstruct_traits == UndefinedPixelTrait) ||
634             ((reconstruct_traits & UpdatePixelTrait) == 0))
635           continue;
636         distance=QuantumScale*fabs(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
637           channel,q));
638         channel_distortion[i]+=distance;
639         channel_distortion[CompositePixelChannel]+=distance;
640       }
641       p+=GetPixelChannels(image);
642       q+=GetPixelChannels(reconstruct_image);
643     }
644 #if defined(MAGICKCORE_OPENMP_SUPPORT)
645     #pragma omp critical (MagickCore_GetMeanAbsoluteError)
646 #endif
647     for (i=0; i <= MaxPixelChannels; i++)
648       distortion[i]+=channel_distortion[i];
649   }
650   reconstruct_view=DestroyCacheView(reconstruct_view);
651   image_view=DestroyCacheView(image_view);
652   for (i=0; i <= MaxPixelChannels; i++)
653     distortion[i]/=((double) image->columns*image->rows);
654   distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
655   return(status);
656 }
657
658 static MagickBooleanType GetMeanErrorPerPixel(Image *image,
659   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
660 {
661   CacheView
662     *image_view,
663     *reconstruct_view;
664
665   MagickBooleanType
666     status;
667
668   double
669     area,
670     maximum_error,
671     mean_error;
672
673   ssize_t
674     y;
675
676   status=MagickTrue;
677   area=0.0;
678   maximum_error=0.0;
679   mean_error=0.0;
680   image_view=AcquireVirtualCacheView(image,exception);
681   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
682   for (y=0; y < (ssize_t) image->rows; y++)
683   {
684     register const Quantum
685       *restrict p,
686       *restrict q;
687
688     register ssize_t
689       x;
690
691     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
692     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
693       1,exception);
694     if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
695       {
696         status=MagickFalse;
697         break;
698       }
699     for (x=0; x < (ssize_t) image->columns; x++)
700     {
701       double
702         Da,
703         Sa;
704
705       register ssize_t
706         i;
707
708       if (GetPixelReadMask(image,p) == 0)
709         {
710           p+=GetPixelChannels(image);
711           q+=GetPixelChannels(reconstruct_image);
712           continue;
713         }
714       Sa=QuantumScale*GetPixelAlpha(image,p);
715       Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
716       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
717       {
718         double
719           distance;
720
721         PixelChannel channel=GetPixelChannelChannel(image,i);
722         PixelTrait traits=GetPixelChannelTraits(image,channel);
723         PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
724           channel);
725         if ((traits == UndefinedPixelTrait) ||
726             (reconstruct_traits == UndefinedPixelTrait) ||
727             ((reconstruct_traits & UpdatePixelTrait) == 0))
728           continue;
729         distance=fabs((double) (Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
730           channel,q)));
731         distortion[i]+=distance;
732         distortion[CompositePixelChannel]+=distance;
733         mean_error+=distance*distance;
734         if (distance > maximum_error)
735           maximum_error=distance;
736         area++;
737       }
738       p+=GetPixelChannels(image);
739       q+=GetPixelChannels(reconstruct_image);
740     }
741   }
742   reconstruct_view=DestroyCacheView(reconstruct_view);
743   image_view=DestroyCacheView(image_view);
744   image->error.mean_error_per_pixel=distortion[CompositePixelChannel]/area;
745   image->error.normalized_mean_error=QuantumScale*QuantumScale*mean_error/area;
746   image->error.normalized_maximum_error=QuantumScale*maximum_error;
747   return(status);
748 }
749
750 static MagickBooleanType GetMeanSquaredDistortion(const Image *image,
751   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
752 {
753   CacheView
754     *image_view,
755     *reconstruct_view;
756
757   MagickBooleanType
758     status;
759
760   register ssize_t
761     i;
762
763   ssize_t
764     y;
765
766   status=MagickTrue;
767   image_view=AcquireVirtualCacheView(image,exception);
768   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
769 #if defined(MAGICKCORE_OPENMP_SUPPORT)
770   #pragma omp parallel for schedule(static,4) shared(status) \
771     magick_threads(image,image,image->rows,1)
772 #endif
773   for (y=0; y < (ssize_t) image->rows; y++)
774   {
775     double
776       channel_distortion[MaxPixelChannels+1];
777
778     register const Quantum
779       *restrict p,
780       *restrict q;
781
782     register ssize_t
783       i,
784       x;
785
786     if (status == MagickFalse)
787       continue;
788     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
789     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
790       1,exception);
791     if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
792       {
793         status=MagickFalse;
794         continue;
795       }
796     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
797     for (x=0; x < (ssize_t) image->columns; x++)
798     {
799       double
800         Da,
801         Sa;
802
803       register ssize_t
804         i;
805
806       if (GetPixelReadMask(image,p) == 0)
807         {
808           p+=GetPixelChannels(image);
809           q+=GetPixelChannels(reconstruct_image);
810           continue;
811         }
812       Sa=QuantumScale*GetPixelAlpha(image,p);
813       Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
814       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
815       {
816         double
817           distance;
818
819         PixelChannel channel=GetPixelChannelChannel(image,i);
820         PixelTrait traits=GetPixelChannelTraits(image,channel);
821         PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
822           channel);
823         if ((traits == UndefinedPixelTrait) ||
824             (reconstruct_traits == UndefinedPixelTrait) ||
825             ((reconstruct_traits & UpdatePixelTrait) == 0))
826           continue;
827         distance=QuantumScale*(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
828           channel,q));
829         channel_distortion[i]+=distance*distance;
830         channel_distortion[CompositePixelChannel]+=distance*distance;
831       }
832       p+=GetPixelChannels(image);
833       q+=GetPixelChannels(reconstruct_image);
834     }
835 #if defined(MAGICKCORE_OPENMP_SUPPORT)
836     #pragma omp critical (MagickCore_GetMeanSquaredError)
837 #endif
838     for (i=0; i <= MaxPixelChannels; i++)
839       distortion[i]+=channel_distortion[i];
840   }
841   reconstruct_view=DestroyCacheView(reconstruct_view);
842   image_view=DestroyCacheView(image_view);
843   for (i=0; i <= MaxPixelChannels; i++)
844     distortion[i]/=((double) image->columns*image->rows);
845   distortion[CompositePixelChannel]/=GetImageChannels(image);
846   return(status);
847 }
848
849 static MagickBooleanType GetNormalizedCrossCorrelationDistortion(
850   const Image *image,const Image *reconstruct_image,double *distortion,
851   ExceptionInfo *exception)
852 {
853 #define SimilarityImageTag  "Similarity/Image"
854
855   CacheView
856     *image_view,
857     *reconstruct_view;
858
859   ChannelStatistics
860     *image_statistics,
861     *reconstruct_statistics;
862
863   double
864     area;
865
866   MagickBooleanType
867     status;
868
869   MagickOffsetType
870     progress;
871
872   register ssize_t
873     i;
874
875   ssize_t
876     y;
877
878   /*
879     Normalize to account for variation due to lighting and exposure condition.
880   */
881   image_statistics=GetImageStatistics(image,exception);
882   reconstruct_statistics=GetImageStatistics(reconstruct_image,exception);
883   if ((image_statistics == (ChannelStatistics *) NULL) ||
884       (reconstruct_statistics == (ChannelStatistics *) NULL))
885     {
886       if (image_statistics != (ChannelStatistics *) NULL)
887         image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
888           image_statistics);
889       if (reconstruct_statistics != (ChannelStatistics *) NULL)
890         reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
891           reconstruct_statistics);
892       return(MagickFalse);
893     }
894   status=MagickTrue;
895   progress=0;
896   for (i=0; i <= MaxPixelChannels; i++)
897     distortion[i]=0.0;
898   area=1.0/((double) image->columns*image->rows);
899   image_view=AcquireVirtualCacheView(image,exception);
900   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
901   for (y=0; y < (ssize_t) image->rows; y++)
902   {
903     register const Quantum
904       *restrict p,
905       *restrict q;
906
907     register ssize_t
908       x;
909
910     if (status == MagickFalse)
911       continue;
912     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
913     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
914       1,exception);
915     if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
916       {
917         status=MagickFalse;
918         continue;
919       }
920     for (x=0; x < (ssize_t) image->columns; x++)
921     {
922       double
923         Da,
924         Sa;
925
926       register ssize_t
927         i;
928
929       if (GetPixelReadMask(image,p) == 0)
930         {
931           p+=GetPixelChannels(image);
932           q+=GetPixelChannels(reconstruct_image);
933           continue;
934         }
935       Sa=QuantumScale*GetPixelAlpha(image,p);
936       Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
937       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
938       {
939         PixelChannel channel=GetPixelChannelChannel(image,i);
940         PixelTrait traits=GetPixelChannelTraits(image,channel);
941         PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
942           channel);
943         if ((traits == UndefinedPixelTrait) ||
944             (reconstruct_traits == UndefinedPixelTrait) ||
945             ((reconstruct_traits & UpdatePixelTrait) == 0))
946           continue;
947         distortion[i]+=area*QuantumScale*(Sa*p[i]-image_statistics[i].mean)*
948           (Da*GetPixelChannel(reconstruct_image,channel,q)-
949           reconstruct_statistics[channel].mean);
950       }
951       p+=GetPixelChannels(image);
952       q+=GetPixelChannels(reconstruct_image);
953     }
954     if (image->progress_monitor != (MagickProgressMonitor) NULL)
955       {
956         MagickBooleanType
957           proceed;
958
959         proceed=SetImageProgress(image,SimilarityImageTag,progress++,
960           image->rows);
961         if (proceed == MagickFalse)
962           status=MagickFalse;
963       }
964   }
965   reconstruct_view=DestroyCacheView(reconstruct_view);
966   image_view=DestroyCacheView(image_view);
967   /*
968     Divide by the standard deviation.
969   */
970   distortion[CompositePixelChannel]=0.0;
971   for (i=0; i < MaxPixelChannels; i++)
972   {
973     double
974       gamma;
975
976     PixelChannel channel=GetPixelChannelChannel(image,i);
977     gamma=image_statistics[i].standard_deviation*
978       reconstruct_statistics[channel].standard_deviation;
979     gamma=PerceptibleReciprocal(gamma);
980     distortion[i]=QuantumRange*gamma*distortion[i];
981     distortion[CompositePixelChannel]+=distortion[i]*distortion[i];
982   }
983   distortion[CompositePixelChannel]=sqrt(distortion[CompositePixelChannel]/
984     GetImageChannels(image));
985   /*
986     Free resources.
987   */
988   reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
989     reconstruct_statistics);
990   image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
991     image_statistics);
992   return(status);
993 }
994
995 static MagickBooleanType GetPeakAbsoluteDistortion(const Image *image,
996   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
997 {
998   CacheView
999     *image_view,
1000     *reconstruct_view;
1001
1002   MagickBooleanType
1003     status;
1004
1005   ssize_t
1006     y;
1007
1008   status=MagickTrue;
1009   image_view=AcquireVirtualCacheView(image,exception);
1010   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1011 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1012   #pragma omp parallel for schedule(static,4) shared(status) \
1013     magick_threads(image,image,image->rows,1)
1014 #endif
1015   for (y=0; y < (ssize_t) image->rows; y++)
1016   {
1017     double
1018       channel_distortion[MaxPixelChannels+1];
1019
1020     register const Quantum
1021       *restrict p,
1022       *restrict q;
1023
1024     register ssize_t
1025       i,
1026       x;
1027
1028     if (status == MagickFalse)
1029       continue;
1030     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1031     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
1032       1,exception);
1033     if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1034       {
1035         status=MagickFalse;
1036         continue;
1037       }
1038     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
1039     for (x=0; x < (ssize_t) image->columns; x++)
1040     {
1041       double
1042         Da,
1043         Sa;
1044
1045       register ssize_t
1046         i;
1047
1048       if (GetPixelReadMask(image,p) == 0)
1049         {
1050           p+=GetPixelChannels(image);
1051           q+=GetPixelChannels(reconstruct_image);
1052           continue;
1053         }
1054       Sa=QuantumScale*GetPixelAlpha(image,p);
1055       Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
1056       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1057       {
1058         double
1059           distance;
1060
1061         PixelChannel channel=GetPixelChannelChannel(image,i);
1062         PixelTrait traits=GetPixelChannelTraits(image,channel);
1063         PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
1064           channel);
1065         if ((traits == UndefinedPixelTrait) ||
1066             (reconstruct_traits == UndefinedPixelTrait) ||
1067             ((reconstruct_traits & UpdatePixelTrait) == 0))
1068           continue;
1069         distance=QuantumScale*fabs(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
1070           channel,q));
1071         if (distance > channel_distortion[i])
1072           channel_distortion[i]=distance;
1073         if (distance > channel_distortion[CompositePixelChannel])
1074           channel_distortion[CompositePixelChannel]=distance;
1075       }
1076       p+=GetPixelChannels(image);
1077       q+=GetPixelChannels(reconstruct_image);
1078     }
1079 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1080     #pragma omp critical (MagickCore_GetPeakAbsoluteError)
1081 #endif
1082     for (i=0; i <= MaxPixelChannels; i++)
1083       if (channel_distortion[i] > distortion[i])
1084         distortion[i]=channel_distortion[i];
1085   }
1086   reconstruct_view=DestroyCacheView(reconstruct_view);
1087   image_view=DestroyCacheView(image_view);
1088   return(status);
1089 }
1090
1091 static inline double MagickLog10(const double x)
1092 {
1093 #define Log10Epsilon  (1.0e-11)
1094
1095  if (fabs(x) < Log10Epsilon)
1096    return(log10(Log10Epsilon));
1097  return(log10(fabs(x)));
1098 }
1099
1100 static MagickBooleanType GetPeakSignalToNoiseRatio(const Image *image,
1101   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1102 {
1103   MagickBooleanType
1104     status;
1105
1106   register ssize_t
1107     i;
1108
1109   status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1110   for (i=0; i <= MaxPixelChannels; i++)
1111     distortion[i]=20.0*MagickLog10((double) 1.0/sqrt(distortion[i]));
1112   return(status);
1113 }
1114
1115 static MagickBooleanType GetPerceptualHashDistortion(const Image *image,
1116   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1117 {
1118   ChannelPerceptualHash
1119     *image_phash,
1120     *reconstruct_phash;
1121
1122   ssize_t
1123     channel;
1124
1125   /*
1126     Compute perceptual hash in the sRGB colorspace.
1127   */
1128   image_phash=GetImagePerceptualHash(image,exception);
1129   if (image_phash == (ChannelPerceptualHash *) NULL)
1130     return(MagickFalse);
1131   reconstruct_phash=GetImagePerceptualHash(reconstruct_image,exception);
1132   if (image_phash == (ChannelPerceptualHash *) NULL)
1133     {
1134       image_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(image_phash);
1135       return(MagickFalse);
1136     }
1137 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1138   #pragma omp parallel for schedule(static,4)
1139 #endif
1140   for (channel=0; channel < MaxPixelChannels; channel++)
1141   {
1142     double
1143       difference;
1144
1145     register ssize_t
1146       i;
1147
1148     difference=0.0;
1149     for (i=0; i < 7; i++)
1150     {
1151       double
1152         alpha,
1153         beta;
1154
1155       alpha=image_phash[channel].P[i];
1156       beta=reconstruct_phash[channel].P[i];
1157       difference+=(beta-alpha)*(beta-alpha);
1158     }
1159     distortion[channel]+=difference;
1160 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1161     #pragma omp critical (MagickCore_GetPerceptualHashDistortion)
1162 #endif
1163     distortion[CompositePixelChannel]+=difference;
1164   }
1165   /*
1166     Compute perceptual hash in the HCLP colorspace.
1167   */
1168 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1169   #pragma omp parallel for schedule(static,4)
1170 #endif
1171   for (channel=0; channel < MaxPixelChannels; channel++)
1172   {
1173     double
1174       difference;
1175
1176     register ssize_t
1177       i;
1178
1179     difference=0.0;
1180     for (i=0; i < 7; i++)
1181     {
1182       double
1183         alpha,
1184         beta;
1185
1186       alpha=image_phash[channel].Q[i];
1187       beta=reconstruct_phash[channel].Q[i];
1188       difference+=(beta-alpha)*(beta-alpha);
1189     }
1190     distortion[channel]+=difference;
1191 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1192     #pragma omp critical (MagickCore_GetPerceptualHashDistortion)
1193 #endif
1194     distortion[CompositePixelChannel]+=difference;
1195   }
1196   /*
1197     Free resources.
1198   */
1199   reconstruct_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1200     reconstruct_phash);
1201   image_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(image_phash);
1202   return(MagickTrue);
1203 }
1204
1205
1206
1207 static MagickBooleanType GetRootMeanSquaredDistortion(const Image *image,
1208   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1209 {
1210   MagickBooleanType
1211     status;
1212
1213   register ssize_t
1214     i;
1215
1216   status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1217   for (i=0; i <= MaxPixelChannels; i++)
1218     distortion[i]=sqrt(distortion[i]);
1219   return(status);
1220 }
1221
1222 MagickExport MagickBooleanType GetImageDistortion(Image *image,
1223   const Image *reconstruct_image,const MetricType metric,double *distortion,
1224   ExceptionInfo *exception)
1225 {
1226   double
1227     *channel_distortion;
1228
1229   MagickBooleanType
1230     status;
1231
1232   size_t
1233     length;
1234
1235   assert(image != (Image *) NULL);
1236   assert(image->signature == MagickSignature);
1237   if (image->debug != MagickFalse)
1238     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1239   assert(reconstruct_image != (const Image *) NULL);
1240   assert(reconstruct_image->signature == MagickSignature);
1241   assert(distortion != (double *) NULL);
1242   *distortion=0.0;
1243   if (image->debug != MagickFalse)
1244     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1245   if (metric != PerceptualHashErrorMetric)
1246     if ((reconstruct_image->columns != image->columns) ||
1247         (reconstruct_image->rows != image->rows))
1248       ThrowBinaryException(ImageError,"ImageSizeDiffers",image->filename);
1249   /*
1250     Get image distortion.
1251   */
1252   length=MaxPixelChannels+1;
1253   channel_distortion=(double *) AcquireQuantumMemory(length,
1254     sizeof(*channel_distortion));
1255   if (channel_distortion == (double *) NULL)
1256     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1257   (void) ResetMagickMemory(channel_distortion,0,length*
1258     sizeof(*channel_distortion));
1259   switch (metric)
1260   {
1261     case AbsoluteErrorMetric:
1262     {
1263       status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1264         exception);
1265       break;
1266     }
1267     case FuzzErrorMetric:
1268     {
1269       status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1270         exception);
1271       break;
1272     }
1273     case MeanAbsoluteErrorMetric:
1274     {
1275       status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1276         channel_distortion,exception);
1277       break;
1278     }
1279     case MeanErrorPerPixelErrorMetric:
1280     {
1281       status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1282         exception);
1283       break;
1284     }
1285     case MeanSquaredErrorMetric:
1286     {
1287       status=GetMeanSquaredDistortion(image,reconstruct_image,
1288         channel_distortion,exception);
1289       break;
1290     }
1291     case NormalizedCrossCorrelationErrorMetric:
1292     default:
1293     {
1294       status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1295         channel_distortion,exception);
1296       break;
1297     }
1298     case PeakAbsoluteErrorMetric:
1299     {
1300       status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1301         channel_distortion,exception);
1302       break;
1303     }
1304     case PeakSignalToNoiseRatioErrorMetric:
1305     {
1306       status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1307         channel_distortion,exception);
1308       break;
1309     }
1310     case PerceptualHashErrorMetric:
1311     {
1312       status=GetPerceptualHashDistortion(image,reconstruct_image,
1313         channel_distortion,exception);
1314       break;
1315     }
1316     case RootMeanSquaredErrorMetric:
1317     {
1318       status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1319         channel_distortion,exception);
1320       break;
1321     }
1322   }
1323   *distortion=channel_distortion[CompositePixelChannel];
1324   channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1325   (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
1326     *distortion);
1327   return(status);
1328 }
1329 \f
1330 /*
1331 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1332 %                                                                             %
1333 %                                                                             %
1334 %                                                                             %
1335 %   G e t I m a g e D i s t o r t i o n s                                     %
1336 %                                                                             %
1337 %                                                                             %
1338 %                                                                             %
1339 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1340 %
1341 %  GetImageDistortions() compares the pixel channels of an image to a
1342 %  reconstructed image and returns the specified distortion metric for each
1343 %  channel.
1344 %
1345 %  The format of the GetImageDistortions method is:
1346 %
1347 %      double *GetImageDistortions(const Image *image,
1348 %        const Image *reconstruct_image,const MetricType metric,
1349 %        ExceptionInfo *exception)
1350 %
1351 %  A description of each parameter follows:
1352 %
1353 %    o image: the image.
1354 %
1355 %    o reconstruct_image: the reconstruct image.
1356 %
1357 %    o metric: the metric.
1358 %
1359 %    o exception: return any errors or warnings in this structure.
1360 %
1361 */
1362 MagickExport double *GetImageDistortions(Image *image,
1363   const Image *reconstruct_image,const MetricType metric,
1364   ExceptionInfo *exception)
1365 {
1366   double
1367     *channel_distortion;
1368
1369   MagickBooleanType
1370     status;
1371
1372   size_t
1373     length;
1374
1375   assert(image != (Image *) NULL);
1376   assert(image->signature == MagickSignature);
1377   if (image->debug != MagickFalse)
1378     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1379   assert(reconstruct_image != (const Image *) NULL);
1380   assert(reconstruct_image->signature == MagickSignature);
1381   if (image->debug != MagickFalse)
1382     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1383   if (metric != PerceptualHashErrorMetric)
1384     if ((reconstruct_image->columns != image->columns) ||
1385         (reconstruct_image->rows != image->rows))
1386       {
1387         (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1388           "ImageSizeDiffers","`%s'",image->filename);
1389         return((double *) NULL);
1390       }
1391   /*
1392     Get image distortion.
1393   */
1394   length=MaxPixelChannels+1UL;
1395   channel_distortion=(double *) AcquireQuantumMemory(length,
1396     sizeof(*channel_distortion));
1397   if (channel_distortion == (double *) NULL)
1398     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1399   (void) ResetMagickMemory(channel_distortion,0,length*
1400     sizeof(*channel_distortion));
1401   status=MagickTrue;
1402   switch (metric)
1403   {
1404     case AbsoluteErrorMetric:
1405     {
1406       status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1407         exception);
1408       break;
1409     }
1410     case FuzzErrorMetric:
1411     {
1412       status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1413         exception);
1414       break;
1415     }
1416     case MeanAbsoluteErrorMetric:
1417     {
1418       status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1419         channel_distortion,exception);
1420       break;
1421     }
1422     case MeanErrorPerPixelErrorMetric:
1423     {
1424       status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1425         exception);
1426       break;
1427     }
1428     case MeanSquaredErrorMetric:
1429     {
1430       status=GetMeanSquaredDistortion(image,reconstruct_image,
1431         channel_distortion,exception);
1432       break;
1433     }
1434     case NormalizedCrossCorrelationErrorMetric:
1435     default:
1436     {
1437       status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1438         channel_distortion,exception);
1439       break;
1440     }
1441     case PeakAbsoluteErrorMetric:
1442     {
1443       status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1444         channel_distortion,exception);
1445       break;
1446     }
1447     case PeakSignalToNoiseRatioErrorMetric:
1448     {
1449       status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1450         channel_distortion,exception);
1451       break;
1452     }
1453     case PerceptualHashErrorMetric:
1454     {
1455       status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1456         channel_distortion,exception);
1457       break;
1458     }
1459     case RootMeanSquaredErrorMetric:
1460     {
1461       status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1462         channel_distortion,exception);
1463       break;
1464     }
1465   }
1466   if (status == MagickFalse)
1467     {
1468       channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1469       return((double *) NULL);
1470     }
1471   return(channel_distortion);
1472 }
1473 \f
1474 /*
1475 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1476 %                                                                             %
1477 %                                                                             %
1478 %                                                                             %
1479 %  I s I m a g e s E q u a l                                                  %
1480 %                                                                             %
1481 %                                                                             %
1482 %                                                                             %
1483 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1484 %
1485 %  IsImagesEqual() measures the difference between colors at each pixel
1486 %  location of two images.  A value other than 0 means the colors match
1487 %  exactly.  Otherwise an error measure is computed by summing over all
1488 %  pixels in an image the distance squared in RGB space between each image
1489 %  pixel and its corresponding pixel in the reconstruct image.  The error
1490 %  measure is assigned to these image members:
1491 %
1492 %    o mean_error_per_pixel:  The mean error for any single pixel in
1493 %      the image.
1494 %
1495 %    o normalized_mean_error:  The normalized mean quantization error for
1496 %      any single pixel in the image.  This distance measure is normalized to
1497 %      a range between 0 and 1.  It is independent of the range of red, green,
1498 %      and blue values in the image.
1499 %
1500 %    o normalized_maximum_error:  The normalized maximum quantization
1501 %      error for any single pixel in the image.  This distance measure is
1502 %      normalized to a range between 0 and 1.  It is independent of the range
1503 %      of red, green, and blue values in your image.
1504 %
1505 %  A small normalized mean square error, accessed as
1506 %  image->normalized_mean_error, suggests the images are very similar in
1507 %  spatial layout and color.
1508 %
1509 %  The format of the IsImagesEqual method is:
1510 %
1511 %      MagickBooleanType IsImagesEqual(Image *image,
1512 %        const Image *reconstruct_image,ExceptionInfo *exception)
1513 %
1514 %  A description of each parameter follows.
1515 %
1516 %    o image: the image.
1517 %
1518 %    o reconstruct_image: the reconstruct image.
1519 %
1520 %    o exception: return any errors or warnings in this structure.
1521 %
1522 */
1523 MagickExport MagickBooleanType IsImagesEqual(Image *image,
1524   const Image *reconstruct_image,ExceptionInfo *exception)
1525 {
1526   CacheView
1527     *image_view,
1528     *reconstruct_view;
1529
1530   MagickBooleanType
1531     status;
1532
1533   double
1534     area,
1535     maximum_error,
1536     mean_error,
1537     mean_error_per_pixel;
1538
1539   ssize_t
1540     y;
1541
1542   assert(image != (Image *) NULL);
1543   assert(image->signature == MagickSignature);
1544   assert(reconstruct_image != (const Image *) NULL);
1545   assert(reconstruct_image->signature == MagickSignature);
1546   if ((reconstruct_image->columns != image->columns) ||
1547       (reconstruct_image->rows != image->rows))
1548     ThrowBinaryException(ImageError,"ImageSizeDiffers",image->filename);
1549   area=0.0;
1550   maximum_error=0.0;
1551   mean_error_per_pixel=0.0;
1552   mean_error=0.0;
1553   image_view=AcquireVirtualCacheView(image,exception);
1554   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1555   for (y=0; y < (ssize_t) image->rows; y++)
1556   {
1557     register const Quantum
1558       *restrict p,
1559       *restrict q;
1560
1561     register ssize_t
1562       x;
1563
1564     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1565     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
1566       1,exception);
1567     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
1568       break;
1569     for (x=0; x < (ssize_t) image->columns; x++)
1570     {
1571       register ssize_t
1572         i;
1573
1574       if (GetPixelReadMask(image,p) == 0)
1575         {
1576           p+=GetPixelChannels(image);
1577           q+=GetPixelChannels(reconstruct_image);
1578           continue;
1579         }
1580       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1581       {
1582         double
1583           distance;
1584
1585         PixelChannel channel=GetPixelChannelChannel(image,i);
1586         PixelTrait traits=GetPixelChannelTraits(image,channel);
1587         PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
1588           channel);
1589         if ((traits == UndefinedPixelTrait) ||
1590             (reconstruct_traits == UndefinedPixelTrait) ||
1591             ((reconstruct_traits & UpdatePixelTrait) == 0))
1592           continue;
1593         distance=fabs(p[i]-(double) GetPixelChannel(reconstruct_image,
1594           channel,q));
1595         mean_error_per_pixel+=distance;
1596         mean_error+=distance*distance;
1597         if (distance > maximum_error)
1598           maximum_error=distance;
1599         area++;
1600       }
1601       p+=GetPixelChannels(image);
1602       q+=GetPixelChannels(reconstruct_image);
1603     }
1604   }
1605   reconstruct_view=DestroyCacheView(reconstruct_view);
1606   image_view=DestroyCacheView(image_view);
1607   image->error.mean_error_per_pixel=(double) (mean_error_per_pixel/area);
1608   image->error.normalized_mean_error=(double) (QuantumScale*QuantumScale*
1609     mean_error/area);
1610   image->error.normalized_maximum_error=(double) (QuantumScale*maximum_error);
1611   status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
1612   return(status);
1613 }
1614 \f
1615 /*
1616 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1617 %                                                                             %
1618 %                                                                             %
1619 %                                                                             %
1620 %   S i m i l a r i t y I m a g e                                             %
1621 %                                                                             %
1622 %                                                                             %
1623 %                                                                             %
1624 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1625 %
1626 %  SimilarityImage() compares the reference image of the image and returns the
1627 %  best match offset.  In addition, it returns a similarity image such that an
1628 %  exact match location is completely white and if none of the pixels match,
1629 %  black, otherwise some gray level in-between.
1630 %
1631 %  The format of the SimilarityImageImage method is:
1632 %
1633 %      Image *SimilarityImage(const Image *image,const Image *reference,
1634 %        const MetricType metric,const double similarity_threshold,
1635 %        RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
1636 %
1637 %  A description of each parameter follows:
1638 %
1639 %    o image: the image.
1640 %
1641 %    o reference: find an area of the image that closely resembles this image.
1642 %
1643 %    o metric: the metric.
1644 %
1645 %    o similarity_threshold: minimum distortion for (sub)image match.
1646 %
1647 %    o offset: the best match offset of the reference image within the image.
1648 %
1649 %    o similarity: the computed similarity between the images.
1650 %
1651 %    o exception: return any errors or warnings in this structure.
1652 %
1653 */
1654
1655 static double GetSimilarityMetric(const Image *image,const Image *reference,
1656   const MetricType metric,const ssize_t x_offset,const ssize_t y_offset,
1657   ExceptionInfo *exception)
1658 {
1659   double
1660     distortion;
1661
1662   Image
1663     *similarity_image;
1664
1665   MagickBooleanType
1666     status;
1667
1668   RectangleInfo
1669     geometry;
1670
1671   SetGeometry(reference,&geometry);
1672   geometry.x=x_offset;
1673   geometry.y=y_offset;
1674   similarity_image=CropImage(image,&geometry,exception);
1675   if (similarity_image == (Image *) NULL)
1676     return(0.0);
1677   distortion=0.0;
1678   status=GetImageDistortion(similarity_image,reference,metric,&distortion,
1679     exception);
1680   similarity_image=DestroyImage(similarity_image);
1681   if (status == MagickFalse)
1682     return(0.0);
1683   return(distortion);
1684 }
1685
1686 static inline double MagickMin(const double x,const double y)
1687 {
1688   if (x < y)
1689     return(x);
1690   return(y);
1691 }
1692
1693 MagickExport Image *SimilarityImage(Image *image,const Image *reference,
1694   const MetricType metric,const double similarity_threshold,
1695   RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
1696 {
1697 #define SimilarityImageTag  "Similarity/Image"
1698
1699   CacheView
1700     *similarity_view;
1701
1702   Image
1703     *similarity_image;
1704
1705   MagickBooleanType
1706     status;
1707
1708   MagickOffsetType
1709     progress;
1710
1711   ssize_t
1712     y;
1713
1714   assert(image != (const Image *) NULL);
1715   assert(image->signature == MagickSignature);
1716   if (image->debug != MagickFalse)
1717     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1718   assert(exception != (ExceptionInfo *) NULL);
1719   assert(exception->signature == MagickSignature);
1720   assert(offset != (RectangleInfo *) NULL);
1721   SetGeometry(reference,offset);
1722   *similarity_metric=MagickMaximumValue;
1723   if ((reference->columns > image->columns) || (reference->rows > image->rows))
1724     ThrowImageException(ImageError,"ImageSizeDiffers");
1725   similarity_image=CloneImage(image,image->columns-reference->columns+1,
1726     image->rows-reference->rows+1,MagickTrue,exception);
1727   if (similarity_image == (Image *) NULL)
1728     return((Image *) NULL);
1729   status=SetImageStorageClass(similarity_image,DirectClass,exception);
1730   if (status == MagickFalse)
1731     {
1732       similarity_image=DestroyImage(similarity_image);
1733       return((Image *) NULL);
1734     }
1735   (void) SetImageAlphaChannel(similarity_image,DeactivateAlphaChannel,
1736     exception);
1737   /*
1738     Measure similarity of reference image against image.
1739   */
1740   status=MagickTrue;
1741   progress=0;
1742   similarity_view=AcquireAuthenticCacheView(similarity_image,exception);
1743 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1744   #pragma omp parallel for schedule(static,4) \
1745     shared(progress,status,similarity_metric) \
1746     magick_threads(image,image,image->rows,1)
1747 #endif
1748   for (y=0; y < (ssize_t) (image->rows-reference->rows+1); y++)
1749   {
1750     double
1751       similarity;
1752
1753     register Quantum
1754       *restrict q;
1755
1756     register ssize_t
1757       x;
1758
1759     if (status == MagickFalse)
1760       continue;
1761 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1762     #pragma omp flush(similarity_metric)
1763 #endif
1764     if (*similarity_metric <= similarity_threshold)
1765       continue;
1766     q=GetCacheViewAuthenticPixels(similarity_view,0,y,similarity_image->columns,
1767       1,exception);
1768     if (q == (Quantum *) NULL)
1769       {
1770         status=MagickFalse;
1771         continue;
1772       }
1773     for (x=0; x < (ssize_t) (image->columns-reference->columns+1); x++)
1774     {
1775       register ssize_t
1776         i;
1777
1778 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1779       #pragma omp flush(similarity_metric)
1780 #endif
1781       if (*similarity_metric <= similarity_threshold)
1782         break;
1783       similarity=GetSimilarityMetric(image,reference,metric,x,y,exception);
1784 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1785       #pragma omp critical (MagickCore_SimilarityImage)
1786 #endif
1787       if (similarity < *similarity_metric)
1788         {
1789           offset->x=x;
1790           offset->y=y;
1791           *similarity_metric=similarity;
1792         }
1793       if (metric == PerceptualHashErrorMetric)
1794         similarity=MagickMin(0.01*similarity,1.0);
1795       if (GetPixelReadMask(similarity_image,q) == 0)
1796         {
1797           SetPixelBackgoundColor(similarity_image,q);
1798           q+=GetPixelChannels(similarity_image);
1799           continue;
1800         }
1801       for (i=0; i < (ssize_t) GetPixelChannels(similarity_image); i++)
1802       {
1803         PixelChannel channel=GetPixelChannelChannel(image,i);
1804         PixelTrait traits=GetPixelChannelTraits(image,channel);
1805         PixelTrait similarity_traits=GetPixelChannelTraits(similarity_image,
1806           channel);
1807         if ((traits == UndefinedPixelTrait) ||
1808             (similarity_traits == UndefinedPixelTrait) ||
1809             ((similarity_traits & UpdatePixelTrait) == 0))
1810           continue;
1811         SetPixelChannel(similarity_image,channel,ClampToQuantum(QuantumRange-
1812           QuantumRange*similarity),q);
1813       }
1814       q+=GetPixelChannels(similarity_image);
1815     }
1816     if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
1817       status=MagickFalse;
1818     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1819       {
1820         MagickBooleanType
1821           proceed;
1822
1823 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1824         #pragma omp critical (MagickCore_SimilarityImage)
1825 #endif
1826         proceed=SetImageProgress(image,SimilarityImageTag,progress++,
1827           image->rows);
1828         if (proceed == MagickFalse)
1829           status=MagickFalse;
1830       }
1831   }
1832   similarity_view=DestroyCacheView(similarity_view);
1833   if (status == MagickFalse)
1834     similarity_image=DestroyImage(similarity_image);
1835   return(similarity_image);
1836 }