MagickCore  7.0.10
compare.c
Go to the documentation of this file.
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-2020 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 % https://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"
51 #include "MagickCore/colorspace.h"
53 #include "MagickCore/compare.h"
55 #include "MagickCore/constitute.h"
57 #include "MagickCore/geometry.h"
59 #include "MagickCore/list.h"
60 #include "MagickCore/log.h"
61 #include "MagickCore/memory_.h"
62 #include "MagickCore/monitor.h"
64 #include "MagickCore/option.h"
66 #include "MagickCore/property.h"
67 #include "MagickCore/resource_.h"
68 #include "MagickCore/string_.h"
69 #include "MagickCore/statistic.h"
72 #include "MagickCore/transform.h"
73 #include "MagickCore/utility.h"
74 #include "MagickCore/version.h"
75 
76 /*
77 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
78 % %
79 % %
80 % %
81 % C o m p a r e I m a g e s %
82 % %
83 % %
84 % %
85 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
86 %
87 % CompareImages() compares one or more pixel channels of an image to a
88 % reconstructed image and returns the difference image.
89 %
90 % The format of the CompareImages method is:
91 %
92 % Image *CompareImages(const Image *image,const Image *reconstruct_image,
93 % const MetricType metric,double *distortion,ExceptionInfo *exception)
94 %
95 % A description of each parameter follows:
96 %
97 % o image: the image.
98 %
99 % o reconstruct_image: the reconstruct image.
100 %
101 % o metric: the metric.
102 %
103 % o distortion: the computed distortion between the images.
104 %
105 % o exception: return any errors or warnings in this structure.
106 %
107 */
108 
109 static size_t GetImageChannels(const Image *image)
110 {
111  register ssize_t
112  i;
113 
114  size_t
115  channels;
116 
117  channels=0;
118  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
119  {
120  PixelChannel channel = GetPixelChannelChannel(image,i);
121  PixelTrait traits = GetPixelChannelTraits(image,channel);
122  if ((traits & UpdatePixelTrait) != 0)
123  channels++;
124  }
125  return(channels == 0 ? (size_t) 1 : channels);
126 }
127 
128 MagickExport Image *CompareImages(Image *image,const Image *reconstruct_image,
129  const MetricType metric,double *distortion,ExceptionInfo *exception)
130 {
131  CacheView
132  *highlight_view,
133  *image_view,
134  *reconstruct_view;
135 
136  const char
137  *artifact;
138 
139  double
140  fuzz;
141 
142  Image
143  *clone_image,
144  *difference_image,
145  *highlight_image;
146 
148  status;
149 
150  PixelInfo
151  highlight,
152  lowlight,
153  masklight;
154 
156  geometry;
157 
158  size_t
159  columns,
160  rows;
161 
162  ssize_t
163  y;
164 
165  assert(image != (Image *) NULL);
166  assert(image->signature == MagickCoreSignature);
167  if (image->debug != MagickFalse)
168  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
169  assert(reconstruct_image != (const Image *) NULL);
170  assert(reconstruct_image->signature == MagickCoreSignature);
171  assert(distortion != (double *) NULL);
172  *distortion=0.0;
173  if (image->debug != MagickFalse)
174  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
175  status=GetImageDistortion(image,reconstruct_image,metric,distortion,
176  exception);
177  if (status == MagickFalse)
178  return((Image *) NULL);
179  columns=MagickMax(image->columns,reconstruct_image->columns);
180  rows=MagickMax(image->rows,reconstruct_image->rows);
181  SetGeometry(image,&geometry);
182  geometry.width=columns;
183  geometry.height=rows;
184  clone_image=CloneImage(image,0,0,MagickTrue,exception);
185  if (clone_image == (Image *) NULL)
186  return((Image *) NULL);
187  (void) SetImageMask(clone_image,ReadPixelMask,(Image *) NULL,exception);
188  difference_image=ExtentImage(clone_image,&geometry,exception);
189  clone_image=DestroyImage(clone_image);
190  if (difference_image == (Image *) NULL)
191  return((Image *) NULL);
192  (void) SetImageAlphaChannel(difference_image,OpaqueAlphaChannel,exception);
193  highlight_image=CloneImage(image,columns,rows,MagickTrue,exception);
194  if (highlight_image == (Image *) NULL)
195  {
196  difference_image=DestroyImage(difference_image);
197  return((Image *) NULL);
198  }
199  status=SetImageStorageClass(highlight_image,DirectClass,exception);
200  if (status == MagickFalse)
201  {
202  difference_image=DestroyImage(difference_image);
203  highlight_image=DestroyImage(highlight_image);
204  return((Image *) NULL);
205  }
206  (void) SetImageMask(highlight_image,ReadPixelMask,(Image *) NULL,exception);
207  (void) SetImageAlphaChannel(highlight_image,OpaqueAlphaChannel,exception);
208  (void) QueryColorCompliance("#f1001ecc",AllCompliance,&highlight,exception);
209  artifact=GetImageArtifact(image,"compare:highlight-color");
210  if (artifact != (const char *) NULL)
211  (void) QueryColorCompliance(artifact,AllCompliance,&highlight,exception);
212  (void) QueryColorCompliance("#ffffffcc",AllCompliance,&lowlight,exception);
213  artifact=GetImageArtifact(image,"compare:lowlight-color");
214  if (artifact != (const char *) NULL)
215  (void) QueryColorCompliance(artifact,AllCompliance,&lowlight,exception);
216  (void) QueryColorCompliance("#888888cc",AllCompliance,&masklight,exception);
217  artifact=GetImageArtifact(image,"compare:masklight-color");
218  if (artifact != (const char *) NULL)
219  (void) QueryColorCompliance(artifact,AllCompliance,&masklight,exception);
220  /*
221  Generate difference image.
222  */
223  status=MagickTrue;
224  fuzz=(double) MagickMin(GetPixelChannels(image),
225  GetPixelChannels(reconstruct_image))*
226  GetFuzzyColorDistance(image,reconstruct_image);
227  image_view=AcquireVirtualCacheView(image,exception);
228  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
229  highlight_view=AcquireAuthenticCacheView(highlight_image,exception);
230 #if defined(MAGICKCORE_OPENMP_SUPPORT)
231  #pragma omp parallel for schedule(static) shared(status) \
232  magick_number_threads(image,highlight_image,rows,1)
233 #endif
234  for (y=0; y < (ssize_t) rows; y++)
235  {
237  sync;
238 
239  register const Quantum
240  *magick_restrict p,
241  *magick_restrict q;
242 
243  register Quantum
244  *magick_restrict r;
245 
246  register ssize_t
247  x;
248 
249  if (status == MagickFalse)
250  continue;
251  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
252  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
253  r=QueueCacheViewAuthenticPixels(highlight_view,0,y,columns,1,exception);
254  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL) ||
255  (r == (Quantum *) NULL))
256  {
257  status=MagickFalse;
258  continue;
259  }
260  for (x=0; x < (ssize_t) columns; x++)
261  {
262  double
263  Da,
264  distance,
265  Sa;
266 
268  difference;
269 
270  register ssize_t
271  i;
272 
273  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
274  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
275  {
276  SetPixelViaPixelInfo(highlight_image,&masklight,r);
277  p+=GetPixelChannels(image);
278  q+=GetPixelChannels(reconstruct_image);
279  r+=GetPixelChannels(highlight_image);
280  continue;
281  }
282  difference=MagickFalse;
283  distance=0.0;
284  Sa=QuantumScale*GetPixelAlpha(image,p);
285  Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
286  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
287  {
288  double
289  pixel;
290 
291  PixelChannel channel = GetPixelChannelChannel(image,i);
292  PixelTrait traits = GetPixelChannelTraits(image,channel);
293  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
294  channel);
295  if ((traits == UndefinedPixelTrait) ||
296  (reconstruct_traits == UndefinedPixelTrait) ||
297  ((reconstruct_traits & UpdatePixelTrait) == 0))
298  continue;
299  if (channel == AlphaPixelChannel)
300  pixel=(double) p[i]-GetPixelChannel(reconstruct_image,channel,q);
301  else
302  pixel=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
303  distance+=pixel*pixel;
304  if (distance >= fuzz)
305  {
306  difference=MagickTrue;
307  break;
308  }
309  }
310  if (difference == MagickFalse)
311  SetPixelViaPixelInfo(highlight_image,&lowlight,r);
312  else
313  SetPixelViaPixelInfo(highlight_image,&highlight,r);
314  p+=GetPixelChannels(image);
315  q+=GetPixelChannels(reconstruct_image);
316  r+=GetPixelChannels(highlight_image);
317  }
318  sync=SyncCacheViewAuthenticPixels(highlight_view,exception);
319  if (sync == MagickFalse)
320  status=MagickFalse;
321  }
322  highlight_view=DestroyCacheView(highlight_view);
323  reconstruct_view=DestroyCacheView(reconstruct_view);
324  image_view=DestroyCacheView(image_view);
325  (void) CompositeImage(difference_image,highlight_image,image->compose,
326  MagickTrue,0,0,exception);
327  highlight_image=DestroyImage(highlight_image);
328  if (status == MagickFalse)
329  difference_image=DestroyImage(difference_image);
330  return(difference_image);
331 }
332 
333 /*
334 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
335 % %
336 % %
337 % %
338 % G e t I m a g e D i s t o r t i o n %
339 % %
340 % %
341 % %
342 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
343 %
344 % GetImageDistortion() compares one or more pixel channels of an image to a
345 % reconstructed image and returns the specified distortion metric.
346 %
347 % The format of the GetImageDistortion method is:
348 %
349 % MagickBooleanType GetImageDistortion(const Image *image,
350 % const Image *reconstruct_image,const MetricType metric,
351 % double *distortion,ExceptionInfo *exception)
352 %
353 % A description of each parameter follows:
354 %
355 % o image: the image.
356 %
357 % o reconstruct_image: the reconstruct image.
358 %
359 % o metric: the metric.
360 %
361 % o distortion: the computed distortion between the images.
362 %
363 % o exception: return any errors or warnings in this structure.
364 %
365 */
366 
368  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
369 {
370  CacheView
371  *image_view,
372  *reconstruct_view;
373 
374  double
375  fuzz;
376 
378  status;
379 
380  size_t
381  columns,
382  rows;
383 
384  ssize_t
385  y;
386 
387  /*
388  Compute the absolute difference in pixels between two images.
389  */
390  status=MagickTrue;
391  fuzz=(double) MagickMin(GetPixelChannels(image),
392  GetPixelChannels(reconstruct_image))*
393  GetFuzzyColorDistance(image,reconstruct_image);
394  rows=MagickMax(image->rows,reconstruct_image->rows);
395  columns=MagickMax(image->columns,reconstruct_image->columns);
396  image_view=AcquireVirtualCacheView(image,exception);
397  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
398 #if defined(MAGICKCORE_OPENMP_SUPPORT)
399  #pragma omp parallel for schedule(static) shared(status) \
400  magick_number_threads(image,image,rows,1)
401 #endif
402  for (y=0; y < (ssize_t) rows; y++)
403  {
404  double
405  channel_distortion[MaxPixelChannels+1];
406 
407  register const Quantum
408  *magick_restrict p,
409  *magick_restrict q;
410 
411  register ssize_t
412  j,
413  x;
414 
415  if (status == MagickFalse)
416  continue;
417  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
418  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
419  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
420  {
421  status=MagickFalse;
422  continue;
423  }
424  (void) memset(channel_distortion,0,sizeof(channel_distortion));
425  for (x=0; x < (ssize_t) columns; x++)
426  {
427  double
428  Da,
429  distance,
430  Sa;
431 
433  difference;
434 
435  register ssize_t
436  i;
437 
438  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
439  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
440  {
441  p+=GetPixelChannels(image);
442  q+=GetPixelChannels(reconstruct_image);
443  continue;
444  }
445  difference=MagickFalse;
446  distance=0.0;
447  Sa=QuantumScale*GetPixelAlpha(image,p);
448  Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
449  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
450  {
451  double
452  pixel;
453 
454  PixelChannel channel = GetPixelChannelChannel(image,i);
455  PixelTrait traits = GetPixelChannelTraits(image,channel);
456  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
457  channel);
458  if ((traits == UndefinedPixelTrait) ||
459  (reconstruct_traits == UndefinedPixelTrait) ||
460  ((reconstruct_traits & UpdatePixelTrait) == 0))
461  continue;
462  if (channel == AlphaPixelChannel)
463  pixel=(double) p[i]-GetPixelChannel(reconstruct_image,channel,q);
464  else
465  pixel=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
466  distance+=pixel*pixel;
467  if (distance >= fuzz)
468  {
469  channel_distortion[i]++;
470  difference=MagickTrue;
471  }
472  }
473  if (difference != MagickFalse)
474  channel_distortion[CompositePixelChannel]++;
475  p+=GetPixelChannels(image);
476  q+=GetPixelChannels(reconstruct_image);
477  }
478 #if defined(MAGICKCORE_OPENMP_SUPPORT)
479  #pragma omp critical (MagickCore_GetAbsoluteDistortion)
480 #endif
481  for (j=0; j <= MaxPixelChannels; j++)
482  distortion[j]+=channel_distortion[j];
483  }
484  reconstruct_view=DestroyCacheView(reconstruct_view);
485  image_view=DestroyCacheView(image_view);
486  return(status);
487 }
488 
490  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
491 {
492  CacheView
493  *image_view,
494  *reconstruct_view;
495 
496  double
497  area;
498 
500  status;
501 
502  register ssize_t
503  j;
504 
505  size_t
506  columns,
507  rows;
508 
509  ssize_t
510  y;
511 
512  status=MagickTrue;
513  rows=MagickMax(image->rows,reconstruct_image->rows);
514  columns=MagickMax(image->columns,reconstruct_image->columns);
515  area=0.0;
516  image_view=AcquireVirtualCacheView(image,exception);
517  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
518 #if defined(MAGICKCORE_OPENMP_SUPPORT)
519  #pragma omp parallel for schedule(static) shared(status) \
520  magick_number_threads(image,image,rows,1) reduction(+:area)
521 #endif
522  for (y=0; y < (ssize_t) rows; y++)
523  {
524  double
525  channel_distortion[MaxPixelChannels+1];
526 
527  register const Quantum
528  *magick_restrict p,
529  *magick_restrict q;
530 
531  register ssize_t
532  x;
533 
534  if (status == MagickFalse)
535  continue;
536  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
537  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
538  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
539  {
540  status=MagickFalse;
541  continue;
542  }
543  (void) memset(channel_distortion,0,sizeof(channel_distortion));
544  for (x=0; x < (ssize_t) columns; x++)
545  {
546  double
547  Da,
548  Sa;
549 
550  register ssize_t
551  i;
552 
553  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
554  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
555  {
556  p+=GetPixelChannels(image);
557  q+=GetPixelChannels(reconstruct_image);
558  continue;
559  }
560  Sa=QuantumScale*GetPixelAlpha(image,p);
561  Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
562  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
563  {
564  double
565  distance;
566 
567  PixelChannel channel = GetPixelChannelChannel(image,i);
568  PixelTrait traits = GetPixelChannelTraits(image,channel);
569  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
570  channel);
571  if ((traits == UndefinedPixelTrait) ||
572  (reconstruct_traits == UndefinedPixelTrait) ||
573  ((reconstruct_traits & UpdatePixelTrait) == 0))
574  continue;
575  if (channel == AlphaPixelChannel)
576  distance=QuantumScale*(p[i]-GetPixelChannel(reconstruct_image,
577  channel,q));
578  else
579  distance=QuantumScale*(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
580  channel,q));
581  channel_distortion[i]+=distance*distance;
582  channel_distortion[CompositePixelChannel]+=distance*distance;
583  }
584  area++;
585  p+=GetPixelChannels(image);
586  q+=GetPixelChannels(reconstruct_image);
587  }
588 #if defined(MAGICKCORE_OPENMP_SUPPORT)
589  #pragma omp critical (MagickCore_GetFuzzDistortion)
590 #endif
591  for (j=0; j <= MaxPixelChannels; j++)
592  distortion[j]+=channel_distortion[j];
593  }
594  reconstruct_view=DestroyCacheView(reconstruct_view);
595  image_view=DestroyCacheView(image_view);
596  area=PerceptibleReciprocal(area);
597  for (j=0; j <= MaxPixelChannels; j++)
598  distortion[j]*=area;
599  distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
600  distortion[CompositePixelChannel]=sqrt(distortion[CompositePixelChannel]);
601  return(status);
602 }
603 
605  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
606 {
607  CacheView
608  *image_view,
609  *reconstruct_view;
610 
611  double
612  area;
613 
615  status;
616 
617  register ssize_t
618  j;
619 
620  size_t
621  columns,
622  rows;
623 
624  ssize_t
625  y;
626 
627  status=MagickTrue;
628  rows=MagickMax(image->rows,reconstruct_image->rows);
629  columns=MagickMax(image->columns,reconstruct_image->columns);
630  area=0.0;
631  image_view=AcquireVirtualCacheView(image,exception);
632  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
633 #if defined(MAGICKCORE_OPENMP_SUPPORT)
634  #pragma omp parallel for schedule(static) shared(status) \
635  magick_number_threads(image,image,rows,1) reduction(+:area)
636 #endif
637  for (y=0; y < (ssize_t) rows; y++)
638  {
639  double
640  channel_distortion[MaxPixelChannels+1];
641 
642  register const Quantum
643  *magick_restrict p,
644  *magick_restrict q;
645 
646  register ssize_t
647  x;
648 
649  if (status == MagickFalse)
650  continue;
651  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
652  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
653  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
654  {
655  status=MagickFalse;
656  continue;
657  }
658  (void) memset(channel_distortion,0,sizeof(channel_distortion));
659  for (x=0; x < (ssize_t) columns; x++)
660  {
661  double
662  Da,
663  Sa;
664 
665  register ssize_t
666  i;
667 
668  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
669  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
670  {
671  p+=GetPixelChannels(image);
672  q+=GetPixelChannels(reconstruct_image);
673  continue;
674  }
675  Sa=QuantumScale*GetPixelAlpha(image,p);
676  Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
677  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
678  {
679  double
680  distance;
681 
682  PixelChannel channel = GetPixelChannelChannel(image,i);
683  PixelTrait traits = GetPixelChannelTraits(image,channel);
684  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
685  channel);
686  if ((traits == UndefinedPixelTrait) ||
687  (reconstruct_traits == UndefinedPixelTrait) ||
688  ((reconstruct_traits & UpdatePixelTrait) == 0))
689  continue;
690  if (channel == AlphaPixelChannel)
691  distance=QuantumScale*fabs((double) p[i]-
692  GetPixelChannel(reconstruct_image,channel,q));
693  else
694  distance=QuantumScale*fabs(Sa*p[i]-Da*
695  GetPixelChannel(reconstruct_image,channel,q));
696  channel_distortion[i]+=distance;
697  channel_distortion[CompositePixelChannel]+=distance;
698  }
699  area++;
700  p+=GetPixelChannels(image);
701  q+=GetPixelChannels(reconstruct_image);
702  }
703 #if defined(MAGICKCORE_OPENMP_SUPPORT)
704  #pragma omp critical (MagickCore_GetMeanAbsoluteError)
705 #endif
706  for (j=0; j <= MaxPixelChannels; j++)
707  distortion[j]+=channel_distortion[j];
708  }
709  reconstruct_view=DestroyCacheView(reconstruct_view);
710  image_view=DestroyCacheView(image_view);
711  area=PerceptibleReciprocal(area);
712  for (j=0; j <= MaxPixelChannels; j++)
713  distortion[j]*=area;
714  distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
715  return(status);
716 }
717 
719  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
720 {
721  CacheView
722  *image_view,
723  *reconstruct_view;
724 
726  status;
727 
728  double
729  area,
730  maximum_error,
731  mean_error;
732 
733  size_t
734  columns,
735  rows;
736 
737  ssize_t
738  y;
739 
740  status=MagickTrue;
741  area=0.0;
742  maximum_error=0.0;
743  mean_error=0.0;
744  rows=MagickMax(image->rows,reconstruct_image->rows);
745  columns=MagickMax(image->columns,reconstruct_image->columns);
746  image_view=AcquireVirtualCacheView(image,exception);
747  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
748  for (y=0; y < (ssize_t) rows; y++)
749  {
750  register const Quantum
751  *magick_restrict p,
752  *magick_restrict q;
753 
754  register ssize_t
755  x;
756 
757  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
758  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
759  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
760  {
761  status=MagickFalse;
762  break;
763  }
764  for (x=0; x < (ssize_t) columns; x++)
765  {
766  double
767  Da,
768  Sa;
769 
770  register ssize_t
771  i;
772 
773  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
774  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
775  {
776  p+=GetPixelChannels(image);
777  q+=GetPixelChannels(reconstruct_image);
778  continue;
779  }
780  Sa=QuantumScale*GetPixelAlpha(image,p);
781  Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
782  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
783  {
784  double
785  distance;
786 
787  PixelChannel channel = GetPixelChannelChannel(image,i);
788  PixelTrait traits = GetPixelChannelTraits(image,channel);
789  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
790  channel);
791  if ((traits == UndefinedPixelTrait) ||
792  (reconstruct_traits == UndefinedPixelTrait) ||
793  ((reconstruct_traits & UpdatePixelTrait) == 0))
794  continue;
795  if (channel == AlphaPixelChannel)
796  distance=fabs((double) p[i]-
797  GetPixelChannel(reconstruct_image,channel,q));
798  else
799  distance=fabs(Sa*p[i]-Da*
800  GetPixelChannel(reconstruct_image,channel,q));
801  distortion[i]+=distance;
802  distortion[CompositePixelChannel]+=distance;
803  mean_error+=distance*distance;
804  if (distance > maximum_error)
805  maximum_error=distance;
806  area++;
807  }
808  p+=GetPixelChannels(image);
809  q+=GetPixelChannels(reconstruct_image);
810  }
811  }
812  reconstruct_view=DestroyCacheView(reconstruct_view);
813  image_view=DestroyCacheView(image_view);
814  image->error.mean_error_per_pixel=distortion[CompositePixelChannel]/area;
815  image->error.normalized_mean_error=QuantumScale*QuantumScale*mean_error/area;
816  image->error.normalized_maximum_error=QuantumScale*maximum_error;
817  return(status);
818 }
819 
821  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
822 {
823  CacheView
824  *image_view,
825  *reconstruct_view;
826 
827  double
828  area;
829 
831  status;
832 
833  register ssize_t
834  j;
835 
836  size_t
837  columns,
838  rows;
839 
840  ssize_t
841  y;
842 
843  status=MagickTrue;
844  rows=MagickMax(image->rows,reconstruct_image->rows);
845  columns=MagickMax(image->columns,reconstruct_image->columns);
846  area=0.0;
847  image_view=AcquireVirtualCacheView(image,exception);
848  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
849 #if defined(MAGICKCORE_OPENMP_SUPPORT)
850  #pragma omp parallel for schedule(static) shared(status) \
851  magick_number_threads(image,image,rows,1) reduction(+:area)
852 #endif
853  for (y=0; y < (ssize_t) rows; y++)
854  {
855  double
856  channel_distortion[MaxPixelChannels+1];
857 
858  register const Quantum
859  *magick_restrict p,
860  *magick_restrict q;
861 
862  register ssize_t
863  x;
864 
865  if (status == MagickFalse)
866  continue;
867  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
868  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
869  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
870  {
871  status=MagickFalse;
872  continue;
873  }
874  (void) memset(channel_distortion,0,sizeof(channel_distortion));
875  for (x=0; x < (ssize_t) columns; x++)
876  {
877  double
878  Da,
879  Sa;
880 
881  register ssize_t
882  i;
883 
884  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
885  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
886  {
887  p+=GetPixelChannels(image);
888  q+=GetPixelChannels(reconstruct_image);
889  continue;
890  }
891  Sa=QuantumScale*GetPixelAlpha(image,p);
892  Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
893  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
894  {
895  double
896  distance;
897 
898  PixelChannel channel = GetPixelChannelChannel(image,i);
899  PixelTrait traits = GetPixelChannelTraits(image,channel);
900  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
901  channel);
902  if ((traits == UndefinedPixelTrait) ||
903  (reconstruct_traits == UndefinedPixelTrait) ||
904  ((reconstruct_traits & UpdatePixelTrait) == 0))
905  continue;
906  if (channel == AlphaPixelChannel)
907  distance=QuantumScale*(p[i]-GetPixelChannel(reconstruct_image,
908  channel,q));
909  else
910  distance=QuantumScale*(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
911  channel,q));
912  channel_distortion[i]+=distance*distance;
913  channel_distortion[CompositePixelChannel]+=distance*distance;
914  }
915  area++;
916  p+=GetPixelChannels(image);
917  q+=GetPixelChannels(reconstruct_image);
918  }
919 #if defined(MAGICKCORE_OPENMP_SUPPORT)
920  #pragma omp critical (MagickCore_GetMeanSquaredError)
921 #endif
922  for (j=0; j <= MaxPixelChannels; j++)
923  distortion[j]+=channel_distortion[j];
924  }
925  reconstruct_view=DestroyCacheView(reconstruct_view);
926  image_view=DestroyCacheView(image_view);
927  area=PerceptibleReciprocal(area);
928  for (j=0; j <= MaxPixelChannels; j++)
929  distortion[j]*=area;
930  distortion[CompositePixelChannel]/=GetImageChannels(image);
931  return(status);
932 }
933 
935  const Image *image,const Image *reconstruct_image,double *distortion,
936  ExceptionInfo *exception)
937 {
938 #define SimilarityImageTag "Similarity/Image"
939 
940  CacheView
941  *image_view,
942  *reconstruct_view;
943 
945  *image_statistics,
946  *reconstruct_statistics;
947 
948  double
949  area;
950 
952  status;
953 
955  progress;
956 
957  register ssize_t
958  i;
959 
960  size_t
961  columns,
962  rows;
963 
964  ssize_t
965  y;
966 
967  /*
968  Normalize to account for variation due to lighting and exposure condition.
969  */
970  image_statistics=GetImageStatistics(image,exception);
971  reconstruct_statistics=GetImageStatistics(reconstruct_image,exception);
972  if ((image_statistics == (ChannelStatistics *) NULL) ||
973  (reconstruct_statistics == (ChannelStatistics *) NULL))
974  {
975  if (image_statistics != (ChannelStatistics *) NULL)
976  image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
977  image_statistics);
978  if (reconstruct_statistics != (ChannelStatistics *) NULL)
979  reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
980  reconstruct_statistics);
981  return(MagickFalse);
982  }
983  status=MagickTrue;
984  progress=0;
985  for (i=0; i <= MaxPixelChannels; i++)
986  distortion[i]=0.0;
987  rows=MagickMax(image->rows,reconstruct_image->rows);
988  columns=MagickMax(image->columns,reconstruct_image->columns);
989  area=0.0;
990  image_view=AcquireVirtualCacheView(image,exception);
991  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
992  for (y=0; y < (ssize_t) rows; y++)
993  {
994  register const Quantum
995  *magick_restrict p,
996  *magick_restrict q;
997 
998  register ssize_t
999  x;
1000 
1001  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1002  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1003  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1004  {
1005  status=MagickFalse;
1006  break;
1007  }
1008  for (x=0; x < (ssize_t) columns; x++)
1009  {
1010  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1011  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1012  {
1013  p+=GetPixelChannels(image);
1014  q+=GetPixelChannels(reconstruct_image);
1015  continue;
1016  }
1017  area++;
1018  p+=GetPixelChannels(image);
1019  q+=GetPixelChannels(reconstruct_image);
1020  }
1021  }
1022  area=PerceptibleReciprocal(area);
1023  for (y=0; y < (ssize_t) rows; y++)
1024  {
1025  register const Quantum
1026  *magick_restrict p,
1027  *magick_restrict q;
1028 
1029  register ssize_t
1030  x;
1031 
1032  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1033  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1034  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1035  {
1036  status=MagickFalse;
1037  break;
1038  }
1039  for (x=0; x < (ssize_t) columns; x++)
1040  {
1041  double
1042  Da,
1043  Sa;
1044 
1045  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1046  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1047  {
1048  p+=GetPixelChannels(image);
1049  q+=GetPixelChannels(reconstruct_image);
1050  continue;
1051  }
1052  Sa=QuantumScale*GetPixelAlpha(image,p);
1053  Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
1054  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1055  {
1056  PixelChannel channel = GetPixelChannelChannel(image,i);
1057  PixelTrait traits = GetPixelChannelTraits(image,channel);
1058  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1059  channel);
1060  if ((traits == UndefinedPixelTrait) ||
1061  (reconstruct_traits == UndefinedPixelTrait) ||
1062  ((reconstruct_traits & UpdatePixelTrait) == 0))
1063  continue;
1064  if (channel == AlphaPixelChannel)
1065  {
1066  distortion[i]+=area*QuantumScale*(p[i]-
1067  image_statistics[channel].mean)*(GetPixelChannel(
1068  reconstruct_image,channel,q)-
1069  reconstruct_statistics[channel].mean);
1070  }
1071  else
1072  {
1073  distortion[i]+=area*QuantumScale*(Sa*p[i]-
1074  image_statistics[channel].mean)*(Da*GetPixelChannel(
1075  reconstruct_image,channel,q)-
1076  reconstruct_statistics[channel].mean);
1077  }
1078  }
1079  p+=GetPixelChannels(image);
1080  q+=GetPixelChannels(reconstruct_image);
1081  }
1082  if (image->progress_monitor != (MagickProgressMonitor) NULL)
1083  {
1085  proceed;
1086 
1087 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1088  #pragma omp atomic
1089 #endif
1090  progress++;
1091  proceed=SetImageProgress(image,SimilarityImageTag,progress,rows);
1092  if (proceed == MagickFalse)
1093  {
1094  status=MagickFalse;
1095  break;
1096  }
1097  }
1098  }
1099  reconstruct_view=DestroyCacheView(reconstruct_view);
1100  image_view=DestroyCacheView(image_view);
1101  /*
1102  Divide by the standard deviation.
1103  */
1104  distortion[CompositePixelChannel]=0.0;
1105  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1106  {
1107  double
1108  gamma;
1109 
1110  PixelChannel channel = GetPixelChannelChannel(image,i);
1111  gamma=image_statistics[channel].standard_deviation*
1112  reconstruct_statistics[channel].standard_deviation;
1113  gamma=PerceptibleReciprocal(gamma);
1114  distortion[i]=QuantumRange*gamma*distortion[i];
1115  distortion[CompositePixelChannel]+=distortion[i]*distortion[i];
1116  }
1117  distortion[CompositePixelChannel]=sqrt(distortion[CompositePixelChannel]/
1118  GetImageChannels(image));
1119  /*
1120  Free resources.
1121  */
1122  reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1123  reconstruct_statistics);
1124  image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1125  image_statistics);
1126  return(status);
1127 }
1128 
1130  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1131 {
1132  CacheView
1133  *image_view,
1134  *reconstruct_view;
1135 
1137  status;
1138 
1139  size_t
1140  columns,
1141  rows;
1142 
1143  ssize_t
1144  y;
1145 
1146  status=MagickTrue;
1147  rows=MagickMax(image->rows,reconstruct_image->rows);
1148  columns=MagickMax(image->columns,reconstruct_image->columns);
1149  image_view=AcquireVirtualCacheView(image,exception);
1150  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1151 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1152  #pragma omp parallel for schedule(static) shared(status) \
1153  magick_number_threads(image,image,rows,1)
1154 #endif
1155  for (y=0; y < (ssize_t) rows; y++)
1156  {
1157  double
1158  channel_distortion[MaxPixelChannels+1];
1159 
1160  register const Quantum
1161  *magick_restrict p,
1162  *magick_restrict q;
1163 
1164  register ssize_t
1165  j,
1166  x;
1167 
1168  if (status == MagickFalse)
1169  continue;
1170  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1171  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1172  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1173  {
1174  status=MagickFalse;
1175  continue;
1176  }
1177  (void) memset(channel_distortion,0,sizeof(channel_distortion));
1178  for (x=0; x < (ssize_t) columns; x++)
1179  {
1180  double
1181  Da,
1182  Sa;
1183 
1184  register ssize_t
1185  i;
1186 
1187  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1188  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1189  {
1190  p+=GetPixelChannels(image);
1191  q+=GetPixelChannels(reconstruct_image);
1192  continue;
1193  }
1194  Sa=QuantumScale*GetPixelAlpha(image,p);
1195  Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
1196  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1197  {
1198  double
1199  distance;
1200 
1201  PixelChannel channel = GetPixelChannelChannel(image,i);
1202  PixelTrait traits = GetPixelChannelTraits(image,channel);
1203  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1204  channel);
1205  if ((traits == UndefinedPixelTrait) ||
1206  (reconstruct_traits == UndefinedPixelTrait) ||
1207  ((reconstruct_traits & UpdatePixelTrait) == 0))
1208  continue;
1209  if (channel == AlphaPixelChannel)
1210  distance=QuantumScale*fabs((double) p[i]-
1211  GetPixelChannel(reconstruct_image,channel,q));
1212  else
1213  distance=QuantumScale*fabs(Sa*p[i]-Da*
1214  GetPixelChannel(reconstruct_image,channel,q));
1215  if (distance > channel_distortion[i])
1216  channel_distortion[i]=distance;
1217  if (distance > channel_distortion[CompositePixelChannel])
1218  channel_distortion[CompositePixelChannel]=distance;
1219  }
1220  p+=GetPixelChannels(image);
1221  q+=GetPixelChannels(reconstruct_image);
1222  }
1223 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1224  #pragma omp critical (MagickCore_GetPeakAbsoluteError)
1225 #endif
1226  for (j=0; j <= MaxPixelChannels; j++)
1227  if (channel_distortion[j] > distortion[j])
1228  distortion[j]=channel_distortion[j];
1229  }
1230  reconstruct_view=DestroyCacheView(reconstruct_view);
1231  image_view=DestroyCacheView(image_view);
1232  return(status);
1233 }
1234 
1235 static inline double MagickLog10(const double x)
1236 {
1237 #define Log10Epsilon (1.0e-11)
1238 
1239  if (fabs(x) < Log10Epsilon)
1240  return(log10(Log10Epsilon));
1241  return(log10(fabs(x)));
1242 }
1243 
1245  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1246 {
1248  status;
1249 
1250  register ssize_t
1251  i;
1252 
1253  status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1254  for (i=0; i <= MaxPixelChannels; i++)
1255  if (fabs(distortion[i]) < MagickEpsilon)
1256  distortion[i]=INFINITY;
1257  else
1258  distortion[i]=10.0*MagickLog10(1.0)-10.0*MagickLog10(distortion[i]);
1259  return(status);
1260 }
1261 
1263  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1264 {
1266  *channel_phash,
1267  *reconstruct_phash;
1268 
1269  const char
1270  *artifact;
1271 
1273  normalize;
1274 
1275  ssize_t
1276  channel;
1277 
1278  /*
1279  Compute perceptual hash in the sRGB colorspace.
1280  */
1281  channel_phash=GetImagePerceptualHash(image,exception);
1282  if (channel_phash == (ChannelPerceptualHash *) NULL)
1283  return(MagickFalse);
1284  reconstruct_phash=GetImagePerceptualHash(reconstruct_image,exception);
1285  if (reconstruct_phash == (ChannelPerceptualHash *) NULL)
1286  {
1288  channel_phash);
1289  return(MagickFalse);
1290  }
1291  artifact=GetImageArtifact(image,"phash:normalize");
1292  normalize=(artifact == (const char *) NULL) ||
1293  (IsStringTrue(artifact) == MagickFalse) ? MagickFalse : MagickTrue;
1294 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1295  #pragma omp parallel for schedule(static)
1296 #endif
1297  for (channel=0; channel < MaxPixelChannels; channel++)
1298  {
1299  double
1300  difference;
1301 
1302  register ssize_t
1303  i;
1304 
1305  difference=0.0;
1306  for (i=0; i < MaximumNumberOfImageMoments; i++)
1307  {
1308  double
1309  alpha,
1310  beta;
1311 
1312  register ssize_t
1313  j;
1314 
1315  for (j=0; j < (ssize_t) channel_phash[0].number_colorspaces; j++)
1316  {
1317  alpha=channel_phash[channel].phash[j][i];
1318  beta=reconstruct_phash[channel].phash[j][i];
1319  if (normalize == MagickFalse)
1320  difference+=(beta-alpha)*(beta-alpha);
1321  else
1322  difference=sqrt((beta-alpha)*(beta-alpha)/
1323  channel_phash[0].number_channels);
1324  }
1325  }
1326  distortion[channel]+=difference;
1327 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1328  #pragma omp critical (MagickCore_GetPerceptualHashDistortion)
1329 #endif
1330  distortion[CompositePixelChannel]+=difference;
1331  }
1332  /*
1333  Free resources.
1334  */
1335  reconstruct_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1336  reconstruct_phash);
1337  channel_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(channel_phash);
1338  return(MagickTrue);
1339 }
1340 
1342  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1343 {
1345  status;
1346 
1347  register ssize_t
1348  i;
1349 
1350  status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1351  for (i=0; i <= MaxPixelChannels; i++)
1352  distortion[i]=sqrt(distortion[i]);
1353  return(status);
1354 }
1355 
1357  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1358 {
1359 #define SSIMRadius 5.0
1360 #define SSIMSigma 1.5
1361 #define SSIMBlocksize 8
1362 #define SSIMK1 0.01
1363 #define SSIMK2 0.03
1364 #define SSIML 1.0
1365 
1366  CacheView
1367  *image_view,
1368  *reconstruct_view;
1369 
1370  char
1371  geometry[MagickPathExtent];
1372 
1373  const char
1374  *artifact;
1375 
1376  double
1377  c1,
1378  c2,
1379  radius,
1380  sigma;
1381 
1382  KernelInfo
1383  *kernel_info;
1384 
1386  status;
1387 
1388  register ssize_t
1389  i;
1390 
1391  size_t
1392  columns,
1393  rows;
1394 
1395  ssize_t
1396  y;
1397 
1398  /*
1399  Compute structural similarity index @
1400  https://en.wikipedia.org/wiki/Structural_similarity.
1401  */
1402  radius=SSIMRadius;
1403  artifact=GetImageArtifact(image,"compare:ssim-radius");
1404  if (artifact != (const char *) NULL)
1405  radius=StringToDouble(artifact,(char **) NULL);
1406  sigma=SSIMSigma;
1407  artifact=GetImageArtifact(image,"compare:ssim-sigma");
1408  if (artifact != (const char *) NULL)
1409  sigma=StringToDouble(artifact,(char **) NULL);
1410  (void) FormatLocaleString(geometry,MagickPathExtent,"gaussian:%.20gx%.20g",
1411  radius,sigma);
1412  kernel_info=AcquireKernelInfo(geometry,exception);
1413  if (kernel_info == (KernelInfo *) NULL)
1414  ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
1415  image->filename);
1416  c1=pow(SSIMK1*SSIML,2.0);
1417  artifact=GetImageArtifact(image,"compare:ssim-k1");
1418  if (artifact != (const char *) NULL)
1419  c1=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
1420  c2=pow(SSIMK2*SSIML,2.0);
1421  artifact=GetImageArtifact(image,"compare:ssim-k2");
1422  if (artifact != (const char *) NULL)
1423  c2=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
1424  status=MagickTrue;
1425  rows=MagickMax(image->rows,reconstruct_image->rows);
1426  columns=MagickMax(image->columns,reconstruct_image->columns);
1427  image_view=AcquireVirtualCacheView(image,exception);
1428  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1429 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1430  #pragma omp parallel for schedule(static) shared(status) \
1431  magick_number_threads(image,reconstruct_image,rows,1)
1432 #endif
1433  for (y=0; y < (ssize_t) rows; y++)
1434  {
1435  double
1436  channel_distortion[MaxPixelChannels+1];
1437 
1438  register const Quantum
1439  *magick_restrict p,
1440  *magick_restrict q;
1441 
1442  register ssize_t
1443  i,
1444  x;
1445 
1446  if (status == MagickFalse)
1447  continue;
1448  p=GetCacheViewVirtualPixels(image_view,-((ssize_t) kernel_info->width/2L),y-
1449  ((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
1450  kernel_info->height,exception);
1451  q=GetCacheViewVirtualPixels(reconstruct_view,-((ssize_t) kernel_info->width/
1452  2L),y-((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
1453  kernel_info->height,exception);
1454  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1455  {
1456  status=MagickFalse;
1457  continue;
1458  }
1459  (void) memset(channel_distortion,0,sizeof(channel_distortion));
1460  for (x=0; x < (ssize_t) columns; x++)
1461  {
1462  double
1463  x_pixel_mu[MaxPixelChannels+1],
1464  x_pixel_sigma_squared[MaxPixelChannels+1],
1465  xy_sigma[MaxPixelChannels+1],
1466  y_pixel_mu[MaxPixelChannels+1],
1467  y_pixel_sigma_squared[MaxPixelChannels+1];
1468 
1469  register const Quantum
1470  *magick_restrict reference,
1471  *magick_restrict target;
1472 
1473  register MagickRealType
1474  *k;
1475 
1476  ssize_t
1477  v;
1478 
1479  (void) memset(x_pixel_mu,0,sizeof(x_pixel_mu));
1480  (void) memset(x_pixel_sigma_squared,0,sizeof(x_pixel_sigma_squared));
1481  (void) memset(xy_sigma,0,sizeof(xy_sigma));
1482  (void) memset(x_pixel_sigma_squared,0,sizeof(y_pixel_sigma_squared));
1483  (void) memset(y_pixel_mu,0,sizeof(y_pixel_mu));
1484  (void) memset(y_pixel_sigma_squared,0,sizeof(y_pixel_sigma_squared));
1485  k=kernel_info->values;
1486  reference=p;
1487  target=q;
1488  for (v=0; v < (ssize_t) kernel_info->height; v++)
1489  {
1490  register ssize_t
1491  u;
1492 
1493  for (u=0; u < (ssize_t) kernel_info->width; u++)
1494  {
1495  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1496  {
1497  double
1498  x_pixel,
1499  y_pixel;
1500 
1501  PixelChannel channel = GetPixelChannelChannel(image,i);
1502  PixelTrait traits = GetPixelChannelTraits(image,channel);
1503  PixelTrait reconstruct_traits = GetPixelChannelTraits(
1504  reconstruct_image,channel);
1505  if ((traits == UndefinedPixelTrait) ||
1506  (reconstruct_traits == UndefinedPixelTrait) ||
1507  ((reconstruct_traits & UpdatePixelTrait) == 0))
1508  continue;
1509  x_pixel=QuantumScale*reference[i];
1510  x_pixel_mu[i]+=(*k)*x_pixel;
1511  x_pixel_sigma_squared[i]+=(*k)*x_pixel*x_pixel;
1512  y_pixel=QuantumScale*
1513  GetPixelChannel(reconstruct_image,channel,target);
1514  y_pixel_mu[i]+=(*k)*y_pixel;
1515  y_pixel_sigma_squared[i]+=(*k)*y_pixel*y_pixel;
1516  xy_sigma[i]+=(*k)*x_pixel*y_pixel;
1517  }
1518  k++;
1519  reference+=GetPixelChannels(image);
1520  target+=GetPixelChannels(reconstruct_image);
1521  }
1522  reference+=GetPixelChannels(image)*columns;
1523  target+=GetPixelChannels(reconstruct_image)*columns;
1524  }
1525  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1526  {
1527  double
1528  ssim,
1529  x_pixel_mu_squared,
1530  x_pixel_sigmas_squared,
1531  xy_mu,
1532  xy_sigmas,
1533  y_pixel_mu_squared,
1534  y_pixel_sigmas_squared;
1535 
1536  PixelChannel channel = GetPixelChannelChannel(image,i);
1537  PixelTrait traits = GetPixelChannelTraits(image,channel);
1538  PixelTrait reconstruct_traits = GetPixelChannelTraits(
1539  reconstruct_image,channel);
1540  if ((traits == UndefinedPixelTrait) ||
1541  (reconstruct_traits == UndefinedPixelTrait) ||
1542  ((reconstruct_traits & UpdatePixelTrait) == 0))
1543  continue;
1544  x_pixel_mu_squared=x_pixel_mu[i]*x_pixel_mu[i];
1545  y_pixel_mu_squared=y_pixel_mu[i]*y_pixel_mu[i];
1546  xy_mu=x_pixel_mu[i]*y_pixel_mu[i];
1547  xy_sigmas=xy_sigma[i]-xy_mu;
1548  x_pixel_sigmas_squared=x_pixel_sigma_squared[i]-x_pixel_mu_squared;
1549  y_pixel_sigmas_squared=y_pixel_sigma_squared[i]-y_pixel_mu_squared;
1550  ssim=((2.0*xy_mu+c1)*(2.0*xy_sigmas+c2))/
1551  ((x_pixel_mu_squared+y_pixel_mu_squared+c1)*
1552  (x_pixel_sigmas_squared+y_pixel_sigmas_squared+c2));
1553  channel_distortion[i]+=ssim;
1554  channel_distortion[CompositePixelChannel]+=ssim;
1555  }
1556  p+=GetPixelChannels(image);
1557  q+=GetPixelChannels(reconstruct_image);
1558  }
1559 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1560  #pragma omp critical (MagickCore_GetStructuralSimilarityDistortion)
1561 #endif
1562  for (i=0; i <= MaxPixelChannels; i++)
1563  distortion[i]+=channel_distortion[i];
1564  }
1565  image_view=DestroyCacheView(image_view);
1566  reconstruct_view=DestroyCacheView(reconstruct_view);
1567  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1568  {
1569  PixelChannel channel = GetPixelChannelChannel(image,i);
1570  PixelTrait traits = GetPixelChannelTraits(image,channel);
1571  if ((traits == UndefinedPixelTrait) || ((traits & UpdatePixelTrait) == 0))
1572  continue;
1573  distortion[i]/=((double) columns*rows);
1574  }
1575  distortion[CompositePixelChannel]/=((double) columns*rows);
1576  distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
1577  kernel_info=DestroyKernelInfo(kernel_info);
1578  return(status);
1579 }
1580 
1582  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1583 {
1585  status;
1586 
1587  register ssize_t
1588  i;
1589 
1590  status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1591  distortion,exception);
1592  for (i=0; i <= MaxPixelChannels; i++)
1593  distortion[i]=(1.0-(distortion[i]))/2.0;
1594  return(status);
1595 }
1596 
1598  const Image *reconstruct_image,const MetricType metric,double *distortion,
1599  ExceptionInfo *exception)
1600 {
1601  double
1602  *channel_distortion;
1603 
1605  status;
1606 
1607  size_t
1608  length;
1609 
1610  assert(image != (Image *) NULL);
1611  assert(image->signature == MagickCoreSignature);
1612  if (image->debug != MagickFalse)
1613  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1614  assert(reconstruct_image != (const Image *) NULL);
1615  assert(reconstruct_image->signature == MagickCoreSignature);
1616  assert(distortion != (double *) NULL);
1617  *distortion=0.0;
1618  if (image->debug != MagickFalse)
1619  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1620  /*
1621  Get image distortion.
1622  */
1623  length=MaxPixelChannels+1;
1624  channel_distortion=(double *) AcquireQuantumMemory(length,
1625  sizeof(*channel_distortion));
1626  if (channel_distortion == (double *) NULL)
1627  ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1628  (void) memset(channel_distortion,0,length*
1629  sizeof(*channel_distortion));
1630  switch (metric)
1631  {
1632  case AbsoluteErrorMetric:
1633  {
1634  status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1635  exception);
1636  break;
1637  }
1638  case FuzzErrorMetric:
1639  {
1640  status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1641  exception);
1642  break;
1643  }
1645  {
1646  status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1647  channel_distortion,exception);
1648  break;
1649  }
1651  {
1652  status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1653  exception);
1654  break;
1655  }
1657  {
1658  status=GetMeanSquaredDistortion(image,reconstruct_image,
1659  channel_distortion,exception);
1660  break;
1661  }
1663  default:
1664  {
1665  status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1666  channel_distortion,exception);
1667  break;
1668  }
1670  {
1671  status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1672  channel_distortion,exception);
1673  break;
1674  }
1676  {
1677  status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1678  channel_distortion,exception);
1679  break;
1680  }
1682  {
1683  status=GetPerceptualHashDistortion(image,reconstruct_image,
1684  channel_distortion,exception);
1685  break;
1686  }
1688  {
1689  status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1690  channel_distortion,exception);
1691  break;
1692  }
1694  {
1695  status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1696  channel_distortion,exception);
1697  break;
1698  }
1700  {
1701  status=GetStructuralDisimilarityDistortion(image,reconstruct_image,
1702  channel_distortion,exception);
1703  break;
1704  }
1705  }
1706  *distortion=channel_distortion[CompositePixelChannel];
1707  channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1708  (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
1709  *distortion);
1710  return(status);
1711 }
1712 
1713 /*
1714 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1715 % %
1716 % %
1717 % %
1718 % G e t I m a g e D i s t o r t i o n s %
1719 % %
1720 % %
1721 % %
1722 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1723 %
1724 % GetImageDistortions() compares the pixel channels of an image to a
1725 % reconstructed image and returns the specified distortion metric for each
1726 % channel.
1727 %
1728 % The format of the GetImageDistortions method is:
1729 %
1730 % double *GetImageDistortions(const Image *image,
1731 % const Image *reconstruct_image,const MetricType metric,
1732 % ExceptionInfo *exception)
1733 %
1734 % A description of each parameter follows:
1735 %
1736 % o image: the image.
1737 %
1738 % o reconstruct_image: the reconstruct image.
1739 %
1740 % o metric: the metric.
1741 %
1742 % o exception: return any errors or warnings in this structure.
1743 %
1744 */
1746  const Image *reconstruct_image,const MetricType metric,
1747  ExceptionInfo *exception)
1748 {
1749  double
1750  *channel_distortion;
1751 
1753  status;
1754 
1755  size_t
1756  length;
1757 
1758  assert(image != (Image *) NULL);
1759  assert(image->signature == MagickCoreSignature);
1760  if (image->debug != MagickFalse)
1761  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1762  assert(reconstruct_image != (const Image *) NULL);
1763  assert(reconstruct_image->signature == MagickCoreSignature);
1764  if (image->debug != MagickFalse)
1765  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1766  /*
1767  Get image distortion.
1768  */
1769  length=MaxPixelChannels+1UL;
1770  channel_distortion=(double *) AcquireQuantumMemory(length,
1771  sizeof(*channel_distortion));
1772  if (channel_distortion == (double *) NULL)
1773  ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1774  (void) memset(channel_distortion,0,length*
1775  sizeof(*channel_distortion));
1776  status=MagickTrue;
1777  switch (metric)
1778  {
1779  case AbsoluteErrorMetric:
1780  {
1781  status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1782  exception);
1783  break;
1784  }
1785  case FuzzErrorMetric:
1786  {
1787  status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1788  exception);
1789  break;
1790  }
1792  {
1793  status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1794  channel_distortion,exception);
1795  break;
1796  }
1798  {
1799  status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1800  exception);
1801  break;
1802  }
1804  {
1805  status=GetMeanSquaredDistortion(image,reconstruct_image,
1806  channel_distortion,exception);
1807  break;
1808  }
1810  default:
1811  {
1812  status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1813  channel_distortion,exception);
1814  break;
1815  }
1817  {
1818  status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1819  channel_distortion,exception);
1820  break;
1821  }
1823  {
1824  status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1825  channel_distortion,exception);
1826  break;
1827  }
1829  {
1830  status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1831  channel_distortion,exception);
1832  break;
1833  }
1835  {
1836  status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1837  channel_distortion,exception);
1838  break;
1839  }
1841  {
1842  status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1843  channel_distortion,exception);
1844  break;
1845  }
1847  {
1848  status=GetStructuralDisimilarityDistortion(image,reconstruct_image,
1849  channel_distortion,exception);
1850  break;
1851  }
1852  }
1853  if (status == MagickFalse)
1854  {
1855  channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1856  return((double *) NULL);
1857  }
1858  return(channel_distortion);
1859 }
1860 
1861 /*
1862 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1863 % %
1864 % %
1865 % %
1866 % I s I m a g e s E q u a l %
1867 % %
1868 % %
1869 % %
1870 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1871 %
1872 % IsImagesEqual() compare the pixels of two images and returns immediately
1873 % if any pixel is not identical.
1874 %
1875 % The format of the IsImagesEqual method is:
1876 %
1877 % MagickBooleanType IsImagesEqual(const Image *image,
1878 % const Image *reconstruct_image,ExceptionInfo *exception)
1879 %
1880 % A description of each parameter follows.
1881 %
1882 % o image: the image.
1883 %
1884 % o reconstruct_image: the reconstruct image.
1885 %
1886 % o exception: return any errors or warnings in this structure.
1887 %
1888 */
1890  const Image *reconstruct_image,ExceptionInfo *exception)
1891 {
1892  CacheView
1893  *image_view,
1894  *reconstruct_view;
1895 
1896  size_t
1897  columns,
1898  rows;
1899 
1900  ssize_t
1901  y;
1902 
1903  assert(image != (Image *) NULL);
1904  assert(image->signature == MagickCoreSignature);
1905  assert(reconstruct_image != (const Image *) NULL);
1906  assert(reconstruct_image->signature == MagickCoreSignature);
1907  rows=MagickMax(image->rows,reconstruct_image->rows);
1908  columns=MagickMax(image->columns,reconstruct_image->columns);
1909  image_view=AcquireVirtualCacheView(image,exception);
1910  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1911  for (y=0; y < (ssize_t) rows; y++)
1912  {
1913  register const Quantum
1914  *magick_restrict p,
1915  *magick_restrict q;
1916 
1917  register ssize_t
1918  x;
1919 
1920  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1921  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1922  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
1923  break;
1924  for (x=0; x < (ssize_t) columns; x++)
1925  {
1926  register ssize_t
1927  i;
1928 
1929  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1930  {
1931  double
1932  distance;
1933 
1934  PixelChannel channel = GetPixelChannelChannel(image,i);
1935  PixelTrait traits = GetPixelChannelTraits(image,channel);
1936  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1937  channel);
1938  if ((traits == UndefinedPixelTrait) ||
1939  (reconstruct_traits == UndefinedPixelTrait) ||
1940  ((reconstruct_traits & UpdatePixelTrait) == 0))
1941  continue;
1942  distance=fabs(p[i]-(double) GetPixelChannel(reconstruct_image,
1943  channel,q));
1944  if (distance >= MagickEpsilon)
1945  break;
1946  }
1947  if (i < (ssize_t) GetPixelChannels(image))
1948  break;
1949  p+=GetPixelChannels(image);
1950  q+=GetPixelChannels(reconstruct_image);
1951  }
1952  if (x < (ssize_t) columns)
1953  break;
1954  }
1955  reconstruct_view=DestroyCacheView(reconstruct_view);
1956  image_view=DestroyCacheView(image_view);
1957  return(y < (ssize_t) rows ? MagickFalse : MagickTrue);
1958 }
1959 
1960 /*
1961 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1962 % %
1963 % %
1964 % %
1965 % S e t I m a g e C o l o r M e t r i c %
1966 % %
1967 % %
1968 % %
1969 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1970 %
1971 % SetImageColorMetric() measures the difference between colors at each pixel
1972 % location of two images. A value other than 0 means the colors match
1973 % exactly. Otherwise an error measure is computed by summing over all
1974 % pixels in an image the distance squared in RGB space between each image
1975 % pixel and its corresponding pixel in the reconstruct image. The error
1976 % measure is assigned to these image members:
1977 %
1978 % o mean_error_per_pixel: The mean error for any single pixel in
1979 % the image.
1980 %
1981 % o normalized_mean_error: The normalized mean quantization error for
1982 % any single pixel in the image. This distance measure is normalized to
1983 % a range between 0 and 1. It is independent of the range of red, green,
1984 % and blue values in the image.
1985 %
1986 % o normalized_maximum_error: The normalized maximum quantization
1987 % error for any single pixel in the image. This distance measure is
1988 % normalized to a range between 0 and 1. It is independent of the range
1989 % of red, green, and blue values in your image.
1990 %
1991 % A small normalized mean square error, accessed as
1992 % image->normalized_mean_error, suggests the images are very similar in
1993 % spatial layout and color.
1994 %
1995 % The format of the SetImageColorMetric method is:
1996 %
1997 % MagickBooleanType SetImageColorMetric(Image *image,
1998 % const Image *reconstruct_image,ExceptionInfo *exception)
1999 %
2000 % A description of each parameter follows.
2001 %
2002 % o image: the image.
2003 %
2004 % o reconstruct_image: the reconstruct image.
2005 %
2006 % o exception: return any errors or warnings in this structure.
2007 %
2008 */
2010  const Image *reconstruct_image,ExceptionInfo *exception)
2011 {
2012  CacheView
2013  *image_view,
2014  *reconstruct_view;
2015 
2016  double
2017  area,
2018  maximum_error,
2019  mean_error,
2020  mean_error_per_pixel;
2021 
2023  status;
2024 
2025  size_t
2026  columns,
2027  rows;
2028 
2029  ssize_t
2030  y;
2031 
2032  assert(image != (Image *) NULL);
2033  assert(image->signature == MagickCoreSignature);
2034  assert(reconstruct_image != (const Image *) NULL);
2035  assert(reconstruct_image->signature == MagickCoreSignature);
2036  area=0.0;
2037  maximum_error=0.0;
2038  mean_error_per_pixel=0.0;
2039  mean_error=0.0;
2040  rows=MagickMax(image->rows,reconstruct_image->rows);
2041  columns=MagickMax(image->columns,reconstruct_image->columns);
2042  image_view=AcquireVirtualCacheView(image,exception);
2043  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
2044  for (y=0; y < (ssize_t) rows; y++)
2045  {
2046  register const Quantum
2047  *magick_restrict p,
2048  *magick_restrict q;
2049 
2050  register ssize_t
2051  x;
2052 
2053  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
2054  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
2055  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2056  break;
2057  for (x=0; x < (ssize_t) columns; x++)
2058  {
2059  register ssize_t
2060  i;
2061 
2062  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2063  {
2064  double
2065  distance;
2066 
2067  PixelChannel channel = GetPixelChannelChannel(image,i);
2068  PixelTrait traits = GetPixelChannelTraits(image,channel);
2069  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2070  channel);
2071  if ((traits == UndefinedPixelTrait) ||
2072  (reconstruct_traits == UndefinedPixelTrait) ||
2073  ((reconstruct_traits & UpdatePixelTrait) == 0))
2074  continue;
2075  distance=fabs(p[i]-(double) GetPixelChannel(reconstruct_image,
2076  channel,q));
2077  if (distance >= MagickEpsilon)
2078  {
2079  mean_error_per_pixel+=distance;
2080  mean_error+=distance*distance;
2081  if (distance > maximum_error)
2082  maximum_error=distance;
2083  }
2084  area++;
2085  }
2086  p+=GetPixelChannels(image);
2087  q+=GetPixelChannels(reconstruct_image);
2088  }
2089  }
2090  reconstruct_view=DestroyCacheView(reconstruct_view);
2091  image_view=DestroyCacheView(image_view);
2092  image->error.mean_error_per_pixel=(double) (mean_error_per_pixel/area);
2094  mean_error/area);
2095  image->error.normalized_maximum_error=(double) (QuantumScale*maximum_error);
2096  status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
2097  return(status);
2098 }
2099 
2100 /*
2101 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2102 % %
2103 % %
2104 % %
2105 % S i m i l a r i t y I m a g e %
2106 % %
2107 % %
2108 % %
2109 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2110 %
2111 % SimilarityImage() compares the reference image of the image and returns the
2112 % best match offset. In addition, it returns a similarity image such that an
2113 % exact match location is completely white and if none of the pixels match,
2114 % black, otherwise some gray level in-between.
2115 %
2116 % The format of the SimilarityImageImage method is:
2117 %
2118 % Image *SimilarityImage(const Image *image,const Image *reference,
2119 % const MetricType metric,const double similarity_threshold,
2120 % RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
2121 %
2122 % A description of each parameter follows:
2123 %
2124 % o image: the image.
2125 %
2126 % o reference: find an area of the image that closely resembles this image.
2127 %
2128 % o metric: the metric.
2129 %
2130 % o similarity_threshold: minimum distortion for (sub)image match.
2131 %
2132 % o offset: the best match offset of the reference image within the image.
2133 %
2134 % o similarity: the computed similarity between the images.
2135 %
2136 % o exception: return any errors or warnings in this structure.
2137 %
2138 */
2139 
2140 static double GetSimilarityMetric(const Image *image,const Image *reference,
2141  const MetricType metric,const ssize_t x_offset,const ssize_t y_offset,
2142  ExceptionInfo *exception)
2143 {
2144  double
2145  distortion;
2146 
2147  Image
2148  *similarity_image;
2149 
2151  status;
2152 
2154  geometry;
2155 
2156  SetGeometry(reference,&geometry);
2157  geometry.x=x_offset;
2158  geometry.y=y_offset;
2159  similarity_image=CropImage(image,&geometry,exception);
2160  if (similarity_image == (Image *) NULL)
2161  return(0.0);
2162  distortion=0.0;
2163  status=GetImageDistortion(similarity_image,reference,metric,&distortion,
2164  exception);
2165  similarity_image=DestroyImage(similarity_image);
2166  if (status == MagickFalse)
2167  return(0.0);
2168  return(distortion);
2169 }
2170 
2171 MagickExport Image *SimilarityImage(const Image *image,const Image *reference,
2172  const MetricType metric,const double similarity_threshold,
2173  RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
2174 {
2175 #define SimilarityImageTag "Similarity/Image"
2176 
2177  CacheView
2178  *similarity_view;
2179 
2180  Image
2181  *similarity_image;
2182 
2184  status;
2185 
2187  progress;
2188 
2189  ssize_t
2190  y;
2191 
2192  assert(image != (const Image *) NULL);
2193  assert(image->signature == MagickCoreSignature);
2194  if (image->debug != MagickFalse)
2195  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2196  assert(exception != (ExceptionInfo *) NULL);
2197  assert(exception->signature == MagickCoreSignature);
2198  assert(offset != (RectangleInfo *) NULL);
2199  SetGeometry(reference,offset);
2200  *similarity_metric=MagickMaximumValue;
2201  similarity_image=CloneImage(image,image->columns-reference->columns+1,
2202  image->rows-reference->rows+1,MagickTrue,exception);
2203  if (similarity_image == (Image *) NULL)
2204  return((Image *) NULL);
2205  status=SetImageStorageClass(similarity_image,DirectClass,exception);
2206  if (status == MagickFalse)
2207  {
2208  similarity_image=DestroyImage(similarity_image);
2209  return((Image *) NULL);
2210  }
2211  (void) SetImageAlphaChannel(similarity_image,DeactivateAlphaChannel,
2212  exception);
2213  /*
2214  Measure similarity of reference image against image.
2215  */
2216  status=MagickTrue;
2217  progress=0;
2218  similarity_view=AcquireAuthenticCacheView(similarity_image,exception);
2219 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2220  #pragma omp parallel for schedule(static) \
2221  shared(progress,status,similarity_metric) \
2222  magick_number_threads(image,image,image->rows-reference->rows+1,1)
2223 #endif
2224  for (y=0; y < (ssize_t) (image->rows-reference->rows+1); y++)
2225  {
2226  double
2227  similarity;
2228 
2229  register Quantum
2230  *magick_restrict q;
2231 
2232  register ssize_t
2233  x;
2234 
2235  if (status == MagickFalse)
2236  continue;
2237 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2238  #pragma omp flush(similarity_metric)
2239 #endif
2240  if (*similarity_metric <= similarity_threshold)
2241  continue;
2242  q=GetCacheViewAuthenticPixels(similarity_view,0,y,similarity_image->columns,
2243  1,exception);
2244  if (q == (Quantum *) NULL)
2245  {
2246  status=MagickFalse;
2247  continue;
2248  }
2249  for (x=0; x < (ssize_t) (image->columns-reference->columns+1); x++)
2250  {
2251  register ssize_t
2252  i;
2253 
2254 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2255  #pragma omp flush(similarity_metric)
2256 #endif
2257  if (*similarity_metric <= similarity_threshold)
2258  break;
2259  similarity=GetSimilarityMetric(image,reference,metric,x,y,exception);
2260 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2261  #pragma omp critical (MagickCore_SimilarityImage)
2262 #endif
2263  if ((metric == NormalizedCrossCorrelationErrorMetric) ||
2264  (metric == UndefinedErrorMetric))
2265  similarity=1.0-similarity;
2266  if (similarity < *similarity_metric)
2267  {
2268  offset->x=x;
2269  offset->y=y;
2270  *similarity_metric=similarity;
2271  }
2272  if (metric == PerceptualHashErrorMetric)
2273  similarity=MagickMin(0.01*similarity,1.0);
2274  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2275  {
2276  PixelChannel channel = GetPixelChannelChannel(image,i);
2277  PixelTrait traits = GetPixelChannelTraits(image,channel);
2278  PixelTrait similarity_traits=GetPixelChannelTraits(similarity_image,
2279  channel);
2280  if ((traits == UndefinedPixelTrait) ||
2281  (similarity_traits == UndefinedPixelTrait) ||
2282  ((similarity_traits & UpdatePixelTrait) == 0))
2283  continue;
2284  SetPixelChannel(similarity_image,channel,ClampToQuantum(QuantumRange-
2285  QuantumRange*similarity),q);
2286  }
2287  q+=GetPixelChannels(similarity_image);
2288  }
2289  if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
2290  status=MagickFalse;
2291  if (image->progress_monitor != (MagickProgressMonitor) NULL)
2292  {
2294  proceed;
2295 
2296 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2297  #pragma omp atomic
2298 #endif
2299  progress++;
2300  proceed=SetImageProgress(image,SimilarityImageTag,progress,image->rows);
2301  if (proceed == MagickFalse)
2302  status=MagickFalse;
2303  }
2304  }
2305  similarity_view=DestroyCacheView(similarity_view);
2306  if (status == MagickFalse)
2307  similarity_image=DestroyImage(similarity_image);
2308  return(similarity_image);
2309 }
size_t rows
Definition: image.h:172
#define magick_restrict
Definition: MagickCore.h:41
MagickDoubleType MagickRealType
Definition: magick-type.h:124
MagickExport CacheView * DestroyCacheView(CacheView *cache_view)
Definition: cache-view.c:252
static MagickBooleanType GetNormalizedCrossCorrelationDistortion(const Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c:934
double standard_deviation
Definition: statistic.h:35
static MagickBooleanType GetAbsoluteDistortion(const Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c:367
static MagickBooleanType GetPeakSignalToNoiseRatio(const Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c:1244
MagickExport Image * SimilarityImage(const Image *image, const Image *reference, const MetricType metric, const double similarity_threshold, RectangleInfo *offset, double *similarity_metric, ExceptionInfo *exception)
Definition: compare.c:2171
MagickProgressMonitor progress_monitor
Definition: image.h:303
static Quantum GetPixelAlpha(const Image *magick_restrict image, const Quantum *magick_restrict pixel)
double phash[MaximumNumberOfPerceptualColorspaces+1][MaximumNumberOfImageMoments+1]
Definition: statistic.h:78
size_t height
Definition: morphology.h:108
MagickExport ChannelStatistics * GetImageStatistics(const Image *image, ExceptionInfo *exception)
Definition: statistic.c:1967
MagickExport KernelInfo * DestroyKernelInfo(KernelInfo *kernel)
Definition: morphology.c:2268
static MagickBooleanType GetStructuralDisimilarityDistortion(const Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c:1581
#define ThrowFatalException(severity, tag)
size_t signature
Definition: exception.h:123
static MagickBooleanType GetMeanAbsoluteDistortion(const Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c:604
#define SimilarityImageTag
static Quantum GetPixelChannel(const Image *magick_restrict image, const PixelChannel channel, const Quantum *magick_restrict pixel)
static size_t GetImageChannels(const Image *image)
Definition: compare.c:109
MagickExport const char * GetImageArtifact(const Image *image, const char *artifact)
Definition: artifact.c:273
static double StringToDouble(const char *magick_restrict string, char **magick_restrict sentinal)
double mean_error_per_pixel
Definition: color.h:63
static PixelTrait GetPixelChannelTraits(const Image *magick_restrict image, const PixelChannel channel)
MetricType
Definition: compare.h:27
MagickExport ssize_t FormatLocaleString(char *magick_restrict string, const size_t length, const char *magick_restrict format,...)
Definition: locale.c:499
static void SetPixelViaPixelInfo(const Image *magick_restrict image, const PixelInfo *magick_restrict pixel_info, Quantum *magick_restrict pixel)
static Quantum GetPixelReadMask(const Image *magick_restrict image, const Quantum *magick_restrict pixel)
MagickExport const Quantum * GetCacheViewVirtualPixels(const CacheView *cache_view, const ssize_t x, const ssize_t y, const size_t columns, const size_t rows, ExceptionInfo *exception)
Definition: cache-view.c:651
#define MaximumNumberOfImageMoments
Definition: statistic.h:25
MagickExport MagickBooleanType SetImageColorMetric(Image *image, const Image *reconstruct_image, ExceptionInfo *exception)
Definition: compare.c:2009
#define MagickEpsilon
Definition: magick-type.h:114
MagickExport MagickBooleanType CompositeImage(Image *image, const Image *composite, const CompositeOperator compose, const MagickBooleanType clip_to_self, const ssize_t x_offset, const ssize_t y_offset, ExceptionInfo *exception)
Definition: composite.c:528
size_t width
Definition: geometry.h:130
static double GetFuzzyColorDistance(const Image *p, const Image *q)
Definition: color-private.h:69
#define ThrowBinaryException(severity, tag, context)
Definition: log.h:52
ssize_t MagickOffsetType
Definition: magick-type.h:133
static Quantum ClampToQuantum(const MagickRealType quantum)
Definition: quantum.h:85
Definition: image.h:151
MagickExport Image * CompareImages(Image *image, const Image *reconstruct_image, const MetricType metric, double *distortion, ExceptionInfo *exception)
Definition: compare.c:128
MagickExport Image * CropImage(const Image *image, const RectangleInfo *geometry, ExceptionInfo *exception)
Definition: transform.c:537
MagickExport KernelInfo * AcquireKernelInfo(const char *kernel_string, ExceptionInfo *exception)
Definition: morphology.c:486
#define MagickCoreSignature
double normalized_mean_error
Definition: color.h:63
MagickExport Quantum * GetCacheViewAuthenticPixels(CacheView *cache_view, const ssize_t x, const ssize_t y, const size_t columns, const size_t rows, ExceptionInfo *exception)
Definition: cache-view.c:299
MagickExport MagickBooleanType SetImageMask(Image *image, const PixelMask type, const Image *mask, ExceptionInfo *exception)
Definition: image.c:3172
MagickExport MagickBooleanType SetImageAlphaChannel(Image *image, const AlphaChannelOption alpha_type, ExceptionInfo *exception)
Definition: channel.c:974
MagickBooleanType
Definition: magick-type.h:169
unsigned int MagickStatusType
Definition: magick-type.h:125
static double PerceptibleReciprocal(const double x)
MagickExport void * AcquireQuantumMemory(const size_t count, const size_t quantum)
Definition: memory.c:634
#define MagickPathExtent
static MagickBooleanType GetMeanErrorPerPixel(Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c:718
MagickExport MagickBooleanType IsStringTrue(const char *value)
Definition: string.c:1415
#define INFINITY
Definition: magick-type.h:195
MagickExport int GetMagickPrecision(void)
Definition: magick.c:942
#define MagickMaximumValue
Definition: magick-type.h:115
MagickExport Quantum * QueueCacheViewAuthenticPixels(CacheView *cache_view, const ssize_t x, const ssize_t y, const size_t columns, const size_t rows, ExceptionInfo *exception)
Definition: cache-view.c:977
MagickExport MagickBooleanType LogMagickEvent(const LogEventType type, const char *module, const char *function, const size_t line, const char *format,...)
Definition: log.c:1660
size_t width
Definition: morphology.h:108
size_t signature
Definition: image.h:354
#define QuantumScale
Definition: magick-type.h:119
size_t columns
Definition: image.h:172
ssize_t x
Definition: geometry.h:134
MagickBooleanType(* MagickProgressMonitor)(const char *, const MagickOffsetType, const MagickSizeType, void *)
Definition: monitor.h:26
size_t height
Definition: geometry.h:130
MagickExport MagickBooleanType QueryColorCompliance(const char *name, const ComplianceType compliance, PixelInfo *color, ExceptionInfo *exception)
Definition: color.c:2180
MagickExport MagickBooleanType SetImageStorageClass(Image *image, const ClassType storage_class, ExceptionInfo *exception)
Definition: image.c:2595
static MagickBooleanType GetPeakAbsoluteDistortion(const Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c:1129
static MagickBooleanType GetMeanSquaredDistortion(const Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c:820
PixelChannel
Definition: pixel.h:70
#define MagickMax(x, y)
Definition: image-private.h:36
static MagickBooleanType GetRootMeanSquaredDistortion(const Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c:1341
static size_t GetPixelChannels(const Image *magick_restrict image)
MagickExport MagickBooleanType IsImagesEqual(const Image *image, const Image *reconstruct_image, ExceptionInfo *exception)
Definition: compare.c:1889
#define SSIMRadius
char filename[MagickPathExtent]
Definition: image.h:319
#define GetMagickModule()
Definition: log.h:28
static double GetSimilarityMetric(const Image *image, const Image *reference, const MetricType metric, const ssize_t x_offset, const ssize_t y_offset, ExceptionInfo *exception)
Definition: compare.c:2140
static PixelChannel GetPixelChannelChannel(const Image *magick_restrict image, const ssize_t offset)
MagickExport CacheView * AcquireVirtualCacheView(const Image *image, ExceptionInfo *exception)
Definition: cache-view.c:149
double normalized_maximum_error
Definition: color.h:63
ErrorInfo error
Definition: image.h:297
MagickExport MagickBooleanType GetImageDistortion(Image *image, const Image *reconstruct_image, const MetricType metric, double *distortion, ExceptionInfo *exception)
Definition: compare.c:1597
unsigned short Quantum
Definition: magick-type.h:86
#define SSIML
#define Log10Epsilon
MagickExport Image * ExtentImage(const Image *image, const RectangleInfo *geometry, ExceptionInfo *exception)
Definition: transform.c:1128
#define SSIMK1
static void SetPixelChannel(const Image *magick_restrict image, const PixelChannel channel, const Quantum quantum, Quantum *magick_restrict pixel)
static MagickBooleanType GetPerceptualHashDistortion(const Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c:1262
#define MagickMin(x, y)
Definition: image-private.h:37
MagickExport void SetGeometry(const Image *image, RectangleInfo *geometry)
Definition: geometry.c:1689
MagickExport void * RelinquishMagickMemory(void *memory)
Definition: memory.c:1122
#define MaxPixelChannels
Definition: pixel.h:27
CompositeOperator compose
Definition: image.h:234
MagickExport double * GetImageDistortions(Image *image, const Image *reconstruct_image, const MetricType metric, ExceptionInfo *exception)
Definition: compare.c:1745
#define MagickExport
MagickExport MagickBooleanType SyncCacheViewAuthenticPixels(CacheView *magick_restrict cache_view, ExceptionInfo *exception)
Definition: cache-view.c:1100
ssize_t y
Definition: geometry.h:134
MagickExport MagickBooleanType FormatImageProperty(Image *image, const char *property, const char *format,...)
Definition: property.c:354
MagickExport CacheView * AcquireAuthenticCacheView(const Image *image, ExceptionInfo *exception)
Definition: cache-view.c:112
#define SSIMSigma
static double MagickLog10(const double x)
Definition: compare.c:1235
static MagickBooleanType GetFuzzDistortion(const Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c:489
PixelTrait
Definition: pixel.h:137
MagickExport ChannelPerceptualHash * GetImagePerceptualHash(const Image *image, ExceptionInfo *exception)
Definition: statistic.c:1724
MagickExport Image * DestroyImage(Image *image)
Definition: image.c:1160
MagickExport Image * CloneImage(const Image *image, const size_t columns, const size_t rows, const MagickBooleanType detach, ExceptionInfo *exception)
Definition: image.c:775
#define QuantumRange
Definition: magick-type.h:87
MagickExport MagickBooleanType SetImageProgress(const Image *image, const char *tag, const MagickOffsetType offset, const MagickSizeType extent)
Definition: monitor.c:136
MagickRealType * values
Definition: morphology.h:116
MagickBooleanType debug
Definition: image.h:334
static MagickBooleanType GetStructuralSimilarityDistortion(const Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c:1356
#define SSIMK2