MagickCore  7.0.11
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-2021 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  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  const Quantum
238  *magick_restrict p,
239  *magick_restrict q;
240 
241  Quantum
242  *magick_restrict r;
243 
244  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  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  const Quantum
403  *magick_restrict p,
404  *magick_restrict q;
405 
406  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  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  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  const Quantum
522  *magick_restrict p,
523  *magick_restrict q;
524 
525  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  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  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  const Quantum
637  *magick_restrict p,
638  *magick_restrict q;
639 
640  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  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]-(double)
686  GetPixelChannel(reconstruct_image,channel,q)));
687  else
688  distance=QuantumScale*fabs((double) (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  const Quantum
745  *magick_restrict p,
746  *magick_restrict q;
747 
748  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  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]-(double)
791  GetPixelChannel(reconstruct_image,channel,q)));
792  else
793  distance=fabs((double) (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  area=PerceptibleReciprocal(area);
809  image->error.mean_error_per_pixel=area*distortion[CompositePixelChannel];
810  image->error.normalized_mean_error=area*QuantumScale*QuantumScale*mean_error;
811  image->error.normalized_maximum_error=QuantumScale*maximum_error;
812  return(status);
813 }
814 
816  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
817 {
818  CacheView
819  *image_view,
820  *reconstruct_view;
821 
822  double
823  area;
824 
826  status;
827 
828  ssize_t
829  j;
830 
831  size_t
832  columns,
833  rows;
834 
835  ssize_t
836  y;
837 
838  status=MagickTrue;
839  rows=MagickMax(image->rows,reconstruct_image->rows);
840  columns=MagickMax(image->columns,reconstruct_image->columns);
841  area=0.0;
842  image_view=AcquireVirtualCacheView(image,exception);
843  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
844 #if defined(MAGICKCORE_OPENMP_SUPPORT)
845  #pragma omp parallel for schedule(static) shared(status) \
846  magick_number_threads(image,image,rows,1) reduction(+:area)
847 #endif
848  for (y=0; y < (ssize_t) rows; y++)
849  {
850  double
851  channel_distortion[MaxPixelChannels+1];
852 
853  const Quantum
854  *magick_restrict p,
855  *magick_restrict q;
856 
857  ssize_t
858  x;
859 
860  if (status == MagickFalse)
861  continue;
862  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
863  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
864  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
865  {
866  status=MagickFalse;
867  continue;
868  }
869  (void) memset(channel_distortion,0,sizeof(channel_distortion));
870  for (x=0; x < (ssize_t) columns; x++)
871  {
872  double
873  Da,
874  Sa;
875 
876  ssize_t
877  i;
878 
879  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
880  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
881  {
882  p+=GetPixelChannels(image);
883  q+=GetPixelChannels(reconstruct_image);
884  continue;
885  }
886  Sa=QuantumScale*GetPixelAlpha(image,p);
887  Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
888  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
889  {
890  double
891  distance;
892 
893  PixelChannel channel = GetPixelChannelChannel(image,i);
894  PixelTrait traits = GetPixelChannelTraits(image,channel);
895  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
896  channel);
897  if ((traits == UndefinedPixelTrait) ||
898  (reconstruct_traits == UndefinedPixelTrait) ||
899  ((reconstruct_traits & UpdatePixelTrait) == 0))
900  continue;
901  if (channel == AlphaPixelChannel)
902  distance=QuantumScale*(p[i]-GetPixelChannel(reconstruct_image,
903  channel,q));
904  else
905  distance=QuantumScale*(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
906  channel,q));
907  channel_distortion[i]+=distance*distance;
908  channel_distortion[CompositePixelChannel]+=distance*distance;
909  }
910  area++;
911  p+=GetPixelChannels(image);
912  q+=GetPixelChannels(reconstruct_image);
913  }
914 #if defined(MAGICKCORE_OPENMP_SUPPORT)
915  #pragma omp critical (MagickCore_GetMeanSquaredError)
916 #endif
917  for (j=0; j <= MaxPixelChannels; j++)
918  distortion[j]+=channel_distortion[j];
919  }
920  reconstruct_view=DestroyCacheView(reconstruct_view);
921  image_view=DestroyCacheView(image_view);
922  area=PerceptibleReciprocal(area);
923  for (j=0; j <= MaxPixelChannels; j++)
924  distortion[j]*=area;
925  distortion[CompositePixelChannel]/=GetImageChannels(image);
926  return(status);
927 }
928 
930  const Image *image,const Image *reconstruct_image,double *distortion,
931  ExceptionInfo *exception)
932 {
933 #define SimilarityImageTag "Similarity/Image"
934 
935  CacheView
936  *image_view,
937  *reconstruct_view;
938 
940  *image_statistics,
941  *reconstruct_statistics;
942 
943  double
944  area;
945 
947  status;
948 
950  progress;
951 
952  ssize_t
953  i;
954 
955  size_t
956  columns,
957  rows;
958 
959  ssize_t
960  y;
961 
962  /*
963  Normalize to account for variation due to lighting and exposure condition.
964  */
965  image_statistics=GetImageStatistics(image,exception);
966  reconstruct_statistics=GetImageStatistics(reconstruct_image,exception);
967  if ((image_statistics == (ChannelStatistics *) NULL) ||
968  (reconstruct_statistics == (ChannelStatistics *) NULL))
969  {
970  if (image_statistics != (ChannelStatistics *) NULL)
971  image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
972  image_statistics);
973  if (reconstruct_statistics != (ChannelStatistics *) NULL)
974  reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
975  reconstruct_statistics);
976  return(MagickFalse);
977  }
978  status=MagickTrue;
979  progress=0;
980  for (i=0; i <= MaxPixelChannels; i++)
981  distortion[i]=0.0;
982  rows=MagickMax(image->rows,reconstruct_image->rows);
983  columns=MagickMax(image->columns,reconstruct_image->columns);
984  area=0.0;
985  image_view=AcquireVirtualCacheView(image,exception);
986  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
987  for (y=0; y < (ssize_t) rows; y++)
988  {
989  const Quantum
990  *magick_restrict p,
991  *magick_restrict q;
992 
993  ssize_t
994  x;
995 
996  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
997  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
998  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
999  {
1000  status=MagickFalse;
1001  break;
1002  }
1003  for (x=0; x < (ssize_t) columns; x++)
1004  {
1005  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1006  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1007  {
1008  p+=GetPixelChannels(image);
1009  q+=GetPixelChannels(reconstruct_image);
1010  continue;
1011  }
1012  area++;
1013  p+=GetPixelChannels(image);
1014  q+=GetPixelChannels(reconstruct_image);
1015  }
1016  }
1017  area=PerceptibleReciprocal(area);
1018  for (y=0; y < (ssize_t) rows; y++)
1019  {
1020  const Quantum
1021  *magick_restrict p,
1022  *magick_restrict q;
1023 
1024  ssize_t
1025  x;
1026 
1027  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1028  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1029  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1030  {
1031  status=MagickFalse;
1032  break;
1033  }
1034  for (x=0; x < (ssize_t) columns; x++)
1035  {
1036  double
1037  Da,
1038  Sa;
1039 
1040  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1041  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1042  {
1043  p+=GetPixelChannels(image);
1044  q+=GetPixelChannels(reconstruct_image);
1045  continue;
1046  }
1047  Sa=QuantumScale*GetPixelAlpha(image,p);
1048  Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
1049  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1050  {
1051  PixelChannel channel = GetPixelChannelChannel(image,i);
1052  PixelTrait traits = GetPixelChannelTraits(image,channel);
1053  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1054  channel);
1055  if ((traits == UndefinedPixelTrait) ||
1056  (reconstruct_traits == UndefinedPixelTrait) ||
1057  ((reconstruct_traits & UpdatePixelTrait) == 0))
1058  continue;
1059  if (channel == AlphaPixelChannel)
1060  {
1061  distortion[i]+=area*QuantumScale*(p[i]-
1062  image_statistics[channel].mean)*(GetPixelChannel(
1063  reconstruct_image,channel,q)-
1064  reconstruct_statistics[channel].mean);
1065  }
1066  else
1067  {
1068  distortion[i]+=area*QuantumScale*(Sa*p[i]-
1069  image_statistics[channel].mean)*(Da*GetPixelChannel(
1070  reconstruct_image,channel,q)-
1071  reconstruct_statistics[channel].mean);
1072  }
1073  }
1074  p+=GetPixelChannels(image);
1075  q+=GetPixelChannels(reconstruct_image);
1076  }
1077  if (image->progress_monitor != (MagickProgressMonitor) NULL)
1078  {
1080  proceed;
1081 
1082 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1083  #pragma omp atomic
1084 #endif
1085  progress++;
1086  proceed=SetImageProgress(image,SimilarityImageTag,progress,rows);
1087  if (proceed == MagickFalse)
1088  {
1089  status=MagickFalse;
1090  break;
1091  }
1092  }
1093  }
1094  reconstruct_view=DestroyCacheView(reconstruct_view);
1095  image_view=DestroyCacheView(image_view);
1096  /*
1097  Divide by the standard deviation.
1098  */
1099  distortion[CompositePixelChannel]=0.0;
1100  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1101  {
1102  double
1103  gamma;
1104 
1105  PixelChannel channel = GetPixelChannelChannel(image,i);
1106  gamma=image_statistics[channel].standard_deviation*
1107  reconstruct_statistics[channel].standard_deviation;
1108  gamma=PerceptibleReciprocal(gamma);
1109  distortion[i]=QuantumRange*gamma*distortion[i];
1110  distortion[CompositePixelChannel]+=distortion[i]*distortion[i];
1111  }
1112  distortion[CompositePixelChannel]=sqrt(distortion[CompositePixelChannel]/
1113  GetImageChannels(image));
1114  /*
1115  Free resources.
1116  */
1117  reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1118  reconstruct_statistics);
1119  image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1120  image_statistics);
1121  return(status);
1122 }
1123 
1125  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1126 {
1127  CacheView
1128  *image_view,
1129  *reconstruct_view;
1130 
1132  status;
1133 
1134  size_t
1135  columns,
1136  rows;
1137 
1138  ssize_t
1139  y;
1140 
1141  status=MagickTrue;
1142  rows=MagickMax(image->rows,reconstruct_image->rows);
1143  columns=MagickMax(image->columns,reconstruct_image->columns);
1144  image_view=AcquireVirtualCacheView(image,exception);
1145  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1146 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1147  #pragma omp parallel for schedule(static) shared(status) \
1148  magick_number_threads(image,image,rows,1)
1149 #endif
1150  for (y=0; y < (ssize_t) rows; y++)
1151  {
1152  double
1153  channel_distortion[MaxPixelChannels+1];
1154 
1155  const Quantum
1156  *magick_restrict p,
1157  *magick_restrict q;
1158 
1159  ssize_t
1160  j,
1161  x;
1162 
1163  if (status == MagickFalse)
1164  continue;
1165  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1166  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1167  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1168  {
1169  status=MagickFalse;
1170  continue;
1171  }
1172  (void) memset(channel_distortion,0,sizeof(channel_distortion));
1173  for (x=0; x < (ssize_t) columns; x++)
1174  {
1175  double
1176  Da,
1177  Sa;
1178 
1179  ssize_t
1180  i;
1181 
1182  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1183  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1184  {
1185  p+=GetPixelChannels(image);
1186  q+=GetPixelChannels(reconstruct_image);
1187  continue;
1188  }
1189  Sa=QuantumScale*GetPixelAlpha(image,p);
1190  Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
1191  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1192  {
1193  double
1194  distance;
1195 
1196  PixelChannel channel = GetPixelChannelChannel(image,i);
1197  PixelTrait traits = GetPixelChannelTraits(image,channel);
1198  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1199  channel);
1200  if ((traits == UndefinedPixelTrait) ||
1201  (reconstruct_traits == UndefinedPixelTrait) ||
1202  ((reconstruct_traits & UpdatePixelTrait) == 0))
1203  continue;
1204  if (channel == AlphaPixelChannel)
1205  distance=QuantumScale*fabs((double) (p[i]-(double)
1206  GetPixelChannel(reconstruct_image,channel,q)));
1207  else
1208  distance=QuantumScale*fabs((double) (Sa*p[i]-Da*
1209  GetPixelChannel(reconstruct_image,channel,q)));
1210  if (distance > channel_distortion[i])
1211  channel_distortion[i]=distance;
1212  if (distance > channel_distortion[CompositePixelChannel])
1213  channel_distortion[CompositePixelChannel]=distance;
1214  }
1215  p+=GetPixelChannels(image);
1216  q+=GetPixelChannels(reconstruct_image);
1217  }
1218 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1219  #pragma omp critical (MagickCore_GetPeakAbsoluteError)
1220 #endif
1221  for (j=0; j <= MaxPixelChannels; j++)
1222  if (channel_distortion[j] > distortion[j])
1223  distortion[j]=channel_distortion[j];
1224  }
1225  reconstruct_view=DestroyCacheView(reconstruct_view);
1226  image_view=DestroyCacheView(image_view);
1227  return(status);
1228 }
1229 
1230 static inline double MagickLog10(const double x)
1231 {
1232 #define Log10Epsilon (1.0e-11)
1233 
1234  if (fabs(x) < Log10Epsilon)
1235  return(log10(Log10Epsilon));
1236  return(log10(fabs(x)));
1237 }
1238 
1240  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1241 {
1243  status;
1244 
1245  ssize_t
1246  i;
1247 
1248  status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1249  for (i=0; i <= MaxPixelChannels; i++)
1250  if (fabs(distortion[i]) < MagickEpsilon)
1251  distortion[i]=INFINITY;
1252  else
1253  distortion[i]=10.0*MagickLog10(1.0)-10.0*MagickLog10(distortion[i]);
1254  return(status);
1255 }
1256 
1258  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1259 {
1261  *channel_phash,
1262  *reconstruct_phash;
1263 
1264  const char
1265  *artifact;
1266 
1268  normalize;
1269 
1270  ssize_t
1271  channel;
1272 
1273  /*
1274  Compute perceptual hash in the sRGB colorspace.
1275  */
1276  channel_phash=GetImagePerceptualHash(image,exception);
1277  if (channel_phash == (ChannelPerceptualHash *) NULL)
1278  return(MagickFalse);
1279  reconstruct_phash=GetImagePerceptualHash(reconstruct_image,exception);
1280  if (reconstruct_phash == (ChannelPerceptualHash *) NULL)
1281  {
1283  channel_phash);
1284  return(MagickFalse);
1285  }
1286  artifact=GetImageArtifact(image,"phash:normalize");
1287  normalize=(artifact == (const char *) NULL) ||
1288  (IsStringTrue(artifact) == MagickFalse) ? MagickFalse : MagickTrue;
1289 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1290  #pragma omp parallel for schedule(static)
1291 #endif
1292  for (channel=0; channel < MaxPixelChannels; channel++)
1293  {
1294  double
1295  difference;
1296 
1297  ssize_t
1298  i;
1299 
1300  difference=0.0;
1301  for (i=0; i < MaximumNumberOfImageMoments; i++)
1302  {
1303  double
1304  alpha,
1305  beta;
1306 
1307  ssize_t
1308  j;
1309 
1310  for (j=0; j < (ssize_t) channel_phash[0].number_colorspaces; j++)
1311  {
1312  alpha=channel_phash[channel].phash[j][i];
1313  beta=reconstruct_phash[channel].phash[j][i];
1314  if (normalize == MagickFalse)
1315  difference+=(beta-alpha)*(beta-alpha);
1316  else
1317  difference=sqrt((beta-alpha)*(beta-alpha)/
1318  channel_phash[0].number_channels);
1319  }
1320  }
1321  distortion[channel]+=difference;
1322 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1323  #pragma omp critical (MagickCore_GetPerceptualHashDistortion)
1324 #endif
1325  distortion[CompositePixelChannel]+=difference;
1326  }
1327  /*
1328  Free resources.
1329  */
1330  reconstruct_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1331  reconstruct_phash);
1332  channel_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(channel_phash);
1333  return(MagickTrue);
1334 }
1335 
1337  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1338 {
1340  status;
1341 
1342  ssize_t
1343  i;
1344 
1345  status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1346  for (i=0; i <= MaxPixelChannels; i++)
1347  distortion[i]=sqrt(distortion[i]);
1348  return(status);
1349 }
1350 
1352  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1353 {
1354 #define SSIMRadius 5.0
1355 #define SSIMSigma 1.5
1356 #define SSIMBlocksize 8
1357 #define SSIMK1 0.01
1358 #define SSIMK2 0.03
1359 #define SSIML 1.0
1360 
1361  CacheView
1362  *image_view,
1363  *reconstruct_view;
1364 
1365  char
1366  geometry[MagickPathExtent];
1367 
1368  const char
1369  *artifact;
1370 
1371  double
1372  area,
1373  c1,
1374  c2,
1375  radius,
1376  sigma;
1377 
1378  KernelInfo
1379  *kernel_info;
1380 
1382  status;
1383 
1384  ssize_t
1385  i;
1386 
1387  size_t
1388  columns,
1389  rows;
1390 
1391  ssize_t
1392  y;
1393 
1394  /*
1395  Compute structural similarity index @
1396  https://en.wikipedia.org/wiki/Structural_similarity.
1397  */
1398  radius=SSIMRadius;
1399  artifact=GetImageArtifact(image,"compare:ssim-radius");
1400  if (artifact != (const char *) NULL)
1401  radius=StringToDouble(artifact,(char **) NULL);
1402  sigma=SSIMSigma;
1403  artifact=GetImageArtifact(image,"compare:ssim-sigma");
1404  if (artifact != (const char *) NULL)
1405  sigma=StringToDouble(artifact,(char **) NULL);
1406  (void) FormatLocaleString(geometry,MagickPathExtent,"gaussian:%.20gx%.20g",
1407  radius,sigma);
1408  kernel_info=AcquireKernelInfo(geometry,exception);
1409  if (kernel_info == (KernelInfo *) NULL)
1410  ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
1411  image->filename);
1412  c1=pow(SSIMK1*SSIML,2.0);
1413  artifact=GetImageArtifact(image,"compare:ssim-k1");
1414  if (artifact != (const char *) NULL)
1415  c1=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
1416  c2=pow(SSIMK2*SSIML,2.0);
1417  artifact=GetImageArtifact(image,"compare:ssim-k2");
1418  if (artifact != (const char *) NULL)
1419  c2=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
1420  status=MagickTrue;
1421  area=0.0;
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  const Quantum
1436  *magick_restrict p,
1437  *magick_restrict q;
1438 
1439  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  const Quantum
1467  *magick_restrict reference,
1468  *magick_restrict target;
1469 
1471  *k;
1472 
1473  ssize_t
1474  v;
1475 
1476  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1477  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1478  {
1479  p+=GetPixelChannels(image);
1480  q+=GetPixelChannels(reconstruct_image);
1481  continue;
1482  }
1483  (void) memset(x_pixel_mu,0,sizeof(x_pixel_mu));
1484  (void) memset(x_pixel_sigma_squared,0,sizeof(x_pixel_sigma_squared));
1485  (void) memset(xy_sigma,0,sizeof(xy_sigma));
1486  (void) memset(x_pixel_sigma_squared,0,sizeof(y_pixel_sigma_squared));
1487  (void) memset(y_pixel_mu,0,sizeof(y_pixel_mu));
1488  (void) memset(y_pixel_sigma_squared,0,sizeof(y_pixel_sigma_squared));
1489  k=kernel_info->values;
1490  reference=p;
1491  target=q;
1492  for (v=0; v < (ssize_t) kernel_info->height; v++)
1493  {
1494  ssize_t
1495  u;
1496 
1497  for (u=0; u < (ssize_t) kernel_info->width; u++)
1498  {
1499  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1500  {
1501  double
1502  x_pixel,
1503  y_pixel;
1504 
1505  PixelChannel channel = GetPixelChannelChannel(image,i);
1506  PixelTrait traits = GetPixelChannelTraits(image,channel);
1507  PixelTrait reconstruct_traits = GetPixelChannelTraits(
1508  reconstruct_image,channel);
1509  if ((traits == UndefinedPixelTrait) ||
1510  (reconstruct_traits == UndefinedPixelTrait) ||
1511  ((reconstruct_traits & UpdatePixelTrait) == 0))
1512  continue;
1513  x_pixel=QuantumScale*reference[i];
1514  x_pixel_mu[i]+=(*k)*x_pixel;
1515  x_pixel_sigma_squared[i]+=(*k)*x_pixel*x_pixel;
1516  y_pixel=QuantumScale*
1517  GetPixelChannel(reconstruct_image,channel,target);
1518  y_pixel_mu[i]+=(*k)*y_pixel;
1519  y_pixel_sigma_squared[i]+=(*k)*y_pixel*y_pixel;
1520  xy_sigma[i]+=(*k)*x_pixel*y_pixel;
1521  }
1522  k++;
1523  reference+=GetPixelChannels(image);
1524  target+=GetPixelChannels(reconstruct_image);
1525  }
1526  reference+=GetPixelChannels(image)*columns;
1527  target+=GetPixelChannels(reconstruct_image)*columns;
1528  }
1529  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1530  {
1531  double
1532  ssim,
1533  x_pixel_mu_squared,
1534  x_pixel_sigmas_squared,
1535  xy_mu,
1536  xy_sigmas,
1537  y_pixel_mu_squared,
1538  y_pixel_sigmas_squared;
1539 
1540  PixelChannel channel = GetPixelChannelChannel(image,i);
1541  PixelTrait traits = GetPixelChannelTraits(image,channel);
1542  PixelTrait reconstruct_traits = GetPixelChannelTraits(
1543  reconstruct_image,channel);
1544  if ((traits == UndefinedPixelTrait) ||
1545  (reconstruct_traits == UndefinedPixelTrait) ||
1546  ((reconstruct_traits & UpdatePixelTrait) == 0))
1547  continue;
1548  x_pixel_mu_squared=x_pixel_mu[i]*x_pixel_mu[i];
1549  y_pixel_mu_squared=y_pixel_mu[i]*y_pixel_mu[i];
1550  xy_mu=x_pixel_mu[i]*y_pixel_mu[i];
1551  xy_sigmas=xy_sigma[i]-xy_mu;
1552  x_pixel_sigmas_squared=x_pixel_sigma_squared[i]-x_pixel_mu_squared;
1553  y_pixel_sigmas_squared=y_pixel_sigma_squared[i]-y_pixel_mu_squared;
1554  ssim=((2.0*xy_mu+c1)*(2.0*xy_sigmas+c2))/
1555  ((x_pixel_mu_squared+y_pixel_mu_squared+c1)*
1556  (x_pixel_sigmas_squared+y_pixel_sigmas_squared+c2));
1557  channel_distortion[i]+=ssim;
1558  channel_distortion[CompositePixelChannel]+=ssim;
1559  }
1560  area++;
1561  p+=GetPixelChannels(image);
1562  q+=GetPixelChannels(reconstruct_image);
1563  }
1564 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1565  #pragma omp critical (MagickCore_GetStructuralSimilarityDistortion)
1566 #endif
1567  for (i=0; i <= MaxPixelChannels; i++)
1568  distortion[i]+=channel_distortion[i];
1569  }
1570  image_view=DestroyCacheView(image_view);
1571  reconstruct_view=DestroyCacheView(reconstruct_view);
1572  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1573  {
1574  PixelChannel channel = GetPixelChannelChannel(image,i);
1575  PixelTrait traits = GetPixelChannelTraits(image,channel);
1576  if ((traits == UndefinedPixelTrait) || ((traits & UpdatePixelTrait) == 0))
1577  continue;
1578  distortion[i]/=area;
1579  }
1580  distortion[CompositePixelChannel]/=area;
1581  distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
1582  kernel_info=DestroyKernelInfo(kernel_info);
1583  return(status);
1584 }
1585 
1587  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1588 {
1590  status;
1591 
1592  ssize_t
1593  i;
1594 
1595  status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1596  distortion,exception);
1597  for (i=0; i <= MaxPixelChannels; i++)
1598  distortion[i]=(1.0-(distortion[i]))/2.0;
1599  return(status);
1600 }
1601 
1603  const Image *reconstruct_image,const MetricType metric,double *distortion,
1604  ExceptionInfo *exception)
1605 {
1606  double
1607  *channel_distortion;
1608 
1610  status;
1611 
1612  size_t
1613  length;
1614 
1615  assert(image != (Image *) NULL);
1616  assert(image->signature == MagickCoreSignature);
1617  if (image->debug != MagickFalse)
1618  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1619  assert(reconstruct_image != (const Image *) NULL);
1620  assert(reconstruct_image->signature == MagickCoreSignature);
1621  assert(distortion != (double *) NULL);
1622  *distortion=0.0;
1623  if (image->debug != MagickFalse)
1624  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1625  /*
1626  Get image distortion.
1627  */
1628  length=MaxPixelChannels+1UL;
1629  channel_distortion=(double *) AcquireQuantumMemory(length,
1630  sizeof(*channel_distortion));
1631  if (channel_distortion == (double *) NULL)
1632  ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1633  (void) memset(channel_distortion,0,length*
1634  sizeof(*channel_distortion));
1635  switch (metric)
1636  {
1637  case AbsoluteErrorMetric:
1638  {
1639  status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1640  exception);
1641  break;
1642  }
1643  case FuzzErrorMetric:
1644  {
1645  status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1646  exception);
1647  break;
1648  }
1650  {
1651  status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1652  channel_distortion,exception);
1653  break;
1654  }
1656  {
1657  status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1658  exception);
1659  break;
1660  }
1662  {
1663  status=GetMeanSquaredDistortion(image,reconstruct_image,
1664  channel_distortion,exception);
1665  break;
1666  }
1668  default:
1669  {
1670  status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1671  channel_distortion,exception);
1672  break;
1673  }
1675  {
1676  status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1677  channel_distortion,exception);
1678  break;
1679  }
1681  {
1682  status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1683  channel_distortion,exception);
1684  break;
1685  }
1687  {
1688  status=GetPerceptualHashDistortion(image,reconstruct_image,
1689  channel_distortion,exception);
1690  break;
1691  }
1693  {
1694  status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1695  channel_distortion,exception);
1696  break;
1697  }
1699  {
1700  status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1701  channel_distortion,exception);
1702  break;
1703  }
1705  {
1706  status=GetStructuralDisimilarityDistortion(image,reconstruct_image,
1707  channel_distortion,exception);
1708  break;
1709  }
1710  }
1711  *distortion=channel_distortion[CompositePixelChannel];
1712  channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1713  (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
1714  *distortion);
1715  return(status);
1716 }
1717 
1718 /*
1719 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1720 % %
1721 % %
1722 % %
1723 % G e t I m a g e D i s t o r t i o n s %
1724 % %
1725 % %
1726 % %
1727 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1728 %
1729 % GetImageDistortions() compares the pixel channels of an image to a
1730 % reconstructed image and returns the specified distortion metric for each
1731 % channel.
1732 %
1733 % The format of the GetImageDistortions method is:
1734 %
1735 % double *GetImageDistortions(const Image *image,
1736 % const Image *reconstruct_image,const MetricType metric,
1737 % ExceptionInfo *exception)
1738 %
1739 % A description of each parameter follows:
1740 %
1741 % o image: the image.
1742 %
1743 % o reconstruct_image: the reconstruct image.
1744 %
1745 % o metric: the metric.
1746 %
1747 % o exception: return any errors or warnings in this structure.
1748 %
1749 */
1751  const Image *reconstruct_image,const MetricType metric,
1752  ExceptionInfo *exception)
1753 {
1754  double
1755  *channel_distortion;
1756 
1758  status;
1759 
1760  size_t
1761  length;
1762 
1763  assert(image != (Image *) NULL);
1764  assert(image->signature == MagickCoreSignature);
1765  if (image->debug != MagickFalse)
1766  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1767  assert(reconstruct_image != (const Image *) NULL);
1768  assert(reconstruct_image->signature == MagickCoreSignature);
1769  if (image->debug != MagickFalse)
1770  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1771  /*
1772  Get image distortion.
1773  */
1774  length=MaxPixelChannels+1UL;
1775  channel_distortion=(double *) AcquireQuantumMemory(length,
1776  sizeof(*channel_distortion));
1777  if (channel_distortion == (double *) NULL)
1778  ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1779  (void) memset(channel_distortion,0,length*
1780  sizeof(*channel_distortion));
1781  status=MagickTrue;
1782  switch (metric)
1783  {
1784  case AbsoluteErrorMetric:
1785  {
1786  status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1787  exception);
1788  break;
1789  }
1790  case FuzzErrorMetric:
1791  {
1792  status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1793  exception);
1794  break;
1795  }
1797  {
1798  status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1799  channel_distortion,exception);
1800  break;
1801  }
1803  {
1804  status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1805  exception);
1806  break;
1807  }
1809  {
1810  status=GetMeanSquaredDistortion(image,reconstruct_image,
1811  channel_distortion,exception);
1812  break;
1813  }
1815  default:
1816  {
1817  status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1818  channel_distortion,exception);
1819  break;
1820  }
1822  {
1823  status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1824  channel_distortion,exception);
1825  break;
1826  }
1828  {
1829  status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1830  channel_distortion,exception);
1831  break;
1832  }
1834  {
1835  status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1836  channel_distortion,exception);
1837  break;
1838  }
1840  {
1841  status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1842  channel_distortion,exception);
1843  break;
1844  }
1846  {
1847  status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1848  channel_distortion,exception);
1849  break;
1850  }
1852  {
1853  status=GetStructuralDisimilarityDistortion(image,reconstruct_image,
1854  channel_distortion,exception);
1855  break;
1856  }
1857  }
1858  if (status == MagickFalse)
1859  {
1860  channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1861  return((double *) NULL);
1862  }
1863  return(channel_distortion);
1864 }
1865 
1866 /*
1867 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1868 % %
1869 % %
1870 % %
1871 % I s I m a g e s E q u a l %
1872 % %
1873 % %
1874 % %
1875 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1876 %
1877 % IsImagesEqual() compare the pixels of two images and returns immediately
1878 % if any pixel is not identical.
1879 %
1880 % The format of the IsImagesEqual method is:
1881 %
1882 % MagickBooleanType IsImagesEqual(const Image *image,
1883 % const Image *reconstruct_image,ExceptionInfo *exception)
1884 %
1885 % A description of each parameter follows.
1886 %
1887 % o image: the image.
1888 %
1889 % o reconstruct_image: the reconstruct image.
1890 %
1891 % o exception: return any errors or warnings in this structure.
1892 %
1893 */
1895  const Image *reconstruct_image,ExceptionInfo *exception)
1896 {
1897  CacheView
1898  *image_view,
1899  *reconstruct_view;
1900 
1901  size_t
1902  columns,
1903  rows;
1904 
1905  ssize_t
1906  y;
1907 
1908  assert(image != (Image *) NULL);
1909  assert(image->signature == MagickCoreSignature);
1910  assert(reconstruct_image != (const Image *) NULL);
1911  assert(reconstruct_image->signature == MagickCoreSignature);
1912  rows=MagickMax(image->rows,reconstruct_image->rows);
1913  columns=MagickMax(image->columns,reconstruct_image->columns);
1914  image_view=AcquireVirtualCacheView(image,exception);
1915  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1916  for (y=0; y < (ssize_t) rows; y++)
1917  {
1918  const Quantum
1919  *magick_restrict p,
1920  *magick_restrict q;
1921 
1922  ssize_t
1923  x;
1924 
1925  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1926  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1927  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
1928  break;
1929  for (x=0; x < (ssize_t) columns; x++)
1930  {
1931  ssize_t
1932  i;
1933 
1934  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1935  {
1936  double
1937  distance;
1938 
1939  PixelChannel channel = GetPixelChannelChannel(image,i);
1940  PixelTrait traits = GetPixelChannelTraits(image,channel);
1941  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1942  channel);
1943  if ((traits == UndefinedPixelTrait) ||
1944  (reconstruct_traits == UndefinedPixelTrait) ||
1945  ((reconstruct_traits & UpdatePixelTrait) == 0))
1946  continue;
1947  distance=fabs((double) (p[i]-(double) GetPixelChannel(reconstruct_image,
1948  channel,q)));
1949  if (distance >= MagickEpsilon)
1950  break;
1951  }
1952  if (i < (ssize_t) GetPixelChannels(image))
1953  break;
1954  p+=GetPixelChannels(image);
1955  q+=GetPixelChannels(reconstruct_image);
1956  }
1957  if (x < (ssize_t) columns)
1958  break;
1959  }
1960  reconstruct_view=DestroyCacheView(reconstruct_view);
1961  image_view=DestroyCacheView(image_view);
1962  return(y < (ssize_t) rows ? MagickFalse : MagickTrue);
1963 }
1964 
1965 /*
1966 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1967 % %
1968 % %
1969 % %
1970 % S e t I m a g e C o l o r M e t r i c %
1971 % %
1972 % %
1973 % %
1974 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1975 %
1976 % SetImageColorMetric() measures the difference between colors at each pixel
1977 % location of two images. A value other than 0 means the colors match
1978 % exactly. Otherwise an error measure is computed by summing over all
1979 % pixels in an image the distance squared in RGB space between each image
1980 % pixel and its corresponding pixel in the reconstruct image. The error
1981 % measure is assigned to these image members:
1982 %
1983 % o mean_error_per_pixel: The mean error for any single pixel in
1984 % the image.
1985 %
1986 % o normalized_mean_error: The normalized mean quantization error for
1987 % any single pixel in the image. This distance measure is normalized to
1988 % a range between 0 and 1. It is independent of the range of red, green,
1989 % and blue values in the image.
1990 %
1991 % o normalized_maximum_error: The normalized maximum quantization
1992 % error for any single pixel in the image. This distance measure is
1993 % normalized to a range between 0 and 1. It is independent of the range
1994 % of red, green, and blue values in your image.
1995 %
1996 % A small normalized mean square error, accessed as
1997 % image->normalized_mean_error, suggests the images are very similar in
1998 % spatial layout and color.
1999 %
2000 % The format of the SetImageColorMetric method is:
2001 %
2002 % MagickBooleanType SetImageColorMetric(Image *image,
2003 % const Image *reconstruct_image,ExceptionInfo *exception)
2004 %
2005 % A description of each parameter follows.
2006 %
2007 % o image: the image.
2008 %
2009 % o reconstruct_image: the reconstruct image.
2010 %
2011 % o exception: return any errors or warnings in this structure.
2012 %
2013 */
2015  const Image *reconstruct_image,ExceptionInfo *exception)
2016 {
2017  CacheView
2018  *image_view,
2019  *reconstruct_view;
2020 
2021  double
2022  area,
2023  maximum_error,
2024  mean_error,
2025  mean_error_per_pixel;
2026 
2028  status;
2029 
2030  size_t
2031  columns,
2032  rows;
2033 
2034  ssize_t
2035  y;
2036 
2037  assert(image != (Image *) NULL);
2038  assert(image->signature == MagickCoreSignature);
2039  assert(reconstruct_image != (const Image *) NULL);
2040  assert(reconstruct_image->signature == MagickCoreSignature);
2041  area=0.0;
2042  maximum_error=0.0;
2043  mean_error_per_pixel=0.0;
2044  mean_error=0.0;
2045  rows=MagickMax(image->rows,reconstruct_image->rows);
2046  columns=MagickMax(image->columns,reconstruct_image->columns);
2047  image_view=AcquireVirtualCacheView(image,exception);
2048  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
2049  for (y=0; y < (ssize_t) rows; y++)
2050  {
2051  const Quantum
2052  *magick_restrict p,
2053  *magick_restrict q;
2054 
2055  ssize_t
2056  x;
2057 
2058  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
2059  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
2060  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2061  break;
2062  for (x=0; x < (ssize_t) columns; x++)
2063  {
2064  ssize_t
2065  i;
2066 
2067  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2068  {
2069  double
2070  distance;
2071 
2072  PixelChannel channel = GetPixelChannelChannel(image,i);
2073  PixelTrait traits = GetPixelChannelTraits(image,channel);
2074  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2075  channel);
2076  if ((traits == UndefinedPixelTrait) ||
2077  (reconstruct_traits == UndefinedPixelTrait) ||
2078  ((reconstruct_traits & UpdatePixelTrait) == 0))
2079  continue;
2080  distance=fabs((double) (p[i]-(double) GetPixelChannel(reconstruct_image,
2081  channel,q)));
2082  if (distance >= MagickEpsilon)
2083  {
2084  mean_error_per_pixel+=distance;
2085  mean_error+=distance*distance;
2086  if (distance > maximum_error)
2087  maximum_error=distance;
2088  }
2089  area++;
2090  }
2091  p+=GetPixelChannels(image);
2092  q+=GetPixelChannels(reconstruct_image);
2093  }
2094  }
2095  reconstruct_view=DestroyCacheView(reconstruct_view);
2096  image_view=DestroyCacheView(image_view);
2097  image->error.mean_error_per_pixel=(double) (mean_error_per_pixel/area);
2099  mean_error/area);
2100  image->error.normalized_maximum_error=(double) (QuantumScale*maximum_error);
2101  status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
2102  return(status);
2103 }
2104 
2105 /*
2106 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2107 % %
2108 % %
2109 % %
2110 % S i m i l a r i t y I m a g e %
2111 % %
2112 % %
2113 % %
2114 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2115 %
2116 % SimilarityImage() compares the reference image of the image and returns the
2117 % best match offset. In addition, it returns a similarity image such that an
2118 % exact match location is completely white and if none of the pixels match,
2119 % black, otherwise some gray level in-between.
2120 %
2121 % The format of the SimilarityImageImage method is:
2122 %
2123 % Image *SimilarityImage(const Image *image,const Image *reference,
2124 % const MetricType metric,const double similarity_threshold,
2125 % RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
2126 %
2127 % A description of each parameter follows:
2128 %
2129 % o image: the image.
2130 %
2131 % o reference: find an area of the image that closely resembles this image.
2132 %
2133 % o metric: the metric.
2134 %
2135 % o similarity_threshold: minimum distortion for (sub)image match.
2136 %
2137 % o offset: the best match offset of the reference image within the image.
2138 %
2139 % o similarity: the computed similarity between the images.
2140 %
2141 % o exception: return any errors or warnings in this structure.
2142 %
2143 */
2144 
2145 static double GetSimilarityMetric(const Image *image,const Image *reference,
2146  const MetricType metric,const ssize_t x_offset,const ssize_t y_offset,
2147  ExceptionInfo *exception)
2148 {
2149  double
2150  distortion;
2151 
2152  Image
2153  *similarity_image;
2154 
2156  status;
2157 
2159  geometry;
2160 
2161  SetGeometry(reference,&geometry);
2162  geometry.x=x_offset;
2163  geometry.y=y_offset;
2164  similarity_image=CropImage(image,&geometry,exception);
2165  if (similarity_image == (Image *) NULL)
2166  return(0.0);
2167  distortion=0.0;
2168  status=GetImageDistortion(similarity_image,reference,metric,&distortion,
2169  exception);
2170  similarity_image=DestroyImage(similarity_image);
2171  if (status == MagickFalse)
2172  return(0.0);
2173  return(distortion);
2174 }
2175 
2176 MagickExport Image *SimilarityImage(const Image *image,const Image *reference,
2177  const MetricType metric,const double similarity_threshold,
2178  RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
2179 {
2180 #define SimilarityImageTag "Similarity/Image"
2181 
2182  CacheView
2183  *similarity_view;
2184 
2185  Image
2186  *similarity_image;
2187 
2189  status;
2190 
2192  progress;
2193 
2194  ssize_t
2195  y;
2196 
2197  assert(image != (const Image *) NULL);
2198  assert(image->signature == MagickCoreSignature);
2199  if (image->debug != MagickFalse)
2200  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2201  assert(exception != (ExceptionInfo *) NULL);
2202  assert(exception->signature == MagickCoreSignature);
2203  assert(offset != (RectangleInfo *) NULL);
2204  SetGeometry(reference,offset);
2205  *similarity_metric=MagickMaximumValue;
2206  similarity_image=CloneImage(image,image->columns-reference->columns+1,
2207  image->rows-reference->rows+1,MagickTrue,exception);
2208  if (similarity_image == (Image *) NULL)
2209  return((Image *) NULL);
2210  status=SetImageStorageClass(similarity_image,DirectClass,exception);
2211  if (status == MagickFalse)
2212  {
2213  similarity_image=DestroyImage(similarity_image);
2214  return((Image *) NULL);
2215  }
2216  (void) SetImageAlphaChannel(similarity_image,DeactivateAlphaChannel,
2217  exception);
2218  /*
2219  Measure similarity of reference image against image.
2220  */
2221  status=MagickTrue;
2222  progress=0;
2223  similarity_view=AcquireAuthenticCacheView(similarity_image,exception);
2224 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2225  #pragma omp parallel for schedule(static) \
2226  shared(progress,status,similarity_metric) \
2227  magick_number_threads(image,image,image->rows-reference->rows+1,1)
2228 #endif
2229  for (y=0; y < (ssize_t) (image->rows-reference->rows+1); y++)
2230  {
2231  double
2232  similarity;
2233 
2234  Quantum
2235  *magick_restrict q;
2236 
2237  ssize_t
2238  x;
2239 
2240  if (status == MagickFalse)
2241  continue;
2242 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2243  #pragma omp flush(similarity_metric)
2244 #endif
2245  if (*similarity_metric <= similarity_threshold)
2246  continue;
2247  q=GetCacheViewAuthenticPixels(similarity_view,0,y,similarity_image->columns,
2248  1,exception);
2249  if (q == (Quantum *) NULL)
2250  {
2251  status=MagickFalse;
2252  continue;
2253  }
2254  for (x=0; x < (ssize_t) (image->columns-reference->columns+1); x++)
2255  {
2256  ssize_t
2257  i;
2258 
2259 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2260  #pragma omp flush(similarity_metric)
2261 #endif
2262  if (*similarity_metric <= similarity_threshold)
2263  break;
2264  similarity=GetSimilarityMetric(image,reference,metric,x,y,exception);
2265 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2266  #pragma omp critical (MagickCore_SimilarityImage)
2267 #endif
2268  if ((metric == NormalizedCrossCorrelationErrorMetric) ||
2269  (metric == UndefinedErrorMetric))
2270  similarity=1.0-similarity;
2271  if (similarity < *similarity_metric)
2272  {
2273  offset->x=x;
2274  offset->y=y;
2275  *similarity_metric=similarity;
2276  }
2277  if (metric == PerceptualHashErrorMetric)
2278  similarity=MagickMin(0.01*similarity,1.0);
2279  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2280  {
2281  PixelChannel channel = GetPixelChannelChannel(image,i);
2282  PixelTrait traits = GetPixelChannelTraits(image,channel);
2283  PixelTrait similarity_traits=GetPixelChannelTraits(similarity_image,
2284  channel);
2285  if ((traits == UndefinedPixelTrait) ||
2286  (similarity_traits == UndefinedPixelTrait) ||
2287  ((similarity_traits & UpdatePixelTrait) == 0))
2288  continue;
2289  SetPixelChannel(similarity_image,channel,ClampToQuantum(QuantumRange-
2290  QuantumRange*similarity),q);
2291  }
2292  q+=GetPixelChannels(similarity_image);
2293  }
2294  if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
2295  status=MagickFalse;
2296  if (image->progress_monitor != (MagickProgressMonitor) NULL)
2297  {
2299  proceed;
2300 
2301 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2302  #pragma omp atomic
2303 #endif
2304  progress++;
2305  proceed=SetImageProgress(image,SimilarityImageTag,progress,image->rows);
2306  if (proceed == MagickFalse)
2307  status=MagickFalse;
2308  }
2309  }
2310  similarity_view=DestroyCacheView(similarity_view);
2311  if (status == MagickFalse)
2312  similarity_image=DestroyImage(similarity_image);
2313  return(similarity_image);
2314 }
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:929
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:1239
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:2176
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:2063
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:1586
#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:79
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:2014
#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:131
static double GetFuzzyColorDistance(const Image *p, const Image *q)
Definition: color-private.h:69
#define ThrowBinaryException(severity, tag, context)
Definition: log.h:52
ssize_t MagickOffsetType
Definition: magick-type.h:133
static Quantum ClampToQuantum(const MagickRealType quantum)
Definition: quantum.h:85
Definition: image.h:151
MagickExport Image * CompareImages(Image *image, const Image *reconstruct_image, const MetricType metric, double *distortion, ExceptionInfo *exception)
Definition: compare.c:128
MagickExport Image * CropImage(const Image *image, const RectangleInfo *geometry, ExceptionInfo *exception)
Definition: transform.c:537
MagickExport KernelInfo * AcquireKernelInfo(const char *kernel_string, ExceptionInfo *exception)
Definition: morphology.c:486
#define MagickCoreSignature
double normalized_mean_error
Definition: color.h:79
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:665
#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:195
MagickExport int GetMagickPrecision(void)
Definition: magick.c:942
#define MagickMaximumValue
Definition: magick-type.h:115
MagickExport Quantum * QueueCacheViewAuthenticPixels(CacheView *cache_view, const ssize_t x, const ssize_t y, const size_t columns, const size_t rows, ExceptionInfo *exception)
Definition: cache-view.c:977
MagickExport MagickBooleanType LogMagickEvent(const LogEventType type, const char *module, const char *function, const size_t line, const char *format,...)
Definition: log.c:1660
size_t width
Definition: morphology.h:108
size_t signature
Definition: image.h:354
#define QuantumScale
Definition: magick-type.h:119
size_t columns
Definition: image.h:172
ssize_t x
Definition: geometry.h:135
MagickBooleanType(* MagickProgressMonitor)(const char *, const MagickOffsetType, const MagickSizeType, void *)
Definition: monitor.h:26
size_t height
Definition: geometry.h:131
MagickExport MagickBooleanType QueryColorCompliance(const char *name, const ComplianceType compliance, PixelInfo *color, ExceptionInfo *exception)
Definition: color.c:2265
MagickExport MagickBooleanType SetImageStorageClass(Image *image, const ClassType storage_class, ExceptionInfo *exception)
Definition: image.c:2605
static MagickBooleanType GetPeakAbsoluteDistortion(const Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c:1124
static MagickBooleanType GetMeanSquaredDistortion(const Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c:815
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:1336
static size_t GetPixelChannels(const Image *magick_restrict image)
MagickExport MagickBooleanType IsImagesEqual(const Image *image, const Image *reconstruct_image, ExceptionInfo *exception)
Definition: compare.c:1894
#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:2145
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:79
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:1602
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:1119
#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:1257
#define MagickMin(x, y)
Definition: image-private.h:37
MagickExport void SetGeometry(const Image *image, RectangleInfo *geometry)
Definition: geometry.c:1696
MagickExport void * RelinquishMagickMemory(void *memory)
Definition: memory.c:1162
#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:1750
#define MagickExport
MagickExport MagickBooleanType SyncCacheViewAuthenticPixels(CacheView *magick_restrict cache_view, ExceptionInfo *exception)
Definition: cache-view.c:1100
ssize_t y
Definition: geometry.h:135
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:1230
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:1168
MagickExport Image * CloneImage(const Image *image, const size_t columns, const size_t rows, const MagickBooleanType detach, ExceptionInfo *exception)
Definition: image.c:779
#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:1351
#define SSIMK2