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