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