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