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