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