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