MagickCore  7.1.2-29
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/license/ %
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/compare-private.h"
55 #include "MagickCore/composite-private.h"
56 #include "MagickCore/constitute.h"
57 #include "MagickCore/distort.h"
58 #include "MagickCore/exception-private.h"
59 #include "MagickCore/enhance.h"
60 #include "MagickCore/fourier.h"
61 #include "MagickCore/geometry.h"
62 #include "MagickCore/image-private.h"
63 #include "MagickCore/list.h"
64 #include "MagickCore/log.h"
65 #include "MagickCore/memory_.h"
66 #include "MagickCore/monitor.h"
67 #include "MagickCore/monitor-private.h"
68 #include "MagickCore/option.h"
69 #include "MagickCore/pixel-accessor.h"
70 #include "MagickCore/property.h"
71 #include "MagickCore/registry.h"
72 #include "MagickCore/resource_.h"
73 #include "MagickCore/string_.h"
74 #include "MagickCore/statistic.h"
75 #include "MagickCore/statistic-private.h"
76 #include "MagickCore/string-private.h"
77 #include "MagickCore/thread-private.h"
78 #include "MagickCore/threshold.h"
79 #include "MagickCore/transform.h"
80 #include "MagickCore/utility.h"
81 #include "MagickCore/version.h"
82 ␌
83 /*
84 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
85 % %
86 % %
87 % %
88 % C o m p a r e I m a g e s %
89 % %
90 % %
91 % %
92 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
93 %
94 % CompareImages() compares one or more pixel channels of an image to a
95 % reconstructed image and returns the difference image.
96 %
97 % The format of the CompareImages method is:
98 %
99 % Image *CompareImages(const Image *image,const Image *reconstruct_image,
100 % const MetricType metric,double *distortion,ExceptionInfo *exception)
101 %
102 % A description of each parameter follows:
103 %
104 % o image: the image.
105 %
106 % o reconstruct_image: the reconstruction image.
107 %
108 % o metric: the metric.
109 %
110 % o distortion: the computed distortion between the images.
111 %
112 % o exception: return any errors or warnings in this structure.
113 %
114 */
115 MagickExport Image *CompareImages(Image *image,const Image *reconstruct_image,
116  const MetricType metric,double *distortion,ExceptionInfo *exception)
117 {
118  CacheView
119  *highlight_view,
120  *image_view,
121  *reconstruct_view;
122 
123  const char
124  *artifact;
125 
126  Image
127  *clone_image,
128  *difference_image,
129  *highlight_image;
130 
131  MagickBooleanType
132  status = MagickTrue;
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  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
162  SetGeometry(image,&geometry);
163  geometry.width=columns;
164  geometry.height=rows;
165  clone_image=CloneImage(image,0,0,MagickTrue,exception);
166  if (clone_image == (Image *) NULL)
167  return((Image *) NULL);
168  (void) SetImageMask(clone_image,ReadPixelMask,(Image *) NULL,exception);
169  difference_image=ExtentImage(clone_image,&geometry,exception);
170  clone_image=DestroyImage(clone_image);
171  if (difference_image == (Image *) NULL)
172  return((Image *) NULL);
173  (void) ResetImagePage(difference_image,"0x0+0+0");
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  image_view=AcquireVirtualCacheView(image,exception);
206  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
207  highlight_view=AcquireAuthenticCacheView(highlight_image,exception);
208 #if defined(MAGICKCORE_OPENMP_SUPPORT)
209  #pragma omp parallel for schedule(static) shared(status) \
210  magick_number_threads(image,highlight_image,rows,1)
211 #endif
212  for (y=0; y < (ssize_t) rows; y++)
213  {
214  const Quantum
215  *magick_restrict p,
216  *magick_restrict q;
217 
218  MagickBooleanType
219  sync;
220 
221  Quantum
222  *magick_restrict r;
223 
224  ssize_t
225  x;
226 
227  if (status == MagickFalse)
228  continue;
229  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
230  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
231  r=QueueCacheViewAuthenticPixels(highlight_view,0,y,columns,1,exception);
232  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL) ||
233  (r == (Quantum *) NULL))
234  {
235  status=MagickFalse;
236  continue;
237  }
238  for (x=0; x < (ssize_t) columns; x++)
239  {
240  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
241  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
242  {
243  SetPixelViaPixelInfo(highlight_image,&masklight,r);
244  p+=(ptrdiff_t) GetPixelChannels(image);
245  q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
246  r+=(ptrdiff_t) GetPixelChannels(highlight_image);
247  continue;
248  }
249  if (IsFuzzyEquivalencePixel(image,p,reconstruct_image,q) == MagickFalse)
250  SetPixelViaPixelInfo(highlight_image,&highlight,r);
251  else
252  SetPixelViaPixelInfo(highlight_image,&lowlight,r);
253  p+=(ptrdiff_t) GetPixelChannels(image);
254  q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
255  r+=(ptrdiff_t) GetPixelChannels(highlight_image);
256  }
257  sync=SyncCacheViewAuthenticPixels(highlight_view,exception);
258  if (sync == MagickFalse)
259  status=MagickFalse;
260  }
261  highlight_view=DestroyCacheView(highlight_view);
262  reconstruct_view=DestroyCacheView(reconstruct_view);
263  image_view=DestroyCacheView(image_view);
264  if ((status != MagickFalse) && (difference_image != (Image *) NULL))
265  status=CompositeImage(difference_image,highlight_image,image->compose,
266  MagickTrue,0,0,exception);
267  highlight_image=DestroyImage(highlight_image);
268  if (status == MagickFalse)
269  difference_image=DestroyImage(difference_image);
270  return(difference_image);
271 }
272 ␌
273 /*
274 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
275 % %
276 % %
277 % %
278 % G e t I m a g e D i s t o r t i o n %
279 % %
280 % %
281 % %
282 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
283 %
284 % GetImageDistortion() compares one or more pixel channels of an image to a
285 % reconstructed image and returns the specified distortion metric.
286 %
287 % The format of the GetImageDistortion method is:
288 %
289 % MagickBooleanType GetImageDistortion(const Image *image,
290 % const Image *reconstruct_image,const MetricType metric,
291 % double *distortion,ExceptionInfo *exception)
292 %
293 % A description of each parameter follows:
294 %
295 % o image: the image.
296 %
297 % o reconstruct_image: the reconstruction image.
298 %
299 % o metric: the metric.
300 %
301 % o distortion: the computed distortion between the images.
302 %
303 % o exception: return any errors or warnings in this structure.
304 %
305 */
306 
307 static MagickBooleanType GetAESimilarity(const Image *image,
308  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
309 {
310  CacheView
311  *image_view,
312  *reconstruct_view;
313 
314  double
315  area,
316  fuzz;
317 
318  MagickBooleanType
319  status = MagickTrue;
320 
321  size_t
322  columns,
323  rows;
324 
325  ssize_t
326  channels = 0,
327  k,
328  y;
329 
330  /*
331  Compute the absolute error similarity.
332  */
333  fuzz=GetFuzzyColorDistance(image,reconstruct_image);
334  (void) memset(similarity,0,(MaxPixelChannels+1)*sizeof(*similarity));
335  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
336  image_view=AcquireVirtualCacheView(image,exception);
337  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
338 #if defined(MAGICKCORE_OPENMP_SUPPORT)
339  #pragma omp parallel for schedule(static) shared(similarity,status) \
340  magick_number_threads(image,image,rows,1)
341 #endif
342  for (y=0; y < (ssize_t) rows; y++)
343  {
344  const Quantum
345  *magick_restrict p,
346  *magick_restrict q;
347 
348  double
349  channel_similarity[MaxPixelChannels+1] = { 0.0 };
350 
351  ssize_t
352  x;
353 
354  if (status == MagickFalse)
355  continue;
356  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
357  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
358  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
359  {
360  status=MagickFalse;
361  continue;
362  }
363  for (x=0; x < (ssize_t) columns; x++)
364  {
365  double
366  Da,
367  Sa;
368 
369  ssize_t
370  i;
371 
372  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
373  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
374  {
375  p+=(ptrdiff_t) GetPixelChannels(image);
376  q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
377  continue;
378  }
379  Sa=QuantumScale*(double) GetPixelAlpha(image,p);
380  Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
381  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
382  {
383  double
384  error;
385 
386  PixelChannel channel = GetPixelChannelChannel(image,i);
387  PixelTrait traits = GetPixelChannelTraits(image,channel);
388  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
389  channel);
390  if (((traits & UpdatePixelTrait) == 0) ||
391  ((reconstruct_traits & UpdatePixelTrait) == 0))
392  continue;
393  if (channel == AlphaPixelChannel)
394  error=(double) p[i]-(double) GetPixelChannel(reconstruct_image,
395  channel,q);
396  else
397  error=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
398  if (MagickSafeSignificantError(error*error,fuzz) != MagickFalse)
399  {
400  double ae = fabs(QuantumScale*error);
401  channel_similarity[i]+=ae;
402  channel_similarity[CompositePixelChannel]+=ae;
403  }
404  }
405  p+=(ptrdiff_t) GetPixelChannels(image);
406  q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
407  }
408 #if defined(MAGICKCORE_OPENMP_SUPPORT)
409  #pragma omp critical (MagickCore_GetAESimilarity)
410 #endif
411  {
412  ssize_t
413  j;
414 
415  for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
416  {
417  PixelChannel channel = GetPixelChannelChannel(image,j);
418  PixelTrait traits = GetPixelChannelTraits(image,channel);
419  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
420  channel);
421  if (((traits & UpdatePixelTrait) == 0) ||
422  ((reconstruct_traits & UpdatePixelTrait) == 0))
423  continue;
424  similarity[j]+=channel_similarity[j];
425  }
426  similarity[CompositePixelChannel]+=
427  channel_similarity[CompositePixelChannel];
428  }
429  }
430  reconstruct_view=DestroyCacheView(reconstruct_view);
431  image_view=DestroyCacheView(image_view);
432  area=MagickSafeReciprocal((double) columns*rows);
433  for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
434  {
435  PixelChannel channel = GetPixelChannelChannel(image,k);
436  PixelTrait traits = GetPixelChannelTraits(image,channel);
437  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
438  channel);
439  if (((traits & UpdatePixelTrait) == 0) ||
440  ((reconstruct_traits & UpdatePixelTrait) == 0))
441  continue;
442  similarity[k]*=area;
443  channels++;
444  }
445  similarity[CompositePixelChannel]*=area;
446  if (channels != 0)
447  similarity[CompositePixelChannel]/=(double) channels;
448  return(status);
449 }
450 
451 static MagickBooleanType GetDPCSimilarity(const Image *image,
452  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
453 {
454 #define SimilarityImageTag "Similarity/Image"
455 
456  CacheView
457  *image_view,
458  *reconstruct_view;
459 
460  ChannelStatistics
461  *image_statistics,
462  *reconstruct_statistics;
463 
464  double
465  norm[MaxPixelChannels+1] = { 0.0 },
466  reconstruct_norm[MaxPixelChannels+1] = { 0.0 };
467 
468  MagickBooleanType
469  status = MagickTrue;
470 
471  MagickOffsetType
472  progress = 0;
473 
474  size_t
475  columns,
476  rows;
477 
478  ssize_t
479  k,
480  y;
481 
482  /*
483  Compute the dot product correlation similarity.
484  */
485  image_statistics=GetImageStatistics(image,exception);
486  reconstruct_statistics=GetImageStatistics(reconstruct_image,exception);
487  if ((image_statistics == (ChannelStatistics *) NULL) ||
488  (reconstruct_statistics == (ChannelStatistics *) NULL))
489  {
490  if (image_statistics != (ChannelStatistics *) NULL)
491  image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
492  image_statistics);
493  if (reconstruct_statistics != (ChannelStatistics *) NULL)
494  reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
495  reconstruct_statistics);
496  return(MagickFalse);
497  }
498  (void) memset(similarity,0,(MaxPixelChannels+1)*sizeof(*similarity));
499  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
500  image_view=AcquireVirtualCacheView(image,exception);
501  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
502 #if defined(MAGICKCORE_OPENMP_SUPPORT)
503  #pragma omp parallel for schedule(static) shared(norm,reconstruct_norm,similarity,status) \
504  magick_number_threads(image,image,rows,1)
505 #endif
506  for (y=0; y < (ssize_t) rows; y++)
507  {
508  const Quantum
509  *magick_restrict p,
510  *magick_restrict q;
511 
512  double
513  channel_norm[MaxPixelChannels+1] = { 0.0 },
514  channel_reconstruct_norm[MaxPixelChannels+1] = { 0.0 },
515  channel_similarity[MaxPixelChannels+1] = { 0.0 };
516 
517  ssize_t
518  x;
519 
520  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
521  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
522  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
523  {
524  status=MagickFalse;
525  continue;
526  }
527  for (x=0; x < (ssize_t) columns; x++)
528  {
529  double
530  Da,
531  Sa;
532 
533  ssize_t
534  i;
535 
536  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
537  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
538  {
539  p+=(ptrdiff_t) GetPixelChannels(image);
540  q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
541  continue;
542  }
543  Sa=QuantumScale*(double) GetPixelAlpha(image,p);
544  Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
545  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
546  {
547  double
548  alpha,
549  beta;
550 
551  PixelChannel channel = GetPixelChannelChannel(image,i);
552  PixelTrait traits = GetPixelChannelTraits(image,channel);
553  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
554  channel);
555  if (((traits & UpdatePixelTrait) == 0) ||
556  ((reconstruct_traits & UpdatePixelTrait) == 0))
557  continue;
558  if (channel == AlphaPixelChannel)
559  {
560  alpha=QuantumScale*((double) p[i]-image_statistics[channel].mean);
561  beta=QuantumScale*((double) GetPixelChannel(reconstruct_image,
562  channel,q)-reconstruct_statistics[channel].mean);
563  }
564  else
565  {
566  alpha=QuantumScale*(Sa*p[i]-image_statistics[channel].mean);
567  beta=QuantumScale*(Da*GetPixelChannel(reconstruct_image,channel,
568  q)-reconstruct_statistics[channel].mean);
569  }
570  channel_similarity[i]+=alpha*beta;
571  channel_norm[i]+=alpha*alpha;
572  channel_reconstruct_norm[i]+=beta*beta;
573  }
574  p+=(ptrdiff_t) GetPixelChannels(image);
575  q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
576  }
577 #if defined(MAGICKCORE_OPENMP_SUPPORT)
578  #pragma omp critical (MagickCore_GetDPCSimilarity)
579 #endif
580  {
581  ssize_t
582  j;
583 
584  for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
585  {
586  PixelChannel channel = GetPixelChannelChannel(image,j);
587  PixelTrait traits = GetPixelChannelTraits(image,channel);
588  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
589  channel);
590  if (((traits & UpdatePixelTrait) == 0) ||
591  ((reconstruct_traits & UpdatePixelTrait) == 0))
592  continue;
593  similarity[j]+=channel_similarity[j];
594  similarity[CompositePixelChannel]+=channel_similarity[j];
595  norm[j]+=channel_norm[j];
596  norm[CompositePixelChannel]+=channel_norm[j];
597  reconstruct_norm[j]+=channel_reconstruct_norm[j];
598  reconstruct_norm[CompositePixelChannel]+=channel_reconstruct_norm[j];
599  }
600  }
601  if (image->progress_monitor != (MagickProgressMonitor) NULL)
602  {
603  MagickBooleanType
604  proceed;
605 
606 #if defined(MAGICKCORE_OPENMP_SUPPORT)
607  #pragma omp atomic
608 #endif
609  progress++;
610  proceed=SetImageProgress(image,SimilarityImageTag,progress,rows);
611  if (proceed == MagickFalse)
612  {
613  status=MagickFalse;
614  continue;
615  }
616  }
617  }
618  reconstruct_view=DestroyCacheView(reconstruct_view);
619  image_view=DestroyCacheView(image_view);
620  /*
621  Compute dot product correlation: divide by mean.
622  */
623  for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
624  {
625  PixelChannel channel = GetPixelChannelChannel(image,k);
626  PixelTrait traits = GetPixelChannelTraits(image,channel);
627  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
628  channel);
629  if (((traits & UpdatePixelTrait) == 0) ||
630  ((reconstruct_traits & UpdatePixelTrait) == 0))
631  continue;
632  similarity[k]*=MagickSafeReciprocal(sqrt(norm[k]*reconstruct_norm[k]));
633  }
634  similarity[CompositePixelChannel]*=MagickSafeReciprocal(sqrt(
635  norm[CompositePixelChannel]*reconstruct_norm[CompositePixelChannel]));
636  /*
637  Free resources.
638  */
639  reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
640  reconstruct_statistics);
641  image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
642  image_statistics);
643  return(status);
644 }
645 
646 static MagickBooleanType GetFUZZSimilarity(const Image *image,
647  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
648 {
649  CacheView
650  *image_view,
651  *reconstruct_view;
652 
653  double
654  area = 0.0,
655  fuzz = 0.0;
656 
657  MagickBooleanType
658  status = MagickTrue;
659 
660  size_t
661  columns,
662  rows;
663 
664  ssize_t
665  k,
666  y;
667 
668  /*
669  Compute the MSE similarity within tolerance (fuzz).
670  */
671  fuzz=GetFuzzyColorDistance(image,reconstruct_image);
672  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
673  image_view=AcquireVirtualCacheView(image,exception);
674  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
675 #if defined(MAGICKCORE_OPENMP_SUPPORT)
676  #pragma omp parallel for schedule(static) shared(area,similarity,status) \
677  magick_number_threads(image,image,rows,1)
678 #endif
679  for (y=0; y < (ssize_t) rows; y++)
680  {
681  const Quantum
682  *magick_restrict p,
683  *magick_restrict q;
684 
685  double
686  channel_area = 0.0,
687  channel_similarity[MaxPixelChannels+1] = { 0.0 };
688 
689  ssize_t
690  x;
691 
692  if (status == MagickFalse)
693  continue;
694  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
695  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
696  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
697  {
698  status=MagickFalse;
699  continue;
700  }
701  for (x=0; x < (ssize_t) columns; x++)
702  {
703  double
704  Da,
705  Sa;
706 
707  ssize_t
708  i;
709 
710  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
711  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
712  {
713  p+=(ptrdiff_t) GetPixelChannels(image);
714  q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
715  continue;
716  }
717  Sa=QuantumScale*(double) GetPixelAlpha(image,p);
718  Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
719  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
720  {
721  double
722  error;
723 
724  PixelChannel channel = GetPixelChannelChannel(image,i);
725  PixelTrait traits = GetPixelChannelTraits(image,channel);
726  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
727  channel);
728  if (((traits & UpdatePixelTrait) == 0) ||
729  ((reconstruct_traits & UpdatePixelTrait) == 0))
730  continue;
731  if (channel == AlphaPixelChannel)
732  error=(double) p[i]-(double) GetPixelChannel(reconstruct_image,
733  channel,q);
734  else
735  error=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
736  if (MagickSafeSignificantError(error*error,fuzz) != MagickFalse)
737  {
738  channel_similarity[i]+=QuantumScale*error*QuantumScale*error;
739  channel_similarity[CompositePixelChannel]+=QuantumScale*error*
740  QuantumScale*error;
741  channel_area++;
742  }
743  }
744  p+=(ptrdiff_t) GetPixelChannels(image);
745  q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
746  }
747 #if defined(MAGICKCORE_OPENMP_SUPPORT)
748  #pragma omp critical (MagickCore_GetFUZZSimilarity)
749 #endif
750  {
751  ssize_t
752  j;
753 
754  area+=channel_area;
755  for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
756  {
757  PixelChannel channel = GetPixelChannelChannel(image,j);
758  PixelTrait traits = GetPixelChannelTraits(image,channel);
759  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
760  channel);
761  if (((traits & UpdatePixelTrait) == 0) ||
762  ((reconstruct_traits & UpdatePixelTrait) == 0))
763  continue;
764  similarity[j]+=channel_similarity[j];
765  }
766  similarity[CompositePixelChannel]+=
767  channel_similarity[CompositePixelChannel];
768  }
769  }
770  reconstruct_view=DestroyCacheView(reconstruct_view);
771  image_view=DestroyCacheView(image_view);
772  area=MagickSafeReciprocal(area);
773  for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
774  {
775  PixelChannel channel = GetPixelChannelChannel(image,k);
776  PixelTrait traits = GetPixelChannelTraits(image,channel);
777  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
778  channel);
779  if (((traits & UpdatePixelTrait) == 0) ||
780  ((reconstruct_traits & UpdatePixelTrait) == 0))
781  continue;
782  similarity[k]*=area;
783  }
784  similarity[CompositePixelChannel]*=area;
785  return(status);
786 }
787 
788 static MagickBooleanType GetMAESimilarity(const Image *image,
789  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
790 {
791  CacheView
792  *image_view,
793  *reconstruct_view;
794 
795  double
796  area = 0.0;
797 
798  MagickBooleanType
799  status = MagickTrue;
800 
801  size_t
802  columns,
803  rows;
804 
805  ssize_t
806  channels = 0,
807  k,
808  y;
809 
810  /*
811  Compute the mean absolute error similarity.
812  */
813  (void) memset(similarity,0,(MaxPixelChannels+1)*sizeof(*similarity));
814  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
815  image_view=AcquireVirtualCacheView(image,exception);
816  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
817 #if defined(MAGICKCORE_OPENMP_SUPPORT)
818  #pragma omp parallel for schedule(static) shared(area,similarity,status) \
819  magick_number_threads(image,image,rows,1)
820 #endif
821  for (y=0; y < (ssize_t) rows; y++)
822  {
823  const Quantum
824  *magick_restrict p,
825  *magick_restrict q;
826 
827  double
828  channel_area = 0.0,
829  channel_similarity[MaxPixelChannels+1] = { 0.0 };
830 
831  ssize_t
832  x;
833 
834  if (status == MagickFalse)
835  continue;
836  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
837  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
838  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
839  {
840  status=MagickFalse;
841  continue;
842  }
843  for (x=0; x < (ssize_t) columns; x++)
844  {
845  double
846  Da,
847  Sa;
848 
849  ssize_t
850  i;
851 
852  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
853  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
854  {
855  p+=(ptrdiff_t) GetPixelChannels(image);
856  q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
857  continue;
858  }
859  Sa=QuantumScale*(double) GetPixelAlpha(image,p);
860  Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
861  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
862  {
863  double
864  error;
865 
866  PixelChannel channel = GetPixelChannelChannel(image,i);
867  PixelTrait traits = GetPixelChannelTraits(image,channel);
868  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
869  channel);
870  if (((traits & UpdatePixelTrait) == 0) ||
871  ((reconstruct_traits & UpdatePixelTrait) == 0))
872  continue;
873  if (channel == AlphaPixelChannel)
874  error=QuantumScale*fabs((double) p[i]-(double) GetPixelChannel(
875  reconstruct_image,channel,q));
876  else
877  error=QuantumScale*fabs(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
878  channel,q));
879  channel_similarity[i]+=error;
880  channel_similarity[CompositePixelChannel]+=error;
881  }
882  channel_area++;
883  p+=(ptrdiff_t) GetPixelChannels(image);
884  q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
885  }
886 #if defined(MAGICKCORE_OPENMP_SUPPORT)
887  #pragma omp critical (MagickCore_GetMAESimilarity)
888 #endif
889  {
890  ssize_t
891  j;
892 
893  area+=channel_area;
894  for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
895  {
896  PixelChannel channel = GetPixelChannelChannel(image,j);
897  PixelTrait traits = GetPixelChannelTraits(image,channel);
898  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
899  channel);
900  if (((traits & UpdatePixelTrait) == 0) ||
901  ((reconstruct_traits & UpdatePixelTrait) == 0))
902  continue;
903  similarity[j]+=channel_similarity[j];
904  }
905  similarity[CompositePixelChannel]+=
906  channel_similarity[CompositePixelChannel];
907  }
908  }
909  reconstruct_view=DestroyCacheView(reconstruct_view);
910  image_view=DestroyCacheView(image_view);
911  area=MagickSafeReciprocal(area);
912  for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
913  {
914  PixelChannel channel = GetPixelChannelChannel(image,k);
915  PixelTrait traits = GetPixelChannelTraits(image,channel);
916  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
917  channel);
918  if (((traits & UpdatePixelTrait) == 0) ||
919  ((reconstruct_traits & UpdatePixelTrait) == 0))
920  continue;
921  similarity[k]*=area;
922  channels++;
923  }
924  similarity[CompositePixelChannel]*=area;
925  if (channels != 0)
926  similarity[CompositePixelChannel]/=(double) channels;
927  return(status);
928 }
929 
930 static MagickBooleanType GetMEPPSimilarity(Image *image,
931  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
932 {
933  CacheView
934  *image_view,
935  *reconstruct_view;
936 
937  double
938  area = 0.0,
939  maximum_error = -MagickMaximumValue,
940  mean_error = 0.0;
941 
942  MagickBooleanType
943  status = MagickTrue;
944 
945  size_t
946  columns,
947  rows;
948 
949  ssize_t
950  channels = 0,
951  k,
952  y;
953 
954  /*
955  Compute the mean error per pixel similarity.
956  */
957  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
958  image_view=AcquireVirtualCacheView(image,exception);
959  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
960 #if defined(MAGICKCORE_OPENMP_SUPPORT)
961  #pragma omp parallel for schedule(static) shared(area,similarity,maximum_error,mean_error,status) \
962  magick_number_threads(image,image,rows,1)
963 #endif
964  for (y=0; y < (ssize_t) rows; y++)
965  {
966  const Quantum
967  *magick_restrict p,
968  *magick_restrict q;
969 
970  double
971  channel_area = 0.0,
972  channel_similarity[MaxPixelChannels+1] = { 0.0 },
973  channel_maximum_error = maximum_error,
974  channel_mean_error = 0.0;
975 
976  ssize_t
977  x;
978 
979  if (status == MagickFalse)
980  continue;
981  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
982  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
983  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
984  {
985  status=MagickFalse;
986  continue;
987  }
988  for (x=0; x < (ssize_t) columns; x++)
989  {
990  double
991  Da,
992  Sa;
993 
994  ssize_t
995  i;
996 
997  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
998  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
999  {
1000  p+=(ptrdiff_t) GetPixelChannels(image);
1001  q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1002  continue;
1003  }
1004  Sa=QuantumScale*(double) GetPixelAlpha(image,p);
1005  Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
1006  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1007  {
1008  double
1009  error;
1010 
1011  PixelChannel channel = GetPixelChannelChannel(image,i);
1012  PixelTrait traits = GetPixelChannelTraits(image,channel);
1013  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1014  channel);
1015  if (((traits & UpdatePixelTrait) == 0) ||
1016  ((reconstruct_traits & UpdatePixelTrait) == 0))
1017  continue;
1018  if (channel == AlphaPixelChannel)
1019  error=QuantumScale*fabs((double) p[i]-(double) GetPixelChannel(
1020  reconstruct_image,channel,q));
1021  else
1022  error=QuantumScale*fabs(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
1023  channel,q));
1024  channel_similarity[i]+=error;
1025  channel_similarity[CompositePixelChannel]+=error;
1026  channel_mean_error+=error*error;
1027  if (error > channel_maximum_error)
1028  channel_maximum_error=error;
1029  }
1030  channel_area++;
1031  p+=(ptrdiff_t) GetPixelChannels(image);
1032  q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1033  }
1034 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1035  #pragma omp critical (MagickCore_GetMEPPSimilarity)
1036 #endif
1037  {
1038  ssize_t
1039  j;
1040 
1041  area+=channel_area;
1042  for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
1043  {
1044  PixelChannel channel = GetPixelChannelChannel(image,j);
1045  PixelTrait traits = GetPixelChannelTraits(image,channel);
1046  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1047  channel);
1048  if (((traits & UpdatePixelTrait) == 0) ||
1049  ((reconstruct_traits & UpdatePixelTrait) == 0))
1050  continue;
1051  similarity[j]+=channel_similarity[j];
1052  }
1053  similarity[CompositePixelChannel]+=
1054  channel_similarity[CompositePixelChannel];
1055  mean_error+=channel_mean_error;
1056  if (channel_maximum_error > maximum_error)
1057  maximum_error=channel_maximum_error;
1058  }
1059  }
1060  reconstruct_view=DestroyCacheView(reconstruct_view);
1061  image_view=DestroyCacheView(image_view);
1062  area=MagickSafeReciprocal(area);
1063  for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
1064  {
1065  PixelChannel channel = GetPixelChannelChannel(image,k);
1066  PixelTrait traits = GetPixelChannelTraits(image,channel);
1067  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1068  channel);
1069  if (((traits & UpdatePixelTrait) == 0) ||
1070  ((reconstruct_traits & UpdatePixelTrait) == 0))
1071  continue;
1072  similarity[k]*=area;
1073  channels++;
1074  }
1075  similarity[CompositePixelChannel]*=area;
1076  if (channels != 0)
1077  similarity[CompositePixelChannel]/=(double) channels;
1078  image->error.mean_error_per_pixel=QuantumRange*
1079  similarity[CompositePixelChannel];
1080  image->error.normalized_mean_error=mean_error*area;
1081  image->error.normalized_maximum_error=maximum_error;
1082  return(status);
1083 }
1084 
1085 static MagickBooleanType GetMSESimilarity(const Image *image,
1086  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
1087 {
1088  CacheView
1089  *image_view,
1090  *reconstruct_view;
1091 
1092  double
1093  area = 0.0;
1094 
1095  MagickBooleanType
1096  status = MagickTrue;
1097 
1098  size_t
1099  columns,
1100  rows;
1101 
1102  ssize_t
1103  channels = 0,
1104  k,
1105  y;
1106 
1107  /*
1108  Compute the mean sequared error similarity.
1109  */
1110  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
1111  image_view=AcquireVirtualCacheView(image,exception);
1112  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1113 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1114  #pragma omp parallel for schedule(static) shared(area,similarity,status) \
1115  magick_number_threads(image,image,rows,1)
1116 #endif
1117  for (y=0; y < (ssize_t) rows; y++)
1118  {
1119  const Quantum
1120  *magick_restrict p,
1121  *magick_restrict q;
1122 
1123  double
1124  channel_area = 0.0,
1125  channel_similarity[MaxPixelChannels+1] = { 0.0 };
1126 
1127  ssize_t
1128  x;
1129 
1130  if (status == MagickFalse)
1131  continue;
1132  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1133  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1134  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1135  {
1136  status=MagickFalse;
1137  continue;
1138  }
1139  for (x=0; x < (ssize_t) columns; x++)
1140  {
1141  double
1142  Da,
1143  Sa;
1144 
1145  ssize_t
1146  i;
1147 
1148  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1149  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1150  {
1151  p+=(ptrdiff_t) GetPixelChannels(image);
1152  q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1153  continue;
1154  }
1155  Sa=QuantumScale*(double) GetPixelAlpha(image,p);
1156  Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
1157  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1158  {
1159  double
1160  error;
1161 
1162  PixelChannel channel = GetPixelChannelChannel(image,i);
1163  PixelTrait traits = GetPixelChannelTraits(image,channel);
1164  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1165  channel);
1166  if (((traits & UpdatePixelTrait) == 0) ||
1167  ((reconstruct_traits & UpdatePixelTrait) == 0))
1168  continue;
1169  if (channel == AlphaPixelChannel)
1170  error=QuantumScale*((double) p[i]-(double) GetPixelChannel(
1171  reconstruct_image,channel,q));
1172  else
1173  error=QuantumScale*(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
1174  channel,q));
1175  channel_similarity[i]+=error*error;
1176  channel_similarity[CompositePixelChannel]+=error*error;
1177  }
1178  channel_area++;
1179  p+=(ptrdiff_t) GetPixelChannels(image);
1180  q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1181  }
1182 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1183  #pragma omp critical (MagickCore_GetMSESimilarity)
1184 #endif
1185  {
1186  ssize_t
1187  j;
1188 
1189  area+=channel_area;
1190  for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
1191  {
1192  PixelChannel channel = GetPixelChannelChannel(image,j);
1193  PixelTrait traits = GetPixelChannelTraits(image,channel);
1194  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1195  channel);
1196  if (((traits & UpdatePixelTrait) == 0) ||
1197  ((reconstruct_traits & UpdatePixelTrait) == 0))
1198  continue;
1199  similarity[j]+=channel_similarity[j];
1200  }
1201  similarity[CompositePixelChannel]+=
1202  channel_similarity[CompositePixelChannel];
1203  }
1204  }
1205  reconstruct_view=DestroyCacheView(reconstruct_view);
1206  image_view=DestroyCacheView(image_view);
1207  area=MagickSafeReciprocal(area);
1208  for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
1209  {
1210  PixelChannel channel = GetPixelChannelChannel(image,k);
1211  PixelTrait traits = GetPixelChannelTraits(image,channel);
1212  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1213  channel);
1214  if (((traits & UpdatePixelTrait) == 0) ||
1215  ((reconstruct_traits & UpdatePixelTrait) == 0))
1216  continue;
1217  similarity[k]*=area;
1218  channels++;
1219  }
1220  similarity[CompositePixelChannel]*=area;
1221  if (channels != 0)
1222  similarity[CompositePixelChannel]/=(double) channels;
1223  return(status);
1224 }
1225 
1226 static MagickBooleanType GetNCCSimilarity(const Image *image,
1227  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
1228 {
1229  CacheView
1230  *image_view,
1231  *reconstruct_view;
1232 
1233  ChannelStatistics
1234  *image_statistics,
1235  *reconstruct_statistics;
1236 
1237  double
1238  reconstruct_variance[MaxPixelChannels+1] = { 0.0 },
1239  variance[MaxPixelChannels+1] = { 0.0 };
1240 
1241  MagickBooleanType
1242  status = MagickTrue;
1243 
1244  MagickOffsetType
1245  progress = 0;
1246 
1247  size_t
1248  columns,
1249  rows;
1250 
1251  ssize_t
1252  k,
1253  y;
1254 
1255  /*
1256  Compute the normalized criss-correlation similarity.
1257  */
1258  image_statistics=GetImageStatistics(image,exception);
1259  reconstruct_statistics=GetImageStatistics(reconstruct_image,exception);
1260  if ((image_statistics == (ChannelStatistics *) NULL) ||
1261  (reconstruct_statistics == (ChannelStatistics *) NULL))
1262  {
1263  if (image_statistics != (ChannelStatistics *) NULL)
1264  image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1265  image_statistics);
1266  if (reconstruct_statistics != (ChannelStatistics *) NULL)
1267  reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1268  reconstruct_statistics);
1269  return(MagickFalse);
1270  }
1271  (void) memset(similarity,0,(MaxPixelChannels+1)*sizeof(*similarity));
1272  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
1273  image_view=AcquireVirtualCacheView(image,exception);
1274  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1275 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1276  #pragma omp parallel for schedule(static) shared(variance,reconstruct_variance,similarity,status) \
1277  magick_number_threads(image,image,rows,1)
1278 #endif
1279  for (y=0; y < (ssize_t) rows; y++)
1280  {
1281  const Quantum
1282  *magick_restrict p,
1283  *magick_restrict q;
1284 
1285  double
1286  channel_reconstruct_variance[MaxPixelChannels+1] = { 0.0 },
1287  channel_similarity[MaxPixelChannels+1] = { 0.0 },
1288  channel_variance[MaxPixelChannels+1] = { 0.0 };
1289 
1290  ssize_t
1291  x;
1292 
1293  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1294  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1295  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1296  {
1297  status=MagickFalse;
1298  continue;
1299  }
1300  for (x=0; x < (ssize_t) columns; x++)
1301  {
1302  double
1303  Da,
1304  Sa;
1305 
1306  ssize_t
1307  i;
1308 
1309  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1310  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1311  {
1312  p+=(ptrdiff_t) GetPixelChannels(image);
1313  q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1314  continue;
1315  }
1316  Sa=QuantumScale*(double) GetPixelAlpha(image,p);
1317  Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
1318  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1319  {
1320  double
1321  alpha,
1322  beta;
1323 
1324  PixelChannel channel = GetPixelChannelChannel(image,i);
1325  PixelTrait traits = GetPixelChannelTraits(image,channel);
1326  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1327  channel);
1328  if (((traits & UpdatePixelTrait) == 0) ||
1329  ((reconstruct_traits & UpdatePixelTrait) == 0))
1330  continue;
1331  if (channel == AlphaPixelChannel)
1332  {
1333  alpha=QuantumScale*((double) p[i]-image_statistics[channel].mean);
1334  beta=QuantumScale*((double) GetPixelChannel(reconstruct_image,
1335  channel,q)-reconstruct_statistics[channel].mean);
1336  }
1337  else
1338  {
1339  alpha=QuantumScale*(Sa*p[i]-image_statistics[channel].mean);
1340  beta=QuantumScale*(Da*GetPixelChannel(reconstruct_image,channel,
1341  q)-reconstruct_statistics[channel].mean);
1342  }
1343  channel_similarity[i]+=alpha*beta;
1344  channel_variance[i]+=alpha*alpha;
1345  channel_reconstruct_variance[i]+=beta*beta;
1346  }
1347  p+=(ptrdiff_t) GetPixelChannels(image);
1348  q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1349  }
1350 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1351  #pragma omp critical (MagickCore_GetNCCSimilarity)
1352 #endif
1353  {
1354  ssize_t
1355  j;
1356 
1357  for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
1358  {
1359  PixelChannel channel = GetPixelChannelChannel(image,j);
1360  PixelTrait traits = GetPixelChannelTraits(image,channel);
1361  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1362  channel);
1363  if (((traits & UpdatePixelTrait) == 0) ||
1364  ((reconstruct_traits & UpdatePixelTrait) == 0))
1365  continue;
1366  similarity[j]+=channel_similarity[j];
1367  similarity[CompositePixelChannel]+=channel_similarity[j];
1368  variance[j]+=channel_variance[j];
1369  variance[CompositePixelChannel]+=channel_variance[j];
1370  reconstruct_variance[j]+=channel_reconstruct_variance[j];
1371  reconstruct_variance[CompositePixelChannel]+=
1372  channel_reconstruct_variance[j];
1373  }
1374  }
1375  if (image->progress_monitor != (MagickProgressMonitor) NULL)
1376  {
1377  MagickBooleanType
1378  proceed;
1379 
1380 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1381  #pragma omp atomic
1382 #endif
1383  progress++;
1384  proceed=SetImageProgress(image,SimilarityImageTag,progress,rows);
1385  if (proceed == MagickFalse)
1386  {
1387  status=MagickFalse;
1388  continue;
1389  }
1390  }
1391  }
1392  reconstruct_view=DestroyCacheView(reconstruct_view);
1393  image_view=DestroyCacheView(image_view);
1394  /*
1395  Compute normalized cross correlation: divide by standard deviation.
1396  */
1397  for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
1398  {
1399  PixelChannel channel = GetPixelChannelChannel(image,k);
1400  PixelTrait traits = GetPixelChannelTraits(image,channel);
1401  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1402  channel);
1403  if (((traits & UpdatePixelTrait) == 0) ||
1404  ((reconstruct_traits & UpdatePixelTrait) == 0))
1405  continue;
1406  similarity[k]*=MagickSafeReciprocal(sqrt(variance[k])*
1407  sqrt(reconstruct_variance[k]));
1408  }
1409  similarity[CompositePixelChannel]*=MagickSafeReciprocal(sqrt(
1410  variance[CompositePixelChannel])*sqrt(
1411  reconstruct_variance[CompositePixelChannel]));
1412  /*
1413  Free resources.
1414  */
1415  reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1416  reconstruct_statistics);
1417  image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1418  image_statistics);
1419  return(status);
1420 }
1421 
1422 static MagickBooleanType GetPASimilarity(const Image *image,
1423  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
1424 {
1425  CacheView
1426  *image_view,
1427  *reconstruct_view;
1428 
1429  MagickBooleanType
1430  status = MagickTrue;
1431 
1432  size_t
1433  columns,
1434  rows;
1435 
1436  ssize_t
1437  y;
1438 
1439  /*
1440  Compute the peak absolute similarity.
1441  */
1442  (void) memset(similarity,0,(MaxPixelChannels+1)*sizeof(*similarity));
1443  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
1444  image_view=AcquireVirtualCacheView(image,exception);
1445  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1446 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1447  #pragma omp parallel for schedule(static) shared(similarity,status) \
1448  magick_number_threads(image,image,rows,1)
1449 #endif
1450  for (y=0; y < (ssize_t) rows; y++)
1451  {
1452  const Quantum
1453  *magick_restrict p,
1454  *magick_restrict q;
1455 
1456  double
1457  channel_similarity[MaxPixelChannels+1] = { 0.0 };
1458 
1459  ssize_t
1460  x;
1461 
1462  if (status == MagickFalse)
1463  continue;
1464  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1465  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1466  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1467  {
1468  status=MagickFalse;
1469  continue;
1470  }
1471  for (x=0; x < (ssize_t) columns; x++)
1472  {
1473  double
1474  Da,
1475  Sa;
1476 
1477  ssize_t
1478  i;
1479 
1480  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1481  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1482  {
1483  p+=(ptrdiff_t) GetPixelChannels(image);
1484  q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1485  continue;
1486  }
1487  Sa=QuantumScale*(double) GetPixelAlpha(image,p);
1488  Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
1489  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1490  {
1491  double
1492  distance;
1493 
1494  PixelChannel channel = GetPixelChannelChannel(image,i);
1495  PixelTrait traits = GetPixelChannelTraits(image,channel);
1496  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1497  channel);
1498  if (((traits & UpdatePixelTrait) == 0) ||
1499  ((reconstruct_traits & UpdatePixelTrait) == 0))
1500  continue;
1501  if (channel == AlphaPixelChannel)
1502  distance=QuantumScale*fabs((double) p[i]-(double)
1503  GetPixelChannel(reconstruct_image,channel,q));
1504  else
1505  distance=QuantumScale*fabs(Sa*p[i]-Da*GetPixelChannel(
1506  reconstruct_image,channel,q));
1507  if (distance > channel_similarity[i])
1508  channel_similarity[i]=distance;
1509  if (distance > channel_similarity[CompositePixelChannel])
1510  channel_similarity[CompositePixelChannel]=distance;
1511  }
1512  p+=(ptrdiff_t) GetPixelChannels(image);
1513  q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1514  }
1515 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1516  #pragma omp critical (MagickCore_GetPASimilarity)
1517 #endif
1518  {
1519  ssize_t
1520  j;
1521 
1522  for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
1523  {
1524  PixelChannel channel = GetPixelChannelChannel(image,j);
1525  PixelTrait traits = GetPixelChannelTraits(image,channel);
1526  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1527  channel);
1528  if (((traits & UpdatePixelTrait) == 0) ||
1529  ((reconstruct_traits & UpdatePixelTrait) == 0))
1530  continue;
1531  if (channel_similarity[j] > similarity[j])
1532  similarity[j]=channel_similarity[j];
1533  }
1534  if (channel_similarity[CompositePixelChannel] > similarity[CompositePixelChannel])
1535  similarity[CompositePixelChannel]=
1536  channel_similarity[CompositePixelChannel];
1537  }
1538  }
1539  reconstruct_view=DestroyCacheView(reconstruct_view);
1540  image_view=DestroyCacheView(image_view);
1541  return(status);
1542 }
1543 
1544 static MagickBooleanType GetPDCSimilarity(const Image *image,
1545  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
1546 {
1547  CacheView
1548  *image_view,
1549  *reconstruct_view;
1550 
1551  double
1552  fuzz;
1553 
1554  MagickBooleanType
1555  status = MagickTrue;
1556 
1557  size_t
1558  columns,
1559  rows;
1560 
1561  ssize_t
1562  y;
1563 
1564  /*
1565  Compute the pixel difference count similarity.
1566  */
1567  fuzz=GetFuzzyColorDistance(image,reconstruct_image);
1568  (void) memset(similarity,0,(MaxPixelChannels+1)*sizeof(*similarity));
1569  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
1570  image_view=AcquireVirtualCacheView(image,exception);
1571  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1572 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1573  #pragma omp parallel for schedule(static) shared(similarity,status) \
1574  magick_number_threads(image,image,rows,1)
1575 #endif
1576  for (y=0; y < (ssize_t) rows; y++)
1577  {
1578  const Quantum
1579  *magick_restrict p,
1580  *magick_restrict q;
1581 
1582  double
1583  channel_similarity[MaxPixelChannels+1] = { 0.0 };
1584 
1585  ssize_t
1586  x;
1587 
1588  if (status == MagickFalse)
1589  continue;
1590  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1591  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1592  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1593  {
1594  status=MagickFalse;
1595  continue;
1596  }
1597  for (x=0; x < (ssize_t) columns; x++)
1598  {
1599  double
1600  Da,
1601  Sa;
1602 
1603  size_t
1604  count = 0;
1605 
1606  ssize_t
1607  i;
1608 
1609  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1610  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1611  {
1612  p+=(ptrdiff_t) GetPixelChannels(image);
1613  q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1614  continue;
1615  }
1616  Sa=QuantumScale*(double) GetPixelAlpha(image,p);
1617  Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
1618  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1619  {
1620  double
1621  error;
1622 
1623  PixelChannel channel = GetPixelChannelChannel(image,i);
1624  PixelTrait traits = GetPixelChannelTraits(image,channel);
1625  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1626  channel);
1627  if (((traits & UpdatePixelTrait) == 0) ||
1628  ((reconstruct_traits & UpdatePixelTrait) == 0))
1629  continue;
1630  if (channel == AlphaPixelChannel)
1631  error=(double) p[i]-(double) GetPixelChannel(reconstruct_image,
1632  channel,q);
1633  else
1634  error=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
1635  if (MagickSafeSignificantError(error*error,fuzz) != MagickFalse)
1636  {
1637  channel_similarity[i]++;
1638  count++;
1639  }
1640  }
1641  if (count != 0)
1642  channel_similarity[CompositePixelChannel]++;
1643  p+=(ptrdiff_t) GetPixelChannels(image);
1644  q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1645  }
1646 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1647  #pragma omp critical (MagickCore_GetPDCSimilarity)
1648 #endif
1649  {
1650  ssize_t
1651  j;
1652 
1653  for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
1654  {
1655  PixelChannel channel = GetPixelChannelChannel(image,j);
1656  PixelTrait traits = GetPixelChannelTraits(image,channel);
1657  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1658  channel);
1659  if (((traits & UpdatePixelTrait) == 0) ||
1660  ((reconstruct_traits & UpdatePixelTrait) == 0))
1661  continue;
1662  similarity[j]+=channel_similarity[j];
1663  }
1664  similarity[CompositePixelChannel]+=
1665  channel_similarity[CompositePixelChannel];
1666  }
1667  }
1668  reconstruct_view=DestroyCacheView(reconstruct_view);
1669  image_view=DestroyCacheView(image_view);
1670  return(status);
1671 }
1672 
1673 static Image *GetPHASECorrelationSurface(const Image *image,
1674  ExceptionInfo *exception)
1675 {
1676  Image
1677  *surface;
1678 
1679  KernelInfo
1680  *kernel;
1681 
1682  /*
1683  Build a spatial-domain correlation surface with a 3x3 Laplacian
1684  high-pass convolution.
1685  */
1686  kernel=AcquireKernelInfo("3x3: 0,-1,0 -1,4,-1 0,-1,0",exception);
1687  if (kernel == (KernelInfo *) NULL)
1688  return((Image *) NULL);
1689  surface=MorphologyImage(image,ConvolveMorphology,1,kernel,exception);
1690  kernel=DestroyKernelInfo(kernel);
1691  return(surface);
1692 }
1693 
1694 static MagickBooleanType GetPHASESimilarity(const Image *image,
1695  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
1696 {
1697  CacheView
1698  *image_view,
1699  *reconstruct_view;
1700 
1701  double
1702  area = 0,
1703  correlation[MaxPixelChannels+1] = { 0.0 },
1704  image_sum[MaxPixelChannels+1] = { 0.0 },
1705  image_sum_squared[MaxPixelChannels+1] = { 0.0 },
1706  reconstruct_sum[MaxPixelChannels+1] = { 0.0 },
1707  reconstruct_sum_squared[MaxPixelChannels+1] = { 0.0 };
1708 
1709  Image
1710  *phase_image,
1711  *phase_reconstruct;
1712 
1713  MagickBooleanType
1714  status = MagickTrue;
1715 
1716  size_t
1717  columns = 0,
1718  rows = 0;
1719 
1720  ssize_t
1721  channels = 0,
1722  j,
1723  y;
1724 
1725  /*
1726  Compute the phase congruency similarity from two spatial high-pass
1727  correlation surfaces.
1728  */
1729  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
1730  phase_image=GetPHASECorrelationSurface(image,exception);
1731  phase_reconstruct=GetPHASECorrelationSurface(reconstruct_image,exception);
1732  if ((phase_image == (Image *) NULL) ||
1733  (phase_reconstruct == (Image *) NULL))
1734  {
1735  if (phase_image != (Image *) NULL)
1736  phase_image=DestroyImage(phase_image);
1737  if (phase_reconstruct != (Image *) NULL)
1738  phase_reconstruct=DestroyImage(phase_reconstruct);
1739  return(MagickFalse);
1740  }
1741  image_view=AcquireVirtualCacheView(phase_image,exception);
1742  reconstruct_view=AcquireVirtualCacheView(phase_reconstruct,exception);
1743 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1744  #pragma omp parallel for schedule(static) shared(status) \
1745  magick_number_threads(phase_image,phase_reconstruct,rows,1)
1746 #endif
1747  for (y=0; y < (ssize_t) rows; y++)
1748  {
1749  const Quantum
1750  *magick_restrict p,
1751  *magick_restrict q;
1752 
1753  double
1754  channel_area = 0,
1755  channel_correlation[MaxPixelChannels+1] = { 0.0 },
1756  channel_image_sum[MaxPixelChannels+1] = { 0.0 },
1757  channel_image_sum_squared[MaxPixelChannels+1] = { 0.0 },
1758  channel_reconstruct_sum[MaxPixelChannels+1] = { 0.0 },
1759  channel_reconstruct_sum_squared[MaxPixelChannels+1] = { 0.0 };
1760 
1761  ssize_t
1762  i,
1763  x;
1764 
1765  if (status == MagickFalse)
1766  continue;
1767  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1768  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1769  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1770  {
1771  status=MagickFalse;
1772  continue;
1773  }
1774  for (x=0; x < (ssize_t) columns; x++)
1775  {
1776  if ((GetPixelReadMask(phase_image,p) <= (QuantumRange/2)) ||
1777  (GetPixelReadMask(phase_reconstruct,q) <= (QuantumRange/2)))
1778  {
1779  p+=(ptrdiff_t) GetPixelChannels(phase_image);
1780  q+=(ptrdiff_t) GetPixelChannels(phase_reconstruct);
1781  continue;
1782  }
1783  for (i=0; i < (ssize_t) GetPixelChannels(phase_image); i++)
1784  {
1785  double
1786  alpha,
1787  beta;
1788 
1789  ssize_t
1790  offset;
1791 
1792  PixelChannel
1793  channel;
1794 
1795  PixelTrait
1796  reconstruct_traits,
1797  traits;
1798 
1799  channel=GetPixelChannelChannel(phase_image,i);
1800  traits=GetPixelChannelTraits(phase_image,channel);
1801  reconstruct_traits=GetPixelChannelTraits(phase_reconstruct,channel);
1802  if (((traits & UpdatePixelTrait) == 0) ||
1803  ((reconstruct_traits & UpdatePixelTrait) == 0))
1804  continue;
1805  offset=GetPixelChannelOffset(phase_reconstruct,channel);
1806  if (offset < 0)
1807  continue;
1808  alpha=QuantumScale*(double) p[i];
1809  beta=QuantumScale*(double) q[offset];
1810  channel_image_sum[i]+=alpha;
1811  channel_image_sum_squared[i]+=alpha*alpha;
1812  channel_reconstruct_sum[i]+=beta;
1813  channel_reconstruct_sum_squared[i]+=beta*beta;
1814  channel_correlation[i]+=alpha*beta;
1815  }
1816  channel_area++;
1817  p+=(ptrdiff_t) GetPixelChannels(phase_image);
1818  q+=(ptrdiff_t) GetPixelChannels(phase_reconstruct);
1819  }
1820 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1821  #pragma omp critical (MagickCore_GetPHASESimilarity)
1822 #endif
1823  {
1824  area+=channel_area;
1825  for (i=0; i <= (ssize_t) MaxPixelChannels; i++)
1826  {
1827  correlation[i]+=channel_correlation[i];
1828  image_sum[i]+=channel_image_sum[i];
1829  image_sum_squared[i]+=channel_image_sum_squared[i];
1830  reconstruct_sum[i]+=channel_reconstruct_sum[i];
1831  reconstruct_sum_squared[i]+=channel_reconstruct_sum_squared[i];
1832  }
1833  }
1834  }
1835  reconstruct_view=DestroyCacheView(reconstruct_view);
1836  image_view=DestroyCacheView(image_view);
1837  phase_image=DestroyImage(phase_image);
1838  phase_reconstruct=DestroyImage(phase_reconstruct);
1839  if ((status == MagickFalse) || (area < 1.0))
1840  return(status);
1841  /*
1842  Reduce the accumulated sums to a per-channel Pearson coefficient and
1843  average across channels for the composite value.
1844  */
1845  similarity[CompositePixelChannel]=0.0;
1846  for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
1847  {
1848  double
1849  denominator,
1850  numerator,
1851  pearson;
1852 
1853  PixelChannel
1854  channel;
1855 
1856  PixelTrait
1857  reconstruct_traits,
1858  traits;
1859 
1860  channel=GetPixelChannelChannel(image,j);
1861  traits=GetPixelChannelTraits(image,channel);
1862  reconstruct_traits=GetPixelChannelTraits(reconstruct_image,channel);
1863  if (((traits & UpdatePixelTrait) == 0) ||
1864  ((reconstruct_traits & UpdatePixelTrait) == 0))
1865  continue;
1866  numerator=area*correlation[j]-image_sum[j]*reconstruct_sum[j];
1867  denominator=sqrt(area*image_sum_squared[j]-image_sum[j]*image_sum[j])*
1868  sqrt(area*reconstruct_sum_squared[j]-reconstruct_sum[j]*
1869  reconstruct_sum[j]);
1870  pearson=denominator < MagickEpsilon ? 0.0 : numerator/denominator;
1871  if (pearson < -1.0)
1872  pearson=(-1.0);
1873  if (pearson > 1.0)
1874  pearson=1.0;
1875  similarity[j]=pearson;
1876  similarity[CompositePixelChannel]+=pearson;
1877  channels++;
1878  }
1879  if (channels != 0)
1880  similarity[CompositePixelChannel]/=(double) channels;
1881  return(status);
1882 }
1883 
1884 static MagickBooleanType GetPHASHSimilarity(const Image *image,
1885  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
1886 {
1887  ChannelPerceptualHash
1888  *channel_phash,
1889  *reconstruct_phash;
1890 
1891  const char
1892  *artifact;
1893 
1894  ssize_t
1895  channels = 0,
1896  i;
1897 
1898  /*
1899  Compute the perceptual hash similarity.
1900  */
1901  channel_phash=GetImagePerceptualHash(image,exception);
1902  if (channel_phash == (ChannelPerceptualHash *) NULL)
1903  return(MagickFalse);
1904  reconstruct_phash=GetImagePerceptualHash(reconstruct_image,exception);
1905  if (reconstruct_phash == (ChannelPerceptualHash *) NULL)
1906  {
1907  channel_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1908  channel_phash);
1909  return(MagickFalse);
1910  }
1911  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1912  {
1913  double
1914  difference = 0.0;
1915 
1916  ssize_t
1917  j;
1918 
1919  PixelChannel channel = GetPixelChannelChannel(image,i);
1920  PixelTrait traits = GetPixelChannelTraits(image,channel);
1921  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1922  channel);
1923  if (((traits & UpdatePixelTrait) == 0) ||
1924  ((reconstruct_traits & UpdatePixelTrait) == 0))
1925  continue;
1926  for (j=0; j < (ssize_t) channel_phash[0].number_colorspaces; j++)
1927  {
1928  double
1929  alpha,
1930  beta;
1931 
1932  ssize_t
1933  k;
1934 
1935  for (k=0; k < MaximumNumberOfPerceptualHashes; k++)
1936  {
1937  double
1938  error;
1939 
1940  alpha=channel_phash[i].phash[j][k];
1941  beta=reconstruct_phash[i].phash[j][k];
1942  error=beta-alpha;
1943  if (IsNaN(error) != 0)
1944  error=0.0;
1945  difference+=error*error;
1946  }
1947  }
1948  similarity[i]+=difference;
1949  similarity[CompositePixelChannel]+=difference;
1950  channels++;
1951  }
1952  if (channels != 0)
1953  similarity[CompositePixelChannel]/=(double) channels;
1954  artifact=GetImageArtifact(image,"phash:normalize");
1955  if (IsStringTrue(artifact) != MagickFalse)
1956  {
1957  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1958  {
1959  PixelChannel channel = GetPixelChannelChannel(image,i);
1960  PixelTrait traits = GetPixelChannelTraits(image,channel);
1961  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1962  channel);
1963  if (((traits & UpdatePixelTrait) == 0) ||
1964  ((reconstruct_traits & UpdatePixelTrait) == 0))
1965  continue;
1966  similarity[i]=sqrt(similarity[i]/channel_phash[0].number_colorspaces);
1967  }
1968  similarity[CompositePixelChannel]=sqrt(similarity[CompositePixelChannel]/
1969  channel_phash[0].number_colorspaces);
1970  }
1971  /*
1972  Free resources.
1973  */
1974  reconstruct_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1975  reconstruct_phash);
1976  channel_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(channel_phash);
1977  return(MagickTrue);
1978 }
1979 
1980 static MagickBooleanType GetPSNRSimilarity(const Image *image,
1981  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
1982 {
1983  MagickBooleanType
1984  status = MagickTrue;
1985 
1986  ssize_t
1987  i;
1988 
1989  /*
1990  Compute the peak signal-to-noise ratio similarity.
1991  */
1992  status=GetMSESimilarity(image,reconstruct_image,similarity,exception);
1993  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1994  {
1995  PixelChannel channel = GetPixelChannelChannel(image,i);
1996  PixelTrait traits = GetPixelChannelTraits(image,channel);
1997  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1998  channel);
1999  if (((traits & UpdatePixelTrait) == 0) ||
2000  ((reconstruct_traits & UpdatePixelTrait) == 0))
2001  continue;
2002  similarity[i]=10.0*MagickSafeLog10(MagickSafeReciprocal(
2003  similarity[i]))/MagickSafePSNRRecipicol(10.0);
2004  }
2005  similarity[CompositePixelChannel]=10.0*MagickSafeLog10(
2006  MagickSafeReciprocal(similarity[CompositePixelChannel]))/
2007  MagickSafePSNRRecipicol(10.0);
2008  return(status);
2009 }
2010 
2011 static MagickBooleanType GetRMSESimilarity(const Image *image,
2012  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
2013 {
2014 #define RMSESquareRoot(x) sqrt((x) < 0.0 ? 0.0 : (x))
2015 
2016  MagickBooleanType
2017  status = MagickTrue;
2018 
2019  ssize_t
2020  i;
2021 
2022  /*
2023  Compute the root mean-squared error similarity.
2024  */
2025  status=GetMSESimilarity(image,reconstruct_image,similarity,exception);
2026  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2027  {
2028  PixelChannel channel = GetPixelChannelChannel(image,i);
2029  PixelTrait traits = GetPixelChannelTraits(image,channel);
2030  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2031  channel);
2032  if (((traits & UpdatePixelTrait) == 0) ||
2033  ((reconstruct_traits & UpdatePixelTrait) == 0))
2034  continue;
2035  similarity[i]=RMSESquareRoot(similarity[i]);
2036  }
2037  similarity[CompositePixelChannel]=RMSESquareRoot(
2038  similarity[CompositePixelChannel]);
2039  return(status);
2040 }
2041 
2042 static MagickBooleanType GetSSIMSimularity(const Image *image,
2043  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
2044 {
2045 #define SSIMRadius 5.0
2046 #define SSIMSigma 1.5
2047 #define SSIMK1 0.01
2048 #define SSIMK2 0.03
2049 #define SSIML 1.0
2050 
2051  CacheView
2052  *image_view,
2053  *reconstruct_view;
2054 
2055  char
2056  geometry[MagickPathExtent];
2057 
2058  const char
2059  *artifact;
2060 
2061  double
2062  area = 0.0,
2063  c1,
2064  c2,
2065  radius,
2066  sigma;
2067 
2068  KernelInfo
2069  *kernel_info;
2070 
2071  MagickBooleanType
2072  status = MagickTrue;
2073 
2074  size_t
2075  columns,
2076  rows;
2077 
2078  ssize_t
2079  channels = 0,
2080  l,
2081  y;
2082 
2083  /*
2084  Compute the structual similarity index similarity.
2085  */
2086  radius=SSIMRadius;
2087  artifact=GetImageArtifact(image,"compare:ssim-radius");
2088  if (artifact != (const char *) NULL)
2089  radius=StringToDouble(artifact,(char **) NULL);
2090  sigma=SSIMSigma;
2091  artifact=GetImageArtifact(image,"compare:ssim-sigma");
2092  if (artifact != (const char *) NULL)
2093  sigma=StringToDouble(artifact,(char **) NULL);
2094  (void) FormatLocaleString(geometry,MagickPathExtent,"gaussian:%.17gx%.17g",
2095  radius,sigma);
2096  kernel_info=AcquireKernelInfo(geometry,exception);
2097  if (kernel_info == (KernelInfo *) NULL)
2098  ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
2099  image->filename);
2100  c1=pow(SSIMK1*SSIML,2.0);
2101  artifact=GetImageArtifact(image,"compare:ssim-k1");
2102  if (artifact != (const char *) NULL)
2103  c1=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
2104  c2=pow(SSIMK2*SSIML,2.0);
2105  artifact=GetImageArtifact(image,"compare:ssim-k2");
2106  if (artifact != (const char *) NULL)
2107  c2=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
2108  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
2109  image_view=AcquireVirtualCacheView(image,exception);
2110  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
2111 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2112  #pragma omp parallel for schedule(static) shared(area,similarity,status) \
2113  magick_number_threads(image,reconstruct_image,rows,1)
2114 #endif
2115  for (y=0; y < (ssize_t) rows; y++)
2116  {
2117  const Quantum
2118  *magick_restrict p,
2119  *magick_restrict q;
2120 
2121  double
2122  channel_area = 0.0,
2123  channel_similarity[MaxPixelChannels+1] = { 0.0 };
2124 
2125  ssize_t
2126  i,
2127  x;
2128 
2129  if (status == MagickFalse)
2130  continue;
2131  p=GetCacheViewVirtualPixels(image_view,-((ssize_t) kernel_info->width/2L),y-
2132  ((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
2133  kernel_info->height,exception);
2134  q=GetCacheViewVirtualPixels(reconstruct_view,-((ssize_t) kernel_info->width/
2135  2L),y-((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
2136  kernel_info->height,exception);
2137  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
2138  {
2139  status=MagickFalse;
2140  continue;
2141  }
2142  for (x=0; x < (ssize_t) columns; x++)
2143  {
2144  const Quantum
2145  *magick_restrict reconstruct,
2146  *magick_restrict test;
2147 
2148  double
2149  x_pixel_mu[MaxPixelChannels+1] = { 0.0 },
2150  x_pixel_sigma_squared[MaxPixelChannels+1] = { 0.0 },
2151  xy_sigma[MaxPixelChannels+1] = { 0.0 },
2152  y_pixel_mu[MaxPixelChannels+1] = { 0.0 },
2153  y_pixel_sigma_squared[MaxPixelChannels+1] = { 0.0 };
2154 
2155  MagickRealType
2156  *k;
2157 
2158  ssize_t
2159  v;
2160 
2161  if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
2162  (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
2163  {
2164  p+=(ptrdiff_t) GetPixelChannels(image);
2165  q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
2166  continue;
2167  }
2168  k=kernel_info->values;
2169  test=p;
2170  reconstruct=q;
2171  for (v=0; v < (ssize_t) kernel_info->height; v++)
2172  {
2173  ssize_t
2174  u;
2175 
2176  for (u=0; u < (ssize_t) kernel_info->width; u++)
2177  {
2178  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2179  {
2180  double
2181  x_pixel,
2182  y_pixel;
2183 
2184  PixelChannel channel = GetPixelChannelChannel(image,i);
2185  PixelTrait traits = GetPixelChannelTraits(image,channel);
2186  PixelTrait reconstruct_traits = GetPixelChannelTraits(
2187  reconstruct_image,channel);
2188  if (((traits & UpdatePixelTrait) == 0) ||
2189  ((reconstruct_traits & UpdatePixelTrait) == 0))
2190  continue;
2191  x_pixel=QuantumScale*(double) test[i];
2192  x_pixel_mu[i]+=(*k)*x_pixel;
2193  x_pixel_sigma_squared[i]+=(*k)*x_pixel*x_pixel;
2194  y_pixel=QuantumScale*(double)
2195  GetPixelChannel(reconstruct_image,channel,reconstruct);
2196  y_pixel_mu[i]+=(*k)*y_pixel;
2197  y_pixel_sigma_squared[i]+=(*k)*y_pixel*y_pixel;
2198  xy_sigma[i]+=(*k)*x_pixel*y_pixel;
2199  }
2200  k++;
2201  test+=(ptrdiff_t) GetPixelChannels(image);
2202  reconstruct+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
2203  }
2204  test+=(ptrdiff_t) GetPixelChannels(image)*columns;
2205  reconstruct+=(ptrdiff_t) GetPixelChannels(reconstruct_image)*columns;
2206  }
2207  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2208  {
2209  double
2210  ssim,
2211  x_pixel_mu_squared,
2212  x_pixel_sigmas_squared,
2213  xy_mu,
2214  xy_sigmas,
2215  y_pixel_mu_squared,
2216  y_pixel_sigmas_squared;
2217 
2218  PixelChannel channel = GetPixelChannelChannel(image,i);
2219  PixelTrait traits = GetPixelChannelTraits(image,channel);
2220  PixelTrait reconstruct_traits = GetPixelChannelTraits(
2221  reconstruct_image,channel);
2222  if (((traits & UpdatePixelTrait) == 0) ||
2223  ((reconstruct_traits & UpdatePixelTrait) == 0))
2224  continue;
2225  x_pixel_mu_squared=x_pixel_mu[i]*x_pixel_mu[i];
2226  y_pixel_mu_squared=y_pixel_mu[i]*y_pixel_mu[i];
2227  xy_mu=x_pixel_mu[i]*y_pixel_mu[i];
2228  xy_sigmas=xy_sigma[i]-xy_mu;
2229  x_pixel_sigmas_squared=x_pixel_sigma_squared[i]-x_pixel_mu_squared;
2230  y_pixel_sigmas_squared=y_pixel_sigma_squared[i]-y_pixel_mu_squared;
2231  ssim=((2.0*xy_mu+c1)*(2.0*xy_sigmas+c2))*
2232  MagickSafeReciprocal((x_pixel_mu_squared+y_pixel_mu_squared+c1)*
2233  (x_pixel_sigmas_squared+y_pixel_sigmas_squared+c2));
2234  channel_similarity[i]+=ssim;
2235  channel_similarity[CompositePixelChannel]+=ssim;
2236  }
2237  p+=(ptrdiff_t) GetPixelChannels(image);
2238  q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
2239  channel_area++;
2240  }
2241 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2242  #pragma omp critical (MagickCore_GetSSIMSimularity)
2243 #endif
2244  {
2245  ssize_t
2246  j;
2247 
2248  area+=channel_area;
2249  for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
2250  {
2251  PixelChannel channel = GetPixelChannelChannel(image,j);
2252  PixelTrait traits = GetPixelChannelTraits(image,channel);
2253  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2254  channel);
2255  if (((traits & UpdatePixelTrait) == 0) ||
2256  ((reconstruct_traits & UpdatePixelTrait) == 0))
2257  continue;
2258  similarity[j]+=channel_similarity[j];
2259  }
2260  similarity[CompositePixelChannel]+=
2261  channel_similarity[CompositePixelChannel];
2262  }
2263  }
2264  image_view=DestroyCacheView(image_view);
2265  reconstruct_view=DestroyCacheView(reconstruct_view);
2266  kernel_info=DestroyKernelInfo(kernel_info);
2267  area=MagickSafeReciprocal(area);
2268  for (l=0; l < (ssize_t) GetPixelChannels(image); l++)
2269  {
2270  PixelChannel channel = GetPixelChannelChannel(image,l);
2271  PixelTrait traits = GetPixelChannelTraits(image,channel);
2272  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2273  channel);
2274  if (((traits & UpdatePixelTrait) == 0) ||
2275  ((reconstruct_traits & UpdatePixelTrait) == 0))
2276  continue;
2277  similarity[l]*=area;
2278  channels++;
2279  }
2280  similarity[CompositePixelChannel]*=area;
2281  if (channels != 0)
2282  similarity[CompositePixelChannel]/=(double) channels;
2283  return(status);
2284 }
2285 
2286 static MagickBooleanType GetDSSIMSimilarity(const Image *image,
2287  const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
2288 {
2289  MagickBooleanType
2290  status = MagickTrue;
2291 
2292  ssize_t
2293  i;
2294 
2295  /*
2296  Compute the structual dissimilarity index similarity.
2297  */
2298  status=GetSSIMSimularity(image,reconstruct_image,similarity,exception);
2299  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2300  {
2301  PixelChannel channel = GetPixelChannelChannel(image,i);
2302  PixelTrait traits = GetPixelChannelTraits(image,channel);
2303  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2304  channel);
2305  if (((traits & UpdatePixelTrait) == 0) ||
2306  ((reconstruct_traits & UpdatePixelTrait) == 0))
2307  continue;
2308  similarity[i]=(1.0-similarity[i])/2.0;
2309  }
2310  similarity[CompositePixelChannel]=(1.0-similarity[CompositePixelChannel])/2.0;
2311  return(status);
2312 }
2313 
2314 MagickExport MagickBooleanType GetImageDistortion(Image *image,
2315  const Image *reconstruct_image,const MetricType metric,double *distortion,
2316  ExceptionInfo *exception)
2317 {
2318 #define CompareMetricNotSupportedException "metric not supported"
2319 
2320  double
2321  *channel_similarity;
2322 
2323  MagickBooleanType
2324  status = MagickTrue;
2325 
2326  size_t
2327  length;
2328 
2329  assert(image != (Image *) NULL);
2330  assert(image->signature == MagickCoreSignature);
2331  assert(reconstruct_image != (const Image *) NULL);
2332  assert(reconstruct_image->signature == MagickCoreSignature);
2333  assert(distortion != (double *) NULL);
2334  if (IsEventLogging() != MagickFalse)
2335  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2336  /*
2337  Get image distortion.
2338  */
2339  *distortion=0.0;
2340  length=MaxPixelChannels+1UL;
2341  channel_similarity=(double *) AcquireQuantumMemory(length,
2342  sizeof(*channel_similarity));
2343  if (channel_similarity == (double *) NULL)
2344  ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
2345  (void) memset(channel_similarity,0,length*sizeof(*channel_similarity));
2346  switch (metric)
2347  {
2348  case AbsoluteErrorMetric:
2349  {
2350  status=GetAESimilarity(image,reconstruct_image,channel_similarity,
2351  exception);
2352  break;
2353  }
2354  case DotProductCorrelationErrorMetric:
2355  {
2356  status=GetDPCSimilarity(image,reconstruct_image,channel_similarity,
2357  exception);
2358  break;
2359  }
2360  case FuzzErrorMetric:
2361  {
2362  status=GetFUZZSimilarity(image,reconstruct_image,channel_similarity,
2363  exception);
2364  break;
2365  }
2366  case MeanAbsoluteErrorMetric:
2367  {
2368  status=GetMAESimilarity(image,reconstruct_image,channel_similarity,
2369  exception);
2370  break;
2371  }
2372  case MeanErrorPerPixelErrorMetric:
2373  {
2374  status=GetMEPPSimilarity(image,reconstruct_image,channel_similarity,
2375  exception);
2376  break;
2377  }
2378  case MeanSquaredErrorMetric:
2379  {
2380  status=GetMSESimilarity(image,reconstruct_image,channel_similarity,
2381  exception);
2382  break;
2383  }
2384  case NormalizedCrossCorrelationErrorMetric:
2385  {
2386  status=GetNCCSimilarity(image,reconstruct_image,channel_similarity,
2387  exception);
2388  break;
2389  }
2390  case PeakAbsoluteErrorMetric:
2391  {
2392  status=GetPASimilarity(image,reconstruct_image,channel_similarity,
2393  exception);
2394  break;
2395  }
2396  case PeakSignalToNoiseRatioErrorMetric:
2397  {
2398  status=GetPSNRSimilarity(image,reconstruct_image,channel_similarity,
2399  exception);
2400  break;
2401  }
2402  case PerceptualHashErrorMetric:
2403  {
2404  status=GetPHASHSimilarity(image,reconstruct_image,channel_similarity,
2405  exception);
2406  break;
2407  }
2408  case PhaseCorrelationErrorMetric:
2409  {
2410  status=GetPHASESimilarity(image,reconstruct_image,channel_similarity,
2411  exception);
2412  break;
2413  }
2414  case PixelDifferenceCountErrorMetric:
2415  {
2416  status=GetPDCSimilarity(image,reconstruct_image,channel_similarity,
2417  exception);
2418  break;
2419  }
2420  case RootMeanSquaredErrorMetric:
2421  case UndefinedErrorMetric:
2422  default:
2423  {
2424  status=GetRMSESimilarity(image,reconstruct_image,channel_similarity,
2425  exception);
2426  break;
2427  }
2428  case StructuralDissimilarityErrorMetric:
2429  {
2430  status=GetDSSIMSimilarity(image,reconstruct_image,channel_similarity,
2431  exception);
2432  break;
2433  }
2434  case StructuralSimilarityErrorMetric:
2435  {
2436  status=GetSSIMSimularity(image,reconstruct_image,channel_similarity,
2437  exception);
2438  break;
2439  }
2440  }
2441  *distortion=channel_similarity[CompositePixelChannel];
2442  switch (metric)
2443  {
2444  case DotProductCorrelationErrorMetric:
2445  case NormalizedCrossCorrelationErrorMetric:
2446  case PhaseCorrelationErrorMetric:
2447  case StructuralSimilarityErrorMetric:
2448  {
2449  *distortion=(1.0-(*distortion))/2.0;
2450  break;
2451  }
2452  default: break;
2453  }
2454  channel_similarity=(double *) RelinquishMagickMemory(channel_similarity);
2455  if (fabs(*distortion) < MagickEpsilon)
2456  *distortion=0.0;
2457  (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
2458  *distortion);
2459  return(status);
2460 }
2461 ␌
2462 /*
2463 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2464 % %
2465 % %
2466 % %
2467 % G e t I m a g e D i s t o r t i o n s %
2468 % %
2469 % %
2470 % %
2471 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2472 %
2473 % GetImageDistortions() compares the pixel channels of an image to a
2474 % reconstructed image and returns the specified metric for each channel.
2475 %
2476 % The format of the GetImageDistortions method is:
2477 %
2478 % double *GetImageDistortions(const Image *image,
2479 % const Image *reconstruct_image,const MetricType metric,
2480 % ExceptionInfo *exception)
2481 %
2482 % A description of each parameter follows:
2483 %
2484 % o image: the image.
2485 %
2486 % o reconstruct_image: the reconstruction image.
2487 %
2488 % o metric: the metric.
2489 %
2490 % o exception: return any errors or warnings in this structure.
2491 %
2492 */
2493 MagickExport double *GetImageDistortions(Image *image,
2494  const Image *reconstruct_image,const MetricType metric,
2495  ExceptionInfo *exception)
2496 {
2497  double
2498  *distortion,
2499  *channel_similarity;
2500 
2501  MagickBooleanType
2502  status = MagickTrue;
2503 
2504  size_t
2505  length;
2506 
2507  ssize_t
2508  i;
2509 
2510  assert(image != (Image *) NULL);
2511  assert(image->signature == MagickCoreSignature);
2512  assert(reconstruct_image != (const Image *) NULL);
2513  assert(reconstruct_image->signature == MagickCoreSignature);
2514  if (IsEventLogging() != MagickFalse)
2515  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2516  /*
2517  Get image distortion.
2518  */
2519  length=MaxPixelChannels+1UL;
2520  channel_similarity=(double *) AcquireQuantumMemory(length,
2521  sizeof(*channel_similarity));
2522  if (channel_similarity == (double *) NULL)
2523  ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
2524  (void) memset(channel_similarity,0,length*sizeof(*channel_similarity));
2525  switch (metric)
2526  {
2527  case AbsoluteErrorMetric:
2528  {
2529  status=GetAESimilarity(image,reconstruct_image,channel_similarity,
2530  exception);
2531  break;
2532  }
2533  case DotProductCorrelationErrorMetric:
2534  {
2535  status=GetDPCSimilarity(image,reconstruct_image,channel_similarity,
2536  exception);
2537  break;
2538  }
2539  case FuzzErrorMetric:
2540  {
2541  status=GetFUZZSimilarity(image,reconstruct_image,channel_similarity,
2542  exception);
2543  break;
2544  }
2545  case MeanAbsoluteErrorMetric:
2546  {
2547  status=GetMAESimilarity(image,reconstruct_image,channel_similarity,
2548  exception);
2549  break;
2550  }
2551  case MeanErrorPerPixelErrorMetric:
2552  {
2553  status=GetMEPPSimilarity(image,reconstruct_image,channel_similarity,
2554  exception);
2555  break;
2556  }
2557  case MeanSquaredErrorMetric:
2558  {
2559  status=GetMSESimilarity(image,reconstruct_image,channel_similarity,
2560  exception);
2561  break;
2562  }
2563  case NormalizedCrossCorrelationErrorMetric:
2564  {
2565  status=GetNCCSimilarity(image,reconstruct_image,channel_similarity,
2566  exception);
2567  break;
2568  }
2569  case PeakAbsoluteErrorMetric:
2570  {
2571  status=GetPASimilarity(image,reconstruct_image,channel_similarity,
2572  exception);
2573  break;
2574  }
2575  case PeakSignalToNoiseRatioErrorMetric:
2576  {
2577  status=GetPSNRSimilarity(image,reconstruct_image,channel_similarity,
2578  exception);
2579  break;
2580  }
2581  case PerceptualHashErrorMetric:
2582  {
2583  status=GetPHASHSimilarity(image,reconstruct_image,channel_similarity,
2584  exception);
2585  break;
2586  }
2587  case PhaseCorrelationErrorMetric:
2588  {
2589  status=GetPHASESimilarity(image,reconstruct_image,channel_similarity,
2590  exception);
2591  break;
2592  }
2593  case PixelDifferenceCountErrorMetric:
2594  {
2595  status=GetPDCSimilarity(image,reconstruct_image,channel_similarity,
2596  exception);
2597  break;
2598  }
2599  case RootMeanSquaredErrorMetric:
2600  case UndefinedErrorMetric:
2601  default:
2602  {
2603  status=GetRMSESimilarity(image,reconstruct_image,channel_similarity,
2604  exception);
2605  break;
2606  }
2607  case StructuralDissimilarityErrorMetric:
2608  {
2609  status=GetDSSIMSimilarity(image,reconstruct_image,channel_similarity,
2610  exception);
2611  break;
2612  }
2613  case StructuralSimilarityErrorMetric:
2614  {
2615  status=GetSSIMSimularity(image,reconstruct_image,channel_similarity,
2616  exception);
2617  break;
2618  }
2619  }
2620  if (status == MagickFalse)
2621  {
2622  channel_similarity=(double *) RelinquishMagickMemory(channel_similarity);
2623  return((double *) NULL);
2624  }
2625  distortion=channel_similarity;
2626  switch (metric)
2627  {
2628  case DotProductCorrelationErrorMetric:
2629  case NormalizedCrossCorrelationErrorMetric:
2630  case PhaseCorrelationErrorMetric:
2631  case StructuralSimilarityErrorMetric:
2632  {
2633  for (i=0; i <= MaxPixelChannels; i++)
2634  distortion[i]=(1.0-distortion[i])/2.0;
2635  break;
2636  }
2637  default: break;
2638  }
2639  for (i=0; i <= MaxPixelChannels; i++)
2640  if (fabs(distortion[i]) < MagickEpsilon)
2641  distortion[i]=0.0;
2642  (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
2643  distortion[CompositePixelChannel]);
2644  return(distortion);
2645 }
2646 ␌
2647 /*
2648 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2649 % %
2650 % %
2651 % %
2652 % I s I m a g e s E q u a l %
2653 % %
2654 % %
2655 % %
2656 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2657 %
2658 % IsImagesEqual() compare the pixels of two images and returns immediately
2659 % if any pixel is not identical.
2660 %
2661 % The format of the IsImagesEqual method is:
2662 %
2663 % MagickBooleanType IsImagesEqual(const Image *image,
2664 % const Image *reconstruct_image,ExceptionInfo *exception)
2665 %
2666 % A description of each parameter follows.
2667 %
2668 % o image: the image.
2669 %
2670 % o reconstruct_image: the reconstruction image.
2671 %
2672 % o exception: return any errors or warnings in this structure.
2673 %
2674 */
2675 MagickExport MagickBooleanType IsImagesEqual(const Image *image,
2676  const Image *reconstruct_image,ExceptionInfo *exception)
2677 {
2678  CacheView
2679  *image_view,
2680  *reconstruct_view;
2681 
2682  size_t
2683  columns,
2684  rows;
2685 
2686  ssize_t
2687  y;
2688 
2689  assert(image != (Image *) NULL);
2690  assert(image->signature == MagickCoreSignature);
2691  assert(reconstruct_image != (const Image *) NULL);
2692  assert(reconstruct_image->signature == MagickCoreSignature);
2693  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
2694  image_view=AcquireVirtualCacheView(image,exception);
2695  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
2696  for (y=0; y < (ssize_t) rows; y++)
2697  {
2698  const Quantum
2699  *magick_restrict p,
2700  *magick_restrict q;
2701 
2702  ssize_t
2703  x;
2704 
2705  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
2706  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
2707  if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
2708  break;
2709  for (x=0; x < (ssize_t) columns; x++)
2710  {
2711  ssize_t
2712  i;
2713 
2714  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2715  {
2716  double
2717  distance;
2718 
2719  PixelChannel channel = GetPixelChannelChannel(image,i);
2720  PixelTrait traits = GetPixelChannelTraits(image,channel);
2721  PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2722  channel);
2723  if (((traits & UpdatePixelTrait) == 0) ||
2724  ((reconstruct_traits & UpdatePixelTrait) == 0))
2725  continue;
2726  distance=fabs((double) p[i]-(double) GetPixelChannel(reconstruct_image,
2727  channel,q));
2728  if (distance >= MagickEpsilon)
2729  break;
2730  }
2731  if (i < (ssize_t) GetPixelChannels(image))
2732  break;
2733  p+=(ptrdiff_t) GetPixelChannels(image);
2734  q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
2735  }
2736  if (x < (ssize_t) columns)
2737  break;
2738  }
2739  reconstruct_view=DestroyCacheView(reconstruct_view);
2740  image_view=DestroyCacheView(image_view);
2741  return(y < (ssize_t) rows ? MagickFalse : MagickTrue);
2742 }
2743 ␌
2744 /*
2745 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2746 % %
2747 % %
2748 % %
2749 % S e t I m a g e C o l o r M e t r i c %
2750 % %
2751 % %
2752 % %
2753 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2754 %
2755 % SetImageColorMetric() measures the difference between colors at each pixel
2756 % location of two images. A value other than 0 means the colors match
2757 % exactly. Otherwise an error measure is computed by summing over all
2758 % pixels in an image the distance squared in RGB space between each image
2759 % pixel and its corresponding pixel in the reconstruction image. The error
2760 % measure is assigned to these image members:
2761 %
2762 % o mean_error_per_pixel: The mean error for any single pixel in
2763 % the image.
2764 %
2765 % o normalized_mean_error: The normalized mean quantization error for
2766 % any single pixel in the image. This distance measure is normalized to
2767 % a range between 0 and 1. It is independent of the range of red, green,
2768 % and blue values in the image.
2769 %
2770 % o normalized_maximum_error: The normalized maximum quantization
2771 % error for any single pixel in the image. This distance measure is
2772 % normalized to a range between 0 and 1. It is independent of the range
2773 % of red, green, and blue values in your image.
2774 %
2775 % A small normalized mean square error, accessed as
2776 % image->normalized_mean_error, suggests the images are very similar in
2777 % spatial layout and color.
2778 %
2779 % The format of the SetImageColorMetric method is:
2780 %
2781 % MagickBooleanType SetImageColorMetric(Image *image,
2782 % const Image *reconstruct_image,ExceptionInfo *exception)
2783 %
2784 % A description of each parameter follows.
2785 %
2786 % o image: the image.
2787 %
2788 % o reconstruct_image: the reconstruction image.
2789 %
2790 % o exception: return any errors or warnings in this structure.
2791 %
2792 */
2793 MagickExport MagickBooleanType SetImageColorMetric(Image *image,
2794  const Image *reconstruct_image,ExceptionInfo *exception)
2795 {
2796  double
2797  channel_similarity[MaxPixelChannels+1] = { 0.0 };
2798 
2799  MagickBooleanType
2800  status;
2801 
2802  status=GetMEPPSimilarity(image,reconstruct_image,channel_similarity,
2803  exception);
2804  if (status == MagickFalse)
2805  return(MagickFalse);
2806  status=fabs(image->error.mean_error_per_pixel) < MagickEpsilon ?
2807  MagickTrue : MagickFalse;
2808  return(status);
2809 }
2810 ␌
2811 /*
2812 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2813 % %
2814 % %
2815 % %
2816 % S i m i l a r i t y I m a g e %
2817 % %
2818 % %
2819 % %
2820 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2821 %
2822 % SimilarityImage() compares the reconstruction of the image and returns the
2823 % best match offset. In addition, it returns a similarity image such that an
2824 % exact match location is completely white and if none of the pixels match,
2825 % black, otherwise some gray level in-between.
2826 %
2827 % Contributed by Fred Weinhaus.
2828 %
2829 % The format of the SimilarityImageImage method is:
2830 %
2831 % Image *SimilarityImage(const Image *image,const Image *reconstruct,
2832 % const MetricType metric,const double similarity_threshold,
2833 % RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
2834 %
2835 % A description of each parameter follows:
2836 %
2837 % o image: the image.
2838 %
2839 % o reconstruct: find an area of the image that closely resembles this image.
2840 %
2841 % o metric: the metric.
2842 %
2843 % o similarity_threshold: minimum similarity for (sub)image match.
2844 %
2845 % o offset: the best match offset of the reconstruction image within the
2846 % image.
2847 %
2848 % o similarity: the computed similarity between the images.
2849 %
2850 % o exception: return any errors or warnings in this structure.
2851 %
2852 */
2853 
2854 #if defined(MAGICKCORE_HDRI_SUPPORT) && defined(MAGICKCORE_FFTW_DELEGATE)
2855 static Image *SIMCrossCorrelationImage(const Image *alpha_image,
2856  const Image *beta_image,ExceptionInfo *exception)
2857 {
2858  Image
2859  *alpha_fft = (Image *) NULL,
2860  *beta_fft = (Image *) NULL,
2861  *complex_conjugate = (Image *) NULL,
2862  *complex_multiplication = (Image *) NULL,
2863  *cross_correlation = (Image *) NULL,
2864  *temp_image = (Image *) NULL;
2865 
2866  /*
2867  Take the FFT of beta (reconstruction) image.
2868  */
2869  temp_image=CloneImage(beta_image,0,0,MagickTrue,exception);
2870  if (temp_image == (Image *) NULL)
2871  return((Image *) NULL);
2872  (void) SetImageArtifact(temp_image,"fourier:normalize","inverse");
2873  beta_fft=ForwardFourierTransformImage(temp_image,MagickFalse,exception);
2874  temp_image=DestroyImageList(temp_image);
2875  if (beta_fft == (Image *) NULL)
2876  return((Image *) NULL);
2877  /*
2878  Take the complex conjugate of beta_fft.
2879  */
2880  complex_conjugate=ComplexImages(beta_fft,ConjugateComplexOperator,exception);
2881  beta_fft=DestroyImageList(beta_fft);
2882  if (complex_conjugate == (Image *) NULL)
2883  return((Image *) NULL);
2884  /*
2885  Take the FFT of the alpha (test) image.
2886  */
2887  temp_image=CloneImage(alpha_image,0,0,MagickTrue,exception);
2888  if (temp_image == (Image *) NULL)
2889  {
2890  complex_conjugate=DestroyImageList(complex_conjugate);
2891  return((Image *) NULL);
2892  }
2893  (void) SetImageArtifact(temp_image,"fourier:normalize","inverse");
2894  alpha_fft=ForwardFourierTransformImage(temp_image,MagickFalse,exception);
2895  temp_image=DestroyImageList(temp_image);
2896  if (alpha_fft == (Image *) NULL)
2897  {
2898  complex_conjugate=DestroyImageList(complex_conjugate);
2899  return((Image *) NULL);
2900  }
2901  /*
2902  Do complex multiplication.
2903  */
2904  DisableCompositeClampUnlessSpecified(complex_conjugate);
2905  DisableCompositeClampUnlessSpecified(complex_conjugate->next);
2906  AppendImageToList(&complex_conjugate,alpha_fft);
2907  complex_multiplication=ComplexImages(complex_conjugate,
2908  MultiplyComplexOperator,exception);
2909  complex_conjugate=DestroyImageList(complex_conjugate);
2910  if (complex_multiplication == (Image *) NULL)
2911  return((Image *) NULL);
2912  /*
2913  Do the IFT and return the cross-correlation result.
2914  */
2915  cross_correlation=InverseFourierTransformImage(complex_multiplication,
2916  complex_multiplication->next,MagickFalse,exception);
2917  complex_multiplication=DestroyImageList(complex_multiplication);
2918  return(cross_correlation);
2919 }
2920 
2921 static Image *SIMDerivativeImage(const Image *image,const char *kernel,
2922  ExceptionInfo *exception)
2923 {
2924  Image
2925  *derivative_image;
2926 
2927  KernelInfo
2928  *kernel_info;
2929 
2930  kernel_info=AcquireKernelInfo(kernel,exception);
2931  if (kernel_info == (KernelInfo *) NULL)
2932  return((Image *) NULL);
2933  derivative_image=MorphologyImage(image,ConvolveMorphology,1,kernel_info,
2934  exception);
2935  kernel_info=DestroyKernelInfo(kernel_info);
2936  return(derivative_image);
2937 }
2938 
2939 static Image *SIMDivideImage(const Image *numerator_image,
2940  const Image *denominator_image,ExceptionInfo *exception)
2941 {
2942  CacheView
2943  *denominator_view,
2944  *numerator_view;
2945 
2946  Image
2947  *divide_image;
2948 
2949  MagickBooleanType
2950  status = MagickTrue;
2951 
2952  ssize_t
2953  y;
2954 
2955  /*
2956  Divide one image into another.
2957  */
2958  divide_image=CloneImage(numerator_image,0,0,MagickTrue,exception);
2959  if (divide_image == (Image *) NULL)
2960  return(divide_image);
2961  numerator_view=AcquireAuthenticCacheView(divide_image,exception);
2962  denominator_view=AcquireVirtualCacheView(denominator_image,exception);
2963 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2964  #pragma omp parallel for schedule(static) shared(status) \
2965  magick_number_threads(denominator_image,divide_image,divide_image->rows,1)
2966 #endif
2967  for (y=0; y < (ssize_t) divide_image->rows; y++)
2968  {
2969  const Quantum
2970  *magick_restrict p;
2971 
2972  Quantum
2973  *magick_restrict q;
2974 
2975  ssize_t
2976  x;
2977 
2978  if (status == MagickFalse)
2979  continue;
2980  p=GetCacheViewVirtualPixels(denominator_view,0,y,
2981  denominator_image->columns,1,exception);
2982  q=GetCacheViewAuthenticPixels(numerator_view,0,y,divide_image->columns,1,
2983  exception);
2984  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2985  {
2986  status=MagickFalse;
2987  continue;
2988  }
2989  for (x=0; x < (ssize_t) divide_image->columns; x++)
2990  {
2991  ssize_t
2992  i;
2993 
2994  for (i=0; i < (ssize_t) GetPixelChannels(divide_image); i++)
2995  {
2996  PixelChannel channel = GetPixelChannelChannel(divide_image,i);
2997  PixelTrait traits = GetPixelChannelTraits(divide_image,channel);
2998  PixelTrait denominator_traits = GetPixelChannelTraits(denominator_image,
2999  channel);
3000  if (((traits & UpdatePixelTrait) == 0) ||
3001  ((denominator_traits & UpdatePixelTrait) == 0))
3002  continue;
3003  q[i]=(Quantum) ((double) q[i]*MagickSafeReciprocal(QuantumScale*
3004  (double) GetPixelChannel(denominator_image,channel,p)));
3005  }
3006  p+=(ptrdiff_t) GetPixelChannels(denominator_image);
3007  q+=(ptrdiff_t) GetPixelChannels(divide_image);
3008  }
3009  if (SyncCacheViewAuthenticPixels(numerator_view,exception) == MagickFalse)
3010  status=MagickFalse;
3011  }
3012  denominator_view=DestroyCacheView(denominator_view);
3013  numerator_view=DestroyCacheView(numerator_view);
3014  if (status == MagickFalse)
3015  divide_image=DestroyImage(divide_image);
3016  return(divide_image);
3017 }
3018 
3019 static Image *SIMDivideByMagnitude(Image *image,Image *magnitude_image,
3020  const Image *source_image,ExceptionInfo *exception)
3021 {
3022  Image
3023  *divide_image,
3024  *result_image;
3025 
3026  RectangleInfo
3027  geometry;
3028 
3029  divide_image=SIMDivideImage(image,magnitude_image,exception);
3030  if (divide_image == (Image *) NULL)
3031  return((Image *) NULL);
3032  GetPixelInfoRGBA((Quantum) 0,(Quantum) 0,(Quantum) 0,(Quantum) 0,
3033  &divide_image->background_color);
3034  SetGeometry(source_image,&geometry);
3035  geometry.width=MagickMax(source_image->columns,divide_image->columns);
3036  geometry.height=MagickMax(source_image->rows,divide_image->rows);
3037  result_image=ExtentImage(divide_image,&geometry,exception);
3038  divide_image=DestroyImage(divide_image);
3039  return(result_image);
3040 }
3041 
3042 static MagickBooleanType SIMFilterImageNaNs(Image *image,
3043  ExceptionInfo *exception)
3044 {
3045  CacheView
3046  *image_view;
3047 
3048  MagickBooleanType
3049  status = MagickTrue;
3050 
3051  ssize_t
3052  y;
3053 
3054  /*
3055  Square each pixel in the image.
3056  */
3057  image_view=AcquireAuthenticCacheView(image,exception);
3058 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3059  #pragma omp parallel for schedule(static) shared(status) \
3060  magick_number_threads(image,image,image->rows,1)
3061 #endif
3062  for (y=0; y < (ssize_t) image->rows; y++)
3063  {
3064  Quantum
3065  *magick_restrict q;
3066 
3067  ssize_t
3068  x;
3069 
3070  if (status == MagickFalse)
3071  continue;
3072  q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
3073  if (q == (Quantum *) NULL)
3074  {
3075  status=MagickFalse;
3076  continue;
3077  }
3078  for (x=0; x < (ssize_t) image->columns; x++)
3079  {
3080  ssize_t
3081  i;
3082 
3083  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3084  {
3085  PixelChannel channel = GetPixelChannelChannel(image,i);
3086  PixelTrait traits = GetPixelChannelTraits(image,channel);
3087  if ((traits & UpdatePixelTrait) == 0)
3088  continue;
3089  if (IsNaN((double) q[i]) != 0)
3090  q[i]=(Quantum) 0;
3091  }
3092  q+=(ptrdiff_t) GetPixelChannels(image);
3093  }
3094  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
3095  status=MagickFalse;
3096  }
3097  image_view=DestroyCacheView(image_view);
3098  return(status);
3099 }
3100 
3101 static Image *SIMSquareImage(const Image *image,ExceptionInfo *exception)
3102 {
3103  CacheView
3104  *image_view;
3105 
3106  Image
3107  *square_image;
3108 
3109  MagickBooleanType
3110  status = MagickTrue;
3111 
3112  ssize_t
3113  y;
3114 
3115  /*
3116  Square each pixel in the image.
3117  */
3118  square_image=CloneImage(image,0,0,MagickTrue,exception);
3119  if (square_image == (Image *) NULL)
3120  return(square_image);
3121  image_view=AcquireAuthenticCacheView(square_image,exception);
3122 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3123  #pragma omp parallel for schedule(static) shared(status) \
3124  magick_number_threads(square_image,square_image,square_image->rows,1)
3125 #endif
3126  for (y=0; y < (ssize_t) square_image->rows; y++)
3127  {
3128  Quantum
3129  *magick_restrict q;
3130 
3131  ssize_t
3132  x;
3133 
3134  if (status == MagickFalse)
3135  continue;
3136  q=GetCacheViewAuthenticPixels(image_view,0,y,square_image->columns,1,
3137  exception);
3138  if (q == (Quantum *) NULL)
3139  {
3140  status=MagickFalse;
3141  continue;
3142  }
3143  for (x=0; x < (ssize_t) square_image->columns; x++)
3144  {
3145  ssize_t
3146  i;
3147 
3148  for (i=0; i < (ssize_t) GetPixelChannels(square_image); i++)
3149  {
3150  PixelChannel channel = GetPixelChannelChannel(square_image,i);
3151  PixelTrait traits = GetPixelChannelTraits(square_image,channel);
3152  if ((traits & UpdatePixelTrait) == 0)
3153  continue;
3154  q[i]=(Quantum) (QuantumScale*q[i]*q[i]);
3155  }
3156  q+=(ptrdiff_t) GetPixelChannels(square_image);
3157  }
3158  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
3159  status=MagickFalse;
3160  }
3161  image_view=DestroyCacheView(image_view);
3162  if (status == MagickFalse)
3163  square_image=DestroyImage(square_image);
3164  return(square_image);
3165 }
3166 
3167 static Image *SIMMagnitudeImage(Image *alpha_image,Image *beta_image,
3168  ExceptionInfo *exception)
3169 {
3170  Image
3171  *magnitude_image,
3172  *xsq_image,
3173  *ysq_image;
3174 
3175  MagickBooleanType
3176  status = MagickTrue;
3177 
3178  (void) SetImageArtifact(alpha_image,"compose:clamp","False");
3179  xsq_image=SIMSquareImage(alpha_image,exception);
3180  if (xsq_image == (Image *) NULL)
3181  return((Image *) NULL);
3182  (void) SetImageArtifact(beta_image,"compose:clamp","False");
3183  ysq_image=SIMSquareImage(beta_image,exception);
3184  if (ysq_image == (Image *) NULL)
3185  {
3186  xsq_image=DestroyImage(xsq_image);
3187  return((Image *) NULL);
3188  }
3189  status=CompositeImage(xsq_image,ysq_image,PlusCompositeOp,MagickTrue,0,0,
3190  exception);
3191  magnitude_image=xsq_image;
3192  ysq_image=DestroyImage(ysq_image);
3193  if (status == MagickFalse)
3194  {
3195  magnitude_image=DestroyImage(magnitude_image);
3196  return((Image *) NULL);
3197  }
3198  status=EvaluateImage(magnitude_image,PowEvaluateOperator,0.5,exception);
3199  if (status == MagickFalse)
3200  {
3201  magnitude_image=DestroyImage(magnitude_image);
3202  return (Image *) NULL;
3203  }
3204  return(magnitude_image);
3205 }
3206 
3207 static MagickBooleanType SIMMaximaImage(const Image *image,double *maxima,
3208  RectangleInfo *offset,ExceptionInfo *exception)
3209 {
3210  typedef struct
3211  {
3212  double
3213  maxima;
3214 
3215  ssize_t
3216  x,
3217  y;
3218  } MaximaInfo;
3219 
3220  CacheView
3221  *image_view;
3222 
3223  const Quantum
3224  *magick_restrict q;
3225 
3226  MagickBooleanType
3227  status = MagickTrue;
3228 
3229  MaximaInfo
3230  maxima_info = { -MagickMaximumValue, 0, 0 };
3231 
3232  ssize_t
3233  y;
3234 
3235  /*
3236  Identify the maxima value in the image and its location.
3237  */
3238  image_view=AcquireVirtualCacheView(image,exception);
3239  q=GetCacheViewVirtualPixels(image_view,maxima_info.x,maxima_info.y,1,1,
3240  exception);
3241  if (q != (const Quantum *) NULL)
3242  maxima_info.maxima=IsNaN((double) q[0]) != 0 ? 0.0 : (double) q[0];
3243 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3244  #pragma omp parallel for schedule(static) shared(maxima_info,status) \
3245  magick_number_threads(image,image,image->rows,1)
3246 #endif
3247  for (y=0; y < (ssize_t) image->rows; y++)
3248  {
3249  const Quantum
3250  *magick_restrict p;
3251 
3252  MaximaInfo
3253  channel_maxima = { -MagickMaximumValue, 0, 0 };
3254 
3255  ssize_t
3256  x;
3257 
3258  if (status == MagickFalse)
3259  continue;
3260  p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
3261  if (p == (const Quantum *) NULL)
3262  {
3263  status=MagickFalse;
3264  continue;
3265  }
3266  channel_maxima=maxima_info;
3267  for (x=0; x < (ssize_t) image->columns; x++)
3268  {
3269  ssize_t
3270  i;
3271 
3272  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3273  {
3274  double
3275  pixel;
3276 
3277  PixelChannel channel = GetPixelChannelChannel(image,i);
3278  PixelTrait traits = GetPixelChannelTraits(image,channel);
3279  if ((traits & UpdatePixelTrait) == 0)
3280  continue;
3281  pixel=(double) p[i];
3282  if (IsNaN(pixel) != 0)
3283  pixel=0.0;
3284  if (pixel > channel_maxima.maxima)
3285  {
3286  channel_maxima.maxima=(double) p[i];
3287  channel_maxima.x=x;
3288  channel_maxima.y=y;
3289  }
3290  }
3291  p+=(ptrdiff_t) GetPixelChannels(image);
3292  }
3293 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3294  #pragma omp critical (MagickCore_SIMMaximaImage)
3295 #endif
3296  if (channel_maxima.maxima > maxima_info.maxima)
3297  maxima_info=channel_maxima;
3298  }
3299  image_view=DestroyCacheView(image_view);
3300  *maxima=maxima_info.maxima;
3301  offset->x=maxima_info.x;
3302  offset->y=maxima_info.y;
3303  return(status);
3304 }
3305 
3306 static MagickBooleanType SIMMinimaImage(const Image *image,double *minima,
3307  RectangleInfo *offset,ExceptionInfo *exception)
3308 {
3309  typedef struct
3310  {
3311  double
3312  minima;
3313 
3314  ssize_t
3315  x,
3316  y;
3317  } MinimaInfo;
3318 
3319  CacheView
3320  *image_view;
3321 
3322  const Quantum
3323  *magick_restrict q;
3324 
3325  MagickBooleanType
3326  status = MagickTrue;
3327 
3328  MinimaInfo
3329  minima_info = { MagickMaximumValue, 0, 0 };
3330 
3331  ssize_t
3332  y;
3333 
3334  /*
3335  Identify the minima value in the image and its location.
3336  */
3337  image_view=AcquireVirtualCacheView(image,exception);
3338  q=GetCacheViewVirtualPixels(image_view,minima_info.x,minima_info.y,1,1,
3339  exception);
3340  if (q != (const Quantum *) NULL)
3341  minima_info.minima=IsNaN((double) q[0]) != 0 ? 0.0 : (double) q[0];
3342 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3343  #pragma omp parallel for schedule(static) shared(minima_info,status) \
3344  magick_number_threads(image,image,image->rows,1)
3345 #endif
3346  for (y=0; y < (ssize_t) image->rows; y++)
3347  {
3348  const Quantum
3349  *magick_restrict p;
3350 
3351  MinimaInfo
3352  channel_minima = { MagickMaximumValue, 0, 0 };
3353 
3354  ssize_t
3355  x;
3356 
3357  if (status == MagickFalse)
3358  continue;
3359  p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
3360  if (p == (const Quantum *) NULL)
3361  {
3362  status=MagickFalse;
3363  continue;
3364  }
3365  channel_minima=minima_info;
3366  for (x=0; x < (ssize_t) image->columns; x++)
3367  {
3368  ssize_t
3369  i;
3370 
3371  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3372  {
3373  double
3374  pixel;
3375 
3376  PixelChannel channel = GetPixelChannelChannel(image,i);
3377  PixelTrait traits = GetPixelChannelTraits(image,channel);
3378  if ((traits & UpdatePixelTrait) == 0)
3379  continue;
3380  pixel=(double) p[i];
3381  if (IsNaN(pixel) != 0)
3382  pixel=0.0;
3383  if (pixel < channel_minima.minima)
3384  {
3385  channel_minima.minima=pixel;
3386  channel_minima.x=x;
3387  channel_minima.y=y;
3388  }
3389  }
3390  p+=(ptrdiff_t) GetPixelChannels(image);
3391  }
3392 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3393  #pragma omp critical (MagickCore_SIMMinimaImage)
3394 #endif
3395  if (channel_minima.minima < minima_info.minima)
3396  minima_info=channel_minima;
3397  }
3398  image_view=DestroyCacheView(image_view);
3399  *minima=minima_info.minima;
3400  offset->x=minima_info.x;
3401  offset->y=minima_info.y;
3402  return(status);
3403 }
3404 
3405 static MagickBooleanType SIMMultiplyImage(Image *image,const double factor,
3406  const ChannelStatistics *channel_statistics,ExceptionInfo *exception)
3407 {
3408  CacheView
3409  *image_view;
3410 
3411  MagickBooleanType
3412  status = MagickTrue;
3413 
3414  ssize_t
3415  y;
3416 
3417  /*
3418  Multiply each pixel by a factor.
3419  */
3420  image_view=AcquireAuthenticCacheView(image,exception);
3421 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3422  #pragma omp parallel for schedule(static) shared(status) \
3423  magick_number_threads(image,image,image->rows,1)
3424 #endif
3425  for (y=0; y < (ssize_t) image->rows; y++)
3426  {
3427  Quantum
3428  *magick_restrict q;
3429 
3430  ssize_t
3431  x;
3432 
3433  if (status == MagickFalse)
3434  continue;
3435  q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
3436  if (q == (Quantum *) NULL)
3437  {
3438  status=MagickFalse;
3439  continue;
3440  }
3441  for (x=0; x < (ssize_t) image->columns; x++)
3442  {
3443  ssize_t
3444  i;
3445 
3446  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3447  {
3448  PixelChannel channel = GetPixelChannelChannel(image,i);
3449  PixelTrait traits = GetPixelChannelTraits(image,channel);
3450  if ((traits & UpdatePixelTrait) == 0)
3451  continue;
3452  if (channel_statistics != (const ChannelStatistics *) NULL)
3453  q[i]=(Quantum) (factor*q[i]*QuantumScale*
3454  channel_statistics[channel].standard_deviation);
3455  else
3456  q[i]=(Quantum) (factor*q[i]);
3457  }
3458  q+=(ptrdiff_t) GetPixelChannels(image);
3459  }
3460  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
3461  status=MagickFalse;
3462  }
3463  image_view=DestroyCacheView(image_view);
3464  return(status);
3465 }
3466 
3467 static Image *SIMPhaseCorrelationImage(const Image *alpha_image,
3468  const Image *beta_image,const Image *magnitude_image,ExceptionInfo *exception)
3469 {
3470  Image
3471  *alpha_fft = (Image *) NULL,
3472  *beta_fft = (Image *) NULL,
3473  *complex_multiplication = (Image *) NULL,
3474  *cross_correlation = (Image *) NULL;
3475 
3476  /*
3477  Take the FFT of the beta (reconstruction) image.
3478  */
3479  beta_fft=CloneImage(beta_image,0,0,MagickTrue,exception);
3480  if (beta_fft == NULL)
3481  return((Image *) NULL);
3482  (void) SetImageArtifact(beta_fft,"fourier:normalize","inverse");
3483  beta_fft=ForwardFourierTransformImage(beta_fft,MagickFalse,exception);
3484  if (beta_fft == NULL)
3485  return((Image *) NULL);
3486  /*
3487  Take the FFT of the alpha (test) image.
3488  */
3489  alpha_fft=CloneImage(alpha_image,0,0,MagickTrue,exception);
3490  if (alpha_fft == (Image *) NULL)
3491  {
3492  beta_fft=DestroyImageList(beta_fft);
3493  return((Image *) NULL);
3494  }
3495  (void) SetImageArtifact(alpha_fft,"fourier:normalize","inverse");
3496  alpha_fft=ForwardFourierTransformImage(alpha_fft,MagickFalse,exception);
3497  if (alpha_fft == (Image *) NULL)
3498  {
3499  beta_fft=DestroyImageList(beta_fft);
3500  return((Image *) NULL);
3501  }
3502  /*
3503  Take the complex conjugate of the beta FFT.
3504  */
3505  beta_fft=ComplexImages(beta_fft,ConjugateComplexOperator,exception);
3506  if (beta_fft == (Image *) NULL)
3507  {
3508  alpha_fft=DestroyImageList(alpha_fft);
3509  return((Image *) NULL);
3510  }
3511  /*
3512  Do complex multiplication.
3513  */
3514  AppendImageToList(&beta_fft,alpha_fft);
3515  DisableCompositeClampUnlessSpecified(beta_fft);
3516  DisableCompositeClampUnlessSpecified(beta_fft->next);
3517  complex_multiplication=ComplexImages(beta_fft,MultiplyComplexOperator,
3518  exception);
3519  beta_fft=DestroyImageList(beta_fft);
3520  if (complex_multiplication == (Image *) NULL)
3521  return((Image *) NULL);
3522  /*
3523  Divide the results.
3524  */
3525  CompositeLayers((Image *) magnitude_image,DivideDstCompositeOp,
3526  complex_multiplication,0,0,exception);
3527  /*
3528  Do the IFT and return the cross-correlation result.
3529  */
3530  (void) SetImageArtifact(complex_multiplication,"fourier:normalize","inverse");
3531  cross_correlation=InverseFourierTransformImage(complex_multiplication,
3532  complex_multiplication->next,MagickFalse,exception);
3533  complex_multiplication=DestroyImageList(complex_multiplication);
3534  return(cross_correlation);
3535 }
3536 
3537 static MagickBooleanType SIMSetImageMean(Image *image,
3538  const ChannelStatistics *channel_statistics,ExceptionInfo *exception)
3539 {
3540  CacheView
3541  *image_view;
3542 
3543  MagickBooleanType
3544  status = MagickTrue;
3545 
3546  ssize_t
3547  y;
3548 
3549  /*
3550  Set image mean.
3551  */
3552  image_view=AcquireAuthenticCacheView(image,exception);
3553 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3554  #pragma omp parallel for schedule(static) shared(status) \
3555  magick_number_threads(image,image,image->rows,1)
3556 #endif
3557  for (y=0; y < (ssize_t) image->rows; y++)
3558  {
3559  Quantum
3560  *magick_restrict q;
3561 
3562  ssize_t
3563  x;
3564 
3565  if (status == MagickFalse)
3566  continue;
3567  q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
3568  if (q == (Quantum *) NULL)
3569  {
3570  status=MagickFalse;
3571  continue;
3572  }
3573  for (x=0; x < (ssize_t) image->columns; x++)
3574  {
3575  ssize_t
3576  i;
3577 
3578  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3579  {
3580  PixelChannel channel = GetPixelChannelChannel(image,i);
3581  PixelTrait traits = GetPixelChannelTraits(image,channel);
3582  if ((traits & UpdatePixelTrait) == 0)
3583  continue;
3584  q[i]=(Quantum) channel_statistics[channel].mean;
3585  }
3586  q+=(ptrdiff_t) GetPixelChannels(image);
3587  }
3588  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
3589  status=MagickFalse;
3590  }
3591  image_view=DestroyCacheView(image_view);
3592  return(status);
3593 }
3594 
3595 static Image *SIMSubtractImageMean(const Image *alpha_image,
3596  const Image *beta_image,const ChannelStatistics *channel_statistics,
3597  ExceptionInfo *exception)
3598 {
3599  CacheView
3600  *beta_view,
3601  *image_view;
3602 
3603  Image
3604  *subtract_image;
3605 
3606  MagickBooleanType
3607  status = MagickTrue;
3608 
3609  ssize_t
3610  y;
3611 
3612  /*
3613  Subtract the image mean and pad.
3614  */
3615  subtract_image=CloneImage(beta_image,alpha_image->columns,alpha_image->rows,
3616  MagickTrue,exception);
3617  if (subtract_image == (Image *) NULL)
3618  return(subtract_image);
3619  image_view=AcquireAuthenticCacheView(subtract_image,exception);
3620  beta_view=AcquireVirtualCacheView(beta_image,exception);
3621 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3622  #pragma omp parallel for schedule(static) shared(status) \
3623  magick_number_threads(beta_image,subtract_image,subtract_image->rows,1)
3624 #endif
3625  for (y=0; y < (ssize_t) subtract_image->rows; y++)
3626  {
3627  const Quantum
3628  *magick_restrict p;
3629 
3630  Quantum
3631  *magick_restrict q;
3632 
3633  ssize_t
3634  x;
3635 
3636  if (status == MagickFalse)
3637  continue;
3638  p=GetCacheViewVirtualPixels(beta_view,0,y,beta_image->columns,1,exception);
3639  q=GetCacheViewAuthenticPixels(image_view,0,y,subtract_image->columns,1,
3640  exception);
3641  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3642  {
3643  status=MagickFalse;
3644  continue;
3645  }
3646  for (x=0; x < (ssize_t) subtract_image->columns; x++)
3647  {
3648  ssize_t
3649  i;
3650 
3651  for (i=0; i < (ssize_t) GetPixelChannels(subtract_image); i++)
3652  {
3653  PixelChannel channel = GetPixelChannelChannel(subtract_image,i);
3654  PixelTrait traits = GetPixelChannelTraits(subtract_image,channel);
3655  PixelTrait beta_traits = GetPixelChannelTraits(beta_image,channel);
3656  if (((traits & UpdatePixelTrait) == 0) ||
3657  ((beta_traits & UpdatePixelTrait) == 0))
3658  continue;
3659  if ((x >= (ssize_t) beta_image->columns) ||
3660  (y >= (ssize_t) beta_image->rows))
3661  q[i]=(Quantum) 0;
3662  else
3663  q[i]=(Quantum) ((double) GetPixelChannel(beta_image,channel,p)-
3664  channel_statistics[channel].mean);
3665  }
3666  p+=(ptrdiff_t) GetPixelChannels(beta_image);
3667  q+=(ptrdiff_t) GetPixelChannels(subtract_image);
3668  }
3669  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
3670  status=MagickFalse;
3671  }
3672  beta_view=DestroyCacheView(beta_view);
3673  image_view=DestroyCacheView(image_view);
3674  if (status == MagickFalse)
3675  subtract_image=DestroyImage(subtract_image);
3676  return(subtract_image);
3677 }
3678 
3679 static Image *SIMUnityImage(const Image *alpha_image,const Image *beta_image,
3680  ExceptionInfo *exception)
3681 {
3682  CacheView
3683  *image_view;
3684 
3685  Image
3686  *unity_image;
3687 
3688  MagickBooleanType
3689  status = MagickTrue;
3690 
3691  ssize_t
3692  y;
3693 
3694  /*
3695  Create a padded unity image.
3696  */
3697  unity_image=CloneImage(alpha_image,alpha_image->columns,alpha_image->rows,
3698  MagickTrue,exception);
3699  if (unity_image == (Image *) NULL)
3700  return(unity_image);
3701  if (SetImageStorageClass(unity_image,DirectClass,exception) == MagickFalse)
3702  return(DestroyImage(unity_image));
3703  image_view=AcquireAuthenticCacheView(unity_image,exception);
3704 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3705  #pragma omp parallel for schedule(static) shared(status) \
3706  magick_number_threads(unity_image,unity_image,unity_image->rows,1)
3707 #endif
3708  for (y=0; y < (ssize_t) unity_image->rows; y++)
3709  {
3710  Quantum
3711  *magick_restrict q;
3712 
3713  ssize_t
3714  x;
3715 
3716  if (status == MagickFalse)
3717  continue;
3718  q=GetCacheViewAuthenticPixels(image_view,0,y,unity_image->columns,1,
3719  exception);
3720  if (q == (Quantum *) NULL)
3721  {
3722  status=MagickFalse;
3723  continue;
3724  }
3725  for (x=0; x < (ssize_t) unity_image->columns; x++)
3726  {
3727  ssize_t
3728  i;
3729 
3730  for (i=0; i < (ssize_t) GetPixelChannels(unity_image); i++)
3731  {
3732  PixelChannel channel = GetPixelChannelChannel(unity_image,i);
3733  PixelTrait traits = GetPixelChannelTraits(unity_image,channel);
3734  if ((traits & UpdatePixelTrait) == 0)
3735  continue;
3736  if ((x >= (ssize_t) beta_image->columns) ||
3737  (y >= (ssize_t) beta_image->rows))
3738  q[i]=(Quantum) 0;
3739  else
3740  q[i]=QuantumRange;
3741  }
3742  q+=(ptrdiff_t) GetPixelChannels(unity_image);
3743  }
3744  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
3745  status=MagickFalse;
3746  }
3747  image_view=DestroyCacheView(image_view);
3748  if (status == MagickFalse)
3749  unity_image=DestroyImage(unity_image);
3750  return(unity_image);
3751 }
3752 
3753 static Image *SIMVarianceImage(Image *alpha_image,const Image *beta_image,
3754  ExceptionInfo *exception)
3755 {
3756  CacheView
3757  *beta_view,
3758  *image_view;
3759 
3760  Image
3761  *variance_image;
3762 
3763  MagickBooleanType
3764  status = MagickTrue;
3765 
3766  ssize_t
3767  y;
3768 
3769  /*
3770  Compute the variance of the two images.
3771  */
3772  variance_image=CloneImage(alpha_image,0,0,MagickTrue,exception);
3773  if (variance_image == (Image *) NULL)
3774  return(variance_image);
3775  image_view=AcquireAuthenticCacheView(variance_image,exception);
3776  beta_view=AcquireVirtualCacheView(beta_image,exception);
3777 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3778  #pragma omp parallel for schedule(static) shared(status) \
3779  magick_number_threads(beta_image,variance_image,variance_image->rows,1)
3780 #endif
3781  for (y=0; y < (ssize_t) variance_image->rows; y++)
3782  {
3783  const Quantum
3784  *magick_restrict p;
3785 
3786  Quantum
3787  *magick_restrict q;
3788 
3789  ssize_t
3790  x;
3791 
3792  if (status == MagickFalse)
3793  continue;
3794  p=GetCacheViewVirtualPixels(beta_view,0,y,beta_image->columns,1,
3795  exception);
3796  q=GetCacheViewAuthenticPixels(image_view,0,y,variance_image->columns,1,
3797  exception);
3798  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3799  {
3800  status=MagickFalse;
3801  continue;
3802  }
3803  for (x=0; x < (ssize_t) variance_image->columns; x++)
3804  {
3805  ssize_t
3806  i;
3807 
3808  for (i=0; i < (ssize_t) GetPixelChannels(variance_image); i++)
3809  {
3810  double
3811  error;
3812 
3813  PixelChannel channel = GetPixelChannelChannel(variance_image,i);
3814  PixelTrait traits = GetPixelChannelTraits(variance_image,channel);
3815  PixelTrait beta_traits = GetPixelChannelTraits(beta_image,channel);
3816  if (((traits & UpdatePixelTrait) == 0) ||
3817  ((beta_traits & UpdatePixelTrait) == 0))
3818  continue;
3819  error=(double) q[i]-(double) GetPixelChannel(beta_image,channel,p);
3820  q[i]=(Quantum) ((double) ClampToQuantum((double) QuantumRange*
3821  (sqrt(fabs(QuantumScale*error))/sqrt((double) QuantumRange))));
3822  }
3823  p+=(ptrdiff_t) GetPixelChannels(beta_image);
3824  q+=(ptrdiff_t) GetPixelChannels(variance_image);
3825  }
3826  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
3827  status=MagickFalse;
3828  }
3829  beta_view=DestroyCacheView(beta_view);
3830  image_view=DestroyCacheView(image_view);
3831  if (status == MagickFalse)
3832  variance_image=DestroyImage(variance_image);
3833  return(variance_image);
3834 }
3835 
3836 static Image *DPCSimilarityImage(const Image *image,const Image *reconstruct,
3837  RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
3838 {
3839 #define ThrowDPCSimilarityException() \
3840 { \
3841  if (dot_product_image != (Image *) NULL) \
3842  dot_product_image=DestroyImage(dot_product_image); \
3843  if (magnitude_image != (Image *) NULL) \
3844  magnitude_image=DestroyImage(magnitude_image); \
3845  if (reconstruct_image != (Image *) NULL) \
3846  reconstruct_image=DestroyImage(reconstruct_image); \
3847  if (rx_image != (Image *) NULL) \
3848  rx_image=DestroyImage(rx_image); \
3849  if (ry_image != (Image *) NULL) \
3850  ry_image=DestroyImage(ry_image); \
3851  if (target_image != (Image *) NULL) \
3852  target_image=DestroyImage(target_image); \
3853  if (threshold_image != (Image *) NULL) \
3854  threshold_image=DestroyImage(threshold_image); \
3855  if (trx_image != (Image *) NULL) \
3856  trx_image=DestroyImage(trx_image); \
3857  if (try_image != (Image *) NULL) \
3858  try_image=DestroyImage(try_image); \
3859  if (tx_image != (Image *) NULL) \
3860  tx_image=DestroyImage(tx_image); \
3861  if (ty_image != (Image *) NULL) \
3862  ty_image=DestroyImage(ty_image); \
3863  return((Image *) NULL); \
3864 }
3865 
3866  double
3867  edge_factor = 0.0,
3868  maxima = 0.0,
3869  mean = 0.0,
3870  standard_deviation = 0.0;
3871 
3872  Image
3873  *dot_product_image = (Image *) NULL,
3874  *magnitude_image = (Image *) NULL,
3875  *reconstruct_image = (Image *) NULL,
3876  *rx_image = (Image *) NULL,
3877  *ry_image = (Image *) NULL,
3878  *trx_image = (Image *) NULL,
3879  *target_image = (Image *) NULL,
3880  *threshold_image = (Image *) NULL,
3881  *try_image = (Image *) NULL,
3882  *tx_image = (Image *) NULL,
3883  *ty_image = (Image *) NULL;
3884 
3885  MagickBooleanType
3886  status = MagickTrue;
3887 
3888  RectangleInfo
3889  geometry;
3890 
3891  /*
3892  Dot product correlation-based image similarity using FFT local statistics.
3893  */
3894  target_image=CloneImage(image,0,0,MagickTrue,exception);
3895  if (target_image == (Image *) NULL)
3896  return((Image *) NULL);
3897  /*
3898  Compute the cross correlation of the test and reconstruct magnitudes.
3899  */
3900  reconstruct_image=CloneImage(reconstruct,0,0,MagickTrue,exception);
3901  if (reconstruct_image == (Image *) NULL)
3902  ThrowDPCSimilarityException();
3903  /*
3904  Compute X and Y derivatives of reference image.
3905  */
3906  (void) SetImageVirtualPixelMethod(reconstruct_image,EdgeVirtualPixelMethod,
3907  exception);
3908  rx_image=SIMDerivativeImage(reconstruct_image,"Sobel",exception);
3909  if (rx_image == (Image *) NULL)
3910  ThrowDPCSimilarityException();
3911  ry_image=SIMDerivativeImage(reconstruct_image,"Sobel:90",exception);
3912  reconstruct_image=DestroyImage(reconstruct_image);
3913  if (ry_image == (Image *) NULL)
3914  ThrowDPCSimilarityException();
3915  /*
3916  Compute magnitude of derivatives.
3917  */
3918  magnitude_image=SIMMagnitudeImage(rx_image,ry_image,exception);
3919  if (magnitude_image == (Image *) NULL)
3920  ThrowDPCSimilarityException();
3921  /*
3922  Compute an edge normalization correction.
3923  */
3924  threshold_image=CloneImage(magnitude_image,0,0,MagickTrue,exception);
3925  if (threshold_image == (Image *) NULL)
3926  ThrowDPCSimilarityException();
3927  status=BilevelImage(threshold_image,0.0,exception);
3928  if (status == MagickFalse)
3929  ThrowDPCSimilarityException();
3930  status=GetImageMean(threshold_image,&mean,&standard_deviation,exception);
3931  threshold_image=DestroyImage(threshold_image);
3932  if (status == MagickFalse)
3933  ThrowDPCSimilarityException();
3934  edge_factor=MagickSafeReciprocal(QuantumScale*mean*reconstruct->columns*
3935  reconstruct->rows)+QuantumScale;
3936  /*
3937  Divide X and Y derivitives of reference image by magnitude.
3938  */
3939  trx_image=SIMDivideByMagnitude(rx_image,magnitude_image,image,exception);
3940  rx_image=DestroyImage(rx_image);
3941  if (trx_image == (Image *) NULL)
3942  ThrowDPCSimilarityException();
3943  rx_image=trx_image;
3944  try_image=SIMDivideByMagnitude(ry_image,magnitude_image,image,exception);
3945  magnitude_image=DestroyImage(magnitude_image);
3946  ry_image=DestroyImage(ry_image);
3947  if (try_image == (Image *) NULL)
3948  ThrowDPCSimilarityException();
3949  ry_image=try_image;
3950  /*
3951  Compute X and Y derivatives of image.
3952  */
3953  (void) SetImageVirtualPixelMethod(target_image,EdgeVirtualPixelMethod,
3954  exception);
3955  tx_image=SIMDerivativeImage(target_image,"Sobel",exception);
3956  if (tx_image == (Image *) NULL)
3957  ThrowDPCSimilarityException();
3958  ty_image=SIMDerivativeImage(target_image,"Sobel:90",exception);
3959  target_image=DestroyImage(target_image);
3960  if (ty_image == (Image *) NULL)
3961  ThrowDPCSimilarityException();
3962  /*
3963  Compute magnitude of derivatives.
3964  */
3965  magnitude_image=SIMMagnitudeImage(tx_image,ty_image,exception);
3966  if (magnitude_image == (Image *) NULL)
3967  ThrowDPCSimilarityException();
3968  /*
3969  Divide Lx and Ly by magnitude.
3970  */
3971  trx_image=SIMDivideByMagnitude(tx_image,magnitude_image,image,exception);
3972  tx_image=DestroyImage(tx_image);
3973  if (trx_image == (Image *) NULL)
3974  ThrowDPCSimilarityException();
3975  tx_image=trx_image;
3976  try_image=SIMDivideByMagnitude(ty_image,magnitude_image,image,exception);
3977  ty_image=DestroyImage(ty_image);
3978  magnitude_image=DestroyImage(magnitude_image);
3979  if (try_image == (Image *) NULL)
3980  ThrowDPCSimilarityException();
3981  ty_image=try_image;
3982  /*
3983  Compute the cross correlation of the test and reference images.
3984  */
3985  trx_image=SIMCrossCorrelationImage(tx_image,rx_image,exception);
3986  rx_image=DestroyImage(rx_image);
3987  tx_image=DestroyImage(tx_image);
3988  if (trx_image == (Image *) NULL)
3989  ThrowDPCSimilarityException();
3990  try_image=SIMCrossCorrelationImage(ty_image,ry_image,exception);
3991  ry_image=DestroyImage(ry_image);
3992  ty_image=DestroyImage(ty_image);
3993  if (try_image == (Image *) NULL)
3994  ThrowDPCSimilarityException();
3995  /*
3996  Evaluate dot product correlation image.
3997  */
3998  (void) SetImageArtifact(try_image,"compose:clamp","false");
3999  status=CompositeImage(trx_image,try_image,PlusCompositeOp,MagickTrue,0,0,
4000  exception);
4001  try_image=DestroyImage(try_image);
4002  if (status == MagickFalse)
4003  ThrowDPCSimilarityException();
4004  status=SIMMultiplyImage(trx_image,edge_factor,
4005  (const ChannelStatistics *) NULL,exception);
4006  if (status == MagickFalse)
4007  ThrowDPCSimilarityException();
4008  /*
4009  Crop results.
4010  */
4011  SetGeometry(image,&geometry);
4012  geometry.width=image->columns;
4013  geometry.height=image->rows;
4014  (void) ResetImagePage(trx_image,"0x0+0+0");
4015  dot_product_image=CropImage(trx_image,&geometry,exception);
4016  trx_image=DestroyImage(trx_image);
4017  if (dot_product_image == (Image *) NULL)
4018  ThrowDPCSimilarityException();
4019  (void) ResetImagePage(dot_product_image,"0x0+0+0");
4020  /*
4021  Identify the maxima value in the image and its location.
4022  */
4023  status=GrayscaleImage(dot_product_image,AveragePixelIntensityMethod,
4024  exception);
4025  if (status == MagickFalse)
4026  ThrowDPCSimilarityException();
4027  dot_product_image->depth=32;
4028  dot_product_image->colorspace=GRAYColorspace;
4029  dot_product_image->alpha_trait=UndefinedPixelTrait;
4030  status=SIMFilterImageNaNs(dot_product_image,exception);
4031  if (status == MagickFalse)
4032  ThrowDPCSimilarityException();
4033  status=SIMMaximaImage(dot_product_image,&maxima,offset,exception);
4034  if (status == MagickFalse)
4035  ThrowDPCSimilarityException();
4036  if ((QuantumScale*maxima) > 1.0)
4037  {
4038  status=SIMMultiplyImage(dot_product_image,1.0/(QuantumScale*maxima),
4039  (const ChannelStatistics *) NULL,exception);
4040  maxima=(double) QuantumRange;
4041  }
4042  *similarity_metric=QuantumScale*maxima;
4043  return(dot_product_image);
4044 }
4045 
4046 static Image *MSESimilarityImage(const Image *image,const Image *reconstruct,
4047  RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
4048 {
4049 #define ThrowMSESimilarityException() \
4050 { \
4051  if (alpha_image != (Image *) NULL) \
4052  alpha_image=DestroyImage(alpha_image); \
4053  if (beta_image != (Image *) NULL) \
4054  beta_image=DestroyImage(beta_image); \
4055  if (channel_statistics != (ChannelStatistics *) NULL) \
4056  channel_statistics=(ChannelStatistics *) \
4057  RelinquishMagickMemory(channel_statistics); \
4058  if (mean_image != (Image *) NULL) \
4059  mean_image=DestroyImage(mean_image); \
4060  if (mse_image != (Image *) NULL) \
4061  mse_image=DestroyImage(mse_image); \
4062  if (reconstruct_image != (Image *) NULL) \
4063  reconstruct_image=DestroyImage(reconstruct_image); \
4064  if (sum_image != (Image *) NULL) \
4065  sum_image=DestroyImage(sum_image); \
4066  if (alpha_image != (Image *) NULL) \
4067  alpha_image=DestroyImage(alpha_image); \
4068  return((Image *) NULL); \
4069 }
4070 
4071  ChannelStatistics
4072  *channel_statistics = (ChannelStatistics *) NULL;
4073 
4074  double
4075  minima = 0.0;
4076 
4077  Image
4078  *alpha_image = (Image *) NULL,
4079  *beta_image = (Image *) NULL,
4080  *mean_image = (Image *) NULL,
4081  *mse_image = (Image *) NULL,
4082  *reconstruct_image = (Image *) NULL,
4083  *sum_image = (Image *) NULL,
4084  *target_image = (Image *) NULL;
4085 
4086  MagickBooleanType
4087  status = MagickTrue;
4088 
4089  RectangleInfo
4090  geometry;
4091 
4092  /*
4093  MSE correlation-based image similarity using FFT local statistics.
4094  */
4095  target_image=SIMSquareImage(image,exception);
4096  if (target_image == (Image *) NULL)
4097  ThrowMSESimilarityException();
4098  reconstruct_image=SIMUnityImage(image,reconstruct,exception);
4099  if (reconstruct_image == (Image *) NULL)
4100  ThrowMSESimilarityException();
4101  /*
4102  Create (U * test)/# pixels.
4103  */
4104  alpha_image=SIMCrossCorrelationImage(target_image,reconstruct_image,
4105  exception);
4106  target_image=DestroyImage(target_image);
4107  if (alpha_image == (Image *) NULL)
4108  ThrowMSESimilarityException();
4109  status=SIMMultiplyImage(alpha_image,1.0/reconstruct->columns/(double)
4110  reconstruct->rows,(const ChannelStatistics *) NULL,exception);
4111  if (status == MagickFalse)
4112  ThrowMSESimilarityException();
4113  /*
4114  Create 2*(test * reconstruction)# pixels.
4115  */
4116  (void) CompositeImage(reconstruct_image,reconstruct,CopyCompositeOp,
4117  MagickTrue,0,0,exception);
4118  beta_image=SIMCrossCorrelationImage(image,reconstruct_image,exception);
4119  if (beta_image == (Image *) NULL)
4120  {
4121  reconstruct_image=DestroyImage(reconstruct_image);
4122  ThrowMSESimilarityException();
4123  }
4124  status=SIMMultiplyImage(beta_image,-2.0/reconstruct->columns/(double)
4125  reconstruct->rows,(const ChannelStatistics *) NULL,exception);
4126  reconstruct_image=DestroyImage(reconstruct_image);
4127  if (status == MagickFalse)
4128  ThrowMSESimilarityException();
4129  /*
4130  Mean of reconstruction squared.
4131  */
4132  sum_image=SIMSquareImage(reconstruct,exception);
4133  if (sum_image == (Image *) NULL)
4134  ThrowMSESimilarityException();
4135  channel_statistics=GetImageStatistics(sum_image,exception);
4136  if (channel_statistics == (ChannelStatistics *) NULL)
4137  ThrowMSESimilarityException();
4138  status=SetImageStorageClass(sum_image,DirectClass,exception);
4139  if (status == MagickFalse)
4140  ThrowMSESimilarityException();
4141  status=SIMSetImageMean(sum_image,channel_statistics,exception);
4142  channel_statistics=(ChannelStatistics *)
4143  RelinquishMagickMemory(channel_statistics);
4144  if (status == MagickFalse)
4145  ThrowMSESimilarityException();
4146  /*
4147  Create mean image.
4148  */
4149  AppendImageToList(&sum_image,alpha_image);
4150  AppendImageToList(&sum_image,beta_image);
4151  mean_image=EvaluateImages(sum_image,SumEvaluateOperator,exception);
4152  if (mean_image == (Image *) NULL)
4153  ThrowMSESimilarityException();
4154  sum_image=DestroyImage(sum_image);
4155  status=GrayscaleImage(mean_image,AveragePixelIntensityMethod,exception);
4156  if (status == MagickFalse)
4157  ThrowMSESimilarityException();
4158  /*
4159  Crop to difference of reconstruction and test images.
4160  */
4161  SetGeometry(image,&geometry);
4162  geometry.width=image->columns;
4163  geometry.height=image->rows;
4164  (void) ResetImagePage(mean_image,"0x0+0+0");
4165  mse_image=CropImage(mean_image,&geometry,exception);
4166  mean_image=DestroyImage(mean_image);
4167  if (mse_image == (Image *) NULL)
4168  ThrowMSESimilarityException();
4169  /*
4170  Identify the minima value in the correlation image and its location.
4171  */
4172  (void) ResetImagePage(mse_image,"0x0+0+0");
4173  (void) ClampImage(mse_image,exception);
4174  mse_image->depth=32;
4175  mse_image->colorspace=GRAYColorspace;
4176  mse_image->alpha_trait=UndefinedPixelTrait;
4177  status=SIMMinimaImage(mse_image,&minima,offset,exception);
4178  if (status == MagickFalse)
4179  ThrowMSESimilarityException();
4180  status=NegateImage(mse_image,MagickFalse,exception);
4181  if (status == MagickFalse)
4182  ThrowMSESimilarityException();
4183  alpha_image=DestroyImage(alpha_image);
4184  beta_image=DestroyImage(beta_image);
4185  if ((QuantumScale*minima) < FLT_EPSILON)
4186  minima=0.0;
4187  *similarity_metric=QuantumScale*minima;
4188  return(mse_image);
4189 }
4190 
4191 static Image *NCCSimilarityImage(const Image *image,const Image *reconstruct,
4192  RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
4193 {
4194 #define ThrowNCCSimilarityException() \
4195 { \
4196  if (alpha_image != (Image *) NULL) \
4197  alpha_image=DestroyImage(alpha_image); \
4198  if (beta_image != (Image *) NULL) \
4199  beta_image=DestroyImage(beta_image); \
4200  if (channel_statistics != (ChannelStatistics *) NULL) \
4201  channel_statistics=(ChannelStatistics *) \
4202  RelinquishMagickMemory(channel_statistics); \
4203  if (correlation_image != (Image *) NULL) \
4204  correlation_image=DestroyImage(correlation_image); \
4205  if (divide_image != (Image *) NULL) \
4206  divide_image=DestroyImage(divide_image); \
4207  if (ncc_image != (Image *) NULL) \
4208  ncc_image=DestroyImage(ncc_image); \
4209  if (normalize_image != (Image *) NULL) \
4210  normalize_image=DestroyImage(normalize_image); \
4211  if (reconstruct_image != (Image *) NULL) \
4212  reconstruct_image=DestroyImage(reconstruct_image); \
4213  if (target_image != (Image *) NULL) \
4214  target_image=DestroyImage(target_image); \
4215  if (variance_image != (Image *) NULL) \
4216  variance_image=DestroyImage(variance_image); \
4217  return((Image *) NULL); \
4218 }
4219 
4220  ChannelStatistics
4221  *channel_statistics = (ChannelStatistics *) NULL;
4222 
4223  double
4224  maxima = 0.0;
4225 
4226  Image
4227  *alpha_image = (Image *) NULL,
4228  *beta_image = (Image *) NULL,
4229  *correlation_image = (Image *) NULL,
4230  *divide_image = (Image *) NULL,
4231  *ncc_image = (Image *) NULL,
4232  *normalize_image = (Image *) NULL,
4233  *reconstruct_image = (Image *) NULL,
4234  *target_image = (Image *) NULL,
4235  *variance_image = (Image *) NULL;
4236 
4237  MagickBooleanType
4238  status = MagickTrue;
4239 
4240  RectangleInfo
4241  geometry;
4242 
4243  /*
4244  NCC correlation-based image similarity with FFT local statistics.
4245  */
4246  target_image=SIMSquareImage(image,exception);
4247  if (target_image == (Image *) NULL)
4248  ThrowNCCSimilarityException();
4249  reconstruct_image=SIMUnityImage(image,reconstruct,exception);
4250  if (reconstruct_image == (Image *) NULL)
4251  ThrowNCCSimilarityException();
4252  /*
4253  Compute the cross correlation of the test and reconstruction images.
4254  */
4255  alpha_image=SIMCrossCorrelationImage(target_image,reconstruct_image,
4256  exception);
4257  target_image=DestroyImage(target_image);
4258  if (alpha_image == (Image *) NULL)
4259  ThrowNCCSimilarityException();
4260  status=SIMMultiplyImage(alpha_image,(double) QuantumRange*
4261  reconstruct->columns*reconstruct->rows,(const ChannelStatistics *) NULL,
4262  exception);
4263  if (status == MagickFalse)
4264  ThrowNCCSimilarityException();
4265  /*
4266  Compute the cross correlation of the source and reconstruction images.
4267  */
4268  beta_image=SIMCrossCorrelationImage(image,reconstruct_image,exception);
4269  reconstruct_image=DestroyImage(reconstruct_image);
4270  if (beta_image == (Image *) NULL)
4271  ThrowNCCSimilarityException();
4272  target_image=SIMSquareImage(beta_image,exception);
4273  beta_image=DestroyImage(beta_image);
4274  if (target_image == (Image *) NULL)
4275  ThrowNCCSimilarityException();
4276  status=SIMMultiplyImage(target_image,(double) QuantumRange,
4277  (const ChannelStatistics *) NULL,exception);
4278  if (status == MagickFalse)
4279  ThrowNCCSimilarityException();
4280  /*
4281  Compute the variance of the two images.
4282  */
4283  variance_image=SIMVarianceImage(alpha_image,target_image,exception);
4284  target_image=DestroyImage(target_image);
4285  alpha_image=DestroyImage(alpha_image);
4286  if (variance_image == (Image *) NULL)
4287  ThrowNCCSimilarityException();
4288  /*
4289  Subtract the image mean.
4290  */
4291  channel_statistics=GetImageStatistics(reconstruct,exception);
4292  if (channel_statistics == (ChannelStatistics *) NULL)
4293  ThrowNCCSimilarityException();
4294  status=SIMMultiplyImage(variance_image,1.0,channel_statistics,exception);
4295  if (status == MagickFalse)
4296  ThrowNCCSimilarityException();
4297  normalize_image=SIMSubtractImageMean(image,reconstruct,channel_statistics,
4298  exception);
4299  channel_statistics=(ChannelStatistics *)
4300  RelinquishMagickMemory(channel_statistics);
4301  if (normalize_image == (Image *) NULL)
4302  ThrowNCCSimilarityException();
4303  correlation_image=SIMCrossCorrelationImage(image,normalize_image,exception);
4304  normalize_image=DestroyImage(normalize_image);
4305  if (correlation_image == (Image *) NULL)
4306  ThrowNCCSimilarityException();
4307  /*
4308  Divide the two images.
4309  */
4310  divide_image=SIMDivideImage(correlation_image,variance_image,exception);
4311  correlation_image=DestroyImage(correlation_image);
4312  variance_image=DestroyImage(variance_image);
4313  if (divide_image == (Image *) NULL)
4314  ThrowNCCSimilarityException();
4315  /*
4316  Crop padding.
4317  */
4318  SetGeometry(image,&geometry);
4319  geometry.width=image->columns;
4320  geometry.height=image->rows;
4321  (void) ResetImagePage(divide_image,"0x0+0+0");
4322  ncc_image=CropImage(divide_image,&geometry,exception);
4323  divide_image=DestroyImage(divide_image);
4324  if (ncc_image == (Image *) NULL)
4325  ThrowNCCSimilarityException();
4326  /*
4327  Identify the maxima value in the image and its location.
4328  */
4329  (void) ResetImagePage(ncc_image,"0x0+0+0");
4330  status=GrayscaleImage(ncc_image,AveragePixelIntensityMethod,exception);
4331  if (status == MagickFalse)
4332  ThrowNCCSimilarityException();
4333  ncc_image->depth=32;
4334  ncc_image->colorspace=GRAYColorspace;
4335  ncc_image->alpha_trait=UndefinedPixelTrait;
4336  status=SIMMaximaImage(ncc_image,&maxima,offset,exception);
4337  if (status == MagickFalse)
4338  ThrowNCCSimilarityException();
4339  if ((QuantumScale*maxima) > 1.0)
4340  {
4341  status=SIMMultiplyImage(ncc_image,1.0/(QuantumScale*maxima),
4342  (const ChannelStatistics *) NULL,exception);
4343  maxima=(double) QuantumRange;
4344  }
4345  *similarity_metric=QuantumScale*maxima;
4346  return(ncc_image);
4347 }
4348 
4349 static Image *PhaseSimilarityImage(const Image *image,const Image *reconstruct,
4350  RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
4351 {
4352 #define ThrowPhaseSimilarityException() \
4353 { \
4354  if (correlation_image != (Image *) NULL) \
4355  correlation_image=DestroyImage(correlation_image); \
4356  if (fft_images != (Image *) NULL) \
4357  fft_images=DestroyImageList(fft_images); \
4358  if (gamma_image != (Image *) NULL) \
4359  gamma_image=DestroyImage(gamma_image); \
4360  if (magnitude_image != (Image *) NULL) \
4361  magnitude_image=DestroyImage(magnitude_image); \
4362  if (phase_image != (Image *) NULL) \
4363  phase_image=DestroyImage(phase_image); \
4364  if (reconstruct_image != (Image *) NULL) \
4365  reconstruct_image=DestroyImage(reconstruct_image); \
4366  if (reconstruct_magnitude != (Image *) NULL) \
4367  reconstruct_magnitude=DestroyImage(reconstruct_magnitude); \
4368  if (target_image != (Image *) NULL) \
4369  target_image=DestroyImage(target_image); \
4370  if (test_magnitude != (Image *) NULL) \
4371  test_magnitude=DestroyImage(test_magnitude); \
4372  return((Image *) NULL); \
4373 }
4374 
4375  double
4376  maxima = 0.0;
4377 
4378  Image
4379  *correlation_image = (Image *) NULL,
4380  *fft_images = (Image *) NULL,
4381  *gamma_image = (Image *) NULL,
4382  *magnitude_image = (Image *) NULL,
4383  *phase_image = (Image *) NULL,
4384  *reconstruct_image = (Image *) NULL,
4385  *reconstruct_magnitude = (Image *) NULL,
4386  *target_image = (Image *) NULL,
4387  *test_magnitude = (Image *) NULL;
4388 
4389  MagickBooleanType
4390  status = MagickTrue;
4391 
4392  RectangleInfo
4393  geometry;
4394 
4395  /*
4396  Phase correlation-based image similarity using FFT local statistics.
4397  */
4398  target_image=CloneImage(image,0,0,MagickTrue,exception);
4399  if (target_image == (Image *) NULL)
4400  ThrowPhaseSimilarityException();
4401  (void) ResetImagePage(target_image,"0x0+0+0");
4402  GetPixelInfoRGBA((Quantum) 0,(Quantum) 0,(Quantum) 0,(Quantum) 0,
4403  &target_image->background_color);
4404  status=SetImageExtent(target_image,2*(size_t) ceil((double) image->columns/
4405  2.0),2*(size_t) ceil((double) image->rows/2.0),exception);
4406  if (status == MagickFalse)
4407  ThrowPhaseSimilarityException();
4408  /*
4409  Compute the cross correlation of the test and reconstruct magnitudes.
4410  */
4411  reconstruct_image=CloneImage(reconstruct,0,0,MagickTrue,exception);
4412  if (reconstruct_image == (Image *) NULL)
4413  ThrowPhaseSimilarityException();
4414  (void) ResetImagePage(reconstruct_image,"0x0+0+0");
4415  GetPixelInfoRGBA((Quantum) 0,(Quantum) 0,(Quantum) 0,(Quantum) 0,
4416  &reconstruct_image->background_color);
4417  status=SetImageExtent(reconstruct_image,2*(size_t) ceil((double)
4418  image->columns/2.0),2*(size_t) ceil((double) image->rows/2.0),exception);
4419  if (status == MagickFalse)
4420  ThrowPhaseSimilarityException();
4421  /*
4422  Evaluate phase coorelation image and divide by the product magnitude.
4423  */
4424  (void) SetImageArtifact(target_image,"fourier:normalize","inverse");
4425  fft_images=ForwardFourierTransformImage(target_image,MagickTrue,exception);
4426  if (fft_images == (Image *) NULL)
4427  ThrowPhaseSimilarityException();
4428  test_magnitude=CloneImage(fft_images,0,0,MagickTrue,exception);
4429  fft_images=DestroyImageList(fft_images);
4430  if (test_magnitude == (Image *) NULL)
4431  ThrowPhaseSimilarityException();
4432  (void) SetImageArtifact(reconstruct_image,"fourier:normalize","inverse");
4433  fft_images=ForwardFourierTransformImage(reconstruct_image,MagickTrue,
4434  exception);
4435  if (fft_images == (Image *) NULL)
4436  ThrowPhaseSimilarityException();
4437  reconstruct_magnitude=CloneImage(fft_images,0,0,MagickTrue,exception);
4438  fft_images=DestroyImageList(fft_images);
4439  if (reconstruct_magnitude == (Image *) NULL)
4440  ThrowPhaseSimilarityException();
4441  magnitude_image=CloneImage(reconstruct_magnitude,0,0,MagickTrue,exception);
4442  if (magnitude_image == (Image *) NULL)
4443  ThrowPhaseSimilarityException();
4444  DisableCompositeClampUnlessSpecified(magnitude_image);
4445  (void) CompositeImage(magnitude_image,test_magnitude,MultiplyCompositeOp,
4446  MagickTrue,0,0,exception);
4447  /*
4448  Compute the cross correlation of the test and reconstruction images.
4449  */
4450  correlation_image=SIMPhaseCorrelationImage(target_image,reconstruct_image,
4451  magnitude_image,exception);
4452  target_image=DestroyImage(target_image);
4453  reconstruct_image=DestroyImage(reconstruct_image);
4454  test_magnitude=DestroyImage(test_magnitude);
4455  reconstruct_magnitude=DestroyImage(reconstruct_magnitude);
4456  if (correlation_image == (Image *) NULL)
4457  ThrowPhaseSimilarityException();
4458  /*
4459  Identify the maxima value in the image and its location.
4460  */
4461  status=SIMMultiplyImage(correlation_image,(double) QuantumScale,
4462  (const ChannelStatistics *) NULL,exception);
4463  if (status == MagickFalse)
4464  ThrowPhaseSimilarityException();
4465  gamma_image=CloneImage(correlation_image,0,0,MagickTrue,exception);
4466  correlation_image=DestroyImage(correlation_image);
4467  if (gamma_image == (Image *) NULL)
4468  ThrowPhaseSimilarityException();
4469  /*
4470  Crop padding.
4471  */
4472  SetGeometry(image,&geometry);
4473  geometry.width=image->columns;
4474  geometry.height=image->rows;
4475  (void) ResetImagePage(gamma_image,"0x0+0+0");
4476  phase_image=CropImage(gamma_image,&geometry,exception);
4477  gamma_image=DestroyImage(gamma_image);
4478  if (phase_image == (Image *) NULL)
4479  ThrowPhaseSimilarityException();
4480  (void) ResetImagePage(phase_image,"0x0+0+0");
4481  /*
4482  Identify the maxima value in the correlation image and its location.
4483  */
4484  status=GrayscaleImage(phase_image,AveragePixelIntensityMethod,exception);
4485  if (status == MagickFalse)
4486  ThrowPhaseSimilarityException();
4487  phase_image->depth=32;
4488  phase_image->colorspace=GRAYColorspace;
4489  phase_image->alpha_trait=UndefinedPixelTrait;
4490  status=SIMFilterImageNaNs(phase_image,exception);
4491  if (status == MagickFalse)
4492  ThrowPhaseSimilarityException();
4493  status=SIMMaximaImage(phase_image,&maxima,offset,exception);
4494  if (status == MagickFalse)
4495  ThrowPhaseSimilarityException();
4496  magnitude_image=DestroyImage(magnitude_image);
4497  *similarity_metric=QuantumScale*maxima;
4498  return(phase_image);
4499 }
4500 
4501 static Image *PSNRSimilarityImage(const Image *image,const Image *reconstruct,
4502  RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
4503 {
4504  Image
4505  *psnr_image = (Image *) NULL;
4506 
4507  psnr_image=MSESimilarityImage(image,reconstruct,offset,similarity_metric,
4508  exception);
4509  if (psnr_image == (Image *) NULL)
4510  return(psnr_image);
4511  *similarity_metric=10.0*MagickSafeLog10(MagickSafeReciprocal(
4512  *similarity_metric))/MagickSafePSNRRecipicol(10.0);
4513  return(psnr_image);
4514 }
4515 
4516 static Image *RMSESimilarityImage(const Image *image,const Image *reconstruct,
4517  RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
4518 {
4519  Image
4520  *rmse_image = (Image *) NULL;
4521 
4522  rmse_image=MSESimilarityImage(image,reconstruct,offset,similarity_metric,
4523  exception);
4524  if (rmse_image == (Image *) NULL)
4525  return(rmse_image);
4526  *similarity_metric=sqrt(*similarity_metric);
4527  return(rmse_image);
4528 }
4529 #endif
4530 
4531 static double GetSimilarityMetric(const Image *image,
4532  const Image *reconstruct_image,const MetricType metric,
4533  const ssize_t x_offset,const ssize_t y_offset,ExceptionInfo *exception)
4534 {
4535  double
4536  *channel_similarity,
4537  similarity = 0.0;
4538 
4539  ExceptionInfo
4540  *sans_exception = AcquireExceptionInfo();
4541 
4542  Image
4543  *similarity_image;
4544 
4545  MagickBooleanType
4546  status = MagickTrue;
4547 
4548  RectangleInfo
4549  geometry;
4550 
4551  size_t
4552  length = MaxPixelChannels+1UL;
4553 
4554  SetGeometry(reconstruct_image,&geometry);
4555  geometry.x=x_offset;
4556  geometry.y=y_offset;
4557  similarity_image=CropImage(image,&geometry,sans_exception);
4558  sans_exception=DestroyExceptionInfo(sans_exception);
4559  if (similarity_image == (Image *) NULL)
4560  return(NAN);
4561  /*
4562  Get image distortion.
4563  */
4564  channel_similarity=(double *) AcquireQuantumMemory(length,
4565  sizeof(*channel_similarity));
4566  if (channel_similarity == (double *) NULL)
4567  ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
4568  (void) memset(channel_similarity,0,length*sizeof(*channel_similarity));
4569  switch (metric)
4570  {
4571  case AbsoluteErrorMetric:
4572  {
4573  status=GetAESimilarity(similarity_image,reconstruct_image,
4574  channel_similarity,exception);
4575  break;
4576  }
4577  case DotProductCorrelationErrorMetric:
4578  {
4579  status=GetDPCSimilarity(similarity_image,reconstruct_image,
4580  channel_similarity,exception);
4581  break;
4582  }
4583  case FuzzErrorMetric:
4584  {
4585  status=GetFUZZSimilarity(similarity_image,reconstruct_image,
4586  channel_similarity,exception);
4587  break;
4588  }
4589  case MeanAbsoluteErrorMetric:
4590  {
4591  status=GetMAESimilarity(similarity_image,reconstruct_image,
4592  channel_similarity,exception);
4593  break;
4594  }
4595  case MeanErrorPerPixelErrorMetric:
4596  {
4597  status=GetMEPPSimilarity(similarity_image,reconstruct_image,
4598  channel_similarity,exception);
4599  break;
4600  }
4601  case MeanSquaredErrorMetric:
4602  {
4603  status=GetMSESimilarity(similarity_image,reconstruct_image,
4604  channel_similarity,exception);
4605  break;
4606  }
4607  case NormalizedCrossCorrelationErrorMetric:
4608  {
4609  status=GetNCCSimilarity(similarity_image,reconstruct_image,
4610  channel_similarity,exception);
4611  break;
4612  }
4613  case PeakAbsoluteErrorMetric:
4614  {
4615  status=GetPASimilarity(similarity_image,reconstruct_image,
4616  channel_similarity,exception);
4617  break;
4618  }
4619  case PeakSignalToNoiseRatioErrorMetric:
4620  {
4621  status=GetPSNRSimilarity(similarity_image,reconstruct_image,
4622  channel_similarity,exception);
4623  break;
4624  }
4625  case PerceptualHashErrorMetric:
4626  {
4627  status=GetPHASHSimilarity(similarity_image,reconstruct_image,
4628  channel_similarity,exception);
4629  break;
4630  }
4631  case PhaseCorrelationErrorMetric:
4632  {
4633  status=GetPHASESimilarity(similarity_image,reconstruct_image,
4634  channel_similarity,exception);
4635  break;
4636  }
4637  case PixelDifferenceCountErrorMetric:
4638  {
4639  status=GetPDCSimilarity(similarity_image,reconstruct_image,
4640  channel_similarity,exception);
4641  break;
4642  }
4643  case RootMeanSquaredErrorMetric:
4644  case UndefinedErrorMetric:
4645  default:
4646  {
4647  status=GetRMSESimilarity(similarity_image,reconstruct_image,
4648  channel_similarity,exception);
4649  break;
4650  }
4651  case StructuralDissimilarityErrorMetric:
4652  {
4653  status=GetDSSIMSimilarity(similarity_image,reconstruct_image,
4654  channel_similarity,exception);
4655  break;
4656  }
4657  case StructuralSimilarityErrorMetric:
4658  {
4659  status=GetSSIMSimularity(similarity_image,reconstruct_image,
4660  channel_similarity,exception);
4661  break;
4662  }
4663  }
4664  similarity_image=DestroyImage(similarity_image);
4665  similarity=channel_similarity[CompositePixelChannel];
4666  channel_similarity=(double *) RelinquishMagickMemory(channel_similarity);
4667  if (status == MagickFalse)
4668  return(NAN);
4669  return(similarity);
4670 }
4671 
4672 MagickExport Image *SimilarityImage(const Image *image,const Image *reconstruct,
4673  const MetricType metric,const double similarity_threshold,
4674  RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
4675 {
4676 #define SimilarityImageTag "Similarity/Image"
4677 
4678  typedef struct
4679  {
4680  double
4681  similarity;
4682 
4683  ssize_t
4684  x,
4685  y;
4686  } SimilarityInfo;
4687 
4688  CacheView
4689  *similarity_view;
4690 
4691  Image
4692  *similarity_image = (Image *) NULL;
4693 
4694  MagickBooleanType
4695  status = MagickTrue;
4696 
4697  MagickOffsetType
4698  progress = 0;
4699 
4700  SimilarityInfo
4701  similarity_info = { 0.0, 0, 0 };
4702 
4703  size_t
4704  columns,
4705  rows;
4706 
4707  ssize_t
4708  y;
4709 
4710  assert(image != (const Image *) NULL);
4711  assert(image->signature == MagickCoreSignature);
4712  assert(exception != (ExceptionInfo *) NULL);
4713  assert(exception->signature == MagickCoreSignature);
4714  assert(offset != (RectangleInfo *) NULL);
4715  if (IsEventLogging() != MagickFalse)
4716  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4717  SetGeometry(reconstruct,offset);
4718  *similarity_metric=0.0;
4719  offset->x=0;
4720  offset->y=0;
4721 #if defined(MAGICKCORE_HDRI_SUPPORT) && defined(MAGICKCORE_FFTW_DELEGATE)
4722 {
4723  const char *artifact = GetImageArtifact(image,"compare:frequency-domain");
4724  if (artifact == (const char *) NULL)
4725  artifact=GetImageArtifact(image,"compare:accelerate-ncc");
4726  if (((artifact == (const char *) NULL) ||
4727  (IsStringTrue(artifact) != MagickFalse)) &&
4728  ((image->channels & ReadMaskChannel) == 0))
4729  switch (metric)
4730  {
4731  case DotProductCorrelationErrorMetric:
4732  {
4733  similarity_image=DPCSimilarityImage(image,reconstruct,offset,
4734  similarity_metric,exception);
4735  return(similarity_image);
4736  }
4737  case MeanSquaredErrorMetric:
4738  {
4739  similarity_image=MSESimilarityImage(image,reconstruct,offset,
4740  similarity_metric,exception);
4741  return(similarity_image);
4742  }
4743  case NormalizedCrossCorrelationErrorMetric:
4744  {
4745  similarity_image=NCCSimilarityImage(image,reconstruct,offset,
4746  similarity_metric,exception);
4747  return(similarity_image);
4748  }
4749  case PeakSignalToNoiseRatioErrorMetric:
4750  {
4751  similarity_image=PSNRSimilarityImage(image,reconstruct,offset,
4752  similarity_metric,exception);
4753  return(similarity_image);
4754  }
4755  case PhaseCorrelationErrorMetric:
4756  {
4757  similarity_image=PhaseSimilarityImage(image,reconstruct,offset,
4758  similarity_metric,exception);
4759  return(similarity_image);
4760  }
4761  case RootMeanSquaredErrorMetric:
4762  case UndefinedErrorMetric:
4763  {
4764  similarity_image=RMSESimilarityImage(image,reconstruct,offset,
4765  similarity_metric,exception);
4766  return(similarity_image);
4767  }
4768  default:
4769  break;
4770  }
4771 }
4772 #endif
4773  if ((image->columns < reconstruct->columns) ||
4774  (image->rows < reconstruct->rows))
4775  {
4776  (void) ThrowMagickException(exception,GetMagickModule(),OptionWarning,
4777  "GeometryDoesNotContainImage","`%s'",image->filename);
4778  return((Image *) NULL);
4779  }
4780  SetImageCompareBounds(image,reconstruct,&columns,&rows);
4781  similarity_image=CloneImage(image,columns,rows,MagickTrue,exception);
4782  if (similarity_image == (Image *) NULL)
4783  return((Image *) NULL);
4784  similarity_image->depth=32;
4785  similarity_image->colorspace=GRAYColorspace;
4786  similarity_image->alpha_trait=UndefinedPixelTrait;
4787  status=SetImageStorageClass(similarity_image,DirectClass,exception);
4788  if (status == MagickFalse)
4789  return(DestroyImage(similarity_image));
4790  /*
4791  Measure similarity of reconstruction image against image.
4792  */
4793  similarity_info.similarity=GetSimilarityMetric(image,reconstruct,metric,
4794  similarity_info.x,similarity_info.y,exception);
4795  similarity_view=AcquireAuthenticCacheView(similarity_image,exception);
4796 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4797  #pragma omp parallel for schedule(static) shared(similarity_info,status) \
4798  magick_number_threads(image,reconstruct,similarity_image->rows,1)
4799 #endif
4800  for (y=0; y < (ssize_t) similarity_image->rows; y++)
4801  {
4802  double
4803  similarity;
4804 
4805  MagickBooleanType
4806  threshold_trigger = MagickFalse;
4807 
4808  Quantum
4809  *magick_restrict q;
4810 
4811  SimilarityInfo
4812  channel_info = similarity_info;
4813 
4814  ssize_t
4815  x;
4816 
4817  if (status == MagickFalse)
4818  continue;
4819  if (threshold_trigger != MagickFalse)
4820  continue;
4821  q=QueueCacheViewAuthenticPixels(similarity_view,0,y,
4822  similarity_image->columns,1,exception);
4823  if (q == (Quantum *) NULL)
4824  {
4825  status=MagickFalse;
4826  continue;
4827  }
4828  for (x=0; x < (ssize_t) similarity_image->columns; x++)
4829  {
4830  ssize_t
4831  i;
4832 
4833  similarity=GetSimilarityMetric((Image *) image,reconstruct,metric,x,y,
4834  exception);
4835  switch (metric)
4836  {
4837  case DotProductCorrelationErrorMetric:
4838  case NormalizedCrossCorrelationErrorMetric:
4839  case PeakSignalToNoiseRatioErrorMetric:
4840  case PhaseCorrelationErrorMetric:
4841  case StructuralSimilarityErrorMetric:
4842  {
4843  if (similarity <= channel_info.similarity)
4844  break;
4845  channel_info.similarity=similarity;
4846  channel_info.x=x;
4847  channel_info.y=y;
4848  break;
4849  }
4850  default:
4851  {
4852  if (similarity >= channel_info.similarity)
4853  break;
4854  channel_info.similarity=similarity;
4855  channel_info.x=x;
4856  channel_info.y=y;
4857  break;
4858  }
4859  }
4860  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
4861  {
4862  PixelChannel channel = GetPixelChannelChannel(image,i);
4863  PixelTrait traits = GetPixelChannelTraits(image,channel);
4864  PixelTrait similarity_traits = GetPixelChannelTraits(similarity_image,
4865  channel);
4866  if (((traits & UpdatePixelTrait) == 0) ||
4867  ((similarity_traits & UpdatePixelTrait) == 0))
4868  continue;
4869  switch (metric)
4870  {
4871  case DotProductCorrelationErrorMetric:
4872  case NormalizedCrossCorrelationErrorMetric:
4873  case PeakSignalToNoiseRatioErrorMetric:
4874  case PhaseCorrelationErrorMetric:
4875  case StructuralSimilarityErrorMetric:
4876  {
4877  SetPixelChannel(similarity_image,channel,ClampToQuantum((double)
4878  QuantumRange*similarity),q);
4879  break;
4880  }
4881  default:
4882  {
4883  SetPixelChannel(similarity_image,channel,ClampToQuantum((double)
4884  QuantumRange*(1.0-similarity)),q);
4885  break;
4886  }
4887  }
4888  }
4889  q+=(ptrdiff_t) GetPixelChannels(similarity_image);
4890  }
4891 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4892  #pragma omp critical (MagickCore_GetSimilarityMetric)
4893 #endif
4894  switch (metric)
4895  {
4896  case DotProductCorrelationErrorMetric:
4897  case NormalizedCrossCorrelationErrorMetric:
4898  case PeakSignalToNoiseRatioErrorMetric:
4899  case PhaseCorrelationErrorMetric:
4900  case StructuralSimilarityErrorMetric:
4901  {
4902  if (similarity_threshold != DefaultSimilarityThreshold)
4903  if (channel_info.similarity >= similarity_threshold)
4904  threshold_trigger=MagickTrue;
4905  if (channel_info.similarity >= similarity_info.similarity)
4906  similarity_info=channel_info;
4907  break;
4908  }
4909  default:
4910  {
4911  if (similarity_threshold != DefaultSimilarityThreshold)
4912  if (channel_info.similarity < similarity_threshold)
4913  threshold_trigger=MagickTrue;
4914  if (channel_info.similarity < similarity_info.similarity)
4915  similarity_info=channel_info;
4916  break;
4917  }
4918  }
4919  if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
4920  status=MagickFalse;
4921  if (image->progress_monitor != (MagickProgressMonitor) NULL)
4922  {
4923  MagickBooleanType
4924  proceed;
4925 
4926  progress++;
4927  proceed=SetImageProgress(image,SimilarityImageTag,progress,image->rows);
4928  if (proceed == MagickFalse)
4929  status=MagickFalse;
4930  }
4931  }
4932  similarity_view=DestroyCacheView(similarity_view);
4933  if (status == MagickFalse)
4934  similarity_image=DestroyImage(similarity_image);
4935  *similarity_metric=similarity_info.similarity;
4936  if (fabs(*similarity_metric) < MagickEpsilon)
4937  *similarity_metric=0.0;
4938  offset->x=similarity_info.x;
4939  offset->y=similarity_info.y;
4940  (void) FormatImageProperty((Image *) image,"similarity","%.*g",
4941  GetMagickPrecision(),*similarity_metric);
4942  (void) FormatImageProperty((Image *) image,"similarity.offset.x","%.*g",
4943  GetMagickPrecision(),(double) offset->x);
4944  (void) FormatImageProperty((Image *) image,"similarity.offset.y","%.*g",
4945  GetMagickPrecision(),(double) offset->y);
4946  return(similarity_image);
4947 }
_CacheView
Definition: cache-view.c:66
_ChannelPerceptualHash
Definition: statistic.h:74
_ChannelStatistics
Definition: statistic.h:30
_ExceptionInfo
Definition: exception.h:102
_Image
Definition: image.h:132
_KernelInfo
Definition: morphology.h:103
_PixelInfo
Definition: pixel.h:182
_RectangleInfo
Definition: geometry.h:130