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