MagickCore  7.0.3
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-2019 ImageMagick Studio LLC, a non-profit organization %
21 % dedicated to making software imaging solutions freely available. %
22 % %
23 % You may not use this file except in compliance with the License. You may %
24 % obtain a copy of the License at %
25 % %
26 % https://imagemagick.org/script/license.php %
27 % %
28 % Unless required by applicable law or agreed to in writing, software %
29 % distributed under the License is distributed on an "AS IS" BASIS, %
30 % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. %
31 % See the License for the specific language governing permissions and %
32 % limitations under the License. %
33 % %
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35 %
36 %
37 %
38 */
39 
40 /*
41  Include declarations.
42 */
43 #include "MagickCore/studio.h"
44 #include "MagickCore/artifact.h"
45 #include "MagickCore/attribute.h"
46 #include "MagickCore/cache-view.h"
47 #include "MagickCore/channel.h"
48 #include "MagickCore/client.h"
49 #include "MagickCore/color.h"
51 #include "MagickCore/colorspace.h"
53 #include "MagickCore/compare.h"
55 #include "MagickCore/constitute.h"
57 #include "MagickCore/geometry.h"
59 #include "MagickCore/list.h"
60 #include "MagickCore/log.h"
61 #include "MagickCore/memory_.h"
62 #include "MagickCore/monitor.h"
64 #include "MagickCore/option.h"
66 #include "MagickCore/property.h"
67 #include "MagickCore/resource_.h"
68 #include "MagickCore/string_.h"
69 #include "MagickCore/statistic.h"
72 #include "MagickCore/transform.h"
73 #include "MagickCore/utility.h"
74 #include "MagickCore/version.h"
75 
76 /*
77 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
78 % %
79 % %
80 % %
81 % C o m p a r e I m a g e %
82 % %
83 % %
84 % %
85 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
86 %
87 % CompareImages() compares one or more pixel channels of an image to a
88 % reconstructed image and returns the difference image.
89 %
90 % The format of the CompareImages method is:
91 %
92 % Image *CompareImages(const Image *image,const Image *reconstruct_image,
93 % const MetricType metric,double *distortion,ExceptionInfo *exception)
94 %
95 % A description of each parameter follows:
96 %
97 % o image: the image.
98 %
99 % o reconstruct_image: the reconstruct image.
100 %
101 % o metric: the metric.
102 %
103 % o distortion: the computed distortion between the images.
104 %
105 % o exception: return any errors or warnings in this structure.
106 %
107 */
108 
109 static size_t GetImageChannels(const Image *image)
110 {
111  register ssize_t
112  i;
113 
114  size_t
115  channels;
116 
117  channels=0;
118  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
119  {
120  PixelChannel channel = GetPixelChannelChannel(image,i);
121  PixelTrait traits = GetPixelChannelTraits(image,channel);
122  if ((traits & UpdatePixelTrait) != 0)
123  channels++;
124  }
125  return(channels == 0 ? (size_t) 1 : channels);
126 }
127 
128 MagickExport Image *CompareImages(Image *image,const Image *reconstruct_image,
129  const MetricType metric,double *distortion,ExceptionInfo *exception)
130 {
131  CacheView
132  *highlight_view,
133  *image_view,
134  *reconstruct_view;
135 
136  const char
137  *artifact;
138 
139  double
140  fuzz;
141 
142  Image
143  *clone_image,
144  *difference_image,
145  *highlight_image;
146 
148  status;
149 
150  PixelInfo
151  highlight,
152  lowlight,
153  masklight;
154 
156  geometry;
157 
158  size_t
159  columns,
160  rows;
161 
162  ssize_t
163  y;
164 
165  assert(image != (Image *) NULL);
166  assert(image->signature == MagickCoreSignature);
167  if (image->debug != MagickFalse)
168  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
169  assert(reconstruct_image != (const Image *) NULL);
170  assert(reconstruct_image->signature == MagickCoreSignature);
171  assert(distortion != (double *) NULL);
172  *distortion=0.0;
173  if (image->debug != MagickFalse)
174  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
175  status=GetImageDistortion(image,reconstruct_image,metric,distortion,
176  exception);
177  if (status == MagickFalse)
178  return((Image *) NULL);
179  columns=MagickMax(image->columns,reconstruct_image->columns);
180  rows=MagickMax(image->rows,reconstruct_image->rows);
181  SetGeometry(image,&geometry);
182  geometry.width=columns;
183  geometry.height=rows;
184  clone_image=CloneImage(image,0,0,MagickTrue,exception);
185  if (clone_image == (Image *) NULL)
186  return((Image *) NULL);
187  (void) SetImageMask(clone_image,ReadPixelMask,(Image *) NULL,exception);
188  difference_image=ExtentImage(clone_image,&geometry,exception);
189  clone_image=DestroyImage(clone_image);
190  if (difference_image == (Image *) NULL)
191  return((Image *) NULL);
192  (void) SetImageAlphaChannel(difference_image,OpaqueAlphaChannel,exception);
193  highlight_image=CloneImage(image,columns,rows,MagickTrue,exception);
194  if (highlight_image == (Image *) NULL)
195  {
196  difference_image=DestroyImage(difference_image);
197  return((Image *) NULL);
198  }
199  status=SetImageStorageClass(highlight_image,DirectClass,exception);
200  if (status == MagickFalse)
201  {
202  difference_image=DestroyImage(difference_image);
203  highlight_image=DestroyImage(highlight_image);
204  return((Image *) NULL);
205  }
206  (void) SetImageMask(highlight_image,ReadPixelMask,(Image *) NULL,exception);
207  (void) SetImageAlphaChannel(highlight_image,OpaqueAlphaChannel,exception);
208  (void) QueryColorCompliance("#f1001ecc",AllCompliance,&highlight,exception);
209  artifact=GetImageArtifact(image,"compare:highlight-color");
210  if (artifact != (const char *) NULL)
211  (void) QueryColorCompliance(artifact,AllCompliance,&highlight,exception);
212  (void) QueryColorCompliance("#ffffffcc",AllCompliance,&lowlight,exception);
213  artifact=GetImageArtifact(image,"compare:lowlight-color");
214  if (artifact != (const char *) NULL)
215  (void) QueryColorCompliance(artifact,AllCompliance,&lowlight,exception);
216  (void) QueryColorCompliance("#888888cc",AllCompliance,&masklight,exception);
217  artifact=GetImageArtifact(image,"compare:masklight-color");
218  if (artifact != (const char *) NULL)
219  (void) QueryColorCompliance(artifact,AllCompliance,&masklight,exception);
220  /*
221  Generate difference image.
222  */
223  status=MagickTrue;
224  fuzz=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 
287  PixelChannel channel = GetPixelChannelChannel(image,i);
288  PixelTrait traits = GetPixelChannelTraits(image,channel);
289  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
290  channel);
291  if ((traits == UndefinedPixelTrait) ||
292  (reconstruct_traits == UndefinedPixelTrait) ||
293  ((reconstruct_traits & UpdatePixelTrait) == 0))
294  continue;
295  if (channel == AlphaPixelChannel)
296  distance=(double) p[i]-GetPixelChannel(reconstruct_image,channel,q);
297  else
298  distance=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
299  if ((distance*distance) > fuzz)
300  {
301  difference=MagickTrue;
302  break;
303  }
304  }
305  if (difference == MagickFalse)
306  SetPixelViaPixelInfo(highlight_image,&lowlight,r);
307  else
308  SetPixelViaPixelInfo(highlight_image,&highlight,r);
309  p+=GetPixelChannels(image);
310  q+=GetPixelChannels(reconstruct_image);
311  r+=GetPixelChannels(highlight_image);
312  }
313  sync=SyncCacheViewAuthenticPixels(highlight_view,exception);
314  if (sync == MagickFalse)
315  status=MagickFalse;
316  }
317  highlight_view=DestroyCacheView(highlight_view);
318  reconstruct_view=DestroyCacheView(reconstruct_view);
319  image_view=DestroyCacheView(image_view);
320  (void) CompositeImage(difference_image,highlight_image,image->compose,
321  MagickTrue,0,0,exception);
322  highlight_image=DestroyImage(highlight_image);
323  if (status == MagickFalse)
324  difference_image=DestroyImage(difference_image);
325  return(difference_image);
326 }
327 
328 /*
329 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
330 % %
331 % %
332 % %
333 % G e t I m a g e D i s t o r t i o n %
334 % %
335 % %
336 % %
337 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
338 %
339 % GetImageDistortion() compares one or more pixel channels of an image to a
340 % reconstructed image and returns the specified distortion metric.
341 %
342 % The format of the GetImageDistortion method is:
343 %
344 % MagickBooleanType GetImageDistortion(const Image *image,
345 % const Image *reconstruct_image,const MetricType metric,
346 % double *distortion,ExceptionInfo *exception)
347 %
348 % A description of each parameter follows:
349 %
350 % o image: the image.
351 %
352 % o reconstruct_image: the reconstruct image.
353 %
354 % o metric: the metric.
355 %
356 % o distortion: the computed distortion between the images.
357 %
358 % o exception: return any errors or warnings in this structure.
359 %
360 */
361 
363  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
364 {
365  CacheView
366  *image_view,
367  *reconstruct_view;
368 
369  double
370  fuzz;
371 
373  status;
374 
375  size_t
376  columns,
377  rows;
378 
379  ssize_t
380  y;
381 
382  /*
383  Compute the absolute difference in pixels between two images.
384  */
385  status=MagickTrue;
386  fuzz=(double) MagickMin(GetPixelChannels(image),
387  GetPixelChannels(reconstruct_image))*
388  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  distance,
425  Sa;
426 
428  difference;
429 
430  register ssize_t
431  i;
432 
433  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
434  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
435  {
436  p+=GetPixelChannels(image);
437  q+=GetPixelChannels(reconstruct_image);
438  continue;
439  }
440  difference=MagickFalse;
441  distance=0.0;
442  Sa=QuantumScale*GetPixelAlpha(image,p);
443  Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
444  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
445  {
446  double
447  pixel;
448 
449  PixelChannel channel = GetPixelChannelChannel(image,i);
450  PixelTrait traits = GetPixelChannelTraits(image,channel);
451  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
452  channel);
453  if ((traits == UndefinedPixelTrait) ||
454  (reconstruct_traits == UndefinedPixelTrait) ||
455  ((reconstruct_traits & UpdatePixelTrait) == 0))
456  continue;
457  if (channel == AlphaPixelChannel)
458  pixel=(double) p[i]-GetPixelChannel(reconstruct_image,channel,q);
459  else
460  pixel=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
461  distance+=pixel*pixel;
462  if (distance > fuzz)
463  {
464  channel_distortion[i]++;
465  difference=MagickTrue;
466  }
467  }
468  if (difference != MagickFalse)
469  channel_distortion[CompositePixelChannel]++;
470  p+=GetPixelChannels(image);
471  q+=GetPixelChannels(reconstruct_image);
472  }
473 #if defined(MAGICKCORE_OPENMP_SUPPORT)
474  #pragma omp critical (MagickCore_GetAbsoluteDistortion)
475 #endif
476  for (j=0; j <= MaxPixelChannels; j++)
477  distortion[j]+=channel_distortion[j];
478  }
479  reconstruct_view=DestroyCacheView(reconstruct_view);
480  image_view=DestroyCacheView(image_view);
481  return(status);
482 }
483 
485  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
486 {
487  CacheView
488  *image_view,
489  *reconstruct_view;
490 
491  double
492  area;
493 
495  status;
496 
497  register ssize_t
498  j;
499 
500  size_t
501  columns,
502  rows;
503 
504  ssize_t
505  y;
506 
507  status=MagickTrue;
508  rows=MagickMax(image->rows,reconstruct_image->rows);
509  columns=MagickMax(image->columns,reconstruct_image->columns);
510  area=0.0;
511  image_view=AcquireVirtualCacheView(image,exception);
512  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
513 #if defined(MAGICKCORE_OPENMP_SUPPORT)
514  #pragma omp parallel for schedule(static) shared(status) \
515  magick_number_threads(image,image,rows,1) reduction(+:area)
516 #endif
517  for (y=0; y < (ssize_t) rows; y++)
518  {
519  double
520  channel_distortion[MaxPixelChannels+1];
521 
522  register const Quantum
523  *magick_restrict p,
524  *magick_restrict q;
525 
526  register ssize_t
527  x;
528 
529  if (status == MagickFalse)
530  continue;
531  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
532  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
533  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
534  {
535  status=MagickFalse;
536  continue;
537  }
538  (void) memset(channel_distortion,0,sizeof(channel_distortion));
539  for (x=0; x < (ssize_t) columns; x++)
540  {
541  double
542  Da,
543  Sa;
544 
545  register ssize_t
546  i;
547 
548  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
549  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
550  {
551  p+=GetPixelChannels(image);
552  q+=GetPixelChannels(reconstruct_image);
553  continue;
554  }
555  Sa=QuantumScale*GetPixelAlpha(image,p);
556  Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
557  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
558  {
559  double
560  distance;
561 
562  PixelChannel channel = GetPixelChannelChannel(image,i);
563  PixelTrait traits = GetPixelChannelTraits(image,channel);
564  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
565  channel);
566  if ((traits == UndefinedPixelTrait) ||
567  (reconstruct_traits == UndefinedPixelTrait) ||
568  ((reconstruct_traits & UpdatePixelTrait) == 0))
569  continue;
570  if (channel == AlphaPixelChannel)
571  distance=QuantumScale*(p[i]-GetPixelChannel(reconstruct_image,
572  channel,q));
573  else
574  distance=QuantumScale*(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
575  channel,q));
576  channel_distortion[i]+=distance*distance;
577  channel_distortion[CompositePixelChannel]+=distance*distance;
578  }
579  area++;
580  p+=GetPixelChannels(image);
581  q+=GetPixelChannels(reconstruct_image);
582  }
583 #if defined(MAGICKCORE_OPENMP_SUPPORT)
584  #pragma omp critical (MagickCore_GetFuzzDistortion)
585 #endif
586  for (j=0; j <= MaxPixelChannels; j++)
587  distortion[j]+=channel_distortion[j];
588  }
589  reconstruct_view=DestroyCacheView(reconstruct_view);
590  image_view=DestroyCacheView(image_view);
591  area=PerceptibleReciprocal(area);
592  for (j=0; j <= MaxPixelChannels; j++)
593  distortion[j]*=area;
594  distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
595  distortion[CompositePixelChannel]=sqrt(distortion[CompositePixelChannel]);
596  return(status);
597 }
598 
600  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
601 {
602  CacheView
603  *image_view,
604  *reconstruct_view;
605 
606  double
607  area;
608 
610  status;
611 
612  register ssize_t
613  j;
614 
615  size_t
616  columns,
617  rows;
618 
619  ssize_t
620  y;
621 
622  status=MagickTrue;
623  rows=MagickMax(image->rows,reconstruct_image->rows);
624  columns=MagickMax(image->columns,reconstruct_image->columns);
625  area=0.0;
626  image_view=AcquireVirtualCacheView(image,exception);
627  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
628 #if defined(MAGICKCORE_OPENMP_SUPPORT)
629  #pragma omp parallel for schedule(static) shared(status) \
630  magick_number_threads(image,image,rows,1) reduction(+:area)
631 #endif
632  for (y=0; y < (ssize_t) rows; y++)
633  {
634  double
635  channel_distortion[MaxPixelChannels+1];
636 
637  register const Quantum
638  *magick_restrict p,
639  *magick_restrict q;
640 
641  register ssize_t
642  x;
643 
644  if (status == MagickFalse)
645  continue;
646  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
647  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
648  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
649  {
650  status=MagickFalse;
651  continue;
652  }
653  (void) memset(channel_distortion,0,sizeof(channel_distortion));
654  for (x=0; x < (ssize_t) columns; x++)
655  {
656  double
657  Da,
658  Sa;
659 
660  register ssize_t
661  i;
662 
663  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
664  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
665  {
666  p+=GetPixelChannels(image);
667  q+=GetPixelChannels(reconstruct_image);
668  continue;
669  }
670  Sa=QuantumScale*GetPixelAlpha(image,p);
671  Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
672  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
673  {
674  double
675  distance;
676 
677  PixelChannel channel = GetPixelChannelChannel(image,i);
678  PixelTrait traits = GetPixelChannelTraits(image,channel);
679  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
680  channel);
681  if ((traits == UndefinedPixelTrait) ||
682  (reconstruct_traits == UndefinedPixelTrait) ||
683  ((reconstruct_traits & UpdatePixelTrait) == 0))
684  continue;
685  if (channel == AlphaPixelChannel)
686  distance=QuantumScale*fabs((double) p[i]-
687  GetPixelChannel(reconstruct_image,channel,q));
688  else
689  distance=QuantumScale*fabs(Sa*p[i]-Da*
690  GetPixelChannel(reconstruct_image,channel,q));
691  channel_distortion[i]+=distance;
692  channel_distortion[CompositePixelChannel]+=distance;
693  }
694  area++;
695  p+=GetPixelChannels(image);
696  q+=GetPixelChannels(reconstruct_image);
697  }
698 #if defined(MAGICKCORE_OPENMP_SUPPORT)
699  #pragma omp critical (MagickCore_GetMeanAbsoluteError)
700 #endif
701  for (j=0; j <= MaxPixelChannels; j++)
702  distortion[j]+=channel_distortion[j];
703  }
704  reconstruct_view=DestroyCacheView(reconstruct_view);
705  image_view=DestroyCacheView(image_view);
706  area=PerceptibleReciprocal(area);
707  for (j=0; j <= MaxPixelChannels; j++)
708  distortion[j]*=area;
709  distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
710  return(status);
711 }
712 
714  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
715 {
716  CacheView
717  *image_view,
718  *reconstruct_view;
719 
721  status;
722 
723  double
724  area,
725  maximum_error,
726  mean_error;
727 
728  size_t
729  columns,
730  rows;
731 
732  ssize_t
733  y;
734 
735  status=MagickTrue;
736  area=0.0;
737  maximum_error=0.0;
738  mean_error=0.0;
739  rows=MagickMax(image->rows,reconstruct_image->rows);
740  columns=MagickMax(image->columns,reconstruct_image->columns);
741  image_view=AcquireVirtualCacheView(image,exception);
742  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
743  for (y=0; y < (ssize_t) rows; y++)
744  {
745  register const Quantum
746  *magick_restrict p,
747  *magick_restrict q;
748 
749  register ssize_t
750  x;
751 
752  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
753  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
754  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
755  {
756  status=MagickFalse;
757  break;
758  }
759  for (x=0; x < (ssize_t) columns; x++)
760  {
761  double
762  Da,
763  Sa;
764 
765  register ssize_t
766  i;
767 
768  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
769  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
770  {
771  p+=GetPixelChannels(image);
772  q+=GetPixelChannels(reconstruct_image);
773  continue;
774  }
775  Sa=QuantumScale*GetPixelAlpha(image,p);
776  Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
777  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
778  {
779  double
780  distance;
781 
782  PixelChannel channel = GetPixelChannelChannel(image,i);
783  PixelTrait traits = GetPixelChannelTraits(image,channel);
784  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
785  channel);
786  if ((traits == UndefinedPixelTrait) ||
787  (reconstruct_traits == UndefinedPixelTrait) ||
788  ((reconstruct_traits & UpdatePixelTrait) == 0))
789  continue;
790  if (channel == AlphaPixelChannel)
791  distance=fabs((double) p[i]-
792  GetPixelChannel(reconstruct_image,channel,q));
793  else
794  distance=fabs(Sa*p[i]-Da*
795  GetPixelChannel(reconstruct_image,channel,q));
796  distortion[i]+=distance;
797  distortion[CompositePixelChannel]+=distance;
798  mean_error+=distance*distance;
799  if (distance > maximum_error)
800  maximum_error=distance;
801  area++;
802  }
803  p+=GetPixelChannels(image);
804  q+=GetPixelChannels(reconstruct_image);
805  }
806  }
807  reconstruct_view=DestroyCacheView(reconstruct_view);
808  image_view=DestroyCacheView(image_view);
809  image->error.mean_error_per_pixel=distortion[CompositePixelChannel]/area;
810  image->error.normalized_mean_error=QuantumScale*QuantumScale*mean_error/area;
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  register 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  register const Quantum
854  *magick_restrict p,
855  *magick_restrict q;
856 
857  register 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  register 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  register 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  register const Quantum
990  *magick_restrict p,
991  *magick_restrict q;
992 
993  register 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  register const Quantum
1021  *magick_restrict p,
1022  *magick_restrict q;
1023 
1024  register 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  register const Quantum
1156  *magick_restrict p,
1157  *magick_restrict q;
1158 
1159  register 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  register 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]-
1206  GetPixelChannel(reconstruct_image,channel,q));
1207  else
1208  distance=QuantumScale*fabs(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  register 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  register 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  register 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  register 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  c1,
1373  c2,
1374  radius,
1375  sigma;
1376 
1377  KernelInfo
1378  *kernel_info;
1379 
1381  status;
1382 
1383  register ssize_t
1384  i;
1385 
1386  size_t
1387  columns,
1388  rows;
1389 
1390  ssize_t
1391  y;
1392 
1393  /*
1394  Compute structural similarity index @
1395  https://en.wikipedia.org/wiki/Structural_similarity.
1396  */
1397  radius=SSIMRadius;
1398  artifact=GetImageArtifact(image,"compare:ssim-radius");
1399  if (artifact != (const char *) NULL)
1400  radius=StringToDouble(artifact,(char **) NULL);
1401  sigma=SSIMSigma;
1402  artifact=GetImageArtifact(image,"compare:ssim-sigma");
1403  if (artifact != (const char *) NULL)
1404  sigma=StringToDouble(artifact,(char **) NULL);
1405  (void) FormatLocaleString(geometry,MagickPathExtent,"gaussian:%.20gx%.20g",
1406  radius,sigma);
1407  kernel_info=AcquireKernelInfo(geometry,exception);
1408  if (kernel_info == (KernelInfo *) NULL)
1409  ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
1410  image->filename);
1411  c1=pow(SSIMK1*SSIML,2.0);
1412  artifact=GetImageArtifact(image,"compare:ssim-k1");
1413  if (artifact != (const char *) NULL)
1414  c1=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
1415  c2=pow(SSIMK2*SSIML,2.0);
1416  artifact=GetImageArtifact(image,"compare:ssim-k2");
1417  if (artifact != (const char *) NULL)
1418  c2=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
1419  status=MagickTrue;
1420  rows=MagickMax(image->rows,reconstruct_image->rows);
1421  columns=MagickMax(image->columns,reconstruct_image->columns);
1422  image_view=AcquireVirtualCacheView(image,exception);
1423  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1424 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1425  #pragma omp parallel for schedule(static) shared(status) \
1426  magick_number_threads(image,reconstruct_image,rows,1)
1427 #endif
1428  for (y=0; y < (ssize_t) rows; y++)
1429  {
1430  double
1431  channel_distortion[MaxPixelChannels+1];
1432 
1433  register const Quantum
1434  *magick_restrict p,
1435  *magick_restrict q;
1436 
1437  register ssize_t
1438  i,
1439  x;
1440 
1441  if (status == MagickFalse)
1442  continue;
1443  p=GetCacheViewVirtualPixels(image_view,-((ssize_t) kernel_info->width/2L),y-
1444  ((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
1445  kernel_info->height,exception);
1446  q=GetCacheViewVirtualPixels(reconstruct_view,-((ssize_t) kernel_info->width/
1447  2L),y-((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
1448  kernel_info->height,exception);
1449  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1450  {
1451  status=MagickFalse;
1452  continue;
1453  }
1454  (void) memset(channel_distortion,0,sizeof(channel_distortion));
1455  for (x=0; x < (ssize_t) columns; x++)
1456  {
1457  double
1458  x_pixel_mu[MaxPixelChannels+1],
1459  x_pixel_sigma_squared[MaxPixelChannels+1],
1460  xy_sigma[MaxPixelChannels+1],
1461  y_pixel_mu[MaxPixelChannels+1],
1462  y_pixel_sigma_squared[MaxPixelChannels+1];
1463 
1464  register const Quantum
1465  *magick_restrict reference,
1466  *magick_restrict target;
1467 
1468  register MagickRealType
1469  *k;
1470 
1471  ssize_t
1472  v;
1473 
1474  (void) memset(x_pixel_mu,0,sizeof(x_pixel_mu));
1475  (void) memset(x_pixel_sigma_squared,0,sizeof(x_pixel_sigma_squared));
1476  (void) memset(xy_sigma,0,sizeof(xy_sigma));
1477  (void) memset(x_pixel_sigma_squared,0,sizeof(y_pixel_sigma_squared));
1478  (void) memset(y_pixel_mu,0,sizeof(y_pixel_mu));
1479  (void) memset(y_pixel_sigma_squared,0,sizeof(y_pixel_sigma_squared));
1480  k=kernel_info->values;
1481  reference=p;
1482  target=q;
1483  for (v=0; v < (ssize_t) kernel_info->height; v++)
1484  {
1485  register ssize_t
1486  u;
1487 
1488  for (u=0; u < (ssize_t) kernel_info->width; u++)
1489  {
1490  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1491  {
1492  double
1493  x_pixel,
1494  y_pixel;
1495 
1496  PixelChannel channel = GetPixelChannelChannel(image,i);
1497  PixelTrait traits = GetPixelChannelTraits(image,channel);
1498  PixelTrait reconstruct_traits = GetPixelChannelTraits(
1499  reconstruct_image,channel);
1500  if ((traits == UndefinedPixelTrait) ||
1501  (reconstruct_traits == UndefinedPixelTrait) ||
1502  ((reconstruct_traits & UpdatePixelTrait) == 0))
1503  continue;
1504  x_pixel=QuantumScale*reference[i];
1505  x_pixel_mu[i]+=(*k)*x_pixel;
1506  x_pixel_sigma_squared[i]+=(*k)*x_pixel*x_pixel;
1507  y_pixel=QuantumScale*
1508  GetPixelChannel(reconstruct_image,channel,target);
1509  y_pixel_mu[i]+=(*k)*y_pixel;
1510  y_pixel_sigma_squared[i]+=(*k)*y_pixel*y_pixel;
1511  xy_sigma[i]+=(*k)*x_pixel*y_pixel;
1512  }
1513  k++;
1514  reference+=GetPixelChannels(image);
1515  target+=GetPixelChannels(reconstruct_image);
1516  }
1517  reference+=GetPixelChannels(image)*columns;
1518  target+=GetPixelChannels(reconstruct_image)*columns;
1519  }
1520  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1521  {
1522  double
1523  ssim,
1524  x_pixel_mu_squared,
1525  x_pixel_sigmas_squared,
1526  xy_mu,
1527  xy_sigmas,
1528  y_pixel_mu_squared,
1529  y_pixel_sigmas_squared;
1530 
1531  PixelChannel channel = GetPixelChannelChannel(image,i);
1532  PixelTrait traits = GetPixelChannelTraits(image,channel);
1533  PixelTrait reconstruct_traits = GetPixelChannelTraits(
1534  reconstruct_image,channel);
1535  if ((traits == UndefinedPixelTrait) ||
1536  (reconstruct_traits == UndefinedPixelTrait) ||
1537  ((reconstruct_traits & UpdatePixelTrait) == 0))
1538  continue;
1539  x_pixel_mu_squared=x_pixel_mu[i]*x_pixel_mu[i];
1540  y_pixel_mu_squared=y_pixel_mu[i]*y_pixel_mu[i];
1541  xy_mu=x_pixel_mu[i]*y_pixel_mu[i];
1542  xy_sigmas=xy_sigma[i]-xy_mu;
1543  x_pixel_sigmas_squared=x_pixel_sigma_squared[i]-x_pixel_mu_squared;
1544  y_pixel_sigmas_squared=y_pixel_sigma_squared[i]-y_pixel_mu_squared;
1545  ssim=((2.0*xy_mu+c1)*(2.0*xy_sigmas+c2))/
1546  ((x_pixel_mu_squared+y_pixel_mu_squared+c1)*
1547  (x_pixel_sigmas_squared+y_pixel_sigmas_squared+c2));
1548  channel_distortion[i]+=ssim;
1549  channel_distortion[CompositePixelChannel]+=ssim;
1550  }
1551  p+=GetPixelChannels(image);
1552  q+=GetPixelChannels(reconstruct_image);
1553  }
1554 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1555  #pragma omp critical (MagickCore_GetStructuralSimilarityDistortion)
1556 #endif
1557  for (i=0; i <= MaxPixelChannels; i++)
1558  distortion[i]+=channel_distortion[i];
1559  }
1560  image_view=DestroyCacheView(image_view);
1561  reconstruct_view=DestroyCacheView(reconstruct_view);
1562  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1563  {
1564  PixelChannel channel = GetPixelChannelChannel(image,i);
1565  PixelTrait traits = GetPixelChannelTraits(image,channel);
1566  if ((traits == UndefinedPixelTrait) || ((traits & UpdatePixelTrait) == 0))
1567  continue;
1568  distortion[i]/=((double) columns*rows);
1569  }
1570  distortion[CompositePixelChannel]/=((double) columns*rows);
1571  distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
1572  kernel_info=DestroyKernelInfo(kernel_info);
1573  return(status);
1574 }
1575 
1577  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1578 {
1580  status;
1581 
1582  register ssize_t
1583  i;
1584 
1585  status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1586  distortion,exception);
1587  for (i=0; i <= MaxPixelChannels; i++)
1588  distortion[i]=(1.0-(distortion[i]))/2.0;
1589  return(status);
1590 }
1591 
1593  const Image *reconstruct_image,const MetricType metric,double *distortion,
1594  ExceptionInfo *exception)
1595 {
1596  double
1597  *channel_distortion;
1598 
1600  status;
1601 
1602  size_t
1603  length;
1604 
1605  assert(image != (Image *) NULL);
1606  assert(image->signature == MagickCoreSignature);
1607  if (image->debug != MagickFalse)
1608  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1609  assert(reconstruct_image != (const Image *) NULL);
1610  assert(reconstruct_image->signature == MagickCoreSignature);
1611  assert(distortion != (double *) NULL);
1612  *distortion=0.0;
1613  if (image->debug != MagickFalse)
1614  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1615  /*
1616  Get image distortion.
1617  */
1618  length=MaxPixelChannels+1;
1619  channel_distortion=(double *) AcquireQuantumMemory(length,
1620  sizeof(*channel_distortion));
1621  if (channel_distortion == (double *) NULL)
1622  ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1623  (void) memset(channel_distortion,0,length*
1624  sizeof(*channel_distortion));
1625  switch (metric)
1626  {
1627  case AbsoluteErrorMetric:
1628  {
1629  status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1630  exception);
1631  break;
1632  }
1633  case FuzzErrorMetric:
1634  {
1635  status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1636  exception);
1637  break;
1638  }
1640  {
1641  status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1642  channel_distortion,exception);
1643  break;
1644  }
1646  {
1647  status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1648  exception);
1649  break;
1650  }
1652  {
1653  status=GetMeanSquaredDistortion(image,reconstruct_image,
1654  channel_distortion,exception);
1655  break;
1656  }
1658  default:
1659  {
1660  status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1661  channel_distortion,exception);
1662  break;
1663  }
1665  {
1666  status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1667  channel_distortion,exception);
1668  break;
1669  }
1671  {
1672  status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1673  channel_distortion,exception);
1674  break;
1675  }
1677  {
1678  status=GetPerceptualHashDistortion(image,reconstruct_image,
1679  channel_distortion,exception);
1680  break;
1681  }
1683  {
1684  status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1685  channel_distortion,exception);
1686  break;
1687  }
1689  {
1690  status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1691  channel_distortion,exception);
1692  break;
1693  }
1695  {
1696  status=GetStructuralDisimilarityDistortion(image,reconstruct_image,
1697  channel_distortion,exception);
1698  break;
1699  }
1700  }
1701  *distortion=channel_distortion[CompositePixelChannel];
1702  channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1703  (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
1704  *distortion);
1705  return(status);
1706 }
1707 
1708 /*
1709 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1710 % %
1711 % %
1712 % %
1713 % G e t I m a g e D i s t o r t i o n s %
1714 % %
1715 % %
1716 % %
1717 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1718 %
1719 % GetImageDistortions() compares the pixel channels of an image to a
1720 % reconstructed image and returns the specified distortion metric for each
1721 % channel.
1722 %
1723 % The format of the GetImageDistortions method is:
1724 %
1725 % double *GetImageDistortions(const Image *image,
1726 % const Image *reconstruct_image,const MetricType metric,
1727 % ExceptionInfo *exception)
1728 %
1729 % A description of each parameter follows:
1730 %
1731 % o image: the image.
1732 %
1733 % o reconstruct_image: the reconstruct image.
1734 %
1735 % o metric: the metric.
1736 %
1737 % o exception: return any errors or warnings in this structure.
1738 %
1739 */
1741  const Image *reconstruct_image,const MetricType metric,
1742  ExceptionInfo *exception)
1743 {
1744  double
1745  *channel_distortion;
1746 
1748  status;
1749 
1750  size_t
1751  length;
1752 
1753  assert(image != (Image *) NULL);
1754  assert(image->signature == MagickCoreSignature);
1755  if (image->debug != MagickFalse)
1756  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1757  assert(reconstruct_image != (const Image *) NULL);
1758  assert(reconstruct_image->signature == MagickCoreSignature);
1759  if (image->debug != MagickFalse)
1760  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1761  /*
1762  Get image distortion.
1763  */
1764  length=MaxPixelChannels+1UL;
1765  channel_distortion=(double *) AcquireQuantumMemory(length,
1766  sizeof(*channel_distortion));
1767  if (channel_distortion == (double *) NULL)
1768  ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1769  (void) memset(channel_distortion,0,length*
1770  sizeof(*channel_distortion));
1771  status=MagickTrue;
1772  switch (metric)
1773  {
1774  case AbsoluteErrorMetric:
1775  {
1776  status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1777  exception);
1778  break;
1779  }
1780  case FuzzErrorMetric:
1781  {
1782  status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1783  exception);
1784  break;
1785  }
1787  {
1788  status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1789  channel_distortion,exception);
1790  break;
1791  }
1793  {
1794  status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1795  exception);
1796  break;
1797  }
1799  {
1800  status=GetMeanSquaredDistortion(image,reconstruct_image,
1801  channel_distortion,exception);
1802  break;
1803  }
1805  default:
1806  {
1807  status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1808  channel_distortion,exception);
1809  break;
1810  }
1812  {
1813  status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1814  channel_distortion,exception);
1815  break;
1816  }
1818  {
1819  status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1820  channel_distortion,exception);
1821  break;
1822  }
1824  {
1825  status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1826  channel_distortion,exception);
1827  break;
1828  }
1830  {
1831  status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1832  channel_distortion,exception);
1833  break;
1834  }
1836  {
1837  status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1838  channel_distortion,exception);
1839  break;
1840  }
1842  {
1843  status=GetStructuralDisimilarityDistortion(image,reconstruct_image,
1844  channel_distortion,exception);
1845  break;
1846  }
1847  }
1848  if (status == MagickFalse)
1849  {
1850  channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1851  return((double *) NULL);
1852  }
1853  return(channel_distortion);
1854 }
1855 
1856 /*
1857 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1858 % %
1859 % %
1860 % %
1861 % I s I m a g e s E q u a l %
1862 % %
1863 % %
1864 % %
1865 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1866 %
1867 % IsImagesEqual() compare the pixels of two images and returns immediately
1868 % if any pixel is not identical.
1869 %
1870 % The format of the IsImagesEqual method is:
1871 %
1872 % MagickBooleanType IsImagesEqual(const Image *image,
1873 % const Image *reconstruct_image,ExceptionInfo *exception)
1874 %
1875 % A description of each parameter follows.
1876 %
1877 % o image: the image.
1878 %
1879 % o reconstruct_image: the reconstruct image.
1880 %
1881 % o exception: return any errors or warnings in this structure.
1882 %
1883 */
1885  const Image *reconstruct_image,ExceptionInfo *exception)
1886 {
1887  CacheView
1888  *image_view,
1889  *reconstruct_view;
1890 
1891  size_t
1892  columns,
1893  rows;
1894 
1895  ssize_t
1896  y;
1897 
1898  assert(image != (Image *) NULL);
1899  assert(image->signature == MagickCoreSignature);
1900  assert(reconstruct_image != (const Image *) NULL);
1901  assert(reconstruct_image->signature == MagickCoreSignature);
1902  rows=MagickMax(image->rows,reconstruct_image->rows);
1903  columns=MagickMax(image->columns,reconstruct_image->columns);
1904  image_view=AcquireVirtualCacheView(image,exception);
1905  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1906  for (y=0; y < (ssize_t) rows; y++)
1907  {
1908  register const Quantum
1909  *magick_restrict p,
1910  *magick_restrict q;
1911 
1912  register ssize_t
1913  x;
1914 
1915  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1916  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1917  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
1918  break;
1919  for (x=0; x < (ssize_t) columns; x++)
1920  {
1921  register ssize_t
1922  i;
1923 
1924  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1925  {
1926  double
1927  distance;
1928 
1929  PixelChannel channel = GetPixelChannelChannel(image,i);
1930  PixelTrait traits = GetPixelChannelTraits(image,channel);
1931  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1932  channel);
1933  if ((traits == UndefinedPixelTrait) ||
1934  (reconstruct_traits == UndefinedPixelTrait) ||
1935  ((reconstruct_traits & UpdatePixelTrait) == 0))
1936  continue;
1937  distance=fabs(p[i]-(double) GetPixelChannel(reconstruct_image,
1938  channel,q));
1939  if (distance >= MagickEpsilon)
1940  break;
1941  }
1942  if (i < (ssize_t) GetPixelChannels(image))
1943  break;
1944  p+=GetPixelChannels(image);
1945  q+=GetPixelChannels(reconstruct_image);
1946  }
1947  if (x < (ssize_t) columns)
1948  break;
1949  }
1950  reconstruct_view=DestroyCacheView(reconstruct_view);
1951  image_view=DestroyCacheView(image_view);
1952  return(y < (ssize_t) rows ? MagickFalse : MagickTrue);
1953 }
1954 
1955 /*
1956 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1957 % %
1958 % %
1959 % %
1960 % S e t I m a g e C o l o r M e t r i c %
1961 % %
1962 % %
1963 % %
1964 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1965 %
1966 % SetImageColorMetric() measures the difference between colors at each pixel
1967 % location of two images. A value other than 0 means the colors match
1968 % exactly. Otherwise an error measure is computed by summing over all
1969 % pixels in an image the distance squared in RGB space between each image
1970 % pixel and its corresponding pixel in the reconstruct image. The error
1971 % measure is assigned to these image members:
1972 %
1973 % o mean_error_per_pixel: The mean error for any single pixel in
1974 % the image.
1975 %
1976 % o normalized_mean_error: The normalized mean quantization error for
1977 % any single pixel in the image. This distance measure is normalized to
1978 % a range between 0 and 1. It is independent of the range of red, green,
1979 % and blue values in the image.
1980 %
1981 % o normalized_maximum_error: The normalized maximum quantization
1982 % error for any single pixel in the image. This distance measure is
1983 % normalized to a range between 0 and 1. It is independent of the range
1984 % of red, green, and blue values in your image.
1985 %
1986 % A small normalized mean square error, accessed as
1987 % image->normalized_mean_error, suggests the images are very similar in
1988 % spatial layout and color.
1989 %
1990 % The format of the SetImageColorMetric method is:
1991 %
1992 % MagickBooleanType SetImageColorMetric(Image *image,
1993 % const Image *reconstruct_image,ExceptionInfo *exception)
1994 %
1995 % A description of each parameter follows.
1996 %
1997 % o image: the image.
1998 %
1999 % o reconstruct_image: the reconstruct image.
2000 %
2001 % o exception: return any errors or warnings in this structure.
2002 %
2003 */
2005  const Image *reconstruct_image,ExceptionInfo *exception)
2006 {
2007  CacheView
2008  *image_view,
2009  *reconstruct_view;
2010 
2011  double
2012  area,
2013  maximum_error,
2014  mean_error,
2015  mean_error_per_pixel;
2016 
2018  status;
2019 
2020  size_t
2021  columns,
2022  rows;
2023 
2024  ssize_t
2025  y;
2026 
2027  assert(image != (Image *) NULL);
2028  assert(image->signature == MagickCoreSignature);
2029  assert(reconstruct_image != (const Image *) NULL);
2030  assert(reconstruct_image->signature == MagickCoreSignature);
2031  area=0.0;
2032  maximum_error=0.0;
2033  mean_error_per_pixel=0.0;
2034  mean_error=0.0;
2035  rows=MagickMax(image->rows,reconstruct_image->rows);
2036  columns=MagickMax(image->columns,reconstruct_image->columns);
2037  image_view=AcquireVirtualCacheView(image,exception);
2038  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
2039  for (y=0; y < (ssize_t) rows; y++)
2040  {
2041  register const Quantum
2042  *magick_restrict p,
2043  *magick_restrict q;
2044 
2045  register ssize_t
2046  x;
2047 
2048  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
2049  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
2050  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2051  break;
2052  for (x=0; x < (ssize_t) columns; x++)
2053  {
2054  register ssize_t
2055  i;
2056 
2057  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2058  {
2059  double
2060  distance;
2061 
2062  PixelChannel channel = GetPixelChannelChannel(image,i);
2063  PixelTrait traits = GetPixelChannelTraits(image,channel);
2064  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2065  channel);
2066  if ((traits == UndefinedPixelTrait) ||
2067  (reconstruct_traits == UndefinedPixelTrait) ||
2068  ((reconstruct_traits & UpdatePixelTrait) == 0))
2069  continue;
2070  distance=fabs(p[i]-(double) GetPixelChannel(reconstruct_image,
2071  channel,q));
2072  if (distance >= MagickEpsilon)
2073  {
2074  mean_error_per_pixel+=distance;
2075  mean_error+=distance*distance;
2076  if (distance > maximum_error)
2077  maximum_error=distance;
2078  }
2079  area++;
2080  }
2081  p+=GetPixelChannels(image);
2082  q+=GetPixelChannels(reconstruct_image);
2083  }
2084  }
2085  reconstruct_view=DestroyCacheView(reconstruct_view);
2086  image_view=DestroyCacheView(image_view);
2087  image->error.mean_error_per_pixel=(double) (mean_error_per_pixel/area);
2089  mean_error/area);
2090  image->error.normalized_maximum_error=(double) (QuantumScale*maximum_error);
2091  status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
2092  return(status);
2093 }
2094 
2095 /*
2096 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2097 % %
2098 % %
2099 % %
2100 % S i m i l a r i t y I m a g e %
2101 % %
2102 % %
2103 % %
2104 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2105 %
2106 % SimilarityImage() compares the reference image of the image and returns the
2107 % best match offset. In addition, it returns a similarity image such that an
2108 % exact match location is completely white and if none of the pixels match,
2109 % black, otherwise some gray level in-between.
2110 %
2111 % The format of the SimilarityImageImage method is:
2112 %
2113 % Image *SimilarityImage(const Image *image,const Image *reference,
2114 % const MetricType metric,const double similarity_threshold,
2115 % RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
2116 %
2117 % A description of each parameter follows:
2118 %
2119 % o image: the image.
2120 %
2121 % o reference: find an area of the image that closely resembles this image.
2122 %
2123 % o metric: the metric.
2124 %
2125 % o similarity_threshold: minimum distortion for (sub)image match.
2126 %
2127 % o offset: the best match offset of the reference image within the image.
2128 %
2129 % o similarity: the computed similarity between the images.
2130 %
2131 % o exception: return any errors or warnings in this structure.
2132 %
2133 */
2134 
2135 static double GetSimilarityMetric(const Image *image,const Image *reference,
2136  const MetricType metric,const ssize_t x_offset,const ssize_t y_offset,
2137  ExceptionInfo *exception)
2138 {
2139  double
2140  distortion;
2141 
2142  Image
2143  *similarity_image;
2144 
2146  status;
2147 
2149  geometry;
2150 
2151  SetGeometry(reference,&geometry);
2152  geometry.x=x_offset;
2153  geometry.y=y_offset;
2154  similarity_image=CropImage(image,&geometry,exception);
2155  if (similarity_image == (Image *) NULL)
2156  return(0.0);
2157  distortion=0.0;
2158  status=GetImageDistortion(similarity_image,reference,metric,&distortion,
2159  exception);
2160  similarity_image=DestroyImage(similarity_image);
2161  if (status == MagickFalse)
2162  return(0.0);
2163  return(distortion);
2164 }
2165 
2166 MagickExport Image *SimilarityImage(const Image *image,const Image *reference,
2167  const MetricType metric,const double similarity_threshold,
2168  RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
2169 {
2170 #define SimilarityImageTag "Similarity/Image"
2171 
2172  CacheView
2173  *similarity_view;
2174 
2175  Image
2176  *similarity_image;
2177 
2179  status;
2180 
2182  progress;
2183 
2184  ssize_t
2185  y;
2186 
2187  assert(image != (const Image *) NULL);
2188  assert(image->signature == MagickCoreSignature);
2189  if (image->debug != MagickFalse)
2190  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2191  assert(exception != (ExceptionInfo *) NULL);
2192  assert(exception->signature == MagickCoreSignature);
2193  assert(offset != (RectangleInfo *) NULL);
2194  SetGeometry(reference,offset);
2195  *similarity_metric=MagickMaximumValue;
2196  similarity_image=CloneImage(image,image->columns-reference->columns+1,
2197  image->rows-reference->rows+1,MagickTrue,exception);
2198  if (similarity_image == (Image *) NULL)
2199  return((Image *) NULL);
2200  status=SetImageStorageClass(similarity_image,DirectClass,exception);
2201  if (status == MagickFalse)
2202  {
2203  similarity_image=DestroyImage(similarity_image);
2204  return((Image *) NULL);
2205  }
2206  (void) SetImageAlphaChannel(similarity_image,DeactivateAlphaChannel,
2207  exception);
2208  /*
2209  Measure similarity of reference image against image.
2210  */
2211  status=MagickTrue;
2212  progress=0;
2213  similarity_view=AcquireAuthenticCacheView(similarity_image,exception);
2214 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2215  #pragma omp parallel for schedule(static) \
2216  shared(progress,status,similarity_metric) \
2217  magick_number_threads(image,image,image->rows-reference->rows+1,1)
2218 #endif
2219  for (y=0; y < (ssize_t) (image->rows-reference->rows+1); y++)
2220  {
2221  double
2222  similarity;
2223 
2224  register Quantum
2225  *magick_restrict q;
2226 
2227  register ssize_t
2228  x;
2229 
2230  if (status == MagickFalse)
2231  continue;
2232 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2233  #pragma omp flush(similarity_metric)
2234 #endif
2235  if (*similarity_metric <= similarity_threshold)
2236  continue;
2237  q=GetCacheViewAuthenticPixels(similarity_view,0,y,similarity_image->columns,
2238  1,exception);
2239  if (q == (Quantum *) NULL)
2240  {
2241  status=MagickFalse;
2242  continue;
2243  }
2244  for (x=0; x < (ssize_t) (image->columns-reference->columns+1); x++)
2245  {
2246  register ssize_t
2247  i;
2248 
2249 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2250  #pragma omp flush(similarity_metric)
2251 #endif
2252  if (*similarity_metric <= similarity_threshold)
2253  break;
2254  similarity=GetSimilarityMetric(image,reference,metric,x,y,exception);
2255 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2256  #pragma omp critical (MagickCore_SimilarityImage)
2257 #endif
2258  if ((metric == NormalizedCrossCorrelationErrorMetric) ||
2259  (metric == UndefinedErrorMetric))
2260  similarity=1.0-similarity;
2261  if (similarity < *similarity_metric)
2262  {
2263  offset->x=x;
2264  offset->y=y;
2265  *similarity_metric=similarity;
2266  }
2267  if (metric == PerceptualHashErrorMetric)
2268  similarity=MagickMin(0.01*similarity,1.0);
2269  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2270  {
2271  PixelChannel channel = GetPixelChannelChannel(image,i);
2272  PixelTrait traits = GetPixelChannelTraits(image,channel);
2273  PixelTrait similarity_traits=GetPixelChannelTraits(similarity_image,
2274  channel);
2275  if ((traits == UndefinedPixelTrait) ||
2276  (similarity_traits == UndefinedPixelTrait) ||
2277  ((similarity_traits & UpdatePixelTrait) == 0))
2278  continue;
2279  SetPixelChannel(similarity_image,channel,ClampToQuantum(QuantumRange-
2280  QuantumRange*similarity),q);
2281  }
2282  q+=GetPixelChannels(similarity_image);
2283  }
2284  if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
2285  status=MagickFalse;
2286  if (image->progress_monitor != (MagickProgressMonitor) NULL)
2287  {
2289  proceed;
2290 
2291 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2292  #pragma omp atomic
2293 #endif
2294  progress++;
2295  proceed=SetImageProgress(image,SimilarityImageTag,progress,image->rows);
2296  if (proceed == MagickFalse)
2297  status=MagickFalse;
2298  }
2299  }
2300  similarity_view=DestroyCacheView(similarity_view);
2301  if (status == MagickFalse)
2302  similarity_image=DestroyImage(similarity_image);
2303  return(similarity_image);
2304 }
size_t rows
Definition: image.h:172
#define magick_restrict
Definition: MagickCore.h:41
MagickDoubleType MagickRealType
Definition: magick-type.h:120
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:362
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:2166
MagickProgressMonitor progress_monitor
Definition: image.h:303
static Quantum GetPixelAlpha(const Image *magick_restrict image, const Quantum *magick_restrict pixel)
double phash[MaximumNumberOfPerceptualColorspaces+1][MaximumNumberOfImageMoments+1]
Definition: statistic.h:78
size_t height
Definition: morphology.h:108
MagickExport ChannelStatistics * GetImageStatistics(const Image *image, ExceptionInfo *exception)
Definition: statistic.c:1930
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:1576
#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:599
#define SimilarityImageTag
static Quantum GetPixelChannel(const Image *magick_restrict image, const PixelChannel channel, const Quantum *magick_restrict pixel)
static size_t GetImageChannels(const Image *image)
Definition: compare.c:109
MagickExport const char * GetImageArtifact(const Image *image, const char *artifact)
Definition: artifact.c:273
static double StringToDouble(const char *magick_restrict string, char **magick_restrict sentinal)
double mean_error_per_pixel
Definition: color.h:63
static PixelTrait GetPixelChannelTraits(const Image *magick_restrict image, const PixelChannel channel)
MetricType
Definition: compare.h:27
MagickExport ssize_t FormatLocaleString(char *magick_restrict string, const size_t length, const char *magick_restrict format,...)
Definition: locale.c:499
static void SetPixelViaPixelInfo(const Image *magick_restrict image, const PixelInfo *magick_restrict pixel_info, Quantum *magick_restrict pixel)
static Quantum GetPixelReadMask(const Image *magick_restrict image, const Quantum *magick_restrict pixel)
MagickExport const Quantum * GetCacheViewVirtualPixels(const CacheView *cache_view, const ssize_t x, const ssize_t y, const size_t columns, const size_t rows, ExceptionInfo *exception)
Definition: cache-view.c:651
#define MaximumNumberOfImageMoments
Definition: statistic.h:25
MagickExport MagickBooleanType SetImageColorMetric(Image *image, const Image *reconstruct_image, ExceptionInfo *exception)
Definition: compare.c:2004
#define MagickEpsilon
Definition: magick-type.h:110
MagickExport MagickBooleanType CompositeImage(Image *image, const Image *composite, const CompositeOperator compose, const MagickBooleanType clip_to_self, const ssize_t x_offset, const ssize_t y_offset, ExceptionInfo *exception)
Definition: composite.c:528
size_t width
Definition: geometry.h:130
static double GetFuzzyColorDistance(const Image *p, const Image *q)
Definition: color-private.h:36
#define ThrowBinaryException(severity, tag, context)
Definition: log.h:52
ssize_t MagickOffsetType
Definition: magick-type.h:129
Definition: image.h:151
MagickExport Image * CompareImages(Image *image, const Image *reconstruct_image, const MetricType metric, double *distortion, ExceptionInfo *exception)
Definition: compare.c:128
MagickExport Image * CropImage(const Image *image, const RectangleInfo *geometry, ExceptionInfo *exception)
Definition: transform.c:535
MagickExport KernelInfo * AcquireKernelInfo(const char *kernel_string, ExceptionInfo *exception)
Definition: morphology.c:486
#define MagickCoreSignature
double normalized_mean_error
Definition: color.h:63
MagickExport Quantum * GetCacheViewAuthenticPixels(CacheView *cache_view, const ssize_t x, const ssize_t y, const size_t columns, const size_t rows, ExceptionInfo *exception)
Definition: cache-view.c:299
MagickExport MagickBooleanType SetImageMask(Image *image, const PixelMask type, const Image *mask, ExceptionInfo *exception)
Definition: image.c:3185
MagickExport MagickBooleanType SetImageAlphaChannel(Image *image, const AlphaChannelOption alpha_type, ExceptionInfo *exception)
Definition: channel.c:974
MagickBooleanType
Definition: magick-type.h:158
unsigned int MagickStatusType
Definition: magick-type.h:121
static double PerceptibleReciprocal(const double x)
MagickExport void * AcquireQuantumMemory(const size_t count, const size_t quantum)
Definition: memory.c:543
#define MagickPathExtent
static MagickBooleanType GetMeanErrorPerPixel(Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c:713
MagickExport MagickBooleanType IsStringTrue(const char *value)
Definition: string.c:1425
#define INFINITY
Definition: magick-type.h:184
MagickExport int GetMagickPrecision(void)
Definition: magick.c:886
#define MagickMaximumValue
Definition: magick-type.h:111
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:1413
size_t width
Definition: morphology.h:108
size_t signature
Definition: image.h:354
#define QuantumScale
Definition: magick-type.h:115
size_t columns
Definition: image.h:172
ssize_t x
Definition: geometry.h:134
size_t height
Definition: geometry.h:130
MagickExport MagickBooleanType QueryColorCompliance(const char *name, const ComplianceType compliance, PixelInfo *color, ExceptionInfo *exception)
Definition: color.c:2216
MagickExport MagickBooleanType SetImageStorageClass(Image *image, const ClassType storage_class, ExceptionInfo *exception)
Definition: image.c:2617
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:67
#define MagickMax(x, y)
Definition: image-private.h:26
static MagickBooleanType GetRootMeanSquaredDistortion(const Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c: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:1884
#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:2135
static Quantum ClampToQuantum(const MagickRealType value)
Definition: quantum.h:84
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:1592
unsigned short Quantum
Definition: magick-type.h:82
#define SSIML
#define Log10Epsilon
MagickExport Image * ExtentImage(const Image *image, const RectangleInfo *geometry, ExceptionInfo *exception)
Definition: transform.c:1115
#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:27
MagickExport void SetGeometry(const Image *image, RectangleInfo *geometry)
Definition: geometry.c:1686
MagickExport void * RelinquishMagickMemory(void *memory)
Definition: memory.c:1069
#define MaxPixelChannels
Definition: pixel.h:27
MagickBooleanType(* MagickProgressMonitor)(const char *, const MagickOffsetType, const MagickSizeType, void *)
Definition: monitor.h:26
CompositeOperator compose
Definition: image.h:234
MagickExport double * GetImageDistortions(Image *image, const Image *reconstruct_image, const MetricType metric, ExceptionInfo *exception)
Definition: compare.c:1740
#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:1230
static MagickBooleanType GetFuzzDistortion(const Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c:484
PixelTrait
Definition: pixel.h:134
MagickExport ChannelPerceptualHash * GetImagePerceptualHash(const Image *image, ExceptionInfo *exception)
Definition: statistic.c:1687
MagickExport Image * DestroyImage(Image *image)
Definition: image.c:1181
MagickExport Image * CloneImage(const Image *image, const size_t columns, const size_t rows, const MagickBooleanType detach, ExceptionInfo *exception)
Definition: image.c:796
#define QuantumRange
Definition: magick-type.h:83
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