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