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-2018 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://www.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  proceed=SetImageProgress(image,SimilarityImageTag,progress++,rows);
1076  if (proceed == MagickFalse)
1077  {
1078  status=MagickFalse;
1079  break;
1080  }
1081  }
1082  }
1083  reconstruct_view=DestroyCacheView(reconstruct_view);
1084  image_view=DestroyCacheView(image_view);
1085  /*
1086  Divide by the standard deviation.
1087  */
1088  distortion[CompositePixelChannel]=0.0;
1089  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1090  {
1091  double
1092  gamma;
1093 
1094  PixelChannel channel = GetPixelChannelChannel(image,i);
1095  gamma=image_statistics[channel].standard_deviation*
1096  reconstruct_statistics[channel].standard_deviation;
1097  gamma=PerceptibleReciprocal(gamma);
1098  distortion[i]=QuantumRange*gamma*distortion[i];
1099  distortion[CompositePixelChannel]+=distortion[i]*distortion[i];
1100  }
1101  distortion[CompositePixelChannel]=sqrt(distortion[CompositePixelChannel]/
1102  GetImageChannels(image));
1103  /*
1104  Free resources.
1105  */
1106  reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1107  reconstruct_statistics);
1108  image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1109  image_statistics);
1110  return(status);
1111 }
1112 
1114  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1115 {
1116  CacheView
1117  *image_view,
1118  *reconstruct_view;
1119 
1121  status;
1122 
1123  size_t
1124  columns,
1125  rows;
1126 
1127  ssize_t
1128  y;
1129 
1130  status=MagickTrue;
1131  rows=MagickMax(image->rows,reconstruct_image->rows);
1132  columns=MagickMax(image->columns,reconstruct_image->columns);
1133  image_view=AcquireVirtualCacheView(image,exception);
1134  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1135 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1136  #pragma omp parallel for schedule(static) shared(status) \
1137  magick_number_threads(image,image,rows,1)
1138 #endif
1139  for (y=0; y < (ssize_t) rows; y++)
1140  {
1141  double
1142  channel_distortion[MaxPixelChannels+1];
1143 
1144  register const Quantum
1145  *magick_restrict p,
1146  *magick_restrict q;
1147 
1148  register ssize_t
1149  j,
1150  x;
1151 
1152  if (status == MagickFalse)
1153  continue;
1154  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1155  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1156  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1157  {
1158  status=MagickFalse;
1159  continue;
1160  }
1161  (void) memset(channel_distortion,0,sizeof(channel_distortion));
1162  for (x=0; x < (ssize_t) columns; x++)
1163  {
1164  double
1165  Da,
1166  Sa;
1167 
1168  register ssize_t
1169  i;
1170 
1171  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1172  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1173  {
1174  p+=GetPixelChannels(image);
1175  q+=GetPixelChannels(reconstruct_image);
1176  continue;
1177  }
1178  Sa=QuantumScale*GetPixelAlpha(image,p);
1179  Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
1180  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1181  {
1182  double
1183  distance;
1184 
1185  PixelChannel channel = GetPixelChannelChannel(image,i);
1186  PixelTrait traits = GetPixelChannelTraits(image,channel);
1187  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1188  channel);
1189  if ((traits == UndefinedPixelTrait) ||
1190  (reconstruct_traits == UndefinedPixelTrait) ||
1191  ((reconstruct_traits & UpdatePixelTrait) == 0))
1192  continue;
1193  if (channel == AlphaPixelChannel)
1194  distance=QuantumScale*fabs((double) p[i]-
1195  GetPixelChannel(reconstruct_image,channel,q));
1196  else
1197  distance=QuantumScale*fabs(Sa*p[i]-Da*
1198  GetPixelChannel(reconstruct_image,channel,q));
1199  if (distance > channel_distortion[i])
1200  channel_distortion[i]=distance;
1201  if (distance > channel_distortion[CompositePixelChannel])
1202  channel_distortion[CompositePixelChannel]=distance;
1203  }
1204  p+=GetPixelChannels(image);
1205  q+=GetPixelChannels(reconstruct_image);
1206  }
1207 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1208  #pragma omp critical (MagickCore_GetPeakAbsoluteError)
1209 #endif
1210  for (j=0; j <= MaxPixelChannels; j++)
1211  if (channel_distortion[j] > distortion[j])
1212  distortion[j]=channel_distortion[j];
1213  }
1214  reconstruct_view=DestroyCacheView(reconstruct_view);
1215  image_view=DestroyCacheView(image_view);
1216  return(status);
1217 }
1218 
1219 static inline double MagickLog10(const double x)
1220 {
1221 #define Log10Epsilon (1.0e-11)
1222 
1223  if (fabs(x) < Log10Epsilon)
1224  return(log10(Log10Epsilon));
1225  return(log10(fabs(x)));
1226 }
1227 
1229  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1230 {
1232  status;
1233 
1234  register ssize_t
1235  i;
1236 
1237  status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1238  for (i=0; i <= MaxPixelChannels; i++)
1239  if (fabs(distortion[i]) < MagickEpsilon)
1240  distortion[i]=INFINITY;
1241  else
1242  distortion[i]=20.0*MagickLog10(1.0/sqrt(distortion[i]));
1243  return(status);
1244 }
1245 
1247  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1248 {
1250  *channel_phash,
1251  *reconstruct_phash;
1252 
1253  const char
1254  *artifact;
1255 
1257  normalize;
1258 
1259  ssize_t
1260  channel;
1261 
1262  /*
1263  Compute perceptual hash in the sRGB colorspace.
1264  */
1265  channel_phash=GetImagePerceptualHash(image,exception);
1266  if (channel_phash == (ChannelPerceptualHash *) NULL)
1267  return(MagickFalse);
1268  reconstruct_phash=GetImagePerceptualHash(reconstruct_image,exception);
1269  if (reconstruct_phash == (ChannelPerceptualHash *) NULL)
1270  {
1272  channel_phash);
1273  return(MagickFalse);
1274  }
1275  artifact=GetImageArtifact(image,"phash:normalize");
1276  normalize=(artifact == (const char *) NULL) ||
1277  (IsStringTrue(artifact) == MagickFalse) ? MagickFalse : MagickTrue;
1278 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1279  #pragma omp parallel for schedule(static)
1280 #endif
1281  for (channel=0; channel < MaxPixelChannels; channel++)
1282  {
1283  double
1284  difference;
1285 
1286  register ssize_t
1287  i;
1288 
1289  difference=0.0;
1290  for (i=0; i < MaximumNumberOfImageMoments; i++)
1291  {
1292  double
1293  alpha,
1294  beta;
1295 
1296  register ssize_t
1297  j;
1298 
1299  for (j=0; j < (ssize_t) channel_phash[0].number_colorspaces; j++)
1300  {
1301  alpha=channel_phash[channel].phash[j][i];
1302  beta=reconstruct_phash[channel].phash[j][i];
1303  if (normalize == MagickFalse)
1304  difference+=(beta-alpha)*(beta-alpha);
1305  else
1306  difference=sqrt((beta-alpha)*(beta-alpha)/
1307  channel_phash[0].number_channels);
1308  }
1309  }
1310  distortion[channel]+=difference;
1311 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1312  #pragma omp critical (MagickCore_GetPerceptualHashDistortion)
1313 #endif
1314  distortion[CompositePixelChannel]+=difference;
1315  }
1316  /*
1317  Free resources.
1318  */
1319  reconstruct_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1320  reconstruct_phash);
1321  channel_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(channel_phash);
1322  return(MagickTrue);
1323 }
1324 
1326  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1327 {
1329  status;
1330 
1331  register ssize_t
1332  i;
1333 
1334  status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1335  for (i=0; i <= MaxPixelChannels; i++)
1336  distortion[i]=sqrt(distortion[i]);
1337  return(status);
1338 }
1339 
1341  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1342 {
1343 #define SSIMRadius 5.0
1344 #define SSIMSigma 1.5
1345 #define SSIMBlocksize 8
1346 #define SSIMK1 0.01
1347 #define SSIMK2 0.03
1348 #define SSIML 1.0
1349 
1350  CacheView
1351  *image_view,
1352  *reconstruct_view;
1353 
1354  char
1355  geometry[MagickPathExtent];
1356 
1357  const char
1358  *artifact;
1359 
1360  double
1361  c1,
1362  c2,
1363  radius,
1364  sigma;
1365 
1366  KernelInfo
1367  *kernel_info;
1368 
1370  status;
1371 
1372  register ssize_t
1373  i;
1374 
1375  size_t
1376  columns,
1377  rows;
1378 
1379  ssize_t
1380  y;
1381 
1382  /*
1383  Compute structural similarity index @
1384  https://en.wikipedia.org/wiki/Structural_similarity.
1385  */
1386  radius=SSIMRadius;
1387  artifact=GetImageArtifact(image,"compare:ssim-radius");
1388  if (artifact != (const char *) NULL)
1389  radius=StringToDouble(artifact,(char **) NULL);
1390  sigma=SSIMSigma;
1391  artifact=GetImageArtifact(image,"compare:ssim-sigma");
1392  if (artifact != (const char *) NULL)
1393  sigma=StringToDouble(artifact,(char **) NULL);
1394  (void) FormatLocaleString(geometry,MagickPathExtent,"gaussian:%.20gx%.20g",
1395  radius,sigma);
1396  kernel_info=AcquireKernelInfo(geometry,exception);
1397  if (kernel_info == (KernelInfo *) NULL)
1398  ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
1399  image->filename);
1400  c1=pow(SSIMK1*SSIML,2.0);
1401  artifact=GetImageArtifact(image,"compare:ssim-k1");
1402  if (artifact != (const char *) NULL)
1403  c1=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
1404  c2=pow(SSIMK2*SSIML,2.0);
1405  artifact=GetImageArtifact(image,"compare:ssim-k2");
1406  if (artifact != (const char *) NULL)
1407  c2=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
1408  status=MagickTrue;
1409  rows=MagickMax(image->rows,reconstruct_image->rows);
1410  columns=MagickMax(image->columns,reconstruct_image->columns);
1411  image_view=AcquireVirtualCacheView(image,exception);
1412  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1413 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1414  #pragma omp parallel for schedule(static) shared(status) \
1415  magick_number_threads(image,reconstruct_image,rows,1)
1416 #endif
1417  for (y=0; y < (ssize_t) rows; y++)
1418  {
1419  double
1420  channel_distortion[MaxPixelChannels+1];
1421 
1422  register const Quantum
1423  *magick_restrict p,
1424  *magick_restrict q;
1425 
1426  register ssize_t
1427  i,
1428  x;
1429 
1430  if (status == MagickFalse)
1431  continue;
1432  p=GetCacheViewVirtualPixels(image_view,-((ssize_t) kernel_info->width/2L),y-
1433  ((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
1434  kernel_info->height,exception);
1435  q=GetCacheViewVirtualPixels(reconstruct_view,-((ssize_t) kernel_info->width/
1436  2L),y-((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
1437  kernel_info->height,exception);
1438  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1439  {
1440  status=MagickFalse;
1441  continue;
1442  }
1443  (void) memset(channel_distortion,0,sizeof(channel_distortion));
1444  for (x=0; x < (ssize_t) columns; x++)
1445  {
1446  double
1447  x_pixel_mu[MaxPixelChannels+1],
1448  x_pixel_sigma_squared[MaxPixelChannels+1],
1449  xy_sigma[MaxPixelChannels+1],
1450  y_pixel_mu[MaxPixelChannels+1],
1451  y_pixel_sigma_squared[MaxPixelChannels+1];
1452 
1453  register const Quantum
1454  *magick_restrict reference,
1455  *magick_restrict target;
1456 
1457  register double
1458  *k;
1459 
1460  ssize_t
1461  v;
1462 
1463  (void) memset(x_pixel_mu,0,sizeof(x_pixel_mu));
1464  (void) memset(x_pixel_sigma_squared,0,
1465  sizeof(x_pixel_sigma_squared));
1466  (void) memset(xy_sigma,0,sizeof(xy_sigma));
1467  (void) memset(x_pixel_sigma_squared,0,
1468  sizeof(y_pixel_sigma_squared));
1469  (void) memset(y_pixel_mu,0,sizeof(y_pixel_mu));
1470  (void) memset(y_pixel_sigma_squared,0,
1471  sizeof(y_pixel_sigma_squared));
1472  k=kernel_info->values;
1473  reference=p;
1474  target=q;
1475  for (v=0; v < (ssize_t) kernel_info->height; v++)
1476  {
1477  register ssize_t
1478  u;
1479 
1480  for (u=0; u < (ssize_t) kernel_info->width; u++)
1481  {
1482  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1483  {
1484  double
1485  x_pixel,
1486  y_pixel;
1487 
1488  PixelChannel channel = GetPixelChannelChannel(image,i);
1489  PixelTrait traits = GetPixelChannelTraits(image,channel);
1490  PixelTrait reconstruct_traits = GetPixelChannelTraits(
1491  reconstruct_image,channel);
1492  if ((traits == UndefinedPixelTrait) ||
1493  (reconstruct_traits == UndefinedPixelTrait) ||
1494  ((reconstruct_traits & UpdatePixelTrait) == 0))
1495  continue;
1496  x_pixel=QuantumScale*reference[i];
1497  x_pixel_mu[i]+=(*k)*x_pixel;
1498  x_pixel_sigma_squared[i]+=(*k)*x_pixel*x_pixel;
1499  y_pixel=QuantumScale*
1500  GetPixelChannel(reconstruct_image,channel,target);
1501  y_pixel_mu[i]+=(*k)*y_pixel;
1502  y_pixel_sigma_squared[i]+=(*k)*y_pixel*y_pixel;
1503  xy_sigma[i]+=(*k)*x_pixel*y_pixel;
1504  }
1505  k++;
1506  reference+=GetPixelChannels(image);
1507  target+=GetPixelChannels(reconstruct_image);
1508  }
1509  reference+=GetPixelChannels(image)*columns;
1510  target+=GetPixelChannels(reconstruct_image)*columns;
1511  }
1512  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1513  {
1514  double
1515  ssim,
1516  x_pixel_mu_squared,
1517  x_pixel_sigmas_squared,
1518  xy_mu,
1519  xy_sigmas,
1520  y_pixel_mu_squared,
1521  y_pixel_sigmas_squared;
1522 
1523  PixelChannel channel = GetPixelChannelChannel(image,i);
1524  PixelTrait traits = GetPixelChannelTraits(image,channel);
1525  PixelTrait reconstruct_traits = GetPixelChannelTraits(
1526  reconstruct_image,channel);
1527  if ((traits == UndefinedPixelTrait) ||
1528  (reconstruct_traits == UndefinedPixelTrait) ||
1529  ((reconstruct_traits & UpdatePixelTrait) == 0))
1530  continue;
1531  x_pixel_mu_squared=x_pixel_mu[i]*x_pixel_mu[i];
1532  y_pixel_mu_squared=y_pixel_mu[i]*y_pixel_mu[i];
1533  xy_mu=x_pixel_mu[i]*y_pixel_mu[i];
1534  xy_sigmas=xy_sigma[i]-xy_mu;
1535  x_pixel_sigmas_squared=x_pixel_sigma_squared[i]-x_pixel_mu_squared;
1536  y_pixel_sigmas_squared=y_pixel_sigma_squared[i]-y_pixel_mu_squared;
1537  ssim=((2.0*xy_mu+c1)*(2.0*xy_sigmas+c2))/
1538  ((x_pixel_mu_squared+y_pixel_mu_squared+c1)*
1539  (x_pixel_sigmas_squared+y_pixel_sigmas_squared+c2));
1540  channel_distortion[i]+=ssim;
1541  channel_distortion[CompositePixelChannel]+=ssim;
1542  }
1543  p+=GetPixelChannels(image);
1544  q+=GetPixelChannels(reconstruct_image);
1545  }
1546 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1547  #pragma omp critical (MagickCore_GetStructuralSimilarityDistortion)
1548 #endif
1549  for (i=0; i <= MaxPixelChannels; i++)
1550  distortion[i]+=channel_distortion[i];
1551  }
1552  image_view=DestroyCacheView(image_view);
1553  reconstruct_view=DestroyCacheView(reconstruct_view);
1554  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1555  {
1556  PixelChannel channel = GetPixelChannelChannel(image,i);
1557  PixelTrait traits = GetPixelChannelTraits(image,channel);
1558  if ((traits == UndefinedPixelTrait) || ((traits & UpdatePixelTrait) == 0))
1559  continue;
1560  distortion[i]/=((double) columns*rows);
1561  }
1562  distortion[CompositePixelChannel]/=((double) columns*rows);
1563  distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
1564  kernel_info=DestroyKernelInfo(kernel_info);
1565  return(status);
1566 }
1567 
1569  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1570 {
1572  status;
1573 
1574  register ssize_t
1575  i;
1576 
1577  status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1578  distortion,exception);
1579  for (i=0; i <= MaxPixelChannels; i++)
1580  distortion[i]=(1.0-(distortion[i]))/2.0;
1581  return(status);
1582 }
1583 
1585  const Image *reconstruct_image,const MetricType metric,double *distortion,
1586  ExceptionInfo *exception)
1587 {
1588  double
1589  *channel_distortion;
1590 
1592  status;
1593 
1594  size_t
1595  length;
1596 
1597  assert(image != (Image *) NULL);
1598  assert(image->signature == MagickCoreSignature);
1599  if (image->debug != MagickFalse)
1600  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1601  assert(reconstruct_image != (const Image *) NULL);
1602  assert(reconstruct_image->signature == MagickCoreSignature);
1603  assert(distortion != (double *) NULL);
1604  *distortion=0.0;
1605  if (image->debug != MagickFalse)
1606  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1607  /*
1608  Get image distortion.
1609  */
1610  length=MaxPixelChannels+1;
1611  channel_distortion=(double *) AcquireQuantumMemory(length,
1612  sizeof(*channel_distortion));
1613  if (channel_distortion == (double *) NULL)
1614  ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1615  (void) memset(channel_distortion,0,length*
1616  sizeof(*channel_distortion));
1617  switch (metric)
1618  {
1619  case AbsoluteErrorMetric:
1620  {
1621  status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1622  exception);
1623  break;
1624  }
1625  case FuzzErrorMetric:
1626  {
1627  status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1628  exception);
1629  break;
1630  }
1632  {
1633  status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1634  channel_distortion,exception);
1635  break;
1636  }
1638  {
1639  status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1640  exception);
1641  break;
1642  }
1644  {
1645  status=GetMeanSquaredDistortion(image,reconstruct_image,
1646  channel_distortion,exception);
1647  break;
1648  }
1650  default:
1651  {
1652  status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1653  channel_distortion,exception);
1654  break;
1655  }
1657  {
1658  status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1659  channel_distortion,exception);
1660  break;
1661  }
1663  {
1664  status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1665  channel_distortion,exception);
1666  break;
1667  }
1669  {
1670  status=GetPerceptualHashDistortion(image,reconstruct_image,
1671  channel_distortion,exception);
1672  break;
1673  }
1675  {
1676  status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1677  channel_distortion,exception);
1678  break;
1679  }
1681  {
1682  status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1683  channel_distortion,exception);
1684  break;
1685  }
1687  {
1688  status=GetStructuralDisimilarityDistortion(image,reconstruct_image,
1689  channel_distortion,exception);
1690  break;
1691  }
1692  }
1693  *distortion=channel_distortion[CompositePixelChannel];
1694  channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1695  (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
1696  *distortion);
1697  return(status);
1698 }
1699 
1700 /*
1701 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1702 % %
1703 % %
1704 % %
1705 % G e t I m a g e D i s t o r t i o n s %
1706 % %
1707 % %
1708 % %
1709 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1710 %
1711 % GetImageDistortions() compares the pixel channels of an image to a
1712 % reconstructed image and returns the specified distortion metric for each
1713 % channel.
1714 %
1715 % The format of the GetImageDistortions method is:
1716 %
1717 % double *GetImageDistortions(const Image *image,
1718 % const Image *reconstruct_image,const MetricType metric,
1719 % ExceptionInfo *exception)
1720 %
1721 % A description of each parameter follows:
1722 %
1723 % o image: the image.
1724 %
1725 % o reconstruct_image: the reconstruct image.
1726 %
1727 % o metric: the metric.
1728 %
1729 % o exception: return any errors or warnings in this structure.
1730 %
1731 */
1733  const Image *reconstruct_image,const MetricType metric,
1734  ExceptionInfo *exception)
1735 {
1736  double
1737  *channel_distortion;
1738 
1740  status;
1741 
1742  size_t
1743  length;
1744 
1745  assert(image != (Image *) NULL);
1746  assert(image->signature == MagickCoreSignature);
1747  if (image->debug != MagickFalse)
1748  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1749  assert(reconstruct_image != (const Image *) NULL);
1750  assert(reconstruct_image->signature == MagickCoreSignature);
1751  if (image->debug != MagickFalse)
1752  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1753  /*
1754  Get image distortion.
1755  */
1756  length=MaxPixelChannels+1UL;
1757  channel_distortion=(double *) AcquireQuantumMemory(length,
1758  sizeof(*channel_distortion));
1759  if (channel_distortion == (double *) NULL)
1760  ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1761  (void) memset(channel_distortion,0,length*
1762  sizeof(*channel_distortion));
1763  status=MagickTrue;
1764  switch (metric)
1765  {
1766  case AbsoluteErrorMetric:
1767  {
1768  status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1769  exception);
1770  break;
1771  }
1772  case FuzzErrorMetric:
1773  {
1774  status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1775  exception);
1776  break;
1777  }
1779  {
1780  status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1781  channel_distortion,exception);
1782  break;
1783  }
1785  {
1786  status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1787  exception);
1788  break;
1789  }
1791  {
1792  status=GetMeanSquaredDistortion(image,reconstruct_image,
1793  channel_distortion,exception);
1794  break;
1795  }
1797  default:
1798  {
1799  status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1800  channel_distortion,exception);
1801  break;
1802  }
1804  {
1805  status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1806  channel_distortion,exception);
1807  break;
1808  }
1810  {
1811  status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1812  channel_distortion,exception);
1813  break;
1814  }
1816  {
1817  status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1818  channel_distortion,exception);
1819  break;
1820  }
1822  {
1823  status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1824  channel_distortion,exception);
1825  break;
1826  }
1828  {
1829  status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1830  channel_distortion,exception);
1831  break;
1832  }
1834  {
1835  status=GetStructuralDisimilarityDistortion(image,reconstruct_image,
1836  channel_distortion,exception);
1837  break;
1838  }
1839  }
1840  if (status == MagickFalse)
1841  {
1842  channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1843  return((double *) NULL);
1844  }
1845  return(channel_distortion);
1846 }
1847 
1848 /*
1849 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1850 % %
1851 % %
1852 % %
1853 % I s I m a g e s E q u a l %
1854 % %
1855 % %
1856 % %
1857 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1858 %
1859 % IsImagesEqual() compare the pixels of two images and returns immediately
1860 % if any pixel is not identical.
1861 %
1862 % The format of the IsImagesEqual method is:
1863 %
1864 % MagickBooleanType IsImagesEqual(const Image *image,
1865 % const Image *reconstruct_image,ExceptionInfo *exception)
1866 %
1867 % A description of each parameter follows.
1868 %
1869 % o image: the image.
1870 %
1871 % o reconstruct_image: the reconstruct image.
1872 %
1873 % o exception: return any errors or warnings in this structure.
1874 %
1875 */
1877  const Image *reconstruct_image,ExceptionInfo *exception)
1878 {
1879  CacheView
1880  *image_view,
1881  *reconstruct_view;
1882 
1883  size_t
1884  columns,
1885  rows;
1886 
1887  ssize_t
1888  y;
1889 
1890  assert(image != (Image *) NULL);
1891  assert(image->signature == MagickCoreSignature);
1892  assert(reconstruct_image != (const Image *) NULL);
1893  assert(reconstruct_image->signature == MagickCoreSignature);
1894  rows=MagickMax(image->rows,reconstruct_image->rows);
1895  columns=MagickMax(image->columns,reconstruct_image->columns);
1896  image_view=AcquireVirtualCacheView(image,exception);
1897  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1898  for (y=0; y < (ssize_t) rows; y++)
1899  {
1900  register const Quantum
1901  *magick_restrict p,
1902  *magick_restrict q;
1903 
1904  register ssize_t
1905  x;
1906 
1907  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1908  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1909  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
1910  break;
1911  for (x=0; x < (ssize_t) columns; x++)
1912  {
1913  register ssize_t
1914  i;
1915 
1916  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1917  {
1918  double
1919  distance;
1920 
1921  PixelChannel channel = GetPixelChannelChannel(image,i);
1922  PixelTrait traits = GetPixelChannelTraits(image,channel);
1923  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1924  channel);
1925  if ((traits == UndefinedPixelTrait) ||
1926  (reconstruct_traits == UndefinedPixelTrait) ||
1927  ((reconstruct_traits & UpdatePixelTrait) == 0))
1928  continue;
1929  distance=fabs(p[i]-(double) GetPixelChannel(reconstruct_image,
1930  channel,q));
1931  if (distance >= MagickEpsilon)
1932  break;
1933  }
1934  if (i < (ssize_t) GetPixelChannels(image))
1935  break;
1936  p+=GetPixelChannels(image);
1937  q+=GetPixelChannels(reconstruct_image);
1938  }
1939  if (x < (ssize_t) columns)
1940  break;
1941  }
1942  reconstruct_view=DestroyCacheView(reconstruct_view);
1943  image_view=DestroyCacheView(image_view);
1944  return(y < (ssize_t) rows ? MagickFalse : MagickTrue);
1945 }
1946 
1947 /*
1948 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1949 % %
1950 % %
1951 % %
1952 % S e t I m a g e C o l o r M e t r i c %
1953 % %
1954 % %
1955 % %
1956 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1957 %
1958 % SetImageColorMetric() measures the difference between colors at each pixel
1959 % location of two images. A value other than 0 means the colors match
1960 % exactly. Otherwise an error measure is computed by summing over all
1961 % pixels in an image the distance squared in RGB space between each image
1962 % pixel and its corresponding pixel in the reconstruct image. The error
1963 % measure is assigned to these image members:
1964 %
1965 % o mean_error_per_pixel: The mean error for any single pixel in
1966 % the image.
1967 %
1968 % o normalized_mean_error: The normalized mean quantization error for
1969 % any single pixel in the image. This distance measure is normalized to
1970 % a range between 0 and 1. It is independent of the range of red, green,
1971 % and blue values in the image.
1972 %
1973 % o normalized_maximum_error: The normalized maximum quantization
1974 % error for any single pixel in the image. This distance measure is
1975 % normalized to a range between 0 and 1. It is independent of the range
1976 % of red, green, and blue values in your image.
1977 %
1978 % A small normalized mean square error, accessed as
1979 % image->normalized_mean_error, suggests the images are very similar in
1980 % spatial layout and color.
1981 %
1982 % The format of the SetImageColorMetric method is:
1983 %
1984 % MagickBooleanType SetImageColorMetric(Image *image,
1985 % const Image *reconstruct_image,ExceptionInfo *exception)
1986 %
1987 % A description of each parameter follows.
1988 %
1989 % o image: the image.
1990 %
1991 % o reconstruct_image: the reconstruct image.
1992 %
1993 % o exception: return any errors or warnings in this structure.
1994 %
1995 */
1997  const Image *reconstruct_image,ExceptionInfo *exception)
1998 {
1999  CacheView
2000  *image_view,
2001  *reconstruct_view;
2002 
2003  double
2004  area,
2005  maximum_error,
2006  mean_error,
2007  mean_error_per_pixel;
2008 
2010  status;
2011 
2012  size_t
2013  columns,
2014  rows;
2015 
2016  ssize_t
2017  y;
2018 
2019  assert(image != (Image *) NULL);
2020  assert(image->signature == MagickCoreSignature);
2021  assert(reconstruct_image != (const Image *) NULL);
2022  assert(reconstruct_image->signature == MagickCoreSignature);
2023  area=0.0;
2024  maximum_error=0.0;
2025  mean_error_per_pixel=0.0;
2026  mean_error=0.0;
2027  rows=MagickMax(image->rows,reconstruct_image->rows);
2028  columns=MagickMax(image->columns,reconstruct_image->columns);
2029  image_view=AcquireVirtualCacheView(image,exception);
2030  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
2031  for (y=0; y < (ssize_t) rows; y++)
2032  {
2033  register const Quantum
2034  *magick_restrict p,
2035  *magick_restrict q;
2036 
2037  register ssize_t
2038  x;
2039 
2040  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
2041  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
2042  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2043  break;
2044  for (x=0; x < (ssize_t) columns; x++)
2045  {
2046  register ssize_t
2047  i;
2048 
2049  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2050  {
2051  double
2052  distance;
2053 
2054  PixelChannel channel = GetPixelChannelChannel(image,i);
2055  PixelTrait traits = GetPixelChannelTraits(image,channel);
2056  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2057  channel);
2058  if ((traits == UndefinedPixelTrait) ||
2059  (reconstruct_traits == UndefinedPixelTrait) ||
2060  ((reconstruct_traits & UpdatePixelTrait) == 0))
2061  continue;
2062  distance=fabs(p[i]-(double) GetPixelChannel(reconstruct_image,
2063  channel,q));
2064  if (distance >= MagickEpsilon)
2065  {
2066  mean_error_per_pixel+=distance;
2067  mean_error+=distance*distance;
2068  if (distance > maximum_error)
2069  maximum_error=distance;
2070  }
2071  area++;
2072  }
2073  p+=GetPixelChannels(image);
2074  q+=GetPixelChannels(reconstruct_image);
2075  }
2076  }
2077  reconstruct_view=DestroyCacheView(reconstruct_view);
2078  image_view=DestroyCacheView(image_view);
2079  image->error.mean_error_per_pixel=(double) (mean_error_per_pixel/area);
2081  mean_error/area);
2082  image->error.normalized_maximum_error=(double) (QuantumScale*maximum_error);
2083  status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
2084  return(status);
2085 }
2086 
2087 /*
2088 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2089 % %
2090 % %
2091 % %
2092 % S i m i l a r i t y I m a g e %
2093 % %
2094 % %
2095 % %
2096 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2097 %
2098 % SimilarityImage() compares the reference image of the image and returns the
2099 % best match offset. In addition, it returns a similarity image such that an
2100 % exact match location is completely white and if none of the pixels match,
2101 % black, otherwise some gray level in-between.
2102 %
2103 % The format of the SimilarityImageImage method is:
2104 %
2105 % Image *SimilarityImage(const Image *image,const Image *reference,
2106 % const MetricType metric,const double similarity_threshold,
2107 % RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
2108 %
2109 % A description of each parameter follows:
2110 %
2111 % o image: the image.
2112 %
2113 % o reference: find an area of the image that closely resembles this image.
2114 %
2115 % o metric: the metric.
2116 %
2117 % o similarity_threshold: minimum distortion for (sub)image match.
2118 %
2119 % o offset: the best match offset of the reference image within the image.
2120 %
2121 % o similarity: the computed similarity between the images.
2122 %
2123 % o exception: return any errors or warnings in this structure.
2124 %
2125 */
2126 
2127 static double GetSimilarityMetric(const Image *image,const Image *reference,
2128  const MetricType metric,const ssize_t x_offset,const ssize_t y_offset,
2129  ExceptionInfo *exception)
2130 {
2131  double
2132  distortion;
2133 
2134  Image
2135  *similarity_image;
2136 
2138  status;
2139 
2141  geometry;
2142 
2143  SetGeometry(reference,&geometry);
2144  geometry.x=x_offset;
2145  geometry.y=y_offset;
2146  similarity_image=CropImage(image,&geometry,exception);
2147  if (similarity_image == (Image *) NULL)
2148  return(0.0);
2149  distortion=0.0;
2150  status=GetImageDistortion(similarity_image,reference,metric,&distortion,
2151  exception);
2152  similarity_image=DestroyImage(similarity_image);
2153  if (status == MagickFalse)
2154  return(0.0);
2155  return(distortion);
2156 }
2157 
2158 MagickExport Image *SimilarityImage(const Image *image,const Image *reference,
2159  const MetricType metric,const double similarity_threshold,
2160  RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
2161 {
2162 #define SimilarityImageTag "Similarity/Image"
2163 
2164  CacheView
2165  *similarity_view;
2166 
2167  Image
2168  *similarity_image;
2169 
2171  status;
2172 
2174  progress;
2175 
2176  ssize_t
2177  y;
2178 
2179  assert(image != (const Image *) NULL);
2180  assert(image->signature == MagickCoreSignature);
2181  if (image->debug != MagickFalse)
2182  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2183  assert(exception != (ExceptionInfo *) NULL);
2184  assert(exception->signature == MagickCoreSignature);
2185  assert(offset != (RectangleInfo *) NULL);
2186  SetGeometry(reference,offset);
2187  *similarity_metric=MagickMaximumValue;
2188  similarity_image=CloneImage(image,image->columns-reference->columns+1,
2189  image->rows-reference->rows+1,MagickTrue,exception);
2190  if (similarity_image == (Image *) NULL)
2191  return((Image *) NULL);
2192  status=SetImageStorageClass(similarity_image,DirectClass,exception);
2193  if (status == MagickFalse)
2194  {
2195  similarity_image=DestroyImage(similarity_image);
2196  return((Image *) NULL);
2197  }
2198  (void) SetImageAlphaChannel(similarity_image,DeactivateAlphaChannel,
2199  exception);
2200  /*
2201  Measure similarity of reference image against image.
2202  */
2203  status=MagickTrue;
2204  progress=0;
2205  similarity_view=AcquireAuthenticCacheView(similarity_image,exception);
2206 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2207  #pragma omp parallel for schedule(static) \
2208  shared(progress,status,similarity_metric) \
2209  magick_number_threads(image,image,image->rows-reference->rows+1,1)
2210 #endif
2211  for (y=0; y < (ssize_t) (image->rows-reference->rows+1); y++)
2212  {
2213  double
2214  similarity;
2215 
2216  register Quantum
2217  *magick_restrict q;
2218 
2219  register ssize_t
2220  x;
2221 
2222  if (status == MagickFalse)
2223  continue;
2224 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2225  #pragma omp flush(similarity_metric)
2226 #endif
2227  if (*similarity_metric <= similarity_threshold)
2228  continue;
2229  q=GetCacheViewAuthenticPixels(similarity_view,0,y,similarity_image->columns,
2230  1,exception);
2231  if (q == (Quantum *) NULL)
2232  {
2233  status=MagickFalse;
2234  continue;
2235  }
2236  for (x=0; x < (ssize_t) (image->columns-reference->columns+1); x++)
2237  {
2238  register ssize_t
2239  i;
2240 
2241 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2242  #pragma omp flush(similarity_metric)
2243 #endif
2244  if (*similarity_metric <= similarity_threshold)
2245  break;
2246  similarity=GetSimilarityMetric(image,reference,metric,x,y,exception);
2247 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2248  #pragma omp critical (MagickCore_SimilarityImage)
2249 #endif
2250  if ((metric == NormalizedCrossCorrelationErrorMetric) ||
2251  (metric == UndefinedErrorMetric))
2252  similarity=1.0-similarity;
2253  if (similarity < *similarity_metric)
2254  {
2255  offset->x=x;
2256  offset->y=y;
2257  *similarity_metric=similarity;
2258  }
2259  if (metric == PerceptualHashErrorMetric)
2260  similarity=MagickMin(0.01*similarity,1.0);
2261  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2262  {
2263  PixelChannel channel = GetPixelChannelChannel(image,i);
2264  PixelTrait traits = GetPixelChannelTraits(image,channel);
2265  PixelTrait similarity_traits=GetPixelChannelTraits(similarity_image,
2266  channel);
2267  if ((traits == UndefinedPixelTrait) ||
2268  (similarity_traits == UndefinedPixelTrait) ||
2269  ((similarity_traits & UpdatePixelTrait) == 0))
2270  continue;
2271  SetPixelChannel(similarity_image,channel,ClampToQuantum(QuantumRange-
2272  QuantumRange*similarity),q);
2273  }
2274  q+=GetPixelChannels(similarity_image);
2275  }
2276  if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
2277  status=MagickFalse;
2278  if (image->progress_monitor != (MagickProgressMonitor) NULL)
2279  {
2281  proceed;
2282 
2283 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2284  #pragma omp critical (MagickCore_SimilarityImage)
2285 #endif
2286  proceed=SetImageProgress(image,SimilarityImageTag,progress++,
2287  image->rows);
2288  if (proceed == MagickFalse)
2289  status=MagickFalse;
2290  }
2291  }
2292  similarity_view=DestroyCacheView(similarity_view);
2293  if (status == MagickFalse)
2294  similarity_image=DestroyImage(similarity_image);
2295  return(similarity_image);
2296 }
size_t rows
Definition: image.h:172
#define magick_restrict
Definition: MagickCore.h:41
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
static MagickBooleanType SetImageProgress(const Image *image, const char *tag, const MagickOffsetType offset, const MagickSizeType extent)
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:1228
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:2158
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:1914
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:1568
#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:1996
#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:127
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:533
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:3178
MagickExport MagickBooleanType SetImageAlphaChannel(Image *image, const AlphaChannelOption alpha_type, ExceptionInfo *exception)
Definition: channel.c:966
MagickBooleanType
Definition: magick-type.h:156
unsigned int MagickStatusType
Definition: magick-type.h:119
static double PerceptibleReciprocal(const double x)
MagickExport void * AcquireQuantumMemory(const size_t count, const size_t quantum)
Definition: memory.c:533
#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:182
MagickExport int GetMagickPrecision(void)
Definition: magick.c:877
#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:1397
size_t width
Definition: morphology.h:108
size_t signature
Definition: image.h:354
#define QuantumScale
Definition: magick-type.h:113
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:2613
static MagickBooleanType GetPeakAbsoluteDistortion(const Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c:1113
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:1325
static size_t GetPixelChannels(const Image *magick_restrict image)
MagickExport MagickBooleanType IsImagesEqual(const Image *image, const Image *reconstruct_image, ExceptionInfo *exception)
Definition: compare.c:1876
#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:2127
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:1584
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:1113
#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:1246
#define MagickMin(x, y)
Definition: image-private.h:27
MagickExport void SetGeometry(const Image *image, RectangleInfo *geometry)
Definition: geometry.c:1658
MagickExport void * RelinquishMagickMemory(void *memory)
Definition: memory.c:1053
#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:1732
#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:363
MagickExport CacheView * AcquireAuthenticCacheView(const Image *image, ExceptionInfo *exception)
Definition: cache-view.c:112
#define SSIMSigma
static double MagickLog10(const double x)
Definition: compare.c:1219
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:1671
MagickExport Image * DestroyImage(Image *image)
Definition: image.c:1177
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
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:1340
#define SSIMK2