MagickCore  7.1.0
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-2021 ImageMagick Studio LLC, a non-profit organization %
21 % dedicated to making software imaging solutions freely available. %
22 % %
23 % You may not use this file except in compliance with the License. You may %
24 % obtain a copy of the License at %
25 % %
26 % https://imagemagick.org/script/license.php %
27 % %
28 % Unless required by applicable law or agreed to in writing, software %
29 % distributed under the License is distributed on an "AS IS" BASIS, %
30 % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. %
31 % See the License for the specific language governing permissions and %
32 % limitations under the License. %
33 % %
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35 %
36 %
37 %
38 */
39 
40 /*
41  Include declarations.
42 */
43 #include "MagickCore/studio.h"
44 #include "MagickCore/artifact.h"
45 #include "MagickCore/attribute.h"
46 #include "MagickCore/cache-view.h"
47 #include "MagickCore/channel.h"
48 #include "MagickCore/client.h"
49 #include "MagickCore/color.h"
51 #include "MagickCore/colorspace.h"
53 #include "MagickCore/compare.h"
55 #include "MagickCore/constitute.h"
57 #include "MagickCore/enhance.h"
58 #include "MagickCore/fourier.h"
59 #include "MagickCore/geometry.h"
61 #include "MagickCore/list.h"
62 #include "MagickCore/log.h"
63 #include "MagickCore/memory_.h"
64 #include "MagickCore/monitor.h"
66 #include "MagickCore/option.h"
68 #include "MagickCore/property.h"
69 #include "MagickCore/registry.h"
70 #include "MagickCore/resource_.h"
71 #include "MagickCore/string_.h"
72 #include "MagickCore/statistic.h"
75 #include "MagickCore/transform.h"
76 #include "MagickCore/utility.h"
77 #include "MagickCore/version.h"
78 
79 /*
80 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
81 % %
82 % %
83 % %
84 % C o m p a r e I m a g e s %
85 % %
86 % %
87 % %
88 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
89 %
90 % CompareImages() compares one or more pixel channels of an image to a
91 % reconstructed image and returns the difference image.
92 %
93 % The format of the CompareImages method is:
94 %
95 % Image *CompareImages(const Image *image,const Image *reconstruct_image,
96 % const MetricType metric,double *distortion,ExceptionInfo *exception)
97 %
98 % A description of each parameter follows:
99 %
100 % o image: the image.
101 %
102 % o reconstruct_image: the reconstruct image.
103 %
104 % o metric: the metric.
105 %
106 % o distortion: the computed distortion between the images.
107 %
108 % o exception: return any errors or warnings in this structure.
109 %
110 */
111 
112 MagickExport Image *CompareImages(Image *image,const Image *reconstruct_image,
113  const MetricType metric,double *distortion,ExceptionInfo *exception)
114 {
115  CacheView
116  *highlight_view,
117  *image_view,
118  *reconstruct_view;
119 
120  const char
121  *artifact;
122 
123  double
124  fuzz;
125 
126  Image
127  *clone_image,
128  *difference_image,
129  *highlight_image;
130 
132  status;
133 
134  PixelInfo
135  highlight,
136  lowlight,
137  masklight;
138 
140  geometry;
141 
142  size_t
143  columns,
144  rows;
145 
146  ssize_t
147  y;
148 
149  assert(image != (Image *) NULL);
150  assert(image->signature == MagickCoreSignature);
151  if (image->debug != MagickFalse)
152  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
153  assert(reconstruct_image != (const Image *) NULL);
154  assert(reconstruct_image->signature == MagickCoreSignature);
155  assert(distortion != (double *) NULL);
156  *distortion=0.0;
157  if (image->debug != MagickFalse)
158  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
159  status=GetImageDistortion(image,reconstruct_image,metric,distortion,
160  exception);
161  if (status == MagickFalse)
162  return((Image *) NULL);
163  columns=MagickMax(image->columns,reconstruct_image->columns);
164  rows=MagickMax(image->rows,reconstruct_image->rows);
165  SetGeometry(image,&geometry);
166  geometry.width=columns;
167  geometry.height=rows;
168  clone_image=CloneImage(image,0,0,MagickTrue,exception);
169  if (clone_image == (Image *) NULL)
170  return((Image *) NULL);
171  (void) SetImageMask(clone_image,ReadPixelMask,(Image *) NULL,exception);
172  difference_image=ExtentImage(clone_image,&geometry,exception);
173  clone_image=DestroyImage(clone_image);
174  if (difference_image == (Image *) NULL)
175  return((Image *) NULL);
176  (void) SetImageAlphaChannel(difference_image,OpaqueAlphaChannel,exception);
177  highlight_image=CloneImage(image,columns,rows,MagickTrue,exception);
178  if (highlight_image == (Image *) NULL)
179  {
180  difference_image=DestroyImage(difference_image);
181  return((Image *) NULL);
182  }
183  status=SetImageStorageClass(highlight_image,DirectClass,exception);
184  if (status == MagickFalse)
185  {
186  difference_image=DestroyImage(difference_image);
187  highlight_image=DestroyImage(highlight_image);
188  return((Image *) NULL);
189  }
190  (void) SetImageMask(highlight_image,ReadPixelMask,(Image *) NULL,exception);
191  (void) SetImageAlphaChannel(highlight_image,OpaqueAlphaChannel,exception);
192  (void) QueryColorCompliance("#f1001ecc",AllCompliance,&highlight,exception);
193  artifact=GetImageArtifact(image,"compare:highlight-color");
194  if (artifact != (const char *) NULL)
195  (void) QueryColorCompliance(artifact,AllCompliance,&highlight,exception);
196  (void) QueryColorCompliance("#ffffffcc",AllCompliance,&lowlight,exception);
197  artifact=GetImageArtifact(image,"compare:lowlight-color");
198  if (artifact != (const char *) NULL)
199  (void) QueryColorCompliance(artifact,AllCompliance,&lowlight,exception);
200  (void) QueryColorCompliance("#888888cc",AllCompliance,&masklight,exception);
201  artifact=GetImageArtifact(image,"compare:masklight-color");
202  if (artifact != (const char *) NULL)
203  (void) QueryColorCompliance(artifact,AllCompliance,&masklight,exception);
204  /*
205  Generate difference image.
206  */
207  status=MagickTrue;
208  fuzz=GetFuzzyColorDistance(image,reconstruct_image);
209  image_view=AcquireVirtualCacheView(image,exception);
210  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
211  highlight_view=AcquireAuthenticCacheView(highlight_image,exception);
212 #if defined(MAGICKCORE_OPENMP_SUPPORT)
213  #pragma omp parallel for schedule(static) shared(status) \
214  magick_number_threads(image,highlight_image,rows,1)
215 #endif
216  for (y=0; y < (ssize_t) rows; y++)
217  {
219  sync;
220 
221  const Quantum
222  *magick_restrict p,
223  *magick_restrict q;
224 
225  Quantum
226  *magick_restrict r;
227 
228  ssize_t
229  x;
230 
231  if (status == MagickFalse)
232  continue;
233  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
234  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
235  r=QueueCacheViewAuthenticPixels(highlight_view,0,y,columns,1,exception);
236  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL) ||
237  (r == (Quantum *) NULL))
238  {
239  status=MagickFalse;
240  continue;
241  }
242  for (x=0; x < (ssize_t) columns; x++)
243  {
244  double
245  Da,
246  Sa;
247 
249  difference;
250 
251  ssize_t
252  i;
253 
254  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
255  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
256  {
257  SetPixelViaPixelInfo(highlight_image,&masklight,r);
258  p+=GetPixelChannels(image);
259  q+=GetPixelChannels(reconstruct_image);
260  r+=GetPixelChannels(highlight_image);
261  continue;
262  }
263  difference=MagickFalse;
264  Sa=QuantumScale*GetPixelAlpha(image,p);
265  Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
266  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
267  {
268  double
269  distance,
270  pixel;
271 
272  PixelChannel channel = GetPixelChannelChannel(image,i);
273  PixelTrait traits = GetPixelChannelTraits(image,channel);
274  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
275  channel);
276  if ((traits == UndefinedPixelTrait) ||
277  (reconstruct_traits == UndefinedPixelTrait) ||
278  ((reconstruct_traits & UpdatePixelTrait) == 0))
279  continue;
280  if (channel == AlphaPixelChannel)
281  pixel=(double) p[i]-GetPixelChannel(reconstruct_image,channel,q);
282  else
283  pixel=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
284  distance=pixel*pixel;
285  if (distance >= fuzz)
286  {
287  difference=MagickTrue;
288  break;
289  }
290  }
291  if (difference == MagickFalse)
292  SetPixelViaPixelInfo(highlight_image,&lowlight,r);
293  else
294  SetPixelViaPixelInfo(highlight_image,&highlight,r);
295  p+=GetPixelChannels(image);
296  q+=GetPixelChannels(reconstruct_image);
297  r+=GetPixelChannels(highlight_image);
298  }
299  sync=SyncCacheViewAuthenticPixels(highlight_view,exception);
300  if (sync == MagickFalse)
301  status=MagickFalse;
302  }
303  highlight_view=DestroyCacheView(highlight_view);
304  reconstruct_view=DestroyCacheView(reconstruct_view);
305  image_view=DestroyCacheView(image_view);
306  (void) CompositeImage(difference_image,highlight_image,image->compose,
307  MagickTrue,0,0,exception);
308  highlight_image=DestroyImage(highlight_image);
309  if (status == MagickFalse)
310  difference_image=DestroyImage(difference_image);
311  return(difference_image);
312 }
313 
314 /*
315 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
316 % %
317 % %
318 % %
319 % G e t I m a g e D i s t o r t i o n %
320 % %
321 % %
322 % %
323 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
324 %
325 % GetImageDistortion() compares one or more pixel channels of an image to a
326 % reconstructed image and returns the specified distortion metric.
327 %
328 % The format of the GetImageDistortion method is:
329 %
330 % MagickBooleanType GetImageDistortion(const Image *image,
331 % const Image *reconstruct_image,const MetricType metric,
332 % double *distortion,ExceptionInfo *exception)
333 %
334 % A description of each parameter follows:
335 %
336 % o image: the image.
337 %
338 % o reconstruct_image: the reconstruct image.
339 %
340 % o metric: the metric.
341 %
342 % o distortion: the computed distortion between the images.
343 %
344 % o exception: return any errors or warnings in this structure.
345 %
346 */
347 
349  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
350 {
351  CacheView
352  *image_view,
353  *reconstruct_view;
354 
355  double
356  fuzz;
357 
359  status;
360 
361  size_t
362  columns,
363  rows;
364 
365  ssize_t
366  y;
367 
368  /*
369  Compute the absolute difference in pixels between two images.
370  */
371  status=MagickTrue;
372  fuzz=GetFuzzyColorDistance(image,reconstruct_image);
373  rows=MagickMax(image->rows,reconstruct_image->rows);
374  columns=MagickMax(image->columns,reconstruct_image->columns);
375  image_view=AcquireVirtualCacheView(image,exception);
376  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
377 #if defined(MAGICKCORE_OPENMP_SUPPORT)
378  #pragma omp parallel for schedule(static) shared(status) \
379  magick_number_threads(image,image,rows,1)
380 #endif
381  for (y=0; y < (ssize_t) rows; y++)
382  {
383  double
384  channel_distortion[MaxPixelChannels+1];
385 
386  const Quantum
387  *magick_restrict p,
388  *magick_restrict q;
389 
390  ssize_t
391  j,
392  x;
393 
394  if (status == MagickFalse)
395  continue;
396  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
397  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
398  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
399  {
400  status=MagickFalse;
401  continue;
402  }
403  (void) memset(channel_distortion,0,sizeof(channel_distortion));
404  for (x=0; x < (ssize_t) columns; x++)
405  {
406  double
407  Da,
408  Sa;
409 
411  difference;
412 
413  ssize_t
414  i;
415 
416  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
417  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
418  {
419  p+=GetPixelChannels(image);
420  q+=GetPixelChannels(reconstruct_image);
421  continue;
422  }
423  difference=MagickFalse;
424  Sa=QuantumScale*GetPixelAlpha(image,p);
425  Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
426  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
427  {
428  double
429  distance,
430  pixel;
431 
432  PixelChannel channel = GetPixelChannelChannel(image,i);
433  PixelTrait traits = GetPixelChannelTraits(image,channel);
434  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
435  channel);
436  if ((traits == UndefinedPixelTrait) ||
437  (reconstruct_traits == UndefinedPixelTrait) ||
438  ((reconstruct_traits & UpdatePixelTrait) == 0))
439  continue;
440  if (channel == AlphaPixelChannel)
441  pixel=(double) p[i]-GetPixelChannel(reconstruct_image,channel,q);
442  else
443  pixel=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
444  distance=pixel*pixel;
445  if (distance >= fuzz)
446  {
447  channel_distortion[i]++;
448  difference=MagickTrue;
449  }
450  }
451  if (difference != MagickFalse)
452  channel_distortion[CompositePixelChannel]++;
453  p+=GetPixelChannels(image);
454  q+=GetPixelChannels(reconstruct_image);
455  }
456 #if defined(MAGICKCORE_OPENMP_SUPPORT)
457  #pragma omp critical (MagickCore_GetAbsoluteDistortion)
458 #endif
459  for (j=0; j <= MaxPixelChannels; j++)
460  distortion[j]+=channel_distortion[j];
461  }
462  reconstruct_view=DestroyCacheView(reconstruct_view);
463  image_view=DestroyCacheView(image_view);
464  return(status);
465 }
466 
468  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
469 {
470  CacheView
471  *image_view,
472  *reconstruct_view;
473 
474  double
475  area;
476 
478  status;
479 
480  ssize_t
481  j;
482 
483  size_t
484  columns,
485  rows;
486 
487  ssize_t
488  y;
489 
490  status=MagickTrue;
491  rows=MagickMax(image->rows,reconstruct_image->rows);
492  columns=MagickMax(image->columns,reconstruct_image->columns);
493  area=0.0;
494  image_view=AcquireVirtualCacheView(image,exception);
495  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
496 #if defined(MAGICKCORE_OPENMP_SUPPORT)
497  #pragma omp parallel for schedule(static) shared(status) \
498  magick_number_threads(image,image,rows,1)
499 #endif
500  for (y=0; y < (ssize_t) rows; y++)
501  {
502  double
503  channel_distortion[MaxPixelChannels+1];
504 
505  const Quantum
506  *magick_restrict p,
507  *magick_restrict q;
508 
509  size_t
510  local_area = 0;
511 
512  ssize_t
513  x;
514 
515  if (status == MagickFalse)
516  continue;
517  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
518  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
519  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
520  {
521  status=MagickFalse;
522  continue;
523  }
524  (void) memset(channel_distortion,0,sizeof(channel_distortion));
525  for (x=0; x < (ssize_t) columns; x++)
526  {
527  double
528  Da,
529  Sa;
530 
531  ssize_t
532  i;
533 
534  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
535  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
536  {
537  p+=GetPixelChannels(image);
538  q+=GetPixelChannels(reconstruct_image);
539  continue;
540  }
541  Sa=QuantumScale*GetPixelAlpha(image,p);
542  Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
543  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
544  {
545  double
546  distance;
547 
548  PixelChannel channel = GetPixelChannelChannel(image,i);
549  PixelTrait traits = GetPixelChannelTraits(image,channel);
550  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
551  channel);
552  if ((traits == UndefinedPixelTrait) ||
553  (reconstruct_traits == UndefinedPixelTrait) ||
554  ((reconstruct_traits & UpdatePixelTrait) == 0))
555  continue;
556  if (channel == AlphaPixelChannel)
557  distance=QuantumScale*(p[i]-GetPixelChannel(reconstruct_image,
558  channel,q));
559  else
560  distance=QuantumScale*(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
561  channel,q));
562  channel_distortion[i]+=distance*distance;
563  channel_distortion[CompositePixelChannel]+=distance*distance;
564  }
565  local_area++;
566  p+=GetPixelChannels(image);
567  q+=GetPixelChannels(reconstruct_image);
568  }
569 #if defined(MAGICKCORE_OPENMP_SUPPORT)
570  #pragma omp critical (MagickCore_GetFuzzDistortion)
571 #endif
572  {
573  area+=local_area;
574  for (j=0; j <= MaxPixelChannels; j++)
575  distortion[j]+=channel_distortion[j];
576  }
577  }
578  reconstruct_view=DestroyCacheView(reconstruct_view);
579  image_view=DestroyCacheView(image_view);
580  area=PerceptibleReciprocal(area);
581  for (j=0; j <= MaxPixelChannels; j++)
582  distortion[j]*=area;
583  distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
584  distortion[CompositePixelChannel]=sqrt(distortion[CompositePixelChannel]);
585  return(status);
586 }
587 
589  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
590 {
591  CacheView
592  *image_view,
593  *reconstruct_view;
594 
595  double
596  area;
597 
599  status;
600 
601  ssize_t
602  j;
603 
604  size_t
605  columns,
606  rows;
607 
608  ssize_t
609  y;
610 
611  status=MagickTrue;
612  rows=MagickMax(image->rows,reconstruct_image->rows);
613  columns=MagickMax(image->columns,reconstruct_image->columns);
614  area=0.0;
615  image_view=AcquireVirtualCacheView(image,exception);
616  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
617 #if defined(MAGICKCORE_OPENMP_SUPPORT)
618  #pragma omp parallel for schedule(static) shared(status) \
619  magick_number_threads(image,image,rows,1)
620 #endif
621  for (y=0; y < (ssize_t) rows; y++)
622  {
623  double
624  channel_distortion[MaxPixelChannels+1];
625 
626  const Quantum
627  *magick_restrict p,
628  *magick_restrict q;
629 
630  size_t
631  local_area = 0;
632 
633  ssize_t
634  x;
635 
636  if (status == MagickFalse)
637  continue;
638  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
639  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
640  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
641  {
642  status=MagickFalse;
643  continue;
644  }
645  (void) memset(channel_distortion,0,sizeof(channel_distortion));
646  for (x=0; x < (ssize_t) columns; x++)
647  {
648  double
649  Da,
650  Sa;
651 
652  ssize_t
653  i;
654 
655  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
656  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
657  {
658  p+=GetPixelChannels(image);
659  q+=GetPixelChannels(reconstruct_image);
660  continue;
661  }
662  Sa=QuantumScale*GetPixelAlpha(image,p);
663  Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
664  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
665  {
666  double
667  distance;
668 
669  PixelChannel channel = GetPixelChannelChannel(image,i);
670  PixelTrait traits = GetPixelChannelTraits(image,channel);
671  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
672  channel);
673  if ((traits == UndefinedPixelTrait) ||
674  (reconstruct_traits == UndefinedPixelTrait) ||
675  ((reconstruct_traits & UpdatePixelTrait) == 0))
676  continue;
677  if (channel == AlphaPixelChannel)
678  distance=QuantumScale*fabs((double) (p[i]-(double)
679  GetPixelChannel(reconstruct_image,channel,q)));
680  else
681  distance=QuantumScale*fabs((double) (Sa*p[i]-Da*
682  GetPixelChannel(reconstruct_image,channel,q)));
683  channel_distortion[i]+=distance;
684  channel_distortion[CompositePixelChannel]+=distance;
685  }
686  local_area++;
687  p+=GetPixelChannels(image);
688  q+=GetPixelChannels(reconstruct_image);
689  }
690 #if defined(MAGICKCORE_OPENMP_SUPPORT)
691  #pragma omp critical (MagickCore_GetMeanAbsoluteError)
692 #endif
693  {
694  area+=local_area;
695  for (j=0; j <= MaxPixelChannels; j++)
696  distortion[j]+=channel_distortion[j];
697  }
698  }
699  reconstruct_view=DestroyCacheView(reconstruct_view);
700  image_view=DestroyCacheView(image_view);
701  area=PerceptibleReciprocal(area);
702  for (j=0; j <= MaxPixelChannels; j++)
703  distortion[j]*=area;
704  distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
705  return(status);
706 }
707 
709  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
710 {
711  CacheView
712  *image_view,
713  *reconstruct_view;
714 
716  status;
717 
718  double
719  area,
720  maximum_error,
721  mean_error;
722 
723  size_t
724  columns,
725  rows;
726 
727  ssize_t
728  y;
729 
730  status=MagickTrue;
731  area=0.0;
732  maximum_error=0.0;
733  mean_error=0.0;
734  rows=MagickMax(image->rows,reconstruct_image->rows);
735  columns=MagickMax(image->columns,reconstruct_image->columns);
736  image_view=AcquireVirtualCacheView(image,exception);
737  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
738  for (y=0; y < (ssize_t) rows; y++)
739  {
740  const Quantum
741  *magick_restrict p,
742  *magick_restrict q;
743 
744  ssize_t
745  x;
746 
747  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
748  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
749  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
750  {
751  status=MagickFalse;
752  break;
753  }
754  for (x=0; x < (ssize_t) columns; x++)
755  {
756  double
757  Da,
758  Sa;
759 
760  ssize_t
761  i;
762 
763  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
764  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
765  {
766  p+=GetPixelChannels(image);
767  q+=GetPixelChannels(reconstruct_image);
768  continue;
769  }
770  Sa=QuantumScale*GetPixelAlpha(image,p);
771  Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
772  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
773  {
774  double
775  distance;
776 
777  PixelChannel channel = GetPixelChannelChannel(image,i);
778  PixelTrait traits = GetPixelChannelTraits(image,channel);
779  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
780  channel);
781  if ((traits == UndefinedPixelTrait) ||
782  (reconstruct_traits == UndefinedPixelTrait) ||
783  ((reconstruct_traits & UpdatePixelTrait) == 0))
784  continue;
785  if (channel == AlphaPixelChannel)
786  distance=fabs((double) (p[i]-(double)
787  GetPixelChannel(reconstruct_image,channel,q)));
788  else
789  distance=fabs((double) (Sa*p[i]-Da*
790  GetPixelChannel(reconstruct_image,channel,q)));
791  distortion[i]+=distance;
792  distortion[CompositePixelChannel]+=distance;
793  mean_error+=distance*distance;
794  if (distance > maximum_error)
795  maximum_error=distance;
796  area++;
797  }
798  p+=GetPixelChannels(image);
799  q+=GetPixelChannels(reconstruct_image);
800  }
801  }
802  reconstruct_view=DestroyCacheView(reconstruct_view);
803  image_view=DestroyCacheView(image_view);
804  area=PerceptibleReciprocal(area);
805  image->error.mean_error_per_pixel=area*distortion[CompositePixelChannel];
806  image->error.normalized_mean_error=area*QuantumScale*QuantumScale*mean_error;
807  image->error.normalized_maximum_error=QuantumScale*maximum_error;
808  return(status);
809 }
810 
812  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
813 {
814  CacheView
815  *image_view,
816  *reconstruct_view;
817 
818  double
819  area;
820 
822  status;
823 
824  ssize_t
825  j;
826 
827  size_t
828  columns,
829  rows;
830 
831  ssize_t
832  y;
833 
834  status=MagickTrue;
835  rows=MagickMax(image->rows,reconstruct_image->rows);
836  columns=MagickMax(image->columns,reconstruct_image->columns);
837  area=0.0;
838  image_view=AcquireVirtualCacheView(image,exception);
839  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
840 #if defined(MAGICKCORE_OPENMP_SUPPORT)
841  #pragma omp parallel for schedule(static) shared(status) \
842  magick_number_threads(image,image,rows,1)
843 #endif
844  for (y=0; y < (ssize_t) rows; y++)
845  {
846  const Quantum
847  *magick_restrict p,
848  *magick_restrict q;
849 
850  double
851  channel_distortion[MaxPixelChannels+1];
852 
853  size_t
854  local_area = 0;
855 
856  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  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  local_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  {
917  area+=local_area;
918  for (j=0; j <= MaxPixelChannels; j++)
919  distortion[j]+=channel_distortion[j];
920  }
921  }
922  reconstruct_view=DestroyCacheView(reconstruct_view);
923  image_view=DestroyCacheView(image_view);
924  area=PerceptibleReciprocal(area);
925  for (j=0; j <= MaxPixelChannels; j++)
926  distortion[j]*=area;
927  distortion[CompositePixelChannel]/=GetImageChannels(image);
928  return(status);
929 }
930 
932  const Image *image,const Image *reconstruct_image,double *distortion,
933  ExceptionInfo *exception)
934 {
935 #define SimilarityImageTag "Similarity/Image"
936 
937  CacheView
938  *image_view,
939  *reconstruct_view;
940 
942  *image_statistics,
943  *reconstruct_statistics;
944 
945  double
946  area;
947 
949  status;
950 
952  progress;
953 
954  ssize_t
955  channels,
956  i;
957 
958  size_t
959  columns,
960  rows;
961 
962  ssize_t
963  y;
964 
965  /*
966  Normalize to account for variation due to lighting and exposure condition.
967  */
968  image_statistics=GetImageStatistics(image,exception);
969  reconstruct_statistics=GetImageStatistics(reconstruct_image,exception);
970  if ((image_statistics == (ChannelStatistics *) NULL) ||
971  (reconstruct_statistics == (ChannelStatistics *) NULL))
972  {
973  if (image_statistics != (ChannelStatistics *) NULL)
974  image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
975  image_statistics);
976  if (reconstruct_statistics != (ChannelStatistics *) NULL)
977  reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
978  reconstruct_statistics);
979  return(MagickFalse);
980  }
981  status=MagickTrue;
982  progress=0;
983  for (i=0; i <= MaxPixelChannels; i++)
984  distortion[i]=0.0;
985  rows=MagickMax(image->rows,reconstruct_image->rows);
986  columns=MagickMax(image->columns,reconstruct_image->columns);
987  area=0.0;
988  image_view=AcquireVirtualCacheView(image,exception);
989  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
990  for (y=0; y < (ssize_t) rows; y++)
991  {
992  const Quantum
993  *magick_restrict p,
994  *magick_restrict q;
995 
996  ssize_t
997  x;
998 
999  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1000  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1001  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1002  {
1003  status=MagickFalse;
1004  break;
1005  }
1006  for (x=0; x < (ssize_t) columns; x++)
1007  {
1008  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1009  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1010  {
1011  p+=GetPixelChannels(image);
1012  q+=GetPixelChannels(reconstruct_image);
1013  continue;
1014  }
1015  area++;
1016  p+=GetPixelChannels(image);
1017  q+=GetPixelChannels(reconstruct_image);
1018  }
1019  }
1020  area=PerceptibleReciprocal(area);
1021  for (y=0; y < (ssize_t) rows; y++)
1022  {
1023  const Quantum
1024  *magick_restrict p,
1025  *magick_restrict q;
1026 
1027  ssize_t
1028  x;
1029 
1030  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1031  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1032  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1033  {
1034  status=MagickFalse;
1035  break;
1036  }
1037  for (x=0; x < (ssize_t) columns; x++)
1038  {
1039  double
1040  Da,
1041  Sa;
1042 
1043  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1044  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1045  {
1046  p+=GetPixelChannels(image);
1047  q+=GetPixelChannels(reconstruct_image);
1048  continue;
1049  }
1050  Sa=QuantumScale*GetPixelAlpha(image,p);
1051  Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
1052  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1053  {
1054  PixelChannel channel = GetPixelChannelChannel(image,i);
1055  PixelTrait traits = GetPixelChannelTraits(image,channel);
1056  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1057  channel);
1058  if ((traits == UndefinedPixelTrait) ||
1059  (reconstruct_traits == UndefinedPixelTrait) ||
1060  ((reconstruct_traits & UpdatePixelTrait) == 0))
1061  continue;
1062  if (channel == AlphaPixelChannel)
1063  distortion[i]+=area*QuantumScale*((double) p[i]-
1064  image_statistics[channel].mean)*(GetPixelChannel(reconstruct_image,
1065  channel,q)-reconstruct_statistics[channel].mean);
1066  else
1067  distortion[i]+=area*QuantumScale*(Sa*p[i]-
1068  image_statistics[channel].mean)*(Da*GetPixelChannel(
1069  reconstruct_image,channel,q)-reconstruct_statistics[channel].mean);
1070  }
1071  p+=GetPixelChannels(image);
1072  q+=GetPixelChannels(reconstruct_image);
1073  }
1074  if (image->progress_monitor != (MagickProgressMonitor) NULL)
1075  {
1077  proceed;
1078 
1079 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1080  #pragma omp atomic
1081 #endif
1082  progress++;
1083  proceed=SetImageProgress(image,SimilarityImageTag,progress,rows);
1084  if (proceed == MagickFalse)
1085  {
1086  status=MagickFalse;
1087  break;
1088  }
1089  }
1090  }
1091  reconstruct_view=DestroyCacheView(reconstruct_view);
1092  image_view=DestroyCacheView(image_view);
1093  /*
1094  Divide by the standard deviation.
1095  */
1096  channels=0;
1097  distortion[CompositePixelChannel]=0.0;
1098  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1099  {
1100  double
1101  gamma;
1102 
1103  PixelChannel channel = GetPixelChannelChannel(image,i);
1104  gamma=image_statistics[channel].standard_deviation*
1105  reconstruct_statistics[channel].standard_deviation;
1106  if (fabs(gamma) >= MagickEpsilon)
1107  {
1108  gamma=PerceptibleReciprocal(gamma);
1109  distortion[i]=QuantumRange*gamma*distortion[i];
1110  distortion[CompositePixelChannel]+=distortion[i]*distortion[i];
1111  channels++;
1112  }
1113  }
1114  distortion[CompositePixelChannel]=sqrt(distortion[CompositePixelChannel]/
1115  channels);
1116  /*
1117  Free resources.
1118  */
1119  reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1120  reconstruct_statistics);
1121  image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1122  image_statistics);
1123  return(status);
1124 }
1125 
1127  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1128 {
1129  CacheView
1130  *image_view,
1131  *reconstruct_view;
1132 
1134  status;
1135 
1136  size_t
1137  columns,
1138  rows;
1139 
1140  ssize_t
1141  y;
1142 
1143  status=MagickTrue;
1144  rows=MagickMax(image->rows,reconstruct_image->rows);
1145  columns=MagickMax(image->columns,reconstruct_image->columns);
1146  image_view=AcquireVirtualCacheView(image,exception);
1147  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1148 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1149  #pragma omp parallel for schedule(static) shared(status) \
1150  magick_number_threads(image,image,rows,1)
1151 #endif
1152  for (y=0; y < (ssize_t) rows; y++)
1153  {
1154  double
1155  channel_distortion[MaxPixelChannels+1];
1156 
1157  const Quantum
1158  *magick_restrict p,
1159  *magick_restrict q;
1160 
1161  ssize_t
1162  j,
1163  x;
1164 
1165  if (status == MagickFalse)
1166  continue;
1167  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1168  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1169  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1170  {
1171  status=MagickFalse;
1172  continue;
1173  }
1174  (void) memset(channel_distortion,0,sizeof(channel_distortion));
1175  for (x=0; x < (ssize_t) columns; x++)
1176  {
1177  double
1178  Da,
1179  Sa;
1180 
1181  ssize_t
1182  i;
1183 
1184  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1185  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1186  {
1187  p+=GetPixelChannels(image);
1188  q+=GetPixelChannels(reconstruct_image);
1189  continue;
1190  }
1191  Sa=QuantumScale*GetPixelAlpha(image,p);
1192  Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
1193  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1194  {
1195  double
1196  distance;
1197 
1198  PixelChannel channel = GetPixelChannelChannel(image,i);
1199  PixelTrait traits = GetPixelChannelTraits(image,channel);
1200  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1201  channel);
1202  if ((traits == UndefinedPixelTrait) ||
1203  (reconstruct_traits == UndefinedPixelTrait) ||
1204  ((reconstruct_traits & UpdatePixelTrait) == 0))
1205  continue;
1206  if (channel == AlphaPixelChannel)
1207  distance=QuantumScale*fabs((double) (p[i]-(double)
1208  GetPixelChannel(reconstruct_image,channel,q)));
1209  else
1210  distance=QuantumScale*fabs((double) (Sa*p[i]-Da*
1211  GetPixelChannel(reconstruct_image,channel,q)));
1212  if (distance > channel_distortion[i])
1213  channel_distortion[i]=distance;
1214  if (distance > channel_distortion[CompositePixelChannel])
1215  channel_distortion[CompositePixelChannel]=distance;
1216  }
1217  p+=GetPixelChannels(image);
1218  q+=GetPixelChannels(reconstruct_image);
1219  }
1220 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1221  #pragma omp critical (MagickCore_GetPeakAbsoluteError)
1222 #endif
1223  for (j=0; j <= MaxPixelChannels; j++)
1224  if (channel_distortion[j] > distortion[j])
1225  distortion[j]=channel_distortion[j];
1226  }
1227  reconstruct_view=DestroyCacheView(reconstruct_view);
1228  image_view=DestroyCacheView(image_view);
1229  return(status);
1230 }
1231 
1232 static inline double MagickLog10(const double x)
1233 {
1234 #define Log10Epsilon (1.0e-11)
1235 
1236  if (fabs(x) < Log10Epsilon)
1237  return(log10(Log10Epsilon));
1238  return(log10(fabs(x)));
1239 }
1240 
1242  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1243 {
1245  status;
1246 
1247  ssize_t
1248  i;
1249 
1250  status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1251  for (i=0; i <= MaxPixelChannels; i++)
1252  if (fabs(distortion[i]) < MagickEpsilon)
1253  distortion[i]=INFINITY;
1254  else
1255  distortion[i]=10.0*MagickLog10(1.0)-10.0*MagickLog10(distortion[i]);
1256  return(status);
1257 }
1258 
1260  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1261 {
1263  *channel_phash,
1264  *reconstruct_phash;
1265 
1266  const char
1267  *artifact;
1268 
1270  normalize;
1271 
1272  ssize_t
1273  channel;
1274 
1275  /*
1276  Compute perceptual hash in the sRGB colorspace.
1277  */
1278  channel_phash=GetImagePerceptualHash(image,exception);
1279  if (channel_phash == (ChannelPerceptualHash *) NULL)
1280  return(MagickFalse);
1281  reconstruct_phash=GetImagePerceptualHash(reconstruct_image,exception);
1282  if (reconstruct_phash == (ChannelPerceptualHash *) NULL)
1283  {
1285  channel_phash);
1286  return(MagickFalse);
1287  }
1288  artifact=GetImageArtifact(image,"phash:normalize");
1289  normalize=(artifact == (const char *) NULL) ||
1290  (IsStringTrue(artifact) == MagickFalse) ? MagickFalse : MagickTrue;
1291 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1292  #pragma omp parallel for schedule(static)
1293 #endif
1294  for (channel=0; channel < MaxPixelChannels; channel++)
1295  {
1296  double
1297  difference;
1298 
1299  ssize_t
1300  i;
1301 
1302  difference=0.0;
1303  for (i=0; i < MaximumNumberOfImageMoments; i++)
1304  {
1305  double
1306  alpha,
1307  beta;
1308 
1309  ssize_t
1310  j;
1311 
1312  for (j=0; j < (ssize_t) channel_phash[0].number_colorspaces; j++)
1313  {
1314  alpha=channel_phash[channel].phash[j][i];
1315  beta=reconstruct_phash[channel].phash[j][i];
1316  if (normalize == MagickFalse)
1317  difference+=(beta-alpha)*(beta-alpha);
1318  else
1319  difference=sqrt((beta-alpha)*(beta-alpha)/
1320  channel_phash[0].number_channels);
1321  }
1322  }
1323  distortion[channel]+=difference;
1324 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1325  #pragma omp critical (MagickCore_GetPerceptualHashDistortion)
1326 #endif
1327  distortion[CompositePixelChannel]+=difference;
1328  }
1329  /*
1330  Free resources.
1331  */
1332  reconstruct_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1333  reconstruct_phash);
1334  channel_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(channel_phash);
1335  return(MagickTrue);
1336 }
1337 
1339  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1340 {
1342  status;
1343 
1344  ssize_t
1345  i;
1346 
1347  status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1348  for (i=0; i <= MaxPixelChannels; i++)
1349  distortion[i]=sqrt(distortion[i]);
1350  return(status);
1351 }
1352 
1354  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1355 {
1356 #define SSIMRadius 5.0
1357 #define SSIMSigma 1.5
1358 #define SSIMBlocksize 8
1359 #define SSIMK1 0.01
1360 #define SSIMK2 0.03
1361 #define SSIML 1.0
1362 
1363  CacheView
1364  *image_view,
1365  *reconstruct_view;
1366 
1367  char
1368  geometry[MagickPathExtent];
1369 
1370  const char
1371  *artifact;
1372 
1373  double
1374  area,
1375  c1,
1376  c2,
1377  radius,
1378  sigma;
1379 
1380  KernelInfo
1381  *kernel_info;
1382 
1384  status;
1385 
1386  ssize_t
1387  j;
1388 
1389  size_t
1390  columns,
1391  rows;
1392 
1393  ssize_t
1394  y;
1395 
1396  /*
1397  Compute structural similarity index @
1398  https://en.wikipedia.org/wiki/Structural_similarity.
1399  */
1400  radius=SSIMRadius;
1401  artifact=GetImageArtifact(image,"compare:ssim-radius");
1402  if (artifact != (const char *) NULL)
1403  radius=StringToDouble(artifact,(char **) NULL);
1404  sigma=SSIMSigma;
1405  artifact=GetImageArtifact(image,"compare:ssim-sigma");
1406  if (artifact != (const char *) NULL)
1407  sigma=StringToDouble(artifact,(char **) NULL);
1408  (void) FormatLocaleString(geometry,MagickPathExtent,"gaussian:%.20gx%.20g",
1409  radius,sigma);
1410  kernel_info=AcquireKernelInfo(geometry,exception);
1411  if (kernel_info == (KernelInfo *) NULL)
1412  ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
1413  image->filename);
1414  c1=pow(SSIMK1*SSIML,2.0);
1415  artifact=GetImageArtifact(image,"compare:ssim-k1");
1416  if (artifact != (const char *) NULL)
1417  c1=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
1418  c2=pow(SSIMK2*SSIML,2.0);
1419  artifact=GetImageArtifact(image,"compare:ssim-k2");
1420  if (artifact != (const char *) NULL)
1421  c2=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
1422  status=MagickTrue;
1423  area=0.0;
1424  rows=MagickMax(image->rows,reconstruct_image->rows);
1425  columns=MagickMax(image->columns,reconstruct_image->columns);
1426  image_view=AcquireVirtualCacheView(image,exception);
1427  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1428 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1429  #pragma omp parallel for schedule(static) shared(status) \
1430  magick_number_threads(image,reconstruct_image,rows,1)
1431 #endif
1432  for (y=0; y < (ssize_t) rows; y++)
1433  {
1434  double
1435  channel_distortion[MaxPixelChannels+1];
1436 
1437  const Quantum
1438  *magick_restrict p,
1439  *magick_restrict q;
1440 
1441  size_t
1442  local_area = 0;
1443 
1444  ssize_t
1445  i,
1446  x;
1447 
1448  if (status == MagickFalse)
1449  continue;
1450  p=GetCacheViewVirtualPixels(image_view,-((ssize_t) kernel_info->width/2L),y-
1451  ((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
1452  kernel_info->height,exception);
1453  q=GetCacheViewVirtualPixels(reconstruct_view,-((ssize_t) kernel_info->width/
1454  2L),y-((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
1455  kernel_info->height,exception);
1456  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1457  {
1458  status=MagickFalse;
1459  continue;
1460  }
1461  (void) memset(channel_distortion,0,sizeof(channel_distortion));
1462  for (x=0; x < (ssize_t) columns; x++)
1463  {
1464  double
1465  x_pixel_mu[MaxPixelChannels+1],
1466  x_pixel_sigma_squared[MaxPixelChannels+1],
1467  xy_sigma[MaxPixelChannels+1],
1468  y_pixel_mu[MaxPixelChannels+1],
1469  y_pixel_sigma_squared[MaxPixelChannels+1];
1470 
1471  const Quantum
1472  *magick_restrict reference,
1473  *magick_restrict target;
1474 
1476  *k;
1477 
1478  ssize_t
1479  v;
1480 
1481  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1482  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1483  {
1484  p+=GetPixelChannels(image);
1485  q+=GetPixelChannels(reconstruct_image);
1486  continue;
1487  }
1488  (void) memset(x_pixel_mu,0,sizeof(x_pixel_mu));
1489  (void) memset(x_pixel_sigma_squared,0,sizeof(x_pixel_sigma_squared));
1490  (void) memset(xy_sigma,0,sizeof(xy_sigma));
1491  (void) memset(x_pixel_sigma_squared,0,sizeof(y_pixel_sigma_squared));
1492  (void) memset(y_pixel_mu,0,sizeof(y_pixel_mu));
1493  (void) memset(y_pixel_sigma_squared,0,sizeof(y_pixel_sigma_squared));
1494  k=kernel_info->values;
1495  reference=p;
1496  target=q;
1497  for (v=0; v < (ssize_t) kernel_info->height; v++)
1498  {
1499  ssize_t
1500  u;
1501 
1502  for (u=0; u < (ssize_t) kernel_info->width; u++)
1503  {
1504  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1505  {
1506  double
1507  x_pixel,
1508  y_pixel;
1509 
1510  PixelChannel channel = GetPixelChannelChannel(image,i);
1511  PixelTrait traits = GetPixelChannelTraits(image,channel);
1512  PixelTrait reconstruct_traits = GetPixelChannelTraits(
1513  reconstruct_image,channel);
1514  if ((traits == UndefinedPixelTrait) ||
1515  (reconstruct_traits == UndefinedPixelTrait) ||
1516  ((reconstruct_traits & UpdatePixelTrait) == 0))
1517  continue;
1518  x_pixel=QuantumScale*reference[i];
1519  x_pixel_mu[i]+=(*k)*x_pixel;
1520  x_pixel_sigma_squared[i]+=(*k)*x_pixel*x_pixel;
1521  y_pixel=QuantumScale*
1522  GetPixelChannel(reconstruct_image,channel,target);
1523  y_pixel_mu[i]+=(*k)*y_pixel;
1524  y_pixel_sigma_squared[i]+=(*k)*y_pixel*y_pixel;
1525  xy_sigma[i]+=(*k)*x_pixel*y_pixel;
1526  }
1527  k++;
1528  reference+=GetPixelChannels(image);
1529  target+=GetPixelChannels(reconstruct_image);
1530  }
1531  reference+=GetPixelChannels(image)*columns;
1532  target+=GetPixelChannels(reconstruct_image)*columns;
1533  }
1534  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1535  {
1536  double
1537  ssim,
1538  x_pixel_mu_squared,
1539  x_pixel_sigmas_squared,
1540  xy_mu,
1541  xy_sigmas,
1542  y_pixel_mu_squared,
1543  y_pixel_sigmas_squared;
1544 
1545  PixelChannel channel = GetPixelChannelChannel(image,i);
1546  PixelTrait traits = GetPixelChannelTraits(image,channel);
1547  PixelTrait reconstruct_traits = GetPixelChannelTraits(
1548  reconstruct_image,channel);
1549  if ((traits == UndefinedPixelTrait) ||
1550  (reconstruct_traits == UndefinedPixelTrait) ||
1551  ((reconstruct_traits & UpdatePixelTrait) == 0))
1552  continue;
1553  x_pixel_mu_squared=x_pixel_mu[i]*x_pixel_mu[i];
1554  y_pixel_mu_squared=y_pixel_mu[i]*y_pixel_mu[i];
1555  xy_mu=x_pixel_mu[i]*y_pixel_mu[i];
1556  xy_sigmas=xy_sigma[i]-xy_mu;
1557  x_pixel_sigmas_squared=x_pixel_sigma_squared[i]-x_pixel_mu_squared;
1558  y_pixel_sigmas_squared=y_pixel_sigma_squared[i]-y_pixel_mu_squared;
1559  ssim=((2.0*xy_mu+c1)*(2.0*xy_sigmas+c2))/
1560  ((x_pixel_mu_squared+y_pixel_mu_squared+c1)*
1561  (x_pixel_sigmas_squared+y_pixel_sigmas_squared+c2));
1562  channel_distortion[i]+=ssim;
1563  channel_distortion[CompositePixelChannel]+=ssim;
1564  }
1565  p+=GetPixelChannels(image);
1566  q+=GetPixelChannels(reconstruct_image);
1567  local_area++;
1568  }
1569 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1570  #pragma omp critical (MagickCore_GetStructuralSimilarityDistortion)
1571 #endif
1572  {
1573  area+=local_area;
1574  for (i=0; i <= MaxPixelChannels; i++)
1575  distortion[i]+=channel_distortion[i];
1576  }
1577  }
1578  image_view=DestroyCacheView(image_view);
1579  reconstruct_view=DestroyCacheView(reconstruct_view);
1580  for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
1581  {
1582  PixelChannel channel = GetPixelChannelChannel(image,j);
1583  PixelTrait traits = GetPixelChannelTraits(image,channel);
1584  if ((traits == UndefinedPixelTrait) || ((traits & UpdatePixelTrait) == 0))
1585  continue;
1586  distortion[j]/=area;
1587  }
1588  distortion[CompositePixelChannel]/=area;
1589  distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
1590  kernel_info=DestroyKernelInfo(kernel_info);
1591  return(status);
1592 }
1593 
1595  const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1596 {
1598  status;
1599 
1600  ssize_t
1601  i;
1602 
1603  status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1604  distortion,exception);
1605  for (i=0; i <= MaxPixelChannels; i++)
1606  distortion[i]=(1.0-(distortion[i]))/2.0;
1607  return(status);
1608 }
1609 
1611  const Image *reconstruct_image,const MetricType metric,double *distortion,
1612  ExceptionInfo *exception)
1613 {
1614  double
1615  *channel_distortion;
1616 
1618  status;
1619 
1620  size_t
1621  length;
1622 
1623  assert(image != (Image *) NULL);
1624  assert(image->signature == MagickCoreSignature);
1625  if (image->debug != MagickFalse)
1626  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1627  assert(reconstruct_image != (const Image *) NULL);
1628  assert(reconstruct_image->signature == MagickCoreSignature);
1629  assert(distortion != (double *) NULL);
1630  *distortion=0.0;
1631  if (image->debug != MagickFalse)
1632  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1633  /*
1634  Get image distortion.
1635  */
1636  length=MaxPixelChannels+1UL;
1637  channel_distortion=(double *) AcquireQuantumMemory(length,
1638  sizeof(*channel_distortion));
1639  if (channel_distortion == (double *) NULL)
1640  ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1641  (void) memset(channel_distortion,0,length*
1642  sizeof(*channel_distortion));
1643  switch (metric)
1644  {
1645  case AbsoluteErrorMetric:
1646  {
1647  status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1648  exception);
1649  break;
1650  }
1651  case FuzzErrorMetric:
1652  {
1653  status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1654  exception);
1655  break;
1656  }
1658  {
1659  status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1660  channel_distortion,exception);
1661  break;
1662  }
1664  {
1665  status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1666  exception);
1667  break;
1668  }
1670  {
1671  status=GetMeanSquaredDistortion(image,reconstruct_image,
1672  channel_distortion,exception);
1673  break;
1674  }
1676  default:
1677  {
1678  status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1679  channel_distortion,exception);
1680  break;
1681  }
1683  {
1684  status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1685  channel_distortion,exception);
1686  break;
1687  }
1689  {
1690  status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1691  channel_distortion,exception);
1692  break;
1693  }
1695  {
1696  status=GetPerceptualHashDistortion(image,reconstruct_image,
1697  channel_distortion,exception);
1698  break;
1699  }
1701  {
1702  status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1703  channel_distortion,exception);
1704  break;
1705  }
1707  {
1708  status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1709  channel_distortion,exception);
1710  break;
1711  }
1713  {
1714  status=GetStructuralDisimilarityDistortion(image,reconstruct_image,
1715  channel_distortion,exception);
1716  break;
1717  }
1718  }
1719  *distortion=channel_distortion[CompositePixelChannel];
1720  channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1721  (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
1722  *distortion);
1723  return(status);
1724 }
1725 
1726 /*
1727 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1728 % %
1729 % %
1730 % %
1731 % G e t I m a g e D i s t o r t i o n s %
1732 % %
1733 % %
1734 % %
1735 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1736 %
1737 % GetImageDistortions() compares the pixel channels of an image to a
1738 % reconstructed image and returns the specified distortion metric for each
1739 % channel.
1740 %
1741 % The format of the GetImageDistortions method is:
1742 %
1743 % double *GetImageDistortions(const Image *image,
1744 % const Image *reconstruct_image,const MetricType metric,
1745 % ExceptionInfo *exception)
1746 %
1747 % A description of each parameter follows:
1748 %
1749 % o image: the image.
1750 %
1751 % o reconstruct_image: the reconstruct image.
1752 %
1753 % o metric: the metric.
1754 %
1755 % o exception: return any errors or warnings in this structure.
1756 %
1757 */
1759  const Image *reconstruct_image,const MetricType metric,
1760  ExceptionInfo *exception)
1761 {
1762  double
1763  *channel_distortion;
1764 
1766  status;
1767 
1768  size_t
1769  length;
1770 
1771  assert(image != (Image *) NULL);
1772  assert(image->signature == MagickCoreSignature);
1773  if (image->debug != MagickFalse)
1774  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1775  assert(reconstruct_image != (const Image *) NULL);
1776  assert(reconstruct_image->signature == MagickCoreSignature);
1777  if (image->debug != MagickFalse)
1778  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1779  /*
1780  Get image distortion.
1781  */
1782  length=MaxPixelChannels+1UL;
1783  channel_distortion=(double *) AcquireQuantumMemory(length,
1784  sizeof(*channel_distortion));
1785  if (channel_distortion == (double *) NULL)
1786  ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1787  (void) memset(channel_distortion,0,length*
1788  sizeof(*channel_distortion));
1789  status=MagickTrue;
1790  switch (metric)
1791  {
1792  case AbsoluteErrorMetric:
1793  {
1794  status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1795  exception);
1796  break;
1797  }
1798  case FuzzErrorMetric:
1799  {
1800  status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1801  exception);
1802  break;
1803  }
1805  {
1806  status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1807  channel_distortion,exception);
1808  break;
1809  }
1811  {
1812  status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1813  exception);
1814  break;
1815  }
1817  {
1818  status=GetMeanSquaredDistortion(image,reconstruct_image,
1819  channel_distortion,exception);
1820  break;
1821  }
1823  default:
1824  {
1825  status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1826  channel_distortion,exception);
1827  break;
1828  }
1830  {
1831  status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1832  channel_distortion,exception);
1833  break;
1834  }
1836  {
1837  status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1838  channel_distortion,exception);
1839  break;
1840  }
1842  {
1843  status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1844  channel_distortion,exception);
1845  break;
1846  }
1848  {
1849  status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1850  channel_distortion,exception);
1851  break;
1852  }
1854  {
1855  status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1856  channel_distortion,exception);
1857  break;
1858  }
1860  {
1861  status=GetStructuralDisimilarityDistortion(image,reconstruct_image,
1862  channel_distortion,exception);
1863  break;
1864  }
1865  }
1866  if (status == MagickFalse)
1867  {
1868  channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1869  return((double *) NULL);
1870  }
1871  return(channel_distortion);
1872 }
1873 
1874 /*
1875 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1876 % %
1877 % %
1878 % %
1879 % I s I m a g e s E q u a l %
1880 % %
1881 % %
1882 % %
1883 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1884 %
1885 % IsImagesEqual() compare the pixels of two images and returns immediately
1886 % if any pixel is not identical.
1887 %
1888 % The format of the IsImagesEqual method is:
1889 %
1890 % MagickBooleanType IsImagesEqual(const Image *image,
1891 % const Image *reconstruct_image,ExceptionInfo *exception)
1892 %
1893 % A description of each parameter follows.
1894 %
1895 % o image: the image.
1896 %
1897 % o reconstruct_image: the reconstruct image.
1898 %
1899 % o exception: return any errors or warnings in this structure.
1900 %
1901 */
1903  const Image *reconstruct_image,ExceptionInfo *exception)
1904 {
1905  CacheView
1906  *image_view,
1907  *reconstruct_view;
1908 
1909  size_t
1910  columns,
1911  rows;
1912 
1913  ssize_t
1914  y;
1915 
1916  assert(image != (Image *) NULL);
1917  assert(image->signature == MagickCoreSignature);
1918  assert(reconstruct_image != (const Image *) NULL);
1919  assert(reconstruct_image->signature == MagickCoreSignature);
1920  rows=MagickMax(image->rows,reconstruct_image->rows);
1921  columns=MagickMax(image->columns,reconstruct_image->columns);
1922  image_view=AcquireVirtualCacheView(image,exception);
1923  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1924  for (y=0; y < (ssize_t) rows; y++)
1925  {
1926  const Quantum
1927  *magick_restrict p,
1928  *magick_restrict q;
1929 
1930  ssize_t
1931  x;
1932 
1933  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1934  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1935  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1936  break;
1937  for (x=0; x < (ssize_t) columns; x++)
1938  {
1939  ssize_t
1940  i;
1941 
1942  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1943  {
1944  double
1945  distance;
1946 
1947  PixelChannel channel = GetPixelChannelChannel(image,i);
1948  PixelTrait traits = GetPixelChannelTraits(image,channel);
1949  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1950  channel);
1951  if ((traits == UndefinedPixelTrait) ||
1952  (reconstruct_traits == UndefinedPixelTrait) ||
1953  ((reconstruct_traits & UpdatePixelTrait) == 0))
1954  continue;
1955  distance=fabs((double) (p[i]-(double) GetPixelChannel(reconstruct_image,
1956  channel,q)));
1957  if (distance >= MagickEpsilon)
1958  break;
1959  }
1960  if (i < (ssize_t) GetPixelChannels(image))
1961  break;
1962  p+=GetPixelChannels(image);
1963  q+=GetPixelChannels(reconstruct_image);
1964  }
1965  if (x < (ssize_t) columns)
1966  break;
1967  }
1968  reconstruct_view=DestroyCacheView(reconstruct_view);
1969  image_view=DestroyCacheView(image_view);
1970  return(y < (ssize_t) rows ? MagickFalse : MagickTrue);
1971 }
1972 
1973 /*
1974 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1975 % %
1976 % %
1977 % %
1978 % S e t I m a g e C o l o r M e t r i c %
1979 % %
1980 % %
1981 % %
1982 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1983 %
1984 % SetImageColorMetric() measures the difference between colors at each pixel
1985 % location of two images. A value other than 0 means the colors match
1986 % exactly. Otherwise an error measure is computed by summing over all
1987 % pixels in an image the distance squared in RGB space between each image
1988 % pixel and its corresponding pixel in the reconstruct image. The error
1989 % measure is assigned to these image members:
1990 %
1991 % o mean_error_per_pixel: The mean error for any single pixel in
1992 % the image.
1993 %
1994 % o normalized_mean_error: The normalized mean quantization error for
1995 % any single pixel in the image. This distance measure is normalized to
1996 % a range between 0 and 1. It is independent of the range of red, green,
1997 % and blue values in the image.
1998 %
1999 % o normalized_maximum_error: The normalized maximum quantization
2000 % error for any single pixel in the image. This distance measure is
2001 % normalized to a range between 0 and 1. It is independent of the range
2002 % of red, green, and blue values in your image.
2003 %
2004 % A small normalized mean square error, accessed as
2005 % image->normalized_mean_error, suggests the images are very similar in
2006 % spatial layout and color.
2007 %
2008 % The format of the SetImageColorMetric method is:
2009 %
2010 % MagickBooleanType SetImageColorMetric(Image *image,
2011 % const Image *reconstruct_image,ExceptionInfo *exception)
2012 %
2013 % A description of each parameter follows.
2014 %
2015 % o image: the image.
2016 %
2017 % o reconstruct_image: the reconstruct image.
2018 %
2019 % o exception: return any errors or warnings in this structure.
2020 %
2021 */
2023  const Image *reconstruct_image,ExceptionInfo *exception)
2024 {
2025  CacheView
2026  *image_view,
2027  *reconstruct_view;
2028 
2029  double
2030  area,
2031  maximum_error,
2032  mean_error,
2033  mean_error_per_pixel;
2034 
2036  status;
2037 
2038  size_t
2039  columns,
2040  rows;
2041 
2042  ssize_t
2043  y;
2044 
2045  assert(image != (Image *) NULL);
2046  assert(image->signature == MagickCoreSignature);
2047  assert(reconstruct_image != (const Image *) NULL);
2048  assert(reconstruct_image->signature == MagickCoreSignature);
2049  area=0.0;
2050  maximum_error=0.0;
2051  mean_error_per_pixel=0.0;
2052  mean_error=0.0;
2053  rows=MagickMax(image->rows,reconstruct_image->rows);
2054  columns=MagickMax(image->columns,reconstruct_image->columns);
2055  image_view=AcquireVirtualCacheView(image,exception);
2056  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
2057  for (y=0; y < (ssize_t) rows; y++)
2058  {
2059  const Quantum
2060  *magick_restrict p,
2061  *magick_restrict q;
2062 
2063  ssize_t
2064  x;
2065 
2066  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
2067  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
2068  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2069  break;
2070  for (x=0; x < (ssize_t) columns; x++)
2071  {
2072  ssize_t
2073  i;
2074 
2075  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2076  {
2077  double
2078  distance;
2079 
2080  PixelChannel channel = GetPixelChannelChannel(image,i);
2081  PixelTrait traits = GetPixelChannelTraits(image,channel);
2082  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2083  channel);
2084  if ((traits == UndefinedPixelTrait) ||
2085  (reconstruct_traits == UndefinedPixelTrait) ||
2086  ((reconstruct_traits & UpdatePixelTrait) == 0))
2087  continue;
2088  distance=fabs((double) (p[i]-(double) GetPixelChannel(reconstruct_image,
2089  channel,q)));
2090  if (distance >= MagickEpsilon)
2091  {
2092  mean_error_per_pixel+=distance;
2093  mean_error+=distance*distance;
2094  if (distance > maximum_error)
2095  maximum_error=distance;
2096  }
2097  area++;
2098  }
2099  p+=GetPixelChannels(image);
2100  q+=GetPixelChannels(reconstruct_image);
2101  }
2102  }
2103  reconstruct_view=DestroyCacheView(reconstruct_view);
2104  image_view=DestroyCacheView(image_view);
2105  image->error.mean_error_per_pixel=(double) (mean_error_per_pixel/area);
2107  mean_error/area);
2108  image->error.normalized_maximum_error=(double) (QuantumScale*maximum_error);
2109  status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
2110  return(status);
2111 }
2112 
2113 /*
2114 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2115 % %
2116 % %
2117 % %
2118 % S i m i l a r i t y I m a g e %
2119 % %
2120 % %
2121 % %
2122 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2123 %
2124 % SimilarityImage() compares the reference image of the image and returns the
2125 % best match offset. In addition, it returns a similarity image such that an
2126 % exact match location is completely white and if none of the pixels match,
2127 % black, otherwise some gray level in-between.
2128 %
2129 % The format of the SimilarityImageImage method is:
2130 %
2131 % Image *SimilarityImage(const Image *image,const Image *reference,
2132 % const MetricType metric,const double similarity_threshold,
2133 % RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
2134 %
2135 % A description of each parameter follows:
2136 %
2137 % o image: the image.
2138 %
2139 % o reference: find an area of the image that closely resembles this image.
2140 %
2141 % o metric: the metric.
2142 %
2143 % o similarity_threshold: minimum distortion for (sub)image match.
2144 %
2145 % o offset: the best match offset of the reference image within the image.
2146 %
2147 % o similarity: the computed similarity between the images.
2148 %
2149 % o exception: return any errors or warnings in this structure.
2150 %
2151 */
2152 
2153 #if defined(MAGICKCORE_HDRI_SUPPORT) && defined(MAGICKCORE_FFTW_DELEGATE)
2154 static Image *CrossCorrelationImage(const Image *alpha_image,
2155  const Image *beta_image,ExceptionInfo *exception)
2156 {
2157  Image
2158  *clone_image,
2159  *complex_conjugate,
2160  *complex_multiplication,
2161  *cross_correlation,
2162  *fft_images;
2163 
2164  /*
2165  Take the FFT of beta image.
2166  */
2167  clone_image=CloneImage(beta_image,0,0,MagickTrue,exception);
2168  if (clone_image == (Image *) NULL)
2169  return(clone_image);
2170  (void) SetImageArtifact(clone_image,"fourier:normalize","inverse");
2171  fft_images=ForwardFourierTransformImage(clone_image,MagickFalse,
2172  exception);
2173  clone_image=DestroyImageList(clone_image);
2174  if (fft_images == (Image *) NULL)
2175  return(fft_images);
2176  /*
2177  Take the complex conjugate of beta image.
2178  */
2179  complex_conjugate=ComplexImages(fft_images,ConjugateComplexOperator,
2180  exception);
2181  fft_images=DestroyImageList(fft_images);
2182  if (complex_conjugate == (Image *) NULL)
2183  return(complex_conjugate);
2184  /*
2185  Take the FFT of the alpha image.
2186  */
2187  clone_image=CloneImage(alpha_image,0,0,MagickTrue,exception);
2188  if (clone_image == (Image *) NULL)
2189  {
2190  complex_conjugate=DestroyImageList(complex_conjugate);
2191  return(clone_image);
2192  }
2193  (void) SetImageArtifact(clone_image,"fourier:normalize","inverse");
2194  fft_images=ForwardFourierTransformImage(clone_image,MagickFalse,exception);
2195  clone_image=DestroyImageList(clone_image);
2196  if (fft_images == (Image *) NULL)
2197  {
2198  complex_conjugate=DestroyImageList(complex_conjugate);
2199  return(fft_images);
2200  }
2201  complex_conjugate->next->next=fft_images;
2202  /*
2203  Do complex multiplication.
2204  */
2205  (void) SetImageArtifact(complex_conjugate,"compose:clamp","false");
2206  complex_multiplication=ComplexImages(complex_conjugate,
2207  MultiplyComplexOperator,exception);
2208  complex_conjugate=DestroyImageList(complex_conjugate);
2209  if (fft_images == (Image *) NULL)
2210  return(fft_images);
2211  /*
2212  Do the IFT and return the cross-correlation result.
2213  */
2214  cross_correlation=InverseFourierTransformImage(complex_multiplication,
2215  complex_multiplication->next,MagickFalse,exception);
2216  complex_multiplication=DestroyImageList(complex_multiplication);
2217  return(cross_correlation);
2218 }
2219 
2220 static Image *NCCDivideImage(const Image *alpha_image,const Image *beta_image,
2221  ExceptionInfo *exception)
2222 {
2223  CacheView
2224  *alpha_view,
2225  *beta_view;
2226 
2227  Image
2228  *divide_image;
2229 
2231  status;
2232 
2233  ssize_t
2234  y;
2235 
2236  /*
2237  Divide one image into another.
2238  */
2239  divide_image=CloneImage(alpha_image,0,0,MagickTrue,exception);
2240  if (divide_image == (Image *) NULL)
2241  return(divide_image);
2242  status=MagickTrue;
2243  alpha_view=AcquireAuthenticCacheView(divide_image,exception);
2244  beta_view=AcquireVirtualCacheView(beta_image,exception);
2245 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2246  #pragma omp parallel for schedule(static) shared(status) \
2247  magick_number_threads(beta_image,divide_image,divide_image->rows,1)
2248 #endif
2249  for (y=0; y < (ssize_t) divide_image->rows; y++)
2250  {
2251  const Quantum
2252  *magick_restrict p;
2253 
2254  Quantum
2255  *magick_restrict q;
2256 
2257  ssize_t
2258  x;
2259 
2260  if (status == MagickFalse)
2261  continue;
2262  p=GetCacheViewVirtualPixels(beta_view,0,y,beta_image->columns,1,
2263  exception);
2264  q=GetCacheViewAuthenticPixels(alpha_view,0,y,divide_image->columns,1,
2265  exception);
2266  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2267  {
2268  status=MagickFalse;
2269  continue;
2270  }
2271  for (x=0; x < (ssize_t) divide_image->columns; x++)
2272  {
2273  ssize_t
2274  i;
2275 
2276  for (i=0; i < (ssize_t) GetPixelChannels(divide_image); i++)
2277  {
2278  PixelChannel channel = GetPixelChannelChannel(divide_image,i);
2279  PixelTrait traits = GetPixelChannelTraits(divide_image,channel);
2280  if ((traits & UpdatePixelTrait) == 0)
2281  continue;
2282  if (fabs(p[i]) >= MagickEpsilon)
2283  q[i]*=PerceptibleReciprocal(QuantumScale*p[i]);
2284  }
2285  p+=GetPixelChannels(beta_image);
2286  q+=GetPixelChannels(divide_image);
2287  }
2288  if (SyncCacheViewAuthenticPixels(alpha_view,exception) == MagickFalse)
2289  status=MagickFalse;
2290  }
2291  beta_view=DestroyCacheView(beta_view);
2292  alpha_view=DestroyCacheView(alpha_view);
2293  if (status == MagickFalse)
2294  divide_image=DestroyImage(divide_image);
2295  return(divide_image);
2296 }
2297 
2298 static MagickBooleanType NCCMaximaImage(const Image *image,double *maxima,
2299  RectangleInfo *offset,ExceptionInfo *exception)
2300 {
2301  CacheView
2302  *image_view;
2303 
2305  status;
2306 
2307  ssize_t
2308  y;
2309 
2310  /*
2311  Identify the maxima value in the image and its location.
2312  */
2313  status=MagickTrue;
2314  *maxima=0.0;
2315  offset->x=0;
2316  offset->y=0;
2317  image_view=AcquireVirtualCacheView(image,exception);
2318  for (y=0; y < (ssize_t) image->rows; y++)
2319  {
2320  const Quantum
2321  *magick_restrict p;
2322 
2323  ssize_t
2324  x;
2325 
2326  if (status == MagickFalse)
2327  continue;
2328  p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2329  if (p == (const Quantum *) NULL)
2330  {
2331  status=MagickFalse;
2332  continue;
2333  }
2334  for (x=0; x < (ssize_t) image->columns; x++)
2335  {
2336  double
2337  sum = 0.0;
2338 
2339  ssize_t
2340  channels = 0,
2341  i;
2342 
2343  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2344  {
2345  PixelChannel channel = GetPixelChannelChannel(image,i);
2346  PixelTrait traits = GetPixelChannelTraits(image,channel);
2347  if ((traits & UpdatePixelTrait) == 0)
2348  continue;
2349  sum+=p[i];
2350  channels++;
2351  }
2352  if ((channels != 0) && ((sum/channels) > *maxima))
2353  {
2354  *maxima=sum/channels;
2355  offset->x=x;
2356  offset->y=y;
2357  }
2358  p+=GetPixelChannels(image);
2359  }
2360  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2361  status=MagickFalse;
2362  }
2363  image_view=DestroyCacheView(image_view);
2364  return(status);
2365 }
2366 
2367 static MagickBooleanType NCCMultiplyImage(Image *image,const double factor,
2368  const ChannelStatistics *channel_statistics,ExceptionInfo *exception)
2369 {
2370  CacheView
2371  *image_view;
2372 
2374  status;
2375 
2376  ssize_t
2377  y;
2378 
2379  /*
2380  Multiply each pixel by a factor.
2381  */
2382  status=MagickTrue;
2383  image_view=AcquireAuthenticCacheView(image,exception);
2384 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2385  #pragma omp parallel for schedule(static) shared(status) \
2386  magick_number_threads(image,image,image->rows,1)
2387 #endif
2388  for (y=0; y < (ssize_t) image->rows; y++)
2389  {
2390  Quantum
2391  *magick_restrict q;
2392 
2393  ssize_t
2394  x;
2395 
2396  if (status == MagickFalse)
2397  continue;
2398  q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
2399  if (q == (Quantum *) NULL)
2400  {
2401  status=MagickFalse;
2402  continue;
2403  }
2404  for (x=0; x < (ssize_t) image->columns; x++)
2405  {
2406  ssize_t
2407  i;
2408 
2409  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2410  {
2411  PixelChannel channel = GetPixelChannelChannel(image,i);
2412  PixelTrait traits = GetPixelChannelTraits(image,channel);
2413  if ((traits & UpdatePixelTrait) == 0)
2414  continue;
2415  if (channel_statistics != (const ChannelStatistics *) NULL)
2416  q[i]*=QuantumScale*channel_statistics[channel].standard_deviation;
2417  q[i]*=factor;
2418  }
2419  q+=GetPixelChannels(image);
2420  }
2421  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2422  status=MagickFalse;
2423  }
2424  image_view=DestroyCacheView(image_view);
2425  return(status);
2426 }
2427 
2428 static Image *NCCSquareImage(const Image *image,ExceptionInfo *exception)
2429 {
2430  CacheView
2431  *image_view;
2432 
2433  Image
2434  *square_image;
2435 
2437  status;
2438 
2439  ssize_t
2440  y;
2441 
2442  /*
2443  Square each pixel in the image.
2444  */
2445  square_image=CloneImage(image,0,0,MagickTrue,exception);
2446  if (square_image == (Image *) NULL)
2447  return(square_image);
2448  status=MagickTrue;
2449  image_view=AcquireAuthenticCacheView(square_image,exception);
2450 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2451  #pragma omp parallel for schedule(static) shared(status) \
2452  magick_number_threads(square_image,square_image,square_image->rows,1)
2453 #endif
2454  for (y=0; y < (ssize_t) square_image->rows; y++)
2455  {
2456  Quantum
2457  *magick_restrict q;
2458 
2459  ssize_t
2460  x;
2461 
2462  if (status == MagickFalse)
2463  continue;
2464  q=GetCacheViewAuthenticPixels(image_view,0,y,square_image->columns,1,
2465  exception);
2466  if (q == (Quantum *) NULL)
2467  {
2468  status=MagickFalse;
2469  continue;
2470  }
2471  for (x=0; x < (ssize_t) square_image->columns; x++)
2472  {
2473  ssize_t
2474  i;
2475 
2476  for (i=0; i < (ssize_t) GetPixelChannels(square_image); i++)
2477  {
2478  PixelChannel channel = GetPixelChannelChannel(square_image,i);
2479  PixelTrait traits = GetPixelChannelTraits(square_image,channel);
2480  if ((traits & UpdatePixelTrait) == 0)
2481  continue;
2482  q[i]*=QuantumScale*q[i];
2483  }
2484  q+=GetPixelChannels(square_image);
2485  }
2486  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2487  status=MagickFalse;
2488  }
2489  image_view=DestroyCacheView(image_view);
2490  if (status == MagickFalse)
2491  square_image=DestroyImage(square_image);
2492  return(square_image);
2493 }
2494 
2495 static Image *NCCSubtractImageMean(const Image *alpha_image,
2496  const Image *beta_image,const ChannelStatistics *channel_statistics,
2497  ExceptionInfo *exception)
2498 {
2499  CacheView
2500  *beta_view,
2501  *image_view;
2502 
2503  Image
2504  *gamma_image;
2505 
2507  status;
2508 
2509  ssize_t
2510  y;
2511 
2512  /*
2513  Subtract the image mean and pad.
2514  */
2515  gamma_image=CloneImage(beta_image,alpha_image->columns,alpha_image->rows,
2516  MagickTrue,exception);
2517  if (gamma_image == (Image *) NULL)
2518  return(gamma_image);
2519  status=MagickTrue;
2520  image_view=AcquireAuthenticCacheView(gamma_image,exception);
2521  beta_view=AcquireVirtualCacheView(beta_image,exception);
2522 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2523  #pragma omp parallel for schedule(static) shared(status) \
2524  magick_number_threads(beta_image,gamma_image,gamma_image->rows,1)
2525 #endif
2526  for (y=0; y < (ssize_t) gamma_image->rows; y++)
2527  {
2528  const Quantum
2529  *magick_restrict p;
2530 
2531  Quantum
2532  *magick_restrict q;
2533 
2534  ssize_t
2535  x;
2536 
2537  if (status == MagickFalse)
2538  continue;
2539  p=GetCacheViewVirtualPixels(beta_view,0,y,beta_image->columns,1,
2540  exception);
2541  q=GetCacheViewAuthenticPixels(image_view,0,y,gamma_image->columns,1,
2542  exception);
2543  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2544  {
2545  status=MagickFalse;
2546  continue;
2547  }
2548  for (x=0; x < (ssize_t) gamma_image->columns; x++)
2549  {
2550  ssize_t
2551  i;
2552 
2553  for (i=0; i < (ssize_t) GetPixelChannels(gamma_image); i++)
2554  {
2555  PixelChannel channel = GetPixelChannelChannel(gamma_image,i);
2556  PixelTrait traits = GetPixelChannelTraits(gamma_image,channel);
2557  if ((traits & UpdatePixelTrait) == 0)
2558  continue;
2559  if ((x >= (ssize_t) beta_image->columns) ||
2560  (y >= (ssize_t) beta_image->rows))
2561  q[i]=(Quantum) 0;
2562  else
2563  q[i]=p[i]-channel_statistics[channel].mean;
2564  }
2565  p+=GetPixelChannels(beta_image);
2566  q+=GetPixelChannels(gamma_image);
2567  }
2568  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2569  status=MagickFalse;
2570  }
2571  beta_view=DestroyCacheView(beta_view);
2572  image_view=DestroyCacheView(image_view);
2573  if (status == MagickFalse)
2574  gamma_image=DestroyImage(gamma_image);
2575  return(gamma_image);
2576 }
2577 
2578 static Image *NCCUnityImage(const Image *alpha_image,const Image *beta_image,
2579  ExceptionInfo *exception)
2580 {
2581  CacheView
2582  *image_view;
2583 
2584  Image
2585  *unity_image;
2586 
2588  status;
2589 
2590  ssize_t
2591  y;
2592 
2593  /*
2594  Create a padded unity image.
2595  */
2596  unity_image=CloneImage(alpha_image,alpha_image->columns,alpha_image->rows,
2597  MagickTrue,exception);
2598  if (unity_image == (Image *) NULL)
2599  return(unity_image);
2600  status=MagickTrue;
2601  image_view=AcquireAuthenticCacheView(unity_image,exception);
2602 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2603  #pragma omp parallel for schedule(static) shared(status) \
2604  magick_number_threads(unity_image,unity_image,unity_image->rows,1)
2605 #endif
2606  for (y=0; y < (ssize_t) unity_image->rows; y++)
2607  {
2608  Quantum
2609  *magick_restrict q;
2610 
2611  ssize_t
2612  x;
2613 
2614  if (status == MagickFalse)
2615  continue;
2616  q=GetCacheViewAuthenticPixels(image_view,0,y,unity_image->columns,1,
2617  exception);
2618  if (q == (Quantum *) NULL)
2619  {
2620  status=MagickFalse;
2621  continue;
2622  }
2623  for (x=0; x < (ssize_t) unity_image->columns; x++)
2624  {
2625  ssize_t
2626  i;
2627 
2628  for (i=0; i < (ssize_t) GetPixelChannels(unity_image); i++)
2629  {
2630  PixelChannel channel = GetPixelChannelChannel(unity_image,i);
2631  PixelTrait traits = GetPixelChannelTraits(unity_image,channel);
2632  if ((traits & UpdatePixelTrait) == 0)
2633  continue;
2634  q[i]=QuantumRange;
2635  if ((x >= (ssize_t) beta_image->columns) ||
2636  (y >= (ssize_t) beta_image->rows))
2637  q[i]=(Quantum) 0;
2638  }
2639  q+=GetPixelChannels(unity_image);
2640  }
2641  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2642  status=MagickFalse;
2643  }
2644  image_view=DestroyCacheView(image_view);
2645  if (status == MagickFalse)
2646  unity_image=DestroyImage(unity_image);
2647  return(unity_image);
2648 }
2649 
2650 static Image *NCCVarianceImage(Image *alpha_image,const Image *beta_image,
2651  ExceptionInfo *exception)
2652 {
2653  CacheView
2654  *beta_view,
2655  *image_view;
2656 
2657  Image
2658  *variance_image;
2659 
2661  status;
2662 
2663  ssize_t
2664  y;
2665 
2666  /*
2667  Compute the variance of the two images.
2668  */
2669  variance_image=CloneImage(alpha_image,0,0,MagickTrue,exception);
2670  if (variance_image == (Image *) NULL)
2671  return(variance_image);
2672  status=MagickTrue;
2673  image_view=AcquireAuthenticCacheView(variance_image,exception);
2674  beta_view=AcquireVirtualCacheView(beta_image,exception);
2675 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2676  #pragma omp parallel for schedule(static) shared(status) \
2677  magick_number_threads(beta_image,variance_image,variance_image->rows,1)
2678 #endif
2679  for (y=0; y < (ssize_t) variance_image->rows; y++)
2680  {
2681  const Quantum
2682  *magick_restrict p;
2683 
2684  Quantum
2685  *magick_restrict q;
2686 
2687  ssize_t
2688  x;
2689 
2690  if (status == MagickFalse)
2691  continue;
2692  p=GetCacheViewVirtualPixels(beta_view,0,y,beta_image->columns,1,
2693  exception);
2694  q=GetCacheViewAuthenticPixels(image_view,0,y,variance_image->columns,1,
2695  exception);
2696  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2697  {
2698  status=MagickFalse;
2699  continue;
2700  }
2701  for (x=0; x < (ssize_t) variance_image->columns; x++)
2702  {
2703  ssize_t
2704  i;
2705 
2706  for (i=0; i < (ssize_t) GetPixelChannels(variance_image); i++)
2707  {
2708  PixelChannel channel = GetPixelChannelChannel(variance_image,i);
2709  PixelTrait traits = GetPixelChannelTraits(variance_image,channel);
2710  if ((traits & UpdatePixelTrait) == 0)
2711  continue;
2712  q[i]=ClampToQuantum((QuantumRange*sqrt(fabs((double) QuantumScale*
2713  (q[i]-p[i])))))/sqrt((double) QuantumRange);
2714  }
2715  p+=GetPixelChannels(beta_image);
2716  q+=GetPixelChannels(variance_image);
2717  }
2718  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2719  status=MagickFalse;
2720  }
2721  beta_view=DestroyCacheView(beta_view);
2722  image_view=DestroyCacheView(image_view);
2723  if (status == MagickFalse)
2724  variance_image=DestroyImage(variance_image);
2725  return(variance_image);
2726 }
2727 
2728 static Image *NCCSimilarityImage(const Image *image,const Image *reference,
2729  const MetricType metric,const double similarity_threshold,
2730  RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
2731 {
2732 #define DestroySimilarityResources() \
2733 { \
2734  if (channel_statistics != (ChannelStatistics *) NULL) \
2735  channel_statistics=(ChannelStatistics *) \
2736  RelinquishMagickMemory(channel_statistics); \
2737  if (beta_image != (Image *) NULL) \
2738  beta_image=DestroyImage(beta_image); \
2739  if (gamma_image != (Image *) NULL) \
2740  gamma_image=DestroyImage(gamma_image); \
2741  if (ncc_image != (Image *) NULL) \
2742  ncc_image=DestroyImage(ncc_image); \
2743  if (normalize_image != (Image *) NULL) \
2744  normalize_image=DestroyImage(normalize_image); \
2745  if (square_image != (Image *) NULL) \
2746  square_image=DestroyImage(square_image); \
2747  if (unity_image != (Image *) NULL) \
2748  unity_image=DestroyImage(unity_image); \
2749 }
2750 #define ThrowSimilarityException() \
2751 { \
2752  DestroySimilarityResources() \
2753  return((Image *) NULL); \
2754 }
2755 
2757  *channel_statistics = (ChannelStatistics *) NULL;
2758 
2759  double
2760  maxima = 0.0;
2761 
2762  Image
2763  *beta_image = (Image *) NULL,
2764  *correlation_image = (Image *) NULL,
2765  *gamma_image = (Image *) NULL,
2766  *ncc_image = (Image *) NULL,
2767  *normalize_image = (Image *) NULL,
2768  *square_image = (Image *) NULL,
2769  *unity_image = (Image *) NULL;
2770 
2772  status;
2773 
2775  geometry;
2776 
2777  /*
2778  Accelerated correlation-based image similary using FFT local statistics.
2779  Contributed by Fred Weinhaus.
2780  */
2781  square_image=NCCSquareImage(image,exception);
2782  if (square_image == (Image *) NULL)
2783  ThrowSimilarityException();
2784  unity_image=NCCUnityImage(image,reference,exception);
2785  if (unity_image == (Image *) NULL)
2786  ThrowSimilarityException();
2787  /*
2788  Compute the cross correlation of the square and unity images.
2789  */
2790  ncc_image=CrossCorrelationImage(square_image,unity_image,exception);
2791  square_image=DestroyImage(square_image); \
2792  if (ncc_image == (Image *) NULL)
2793  ThrowSimilarityException();
2794  status=NCCMultiplyImage(ncc_image,(double) QuantumRange*reference->columns*
2795  reference->rows,(const ChannelStatistics *) NULL,exception);
2796  if (status == MagickFalse)
2797  ThrowSimilarityException();
2798  /*
2799  Compute the cross correlation of the source and unity images.
2800  */
2801  gamma_image=CrossCorrelationImage(image,unity_image,exception);
2802  unity_image=DestroyImage(unity_image);
2803  if (gamma_image == (Image *) NULL)
2804  ThrowSimilarityException();
2805  square_image=NCCSquareImage(gamma_image,exception);
2806  gamma_image=DestroyImage(gamma_image);
2807  status=NCCMultiplyImage(square_image,(double) QuantumRange,
2808  (const ChannelStatistics *) NULL,exception);
2809  if (status == MagickFalse)
2810  ThrowSimilarityException();
2811  /*
2812  Compute the variance of the two images.
2813  */
2814  gamma_image=NCCVarianceImage(ncc_image,square_image,exception);
2815  square_image=DestroyImage(square_image);
2816  ncc_image=DestroyImage(ncc_image);
2817  if (gamma_image == (Image *) NULL)
2818  ThrowSimilarityException();
2819  channel_statistics=GetImageStatistics(reference,exception);
2820  if (channel_statistics == (ChannelStatistics *) NULL)
2821  ThrowSimilarityException();
2822  /*
2823  Subtract the image mean.
2824  */
2825  status=NCCMultiplyImage(gamma_image,1.0,channel_statistics,exception);
2826  if (status == MagickFalse)
2827  ThrowSimilarityException();
2828  normalize_image=NCCSubtractImageMean(image,reference,channel_statistics,
2829  exception);
2830  if (normalize_image == (Image *) NULL)
2831  ThrowSimilarityException();
2832  ncc_image=CrossCorrelationImage(image,normalize_image,exception);
2833  normalize_image=DestroyImage(normalize_image);
2834  if (ncc_image == (Image *) NULL)
2835  ThrowSimilarityException();
2836  /*
2837  Divide the two images.
2838  */
2839  beta_image=NCCDivideImage(ncc_image,gamma_image,exception);
2840  ncc_image=DestroyImage(ncc_image);
2841  gamma_image=DestroyImage(gamma_image);
2842  if (beta_image == (Image *) NULL)
2843  ThrowSimilarityException();
2844  (void) ResetImagePage(beta_image,"0x0+0+0");
2845  SetGeometry(image,&geometry);
2846  geometry.width=image->columns-reference->columns;
2847  geometry.height=image->rows-reference->rows;
2848  /*
2849  Crop padding.
2850  */
2851  correlation_image=CropImage(beta_image,&geometry,exception);
2852  beta_image=DestroyImage(beta_image);
2853  if (correlation_image == (Image *) NULL)
2854  ThrowSimilarityException();
2855  (void) ResetImagePage(correlation_image,"0x0+0+0");
2856  /*
2857  Identify the maxima value in the image and its location.
2858  */
2859  status=GrayscaleImage(correlation_image,AveragePixelIntensityMethod,
2860  exception);
2861  if (status == MagickFalse)
2862  ThrowSimilarityException();
2863  status=NCCMaximaImage(correlation_image,&maxima,offset,exception);
2864  if (status == MagickFalse)
2865  {
2866  correlation_image=DestroyImage(correlation_image);
2867  ThrowSimilarityException();
2868  }
2869  *similarity_metric=1.0-QuantumScale*maxima;
2870  DestroySimilarityResources();
2871  return(correlation_image);
2872 }
2873 #endif
2874 
2875 static double GetSimilarityMetric(const Image *image,const Image *reference,
2876  const MetricType metric,const ssize_t x_offset,const ssize_t y_offset,
2877  ExceptionInfo *exception)
2878 {
2879  double
2880  distortion;
2881 
2882  Image
2883  *similarity_image;
2884 
2886  status;
2887 
2889  geometry;
2890 
2891  SetGeometry(reference,&geometry);
2892  geometry.x=x_offset;
2893  geometry.y=y_offset;
2894  similarity_image=CropImage(image,&geometry,exception);
2895  if (similarity_image == (Image *) NULL)
2896  return(0.0);
2897  distortion=0.0;
2898  status=GetImageDistortion(similarity_image,reference,metric,&distortion,
2899  exception);
2900  similarity_image=DestroyImage(similarity_image);
2901  if (status == MagickFalse)
2902  return(0.0);
2903  return(distortion);
2904 }
2905 
2906 MagickExport Image *SimilarityImage(const Image *image,const Image *reference,
2907  const MetricType metric,const double similarity_threshold,
2908  RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
2909 {
2910 #define SimilarityImageTag "Similarity/Image"
2911 
2912  CacheView
2913  *similarity_view;
2914 
2915  Image
2916  *similarity_image;
2917 
2919  status;
2920 
2922  progress;
2923 
2924  ssize_t
2925  y;
2926 
2927  assert(image != (const Image *) NULL);
2928  assert(image->signature == MagickCoreSignature);
2929  if (image->debug != MagickFalse)
2930  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2931  assert(exception != (ExceptionInfo *) NULL);
2932  assert(exception->signature == MagickCoreSignature);
2933  assert(offset != (RectangleInfo *) NULL);
2934  SetGeometry(reference,offset);
2935  *similarity_metric=MagickMaximumValue;
2936 #if defined(MAGICKCORE_HDRI_SUPPORT) && defined(MAGICKCORE_FFTW_DELEGATE)
2937  {
2938  const char *artifact = GetImageArtifact(image,"compare:accelerate-ncc");
2939  MagickBooleanType accelerate = (artifact != (const char *) NULL) &&
2940  (IsStringTrue(artifact) == MagickFalse) ? MagickFalse : MagickTrue;
2941  if ((accelerate != MagickFalse) &&
2943  {
2944  similarity_image=NCCSimilarityImage(image,reference,metric,
2945  similarity_threshold,offset,similarity_metric,exception);
2946  return(similarity_image);
2947  }
2948  }
2949 #endif
2950  similarity_image=CloneImage(image,image->columns-reference->columns+1,
2951  image->rows-reference->rows+1,MagickTrue,exception);
2952  if (similarity_image == (Image *) NULL)
2953  return((Image *) NULL);
2954  status=SetImageStorageClass(similarity_image,DirectClass,exception);
2955  if (status == MagickFalse)
2956  {
2957  similarity_image=DestroyImage(similarity_image);
2958  return((Image *) NULL);
2959  }
2960  (void) SetImageAlphaChannel(similarity_image,DeactivateAlphaChannel,
2961  exception);
2962  /*
2963  Measure similarity of reference image against image.
2964  */
2965  status=MagickTrue;
2966  progress=0;
2967  similarity_view=AcquireAuthenticCacheView(similarity_image,exception);
2968 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2969  #pragma omp parallel for schedule(static) \
2970  shared(progress,status,similarity_metric) \
2971  magick_number_threads(image,image,image->rows-reference->rows+1,1)
2972 #endif
2973  for (y=0; y < (ssize_t) (image->rows-reference->rows+1); y++)
2974  {
2975  double
2976  similarity;
2977 
2978  Quantum
2979  *magick_restrict q;
2980 
2981  ssize_t
2982  x;
2983 
2984  if (status == MagickFalse)
2985  continue;
2986 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2987  #pragma omp flush(similarity_metric)
2988 #endif
2989  if (*similarity_metric <= similarity_threshold)
2990  continue;
2991  q=GetCacheViewAuthenticPixels(similarity_view,0,y,similarity_image->columns,
2992  1,exception);
2993  if (q == (Quantum *) NULL)
2994  {
2995  status=MagickFalse;
2996  continue;
2997  }
2998  for (x=0; x < (ssize_t) (image->columns-reference->columns+1); x++)
2999  {
3000  ssize_t
3001  i;
3002 
3003 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3004  #pragma omp flush(similarity_metric)
3005 #endif
3006  if (*similarity_metric <= similarity_threshold)
3007  break;
3008  similarity=GetSimilarityMetric(image,reference,metric,x,y,exception);
3009  if ((metric == NormalizedCrossCorrelationErrorMetric) ||
3010  (metric == UndefinedErrorMetric))
3011  similarity=1.0-similarity;
3012 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3013  #pragma omp critical (MagickCore_SimilarityImage)
3014 #endif
3015  if (similarity < *similarity_metric)
3016  {
3017  offset->x=x;
3018  offset->y=y;
3019  *similarity_metric=similarity;
3020  }
3021  if (metric == PerceptualHashErrorMetric)
3022  similarity=MagickMin(0.01*similarity,1.0);
3023  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3024  {
3025  PixelChannel channel = GetPixelChannelChannel(image,i);
3026  PixelTrait traits = GetPixelChannelTraits(image,channel);
3027  PixelTrait similarity_traits=GetPixelChannelTraits(similarity_image,
3028  channel);
3029  if ((traits == UndefinedPixelTrait) ||
3030  (similarity_traits == UndefinedPixelTrait) ||
3031  ((similarity_traits & UpdatePixelTrait) == 0))
3032  continue;
3033  SetPixelChannel(similarity_image,channel,ClampToQuantum(QuantumRange-
3034  QuantumRange*similarity),q);
3035  }
3036  q+=GetPixelChannels(similarity_image);
3037  }
3038  if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
3039  status=MagickFalse;
3040  if (image->progress_monitor != (MagickProgressMonitor) NULL)
3041  {
3043  proceed;
3044 
3045 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3046  #pragma omp atomic
3047 #endif
3048  progress++;
3049  proceed=SetImageProgress(image,SimilarityImageTag,progress,image->rows);
3050  if (proceed == MagickFalse)
3051  status=MagickFalse;
3052  }
3053  }
3054  similarity_view=DestroyCacheView(similarity_view);
3055  if (status == MagickFalse)
3056  similarity_image=DestroyImage(similarity_image);
3057  return(similarity_image);
3058 }
size_t rows
Definition: image.h:172
#define magick_restrict
Definition: MagickCore.h:41
MagickDoubleType MagickRealType
Definition: magick-type.h:124
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:931
double standard_deviation
Definition: statistic.h:35
static MagickBooleanType GetAbsoluteDistortion(const Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c:348
static MagickBooleanType GetPeakSignalToNoiseRatio(const Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c:1241
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:2906
MagickProgressMonitor progress_monitor
Definition: image.h:303
static Quantum GetPixelAlpha(const Image *magick_restrict image, const Quantum *magick_restrict pixel)
double phash[MaximumNumberOfPerceptualColorspaces+1][MaximumNumberOfImageMoments+1]
Definition: statistic.h:79
size_t height
Definition: morphology.h:108
MagickExport ChannelStatistics * GetImageStatistics(const Image *image, ExceptionInfo *exception)
Definition: statistic.c:2029
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:1594
#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:588
#define SimilarityImageTag
MagickExport MagickBooleanType SetImageArtifact(Image *image, const char *artifact, const char *value)
Definition: artifact.c:448
static Quantum GetPixelChannel(const Image *magick_restrict image, const PixelChannel channel, const Quantum *magick_restrict pixel)
MagickExport const char * GetImageArtifact(const Image *image, const char *artifact)
Definition: artifact.c:273
double mean_error_per_pixel
Definition: color.h:79
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:465
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
MagickExport MagickBooleanType GrayscaleImage(Image *image, const PixelIntensityMethod method, ExceptionInfo *exception)
Definition: enhance.c:2467
#define MaximumNumberOfImageMoments
Definition: statistic.h:25
MagickExport MagickBooleanType SetImageColorMetric(Image *image, const Image *reconstruct_image, ExceptionInfo *exception)
Definition: compare.c:2022
#define MagickEpsilon
Definition: magick-type.h:114
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:1525
size_t width
Definition: geometry.h:132
static double GetFuzzyColorDistance(const Image *p, const Image *q)
Definition: color-private.h:69
#define ThrowBinaryException(severity, tag, context)
Definition: log.h:52
ssize_t MagickOffsetType
Definition: magick-type.h:133
static Quantum ClampToQuantum(const MagickRealType quantum)
Definition: quantum.h:85
Definition: image.h:151
MagickExport Image * ForwardFourierTransformImage(const Image *image, const MagickBooleanType modulus, ExceptionInfo *exception)
Definition: fourier.c:905
MagickExport Image * CompareImages(Image *image, const Image *reconstruct_image, const MetricType metric, double *distortion, ExceptionInfo *exception)
Definition: compare.c:112
MagickExport Image * CropImage(const Image *image, const RectangleInfo *geometry, ExceptionInfo *exception)
Definition: transform.c:537
MagickExport KernelInfo * AcquireKernelInfo(const char *kernel_string, ExceptionInfo *exception)
Definition: morphology.c:485
#define MagickCoreSignature
double normalized_mean_error
Definition: color.h:79
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:3201
MagickExport MagickBooleanType SetImageAlphaChannel(Image *image, const AlphaChannelOption alpha_type, ExceptionInfo *exception)
Definition: channel.c:973
MagickBooleanType
Definition: magick-type.h:161
unsigned int MagickStatusType
Definition: magick-type.h:125
static double PerceptibleReciprocal(const double x)
MagickExport void * AcquireQuantumMemory(const size_t count, const size_t quantum)
Definition: memory.c:665
#define MagickPathExtent
static MagickBooleanType GetMeanErrorPerPixel(Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c:708
MagickExport MagickBooleanType IsStringTrue(const char *value)
Definition: string.c:1386
#define INFINITY
Definition: magick-type.h:187
MagickExport int GetMagickPrecision(void)
Definition: magick.c:942
#define MagickMaximumValue
Definition: magick-type.h:115
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:1662
size_t width
Definition: morphology.h:108
size_t signature
Definition: image.h:354
#define QuantumScale
Definition: magick-type.h:119
size_t columns
Definition: image.h:172
ssize_t x
Definition: geometry.h:136
MagickBooleanType(* MagickProgressMonitor)(const char *, const MagickOffsetType, const MagickSizeType, void *)
Definition: monitor.h:26
struct _Image * next
Definition: image.h:348
size_t height
Definition: geometry.h:132
MagickExport MagickBooleanType QueryColorCompliance(const char *name, const ComplianceType compliance, PixelInfo *color, ExceptionInfo *exception)
Definition: color.c:2267
MagickExport MagickBooleanType SetImageStorageClass(Image *image, const ClassType storage_class, ExceptionInfo *exception)
Definition: image.c:2616
static MagickBooleanType GetPeakAbsoluteDistortion(const Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c:1126
MagickExport Image * DestroyImageList(Image *images)
Definition: list.c:477
MagickExport MagickBooleanType ResetImagePage(Image *image, const char *page)
Definition: image.c:2170
static MagickBooleanType GetMeanSquaredDistortion(const Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c:811
PixelChannel
Definition: pixel.h:70
#define MagickMax(x, y)
Definition: image-private.h:38
static MagickBooleanType GetRootMeanSquaredDistortion(const Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c:1338
static size_t GetPixelChannels(const Image *magick_restrict image)
MagickExport MagickBooleanType IsImagesEqual(const Image *image, const Image *reconstruct_image, ExceptionInfo *exception)
Definition: compare.c:1902
#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:2875
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:79
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:1610
unsigned short Quantum
Definition: magick-type.h:86
#define SSIML
#define Log10Epsilon
MagickExport Image * ExtentImage(const Image *image, const RectangleInfo *geometry, ExceptionInfo *exception)
Definition: transform.c:1119
#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:1259
#define MagickMin(x, y)
Definition: image-private.h:39
MagickExport void SetGeometry(const Image *image, RectangleInfo *geometry)
Definition: geometry.c:1714
static double StringToDouble(const char *magick_restrict string, char *magick_restrict *sentinal)
static size_t GetImageChannels(const Image *image)
Definition: image-private.h:71
MagickExport void * RelinquishMagickMemory(void *memory)
Definition: memory.c:1162
#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:1758
#define MagickExport
MagickExport MagickBooleanType SyncCacheViewAuthenticPixels(CacheView *magick_restrict cache_view, ExceptionInfo *exception)
Definition: cache-view.c:1100
ssize_t y
Definition: geometry.h:136
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:1232
static MagickBooleanType GetFuzzDistortion(const Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c:467
PixelTrait
Definition: pixel.h:137
MagickExport ChannelPerceptualHash * GetImagePerceptualHash(const Image *image, ExceptionInfo *exception)
Definition: statistic.c:1736
MagickExport Image * DestroyImage(Image *image)
Definition: image.c:1179
MagickExport Image * CloneImage(const Image *image, const size_t columns, const size_t rows, const MagickBooleanType detach, ExceptionInfo *exception)
Definition: image.c:787
#define QuantumRange
Definition: magick-type.h:87
MagickExport Image * ComplexImages(const Image *images, const ComplexOperator op, ExceptionInfo *exception)
Definition: fourier.c:129
MagickExport MagickBooleanType SetImageProgress(const Image *image, const char *tag, const MagickOffsetType offset, const MagickSizeType extent)
Definition: monitor.c:136
MagickRealType * values
Definition: morphology.h:116
MagickBooleanType debug
Definition: image.h:334
static MagickBooleanType GetStructuralSimilarityDistortion(const Image *image, const Image *reconstruct_image, double *distortion, ExceptionInfo *exception)
Definition: compare.c:1353
MagickExport Image * InverseFourierTransformImage(const Image *magnitude_image, const Image *phase_image, const MagickBooleanType modulus, ExceptionInfo *exception)
Definition: fourier.c:1483
#define SSIMK2