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