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