MagickCore  6.9.13-44
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 dedicated %
21 % to making software imaging solutions freely available. %
22 % %
23 % You may not use this file except in compliance with the License. You may %
24 % obtain a copy of the License at %
25 % %
26 % https://imagemagick.org/license/ %
27 % %
28 % Unless required by applicable law or agreed to in writing, software %
29 % distributed under the License is distributed on an "AS IS" BASIS, %
30 % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. %
31 % See the License for the specific language governing permissions and %
32 % limitations under the License. %
33 % %
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35 %
36 %
37 %
38 */
39 
40 /*
41  Include declarations.
42 */
43 #include "magick/studio.h"
44 #include "magick/artifact.h"
45 #include "magick/attribute.h"
46 #include "magick/cache-view.h"
47 #include "magick/channel.h"
48 #include "magick/client.h"
49 #include "magick/color.h"
50 #include "magick/color-private.h"
51 #include "magick/colorspace.h"
52 #include "magick/colorspace-private.h"
53 #include "magick/compare.h"
54 #include "magick/compare-private.h"
55 #include "magick/composite-private.h"
56 #include "magick/constitute.h"
57 #include "magick/exception-private.h"
58 #include "magick/geometry.h"
59 #include "magick/image-private.h"
60 #include "magick/list.h"
61 #include "magick/log.h"
62 #include "magick/memory_.h"
63 #include "magick/monitor.h"
64 #include "magick/monitor-private.h"
65 #include "magick/option.h"
66 #include "magick/pixel-private.h"
67 #include "magick/property.h"
68 #include "magick/resource_.h"
69 #include "magick/statistic-private.h"
70 #include "magick/string_.h"
71 #include "magick/string-private.h"
72 #include "magick/statistic.h"
73 #include "magick/thread-private.h"
74 #include "magick/transform.h"
75 #include "magick/utility.h"
76 #include "magick/version.h"
77 
78 /*
79 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
80 % %
81 % %
82 % %
83 % C o m p a r e I m a g e C h a n n e l s %
84 % %
85 % %
86 % %
87 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
88 %
89 % CompareImageChannels() compares one or more image channels of an image
90 % to a reconstructed image and returns the difference image.
91 %
92 % The format of the CompareImageChannels method is:
93 %
94 % Image *CompareImageChannels(const Image *image,
95 % const Image *reconstruct_image,const ChannelType channel,
96 % const MetricType metric,double *distortion,ExceptionInfo *exception)
97 %
98 % A description of each parameter follows:
99 %
100 % o image: the image.
101 %
102 % o reconstruct_image: the reconstruct image.
103 %
104 % o channel: the channel.
105 %
106 % o metric: the metric.
107 %
108 % o distortion: the computed distortion between the images.
109 %
110 % o exception: return any errors or warnings in this structure.
111 %
112 */
113 
114 MagickExport Image *CompareImages(Image *image,const Image *reconstruct_image,
115  const MetricType metric,double *distortion,ExceptionInfo *exception)
116 {
117  Image
118  *highlight_image;
119 
120  highlight_image=CompareImageChannels(image,reconstruct_image,
121  CompositeChannels,metric,distortion,exception);
122  return(highlight_image);
123 }
124 
125 static size_t GetNumberChannels(const Image *image,const ChannelType channel)
126 {
127  size_t
128  channels;
129 
130  channels=0;
131  if ((channel & RedChannel) != 0)
132  channels++;
133  if ((channel & GreenChannel) != 0)
134  channels++;
135  if ((channel & BlueChannel) != 0)
136  channels++;
137  if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
138  channels++;
139  if (((channel & IndexChannel) != 0) && (image->colorspace == CMYKColorspace))
140  channels++;
141  return(channels == 0 ? 1UL : channels);
142 }
143 
144 static inline MagickBooleanType ValidateImageMorphology(
145  const Image *magick_restrict image,
146  const Image *magick_restrict reconstruct_image)
147 {
148  /*
149  Does the image match the reconstructed image morphology?
150  */
151  if (GetNumberChannels(image,DefaultChannels) !=
152  GetNumberChannels(reconstruct_image,DefaultChannels))
153  return(MagickFalse);
154  return(MagickTrue);
155 }
156 
157 MagickExport Image *CompareImageChannels(Image *image,
158  const Image *reconstruct_image,const ChannelType channel,
159  const MetricType metric,double *distortion,ExceptionInfo *exception)
160 {
161  CacheView
162  *highlight_view,
163  *image_view,
164  *reconstruct_view;
165 
166  const char
167  *artifact;
168 
169  Image
170  *clone_image,
171  *difference_image,
172  *highlight_image;
173 
174  MagickBooleanType
175  status = MagickTrue;
176 
177  MagickPixelPacket
178  highlight,
179  lowlight,
180  zero;
181 
182  size_t
183  columns,
184  rows;
185 
186  ssize_t
187  y;
188 
189  assert(image != (Image *) NULL);
190  assert(image->signature == MagickCoreSignature);
191  assert(reconstruct_image != (const Image *) NULL);
192  assert(reconstruct_image->signature == MagickCoreSignature);
193  assert(distortion != (double *) NULL);
194  if (IsEventLogging() != MagickFalse)
195  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
196  *distortion=0.0;
197  if (metric != PerceptualHashErrorMetric)
198  if (ValidateImageMorphology(image,reconstruct_image) == MagickFalse)
199  ThrowImageException(ImageError,"ImageMorphologyDiffers");
200  status=GetImageChannelDistortion(image,reconstruct_image,channel,metric,
201  distortion,exception);
202  if (status == MagickFalse)
203  return((Image *) NULL);
204  clone_image=CloneImage(image,0,0,MagickTrue,exception);
205  if (clone_image == (Image *) NULL)
206  return((Image *) NULL);
207  (void) SetImageMask(clone_image,(Image *) NULL);
208  difference_image=CloneImage(clone_image,0,0,MagickTrue,exception);
209  clone_image=DestroyImage(clone_image);
210  if (difference_image == (Image *) NULL)
211  return((Image *) NULL);
212  (void) SetImageAlphaChannel(difference_image,OpaqueAlphaChannel);
213  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
214  highlight_image=CloneImage(image,columns,rows,MagickTrue,exception);
215  if (highlight_image == (Image *) NULL)
216  {
217  difference_image=DestroyImage(difference_image);
218  return((Image *) NULL);
219  }
220  if (SetImageStorageClass(highlight_image,DirectClass) == MagickFalse)
221  {
222  InheritException(exception,&highlight_image->exception);
223  difference_image=DestroyImage(difference_image);
224  highlight_image=DestroyImage(highlight_image);
225  return((Image *) NULL);
226  }
227  (void) SetImageMask(highlight_image,(Image *) NULL);
228  (void) SetImageAlphaChannel(highlight_image,OpaqueAlphaChannel);
229  (void) QueryMagickColor("#f1001ecc",&highlight,exception);
230  artifact=GetImageArtifact(image,"compare:highlight-color");
231  if (artifact != (const char *) NULL)
232  (void) QueryMagickColor(artifact,&highlight,exception);
233  (void) QueryMagickColor("#ffffffcc",&lowlight,exception);
234  artifact=GetImageArtifact(image,"compare:lowlight-color");
235  if (artifact != (const char *) NULL)
236  (void) QueryMagickColor(artifact,&lowlight,exception);
237  if (highlight_image->colorspace == CMYKColorspace)
238  {
239  ConvertRGBToCMYK(&highlight);
240  ConvertRGBToCMYK(&lowlight);
241  }
242  /*
243  Generate difference image.
244  */
245  GetMagickPixelPacket(image,&zero);
246  image_view=AcquireVirtualCacheView(image,exception);
247  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
248  highlight_view=AcquireAuthenticCacheView(highlight_image,exception);
249 #if defined(MAGICKCORE_OPENMP_SUPPORT)
250  #pragma omp parallel for schedule(static) shared(status) \
251  magick_number_threads(image,highlight_image,rows,1)
252 #endif
253  for (y=0; y < (ssize_t) rows; y++)
254  {
255  MagickBooleanType
256  sync;
257 
258  MagickPixelPacket
259  pixel,
260  reconstruct_pixel;
261 
262  const IndexPacket
263  *magick_restrict indexes,
264  *magick_restrict reconstruct_indexes;
265 
266  const PixelPacket
267  *magick_restrict p,
268  *magick_restrict q;
269 
270  IndexPacket
271  *magick_restrict highlight_indexes;
272 
273  PixelPacket
274  *magick_restrict r;
275 
276  ssize_t
277  x;
278 
279  if (status == MagickFalse)
280  continue;
281  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
282  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
283  r=QueueCacheViewAuthenticPixels(highlight_view,0,y,columns,1,exception);
284  if ((p == (const PixelPacket *) NULL) ||
285  (q == (const PixelPacket *) NULL) || (r == (PixelPacket *) NULL))
286  {
287  status=MagickFalse;
288  continue;
289  }
290  indexes=GetCacheViewVirtualIndexQueue(image_view);
291  reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
292  highlight_indexes=GetCacheViewAuthenticIndexQueue(highlight_view);
293  pixel=zero;
294  reconstruct_pixel=zero;
295  for (x=0; x < (ssize_t) columns; x++)
296  {
297  SetMagickPixelPacket(image,p,indexes == (IndexPacket *) NULL ? NULL :
298  indexes+x,&pixel);
299  SetMagickPixelPacket(reconstruct_image,q,reconstruct_indexes ==
300  (IndexPacket *) NULL ? NULL : reconstruct_indexes+x,&reconstruct_pixel);
301  if (IsMagickColorSimilar(&pixel,&reconstruct_pixel) == MagickFalse)
302  SetPixelPacket(highlight_image,&highlight,r,highlight_indexes ==
303  (IndexPacket *) NULL ? NULL : highlight_indexes+x);
304  else
305  SetPixelPacket(highlight_image,&lowlight,r,highlight_indexes ==
306  (IndexPacket *) NULL ? NULL : highlight_indexes+x);
307  p++;
308  q++;
309  r++;
310  }
311  sync=SyncCacheViewAuthenticPixels(highlight_view,exception);
312  if (sync == MagickFalse)
313  status=MagickFalse;
314  }
315  highlight_view=DestroyCacheView(highlight_view);
316  reconstruct_view=DestroyCacheView(reconstruct_view);
317  image_view=DestroyCacheView(image_view);
318  (void) CompositeImage(difference_image,image->compose,highlight_image,0,0);
319  highlight_image=DestroyImage(highlight_image);
320  if (status == MagickFalse)
321  difference_image=DestroyImage(difference_image);
322  return(difference_image);
323 }
324 
325 /*
326 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
327 % %
328 % %
329 % %
330 % G e t I m a g e C h a n n e l D i s t o r t i o n %
331 % %
332 % %
333 % %
334 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
335 %
336 % GetImageChannelDistortion() compares one or more image channels of an image
337 % to a reconstructed image and returns the specified distortion metric.
338 %
339 % The format of the GetImageChannelDistortion method is:
340 %
341 % MagickBooleanType GetImageChannelDistortion(const Image *image,
342 % const Image *reconstruct_image,const ChannelType channel,
343 % const MetricType metric,double *distortion,ExceptionInfo *exception)
344 %
345 % A description of each parameter follows:
346 %
347 % o image: the image.
348 %
349 % o reconstruct_image: the reconstruct image.
350 %
351 % o channel: the channel.
352 %
353 % o metric: the metric.
354 %
355 % o distortion: the computed distortion between the images.
356 %
357 % o exception: return any errors or warnings in this structure.
358 %
359 */
360 
361 MagickExport MagickBooleanType GetImageDistortion(Image *image,
362  const Image *reconstruct_image,const MetricType metric,double *distortion,
363  ExceptionInfo *exception)
364 {
365  MagickBooleanType
366  status;
367 
368  status=GetImageChannelDistortion(image,reconstruct_image,CompositeChannels,
369  metric,distortion,exception);
370  return(status);
371 }
372 
373 static MagickBooleanType GetAESimilarity(const Image *image,
374  const Image *reconstruct_image,const ChannelType channel,double *similarity,
375  ExceptionInfo *exception)
376 {
377  CacheView
378  *image_view,
379  *reconstruct_view;
380 
381  double
382  area,
383  fuzz;
384 
385  MagickBooleanType
386  status = MagickTrue;
387 
388  size_t
389  columns,
390  rows;
391 
392  ssize_t
393  j,
394  y;
395 
396  /*
397  Compute the absolute difference in pixels between two images.
398  */
399  fuzz=GetFuzzyColorDistance(image,reconstruct_image);
400  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
401  image_view=AcquireVirtualCacheView(image,exception);
402  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
403 #if defined(MAGICKCORE_OPENMP_SUPPORT)
404  #pragma omp parallel for schedule(static) shared(similarity,status) \
405  magick_number_threads(image,image,rows,1)
406 #endif
407  for (y=0; y < (ssize_t) rows; y++)
408  {
409  const IndexPacket
410  *magick_restrict indexes,
411  *magick_restrict reconstruct_indexes;
412 
413  const PixelPacket
414  *magick_restrict p,
415  *magick_restrict q;
416 
417  double
418  channel_similarity[CompositeChannels+1] = { 0.0 };
419 
420  ssize_t
421  i,
422  x;
423 
424  if (status == MagickFalse)
425  continue;
426  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
427  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
428  if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
429  {
430  status=MagickFalse;
431  continue;
432  }
433  indexes=GetCacheViewVirtualIndexQueue(image_view);
434  reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
435  (void) memset(channel_similarity,0,sizeof(channel_similarity));
436  for (x=0; x < (ssize_t) columns; x++)
437  {
438  double
439  Da,
440  error,
441  Sa;
442 
443  size_t
444  count = 0;
445 
446  Sa=QuantumScale*(image->matte != MagickFalse ? (double) GetPixelAlpha(p) :
447  ((double) QuantumRange-(double) OpaqueOpacity));
448  Da=QuantumScale*(image->matte != MagickFalse ? (double) GetPixelAlpha(q) :
449  ((double) QuantumRange-(double) OpaqueOpacity));
450  if ((channel & RedChannel) != 0)
451  {
452  error=Sa*(double) GetPixelRed(p)-Da*(double)
453  GetPixelRed(q);
454  if (MagickSafeSignificantError(error*error,fuzz) != MagickFalse)
455  {
456  channel_similarity[RedChannel]++;
457  count++;
458  }
459  }
460  if ((channel & GreenChannel) != 0)
461  {
462  error=Sa*(double) GetPixelGreen(p)-Da*(double)
463  GetPixelGreen(q);
464  if (MagickSafeSignificantError(error*error,fuzz) != MagickFalse)
465  {
466  channel_similarity[GreenChannel]++;
467  count++;
468  }
469  }
470  if ((channel & BlueChannel) != 0)
471  {
472  error=Sa*(double) GetPixelBlue(p)-Da*(double)
473  GetPixelBlue(q);
474  if (MagickSafeSignificantError(error*error,fuzz) != MagickFalse)
475  {
476  channel_similarity[BlueChannel]++;
477  count++;
478  }
479  }
480  if (((channel & OpacityChannel) != 0) &&
481  (image->matte != MagickFalse))
482  {
483  error=(double) GetPixelOpacity(p)-(double) GetPixelOpacity(q);
484  if (MagickSafeSignificantError(error*error,fuzz) != MagickFalse)
485  {
486  channel_similarity[OpacityChannel]++;
487  count++;
488  }
489  }
490  if (((channel & IndexChannel) != 0) &&
491  (image->colorspace == CMYKColorspace))
492  {
493  error=Sa*(double) indexes[x]-Da*(double) reconstruct_indexes[x];
494  if (MagickSafeSignificantError(error*error,fuzz) != MagickFalse)
495  {
496  channel_similarity[IndexChannel]++;
497  count++;
498  }
499  }
500  if (count != 0)
501  channel_similarity[CompositeChannels]++;
502  p++;
503  q++;
504  }
505 #if defined(MAGICKCORE_OPENMP_SUPPORT)
506  #pragma omp critical (MagickCore_GetAESimilarity)
507 #endif
508  for (i=0; i <= (ssize_t) CompositeChannels; i++)
509  similarity[i]+=channel_similarity[i];
510  }
511  reconstruct_view=DestroyCacheView(reconstruct_view);
512  image_view=DestroyCacheView(image_view);
513  area=MagickSafeReciprocal((double) columns*rows);
514  for (j=0; j <= CompositeChannels; j++)
515  similarity[j]*=area;
516  return(status);
517 }
518 
519 static MagickBooleanType GetFUZZSimilarity(const Image *image,
520  const Image *reconstruct_image,const ChannelType channel,
521  double *similarity,ExceptionInfo *exception)
522 {
523  CacheView
524  *image_view,
525  *reconstruct_view;
526 
527  double
528  area = 0.0,
529  fuzz;
530 
531  MagickBooleanType
532  status = MagickTrue;
533 
534  size_t
535  columns,
536  rows;
537 
538  ssize_t
539  i,
540  y;
541 
542  fuzz=GetFuzzyColorDistance(image,reconstruct_image);
543  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
544  image_view=AcquireVirtualCacheView(image,exception);
545  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
546 #if defined(MAGICKCORE_OPENMP_SUPPORT)
547  #pragma omp parallel for schedule(static) shared(status) \
548  magick_number_threads(image,image,rows,1)
549 #endif
550  for (y=0; y < (ssize_t) rows; y++)
551  {
552  double
553  channel_area = 0.0,
554  channel_similarity[CompositeChannels+1] = { 0.0 };
555 
556  const IndexPacket
557  *magick_restrict indexes,
558  *magick_restrict reconstruct_indexes;
559 
560  const PixelPacket
561  *magick_restrict p,
562  *magick_restrict q;
563 
564  ssize_t
565  i,
566  x;
567 
568  if (status == MagickFalse)
569  continue;
570  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
571  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
572  if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
573  {
574  status=MagickFalse;
575  continue;
576  }
577  indexes=GetCacheViewVirtualIndexQueue(image_view);
578  reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
579  for (x=0; x < (ssize_t) columns; x++)
580  {
581  MagickRealType
582  Da,
583  error,
584  Sa;
585 
586  Sa=QuantumScale*(image->matte != MagickFalse ? (double)
587  GetPixelAlpha(p) : ((double) QuantumRange-(double) OpaqueOpacity));
588  Da=QuantumScale*(reconstruct_image->matte != MagickFalse ?
589  (double) GetPixelAlpha(q) : ((double) QuantumRange-(double)
590  OpaqueOpacity));
591  if ((channel & RedChannel) != 0)
592  {
593  error=QuantumScale*(Sa*GetPixelRed(p)-Da*GetPixelRed(q));
594  if (MagickSafeSignificantError(error*error,fuzz) != MagickFalse)
595  {
596  channel_similarity[RedChannel]+=error*error;
597  channel_similarity[CompositeChannels]+=error*error;
598  channel_area++;
599  }
600  }
601  if ((channel & GreenChannel) != 0)
602  {
603  error=QuantumScale*(Sa*GetPixelGreen(p)-Da*GetPixelGreen(q));
604  if (MagickSafeSignificantError(error*error,fuzz) != MagickFalse)
605  {
606  channel_similarity[GreenChannel]+=error*error;
607  channel_similarity[CompositeChannels]+=error*error;
608  channel_area++;
609  }
610  }
611  if ((channel & BlueChannel) != 0)
612  {
613  error=QuantumScale*(Sa*GetPixelBlue(p)-Da*GetPixelBlue(q));
614  if (MagickSafeSignificantError(error*error,fuzz) != MagickFalse)
615  {
616  channel_similarity[BlueChannel]+=error*error;
617  channel_similarity[CompositeChannels]+=error*error;
618  channel_area++;
619  }
620  }
621  if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
622  {
623  error=QuantumScale*((double) GetPixelOpacity(p)-GetPixelOpacity(q));
624  if (MagickSafeSignificantError(error*error,fuzz) != MagickFalse)
625  {
626  channel_similarity[OpacityChannel]+=error*error;
627  channel_similarity[CompositeChannels]+=error*error;
628  channel_area++;
629  }
630  }
631  if (((channel & IndexChannel) != 0) &&
632  (image->colorspace == CMYKColorspace))
633  {
634  error=QuantumScale*(Sa*GetPixelIndex(indexes+x)-Da*
635  GetPixelIndex(reconstruct_indexes+x));
636  if (MagickSafeSignificantError(error*error,fuzz) != MagickFalse)
637  {
638  channel_similarity[BlackChannel]+=error*error;
639  channel_similarity[CompositeChannels]+=error*error;
640  channel_area++;
641  }
642  }
643  p++;
644  q++;
645  }
646 #if defined(MAGICKCORE_OPENMP_SUPPORT)
647  #pragma omp critical (MagickCore_GetMeanAbsoluteError)
648 #endif
649  {
650  area+=channel_area;
651  for (i=0; i <= (ssize_t) CompositeChannels; i++)
652  similarity[i]+=channel_similarity[i];
653  }
654  }
655  reconstruct_view=DestroyCacheView(reconstruct_view);
656  image_view=DestroyCacheView(image_view);
657  area=MagickSafeReciprocal(area);
658  for (i=0; i <= (ssize_t) CompositeChannels; i++)
659  similarity[i]*=area;
660  return(status);
661 }
662 
663 static MagickBooleanType GetMAESimilarity(const Image *image,
664  const Image *reconstruct_image,const ChannelType channel,
665  double *similarity,ExceptionInfo *exception)
666 {
667  CacheView
668  *image_view,
669  *reconstruct_view;
670 
671  MagickBooleanType
672  status;
673 
674  size_t
675  columns,
676  rows;
677 
678  ssize_t
679  i,
680  y;
681 
682  status=MagickTrue;
683  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
684  image_view=AcquireVirtualCacheView(image,exception);
685  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
686 #if defined(MAGICKCORE_OPENMP_SUPPORT)
687  #pragma omp parallel for schedule(static) shared(status) \
688  magick_number_threads(image,image,rows,1)
689 #endif
690  for (y=0; y < (ssize_t) rows; y++)
691  {
692  double
693  channel_similarity[CompositeChannels+1];
694 
695  const IndexPacket
696  *magick_restrict indexes,
697  *magick_restrict reconstruct_indexes;
698 
699  const PixelPacket
700  *magick_restrict p,
701  *magick_restrict q;
702 
703  ssize_t
704  i,
705  x;
706 
707  if (status == MagickFalse)
708  continue;
709  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
710  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
711  if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
712  {
713  status=MagickFalse;
714  continue;
715  }
716  indexes=GetCacheViewVirtualIndexQueue(image_view);
717  reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
718  (void) memset(channel_similarity,0,sizeof(channel_similarity));
719  for (x=0; x < (ssize_t) columns; x++)
720  {
721  MagickRealType
722  distance,
723  Da,
724  Sa;
725 
726  Sa=QuantumScale*(image->matte != MagickFalse ? (double) GetPixelAlpha(p) :
727  ((double) QuantumRange-(double) OpaqueOpacity));
728  Da=QuantumScale*(reconstruct_image->matte != MagickFalse ?
729  (double) GetPixelAlpha(q) : ((double) QuantumRange-(double)
730  OpaqueOpacity));
731  if ((channel & RedChannel) != 0)
732  {
733  distance=QuantumScale*fabs(Sa*(double) GetPixelRed(p)-Da*
734  (double) GetPixelRed(q));
735  channel_similarity[RedChannel]+=distance;
736  channel_similarity[CompositeChannels]+=distance;
737  }
738  if ((channel & GreenChannel) != 0)
739  {
740  distance=QuantumScale*fabs(Sa*(double) GetPixelGreen(p)-Da*
741  (double) GetPixelGreen(q));
742  channel_similarity[GreenChannel]+=distance;
743  channel_similarity[CompositeChannels]+=distance;
744  }
745  if ((channel & BlueChannel) != 0)
746  {
747  distance=QuantumScale*fabs(Sa*(double) GetPixelBlue(p)-Da*
748  (double) GetPixelBlue(q));
749  channel_similarity[BlueChannel]+=distance;
750  channel_similarity[CompositeChannels]+=distance;
751  }
752  if (((channel & OpacityChannel) != 0) &&
753  (image->matte != MagickFalse))
754  {
755  distance=QuantumScale*fabs((double) GetPixelOpacity(p)-(double)
756  GetPixelOpacity(q));
757  channel_similarity[OpacityChannel]+=distance;
758  channel_similarity[CompositeChannels]+=distance;
759  }
760  if (((channel & IndexChannel) != 0) &&
761  (image->colorspace == CMYKColorspace))
762  {
763  distance=QuantumScale*fabs(Sa*(double) GetPixelIndex(indexes+x)-Da*
764  (double) GetPixelIndex(reconstruct_indexes+x));
765  channel_similarity[BlackChannel]+=distance;
766  channel_similarity[CompositeChannels]+=distance;
767  }
768  p++;
769  q++;
770  }
771 #if defined(MAGICKCORE_OPENMP_SUPPORT)
772  #pragma omp critical (MagickCore_GetMeanAbsoluteError)
773 #endif
774  for (i=0; i <= (ssize_t) CompositeChannels; i++)
775  similarity[i]+=channel_similarity[i];
776  }
777  reconstruct_view=DestroyCacheView(reconstruct_view);
778  image_view=DestroyCacheView(image_view);
779  for (i=0; i <= (ssize_t) CompositeChannels; i++)
780  similarity[i]/=((double) columns*rows);
781  similarity[CompositeChannels]/=(double) GetNumberChannels(image,channel);
782  return(status);
783 }
784 
785 static MagickBooleanType GetMEPPSimilarity(Image *image,
786  const Image *reconstruct_image,const ChannelType channel,double *similarity,
787  ExceptionInfo *exception)
788 {
789  CacheView
790  *image_view,
791  *reconstruct_view;
792 
793  double
794  maximum_error = -MagickMaximumValue,
795  mean_error = 0.0;
796 
797  MagickBooleanType
798  status;
799 
800  size_t
801  columns,
802  rows;
803 
804  ssize_t
805  i,
806  y;
807 
808  status=MagickTrue;
809  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
810  image_view=AcquireVirtualCacheView(image,exception);
811  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
812 #if defined(MAGICKCORE_OPENMP_SUPPORT)
813  #pragma omp parallel for schedule(static) shared(maximum_error,status) \
814  magick_number_threads(image,image,rows,1)
815 #endif
816  for (y=0; y < (ssize_t) rows; y++)
817  {
818  double
819  channel_similarity[CompositeChannels+1] = { 0.0 },
820  local_maximum = maximum_error,
821  local_mean_error = 0.0;
822 
823  const IndexPacket
824  *magick_restrict indexes,
825  *magick_restrict reconstruct_indexes;
826 
827  const PixelPacket
828  *magick_restrict p,
829  *magick_restrict q;
830 
831  ssize_t
832  i,
833  x;
834 
835  if (status == MagickFalse)
836  continue;
837  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
838  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
839  if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
840  {
841  status=MagickFalse;
842  continue;
843  }
844  indexes=GetCacheViewVirtualIndexQueue(image_view);
845  reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
846  (void) memset(channel_similarity,0,sizeof(channel_similarity));
847  for (x=0; x < (ssize_t) columns; x++)
848  {
849  MagickRealType
850  distance,
851  Da,
852  Sa;
853 
854  Sa=QuantumScale*(image->matte != MagickFalse ? (double) GetPixelAlpha(p) :
855  ((double) QuantumRange-(double) OpaqueOpacity));
856  Da=QuantumScale*(reconstruct_image->matte != MagickFalse ?
857  (double) GetPixelAlpha(q) : ((double) QuantumRange-(double)
858  OpaqueOpacity));
859  if ((channel & RedChannel) != 0)
860  {
861  distance=QuantumScale*fabs(Sa*(double) GetPixelRed(p)-Da*
862  (double) GetPixelRed(q));
863  channel_similarity[RedChannel]+=distance;
864  channel_similarity[CompositeChannels]+=distance;
865  local_mean_error+=distance*distance;
866  if (distance > local_maximum)
867  local_maximum=distance;
868  }
869  if ((channel & GreenChannel) != 0)
870  {
871  distance=QuantumScale*fabs(Sa*(double) GetPixelGreen(p)-Da*
872  (double) GetPixelGreen(q));
873  channel_similarity[GreenChannel]+=distance;
874  channel_similarity[CompositeChannels]+=distance;
875  local_mean_error+=distance*distance;
876  if (distance > local_maximum)
877  local_maximum=distance;
878  }
879  if ((channel & BlueChannel) != 0)
880  {
881  distance=QuantumScale*fabs(Sa*(double) GetPixelBlue(p)-Da*
882  (double) GetPixelBlue(q));
883  channel_similarity[BlueChannel]+=distance;
884  channel_similarity[CompositeChannels]+=distance;
885  local_mean_error+=distance*distance;
886  if (distance > local_maximum)
887  local_maximum=distance;
888  }
889  if (((channel & OpacityChannel) != 0) &&
890  (image->matte != MagickFalse))
891  {
892  distance=QuantumScale*fabs((double) GetPixelOpacity(p)-(double)
893  GetPixelOpacity(q));
894  channel_similarity[OpacityChannel]+=distance;
895  channel_similarity[CompositeChannels]+=distance;
896  local_mean_error+=distance*distance;
897  if (distance > local_maximum)
898  local_maximum=distance;
899  }
900  if (((channel & IndexChannel) != 0) &&
901  (image->colorspace == CMYKColorspace))
902  {
903  distance=QuantumScale*fabs(Sa*(double) GetPixelIndex(indexes+x)-Da*
904  (double) GetPixelIndex(reconstruct_indexes+x));
905  channel_similarity[BlackChannel]+=distance;
906  channel_similarity[CompositeChannels]+=distance;
907  local_mean_error+=distance*distance;
908  if (distance > local_maximum)
909  local_maximum=distance;
910  }
911  p++;
912  q++;
913  }
914 #if defined(MAGICKCORE_OPENMP_SUPPORT)
915  #pragma omp critical (MagickCore_GetMeanAbsoluteError)
916 #endif
917  {
918  for (i=0; i <= (ssize_t) CompositeChannels; i++)
919  similarity[i]+=channel_similarity[i];
920  mean_error+=local_mean_error;
921  if (local_maximum > maximum_error)
922  maximum_error=local_maximum;
923  }
924  }
925  reconstruct_view=DestroyCacheView(reconstruct_view);
926  image_view=DestroyCacheView(image_view);
927  for (i=0; i <= (ssize_t) CompositeChannels; i++)
928  similarity[i]/=((double) columns*rows);
929  similarity[CompositeChannels]/=(double) GetNumberChannels(image,channel);
930  image->error.mean_error_per_pixel=QuantumRange*similarity[CompositeChannels];
931  image->error.normalized_mean_error=mean_error/((double) columns*rows);
932  image->error.normalized_maximum_error=maximum_error;
933  return(status);
934 }
935 
936 static MagickBooleanType GetMSESimilarity(const Image *image,
937  const Image *reconstruct_image,const ChannelType channel,
938  double *similarity,ExceptionInfo *exception)
939 {
940  CacheView
941  *image_view,
942  *reconstruct_view;
943 
944  double
945  area = 0.0;
946 
947  MagickBooleanType
948  status;
949 
950  size_t
951  columns,
952  rows;
953 
954  ssize_t
955  i,
956  y;
957 
958  status=MagickTrue;
959  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
960  image_view=AcquireVirtualCacheView(image,exception);
961  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
962 #if defined(MAGICKCORE_OPENMP_SUPPORT)
963  #pragma omp parallel for schedule(static) shared(similarity,status) \
964  magick_number_threads(image,image,rows,1)
965 #endif
966  for (y=0; y < (ssize_t) rows; y++)
967  {
968  double
969  channel_similarity[CompositeChannels+1] = { 0.0 };
970 
971  const IndexPacket
972  *magick_restrict indexes,
973  *magick_restrict reconstruct_indexes;
974 
975  const PixelPacket
976  *magick_restrict p,
977  *magick_restrict q;
978 
979  ssize_t
980  i,
981  x;
982 
983  if (status == MagickFalse)
984  continue;
985  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
986  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
987  if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
988  {
989  status=MagickFalse;
990  continue;
991  }
992  indexes=GetCacheViewVirtualIndexQueue(image_view);
993  reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
994  for (x=0; x < (ssize_t) columns; x++)
995  {
996  double
997  distance,
998  Da,
999  Sa;
1000 
1001  Sa=QuantumScale*(image->matte != MagickFalse ? (double) GetPixelAlpha(p) :
1002  ((double) QuantumRange-(double) OpaqueOpacity));
1003  Da=QuantumScale*(reconstruct_image->matte != MagickFalse ?
1004  (double) GetPixelAlpha(q) : ((double) QuantumRange-(double)
1005  OpaqueOpacity));
1006  if ((channel & RedChannel) != 0)
1007  {
1008  distance=QuantumScale*(Sa*(double) GetPixelRed(p)-Da*(double)
1009  GetPixelRed(q));
1010  channel_similarity[RedChannel]+=distance*distance;
1011  channel_similarity[CompositeChannels]+=distance*distance;
1012  }
1013  if ((channel & GreenChannel) != 0)
1014  {
1015  distance=QuantumScale*(Sa*(double) GetPixelGreen(p)-Da*(double)
1016  GetPixelGreen(q));
1017  channel_similarity[GreenChannel]+=distance*distance;
1018  channel_similarity[CompositeChannels]+=distance*distance;
1019  }
1020  if ((channel & BlueChannel) != 0)
1021  {
1022  distance=QuantumScale*(Sa*(double) GetPixelBlue(p)-Da*(double)
1023  GetPixelBlue(q));
1024  channel_similarity[BlueChannel]+=distance*distance;
1025  channel_similarity[CompositeChannels]+=distance*distance;
1026  }
1027  if (((channel & OpacityChannel) != 0) &&
1028  (image->matte != MagickFalse))
1029  {
1030  distance=QuantumScale*((double) GetPixelOpacity(p)-(double)
1031  GetPixelOpacity(q));
1032  channel_similarity[OpacityChannel]+=distance*distance;
1033  channel_similarity[CompositeChannels]+=distance*distance;
1034  }
1035  if (((channel & IndexChannel) != 0) &&
1036  (image->colorspace == CMYKColorspace) &&
1037  (reconstruct_image->colorspace == CMYKColorspace))
1038  {
1039  distance=QuantumScale*(Sa*(double) GetPixelIndex(indexes+x)-Da*
1040  (double) GetPixelIndex(reconstruct_indexes+x));
1041  channel_similarity[BlackChannel]+=distance*distance;
1042  channel_similarity[CompositeChannels]+=distance*distance;
1043  }
1044  p++;
1045  q++;
1046  }
1047 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1048  #pragma omp critical (MagickCore_GetMeanSquaredError)
1049 #endif
1050  for (i=0; i <= (ssize_t) CompositeChannels; i++)
1051  similarity[i]+=channel_similarity[i];
1052  }
1053  reconstruct_view=DestroyCacheView(reconstruct_view);
1054  image_view=DestroyCacheView(image_view);
1055  area=MagickSafeReciprocal((double) columns*rows);
1056  for (i=0; i <= (ssize_t) CompositeChannels; i++)
1057  similarity[i]*=area;
1058  similarity[CompositeChannels]/=(double) GetNumberChannels(image,channel);
1059  return(status);
1060 }
1061 
1062 static MagickBooleanType GetNCCSimilarity(const Image *image,
1063  const Image *reconstruct_image,const ChannelType channel,double *similarity,
1064  ExceptionInfo *exception)
1065 {
1066 #define SimilarityImageTag "Similarity/Image"
1067 
1068  CacheView
1069  *image_view,
1070  *reconstruct_view;
1071 
1072  ChannelStatistics
1073  *image_statistics,
1074  *reconstruct_statistics;
1075 
1076  double
1077  alpha_variance[CompositeChannels+1] = { 0.0 },
1078  beta_variance[CompositeChannels+1] = { 0.0 };
1079 
1080  MagickBooleanType
1081  status;
1082 
1083  MagickOffsetType
1084  progress;
1085 
1086  size_t
1087  columns,
1088  rows;
1089 
1090  ssize_t
1091  i,
1092  y;
1093 
1094  /*
1095  Normalize to account for variation due to lighting and exposure condition.
1096  */
1097  image_statistics=GetImageChannelStatistics(image,exception);
1098  reconstruct_statistics=GetImageChannelStatistics(reconstruct_image,exception);
1099  if ((image_statistics == (ChannelStatistics *) NULL) ||
1100  (reconstruct_statistics == (ChannelStatistics *) NULL))
1101  {
1102  if (image_statistics != (ChannelStatistics *) NULL)
1103  image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1104  image_statistics);
1105  if (reconstruct_statistics != (ChannelStatistics *) NULL)
1106  reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1107  reconstruct_statistics);
1108  return(MagickFalse);
1109  }
1110  (void) memset(similarity,0,(CompositeChannels+1)*sizeof(*similarity));
1111  status=MagickTrue;
1112  progress=0;
1113  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
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  const IndexPacket
1123  *magick_restrict indexes,
1124  *magick_restrict reconstruct_indexes;
1125 
1126  const PixelPacket
1127  *magick_restrict p,
1128  *magick_restrict q;
1129 
1130  double
1131  channel_alpha_variance[CompositeChannels+1] = { 0.0 },
1132  channel_beta_variance[CompositeChannels+1] = { 0.0 },
1133  channel_similarity[CompositeChannels+1] = { 0.0 };
1134 
1135  ssize_t
1136  x;
1137 
1138  if (status == MagickFalse)
1139  continue;
1140  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1141  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1142  if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1143  {
1144  status=MagickFalse;
1145  continue;
1146  }
1147  indexes=GetCacheViewVirtualIndexQueue(image_view);
1148  reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
1149  for (x=0; x < (ssize_t) columns; x++)
1150  {
1151  MagickRealType
1152  alpha,
1153  beta,
1154  Da,
1155  Sa;
1156 
1157  Sa=QuantumScale*(image->matte != MagickFalse ? (double) GetPixelAlpha(p) :
1158  (double) QuantumRange);
1159  Da=QuantumScale*(reconstruct_image->matte != MagickFalse ?
1160  (double) GetPixelAlpha(q) : (double) QuantumRange);
1161  if ((channel & RedChannel) != 0)
1162  {
1163  alpha=QuantumScale*(Sa*(double) GetPixelRed(p)-
1164  image_statistics[RedChannel].mean);
1165  beta=QuantumScale*(Da*(double) GetPixelRed(q)-
1166  reconstruct_statistics[RedChannel].mean);
1167  channel_similarity[RedChannel]+=alpha*beta;
1168  channel_similarity[CompositeChannels]+=alpha*beta;
1169  channel_alpha_variance[RedChannel]+=alpha*alpha;
1170  channel_alpha_variance[CompositeChannels]+=alpha*alpha;
1171  channel_beta_variance[RedChannel]+=beta*beta;
1172  channel_beta_variance[CompositeChannels]+=beta*beta;
1173  }
1174  if ((channel & GreenChannel) != 0)
1175  {
1176  alpha=QuantumScale*(Sa*(double) GetPixelGreen(p)-
1177  image_statistics[GreenChannel].mean);
1178  beta=QuantumScale*(Da*(double) GetPixelGreen(q)-
1179  reconstruct_statistics[GreenChannel].mean);
1180  channel_similarity[GreenChannel]+=alpha*beta;
1181  channel_similarity[CompositeChannels]+=alpha*beta;
1182  channel_alpha_variance[GreenChannel]+=alpha*alpha;
1183  channel_alpha_variance[CompositeChannels]+=alpha*alpha;
1184  channel_beta_variance[GreenChannel]+=beta*beta;
1185  channel_beta_variance[CompositeChannels]+=beta*beta;
1186  }
1187  if ((channel & BlueChannel) != 0)
1188  {
1189  alpha=QuantumScale*(Sa*(double) GetPixelBlue(p)-
1190  image_statistics[BlueChannel].mean);
1191  beta=QuantumScale*(Da*(double) GetPixelBlue(q)-
1192  reconstruct_statistics[BlueChannel].mean);
1193  channel_similarity[BlueChannel]+=alpha*beta;
1194  channel_alpha_variance[BlueChannel]+=alpha*alpha;
1195  channel_beta_variance[BlueChannel]+=beta*beta;
1196  }
1197  if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
1198  {
1199  alpha=QuantumScale*((double) GetPixelAlpha(p)-
1200  image_statistics[AlphaChannel].mean);
1201  beta=QuantumScale*((double) GetPixelAlpha(q)-
1202  reconstruct_statistics[AlphaChannel].mean);
1203  channel_similarity[OpacityChannel]+=alpha*beta;
1204  channel_similarity[CompositeChannels]+=alpha*beta;
1205  channel_alpha_variance[OpacityChannel]+=alpha*alpha;
1206  channel_alpha_variance[CompositeChannels]+=alpha*alpha;
1207  channel_beta_variance[OpacityChannel]+=beta*beta;
1208  channel_beta_variance[CompositeChannels]+=beta*beta;
1209  }
1210  if (((channel & IndexChannel) != 0) &&
1211  (image->colorspace == CMYKColorspace) &&
1212  (reconstruct_image->colorspace == CMYKColorspace))
1213  {
1214  alpha=QuantumScale*(Sa*(double) GetPixelIndex(indexes+x)-
1215  image_statistics[BlackChannel].mean);
1216  beta=QuantumScale*(Da*(double) GetPixelIndex(reconstruct_indexes+
1217  x)-reconstruct_statistics[BlackChannel].mean);
1218  channel_similarity[BlackChannel]+=alpha*beta;
1219  channel_similarity[CompositeChannels]+=alpha*beta;
1220  channel_alpha_variance[BlackChannel]+=alpha*alpha;
1221  channel_alpha_variance[CompositeChannels]+=alpha*alpha;
1222  channel_beta_variance[BlackChannel]+=beta*beta;
1223  channel_beta_variance[CompositeChannels]+=beta*beta;
1224  }
1225  p++;
1226  q++;
1227  }
1228 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1229  #pragma omp critical (GetNCCSimilarity)
1230 #endif
1231  {
1232  ssize_t
1233  j;
1234 
1235  for (j=0; j <= (ssize_t) CompositeChannels; j++)
1236  {
1237  similarity[j]+=channel_similarity[j];
1238  alpha_variance[j]+=channel_alpha_variance[j];
1239  beta_variance[j]+=channel_beta_variance[j];
1240  }
1241  }
1242  if (image->progress_monitor != (MagickProgressMonitor) NULL)
1243  {
1244  MagickBooleanType
1245  proceed;
1246 
1247 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1248  #pragma omp atomic
1249 #endif
1250  progress++;
1251  proceed=SetImageProgress(image,SimilarityImageTag,progress,rows);
1252  if (proceed == MagickFalse)
1253  status=MagickFalse;
1254  }
1255  }
1256  reconstruct_view=DestroyCacheView(reconstruct_view);
1257  image_view=DestroyCacheView(image_view);
1258  /*
1259  Divide by the standard deviation.
1260  */
1261  for (i=0; i <= (ssize_t) CompositeChannels; i++)
1262  similarity[i]*=MagickSafeReciprocal(sqrt(alpha_variance[i])*
1263  sqrt(beta_variance[i]));
1264  /*
1265  Free resources.
1266  */
1267  reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1268  reconstruct_statistics);
1269  image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1270  image_statistics);
1271  return(status);
1272 }
1273 
1274 static MagickBooleanType GetPASimilarity(const Image *image,
1275  const Image *reconstruct_image,const ChannelType channel,
1276  double *similarity,ExceptionInfo *exception)
1277 {
1278  CacheView
1279  *image_view,
1280  *reconstruct_view;
1281 
1282  MagickBooleanType
1283  status;
1284 
1285  size_t
1286  columns,
1287  rows;
1288 
1289  ssize_t
1290  y;
1291 
1292  status=MagickTrue;
1293  (void) memset(similarity,0,(CompositeChannels+1)*sizeof(*similarity));
1294  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
1295  image_view=AcquireVirtualCacheView(image,exception);
1296  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1297 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1298  #pragma omp parallel for schedule(static) shared(status) \
1299  magick_number_threads(image,image,rows,1)
1300 #endif
1301  for (y=0; y < (ssize_t) rows; y++)
1302  {
1303  double
1304  channel_similarity[CompositeChannels+1];
1305 
1306  const IndexPacket
1307  *magick_restrict indexes,
1308  *magick_restrict reconstruct_indexes;
1309 
1310  const PixelPacket
1311  *magick_restrict p,
1312  *magick_restrict q;
1313 
1314  ssize_t
1315  i,
1316  x;
1317 
1318  if (status == MagickFalse)
1319  continue;
1320  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1321  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1322  if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1323  {
1324  status=MagickFalse;
1325  continue;
1326  }
1327  indexes=GetCacheViewVirtualIndexQueue(image_view);
1328  reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
1329  (void) memset(channel_similarity,0,(CompositeChannels+1)*
1330  sizeof(*channel_similarity));
1331  for (x=0; x < (ssize_t) columns; x++)
1332  {
1333  MagickRealType
1334  distance,
1335  Da,
1336  Sa;
1337 
1338  Sa=QuantumScale*(image->matte != MagickFalse ? (double) GetPixelAlpha(p) :
1339  ((double) QuantumRange-(double) OpaqueOpacity));
1340  Da=QuantumScale*(reconstruct_image->matte != MagickFalse ?
1341  (double) GetPixelAlpha(q) : ((double) QuantumRange-(double)
1342  OpaqueOpacity));
1343  if ((channel & RedChannel) != 0)
1344  {
1345  distance=QuantumScale*fabs(Sa*(double) GetPixelRed(p)-Da*
1346  (double) GetPixelRed(q));
1347  if (distance > channel_similarity[RedChannel])
1348  channel_similarity[RedChannel]=distance;
1349  if (distance > channel_similarity[CompositeChannels])
1350  channel_similarity[CompositeChannels]=distance;
1351  }
1352  if ((channel & GreenChannel) != 0)
1353  {
1354  distance=QuantumScale*fabs(Sa*(double) GetPixelGreen(p)-Da*
1355  (double) GetPixelGreen(q));
1356  if (distance > channel_similarity[GreenChannel])
1357  channel_similarity[GreenChannel]=distance;
1358  if (distance > channel_similarity[CompositeChannels])
1359  channel_similarity[CompositeChannels]=distance;
1360  }
1361  if ((channel & BlueChannel) != 0)
1362  {
1363  distance=QuantumScale*fabs(Sa*(double) GetPixelBlue(p)-Da*
1364  (double) GetPixelBlue(q));
1365  if (distance > channel_similarity[BlueChannel])
1366  channel_similarity[BlueChannel]=distance;
1367  if (distance > channel_similarity[CompositeChannels])
1368  channel_similarity[CompositeChannels]=distance;
1369  }
1370  if (((channel & OpacityChannel) != 0) &&
1371  (image->matte != MagickFalse))
1372  {
1373  distance=QuantumScale*fabs((double) GetPixelOpacity(p)-(double)
1374  GetPixelOpacity(q));
1375  if (distance > channel_similarity[OpacityChannel])
1376  channel_similarity[OpacityChannel]=distance;
1377  if (distance > channel_similarity[CompositeChannels])
1378  channel_similarity[CompositeChannels]=distance;
1379  }
1380  if (((channel & IndexChannel) != 0) &&
1381  (image->colorspace == CMYKColorspace) &&
1382  (reconstruct_image->colorspace == CMYKColorspace))
1383  {
1384  distance=QuantumScale*fabs(Sa*(double) GetPixelIndex(indexes+x)-Da*
1385  (double) GetPixelIndex(reconstruct_indexes+x));
1386  if (distance > channel_similarity[BlackChannel])
1387  channel_similarity[BlackChannel]=distance;
1388  if (distance > channel_similarity[CompositeChannels])
1389  channel_similarity[CompositeChannels]=distance;
1390  }
1391  p++;
1392  q++;
1393  }
1394 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1395  #pragma omp critical (MagickCore_GetPeakAbsoluteError)
1396 #endif
1397  for (i=0; i <= (ssize_t) CompositeChannels; i++)
1398  if (channel_similarity[i] > similarity[i])
1399  similarity[i]=channel_similarity[i];
1400  }
1401  reconstruct_view=DestroyCacheView(reconstruct_view);
1402  image_view=DestroyCacheView(image_view);
1403  return(status);
1404 }
1405 
1406 static MagickBooleanType GetPSNRSimilarity(const Image *image,
1407  const Image *reconstruct_image,const ChannelType channel,
1408  double *similarity,ExceptionInfo *exception)
1409 {
1410  MagickBooleanType
1411  status;
1412 
1413  status=GetMSESimilarity(image,reconstruct_image,channel,similarity,
1414  exception);
1415  if ((channel & RedChannel) != 0)
1416  similarity[RedChannel]=10.0*MagickSafeLog10(MagickSafeReciprocal(
1417  similarity[RedChannel]))/MagickSafePSNRRecipicol(10.0);
1418  if ((channel & GreenChannel) != 0)
1419  similarity[GreenChannel]=10.0*MagickSafeLog10(MagickSafeReciprocal(
1420  similarity[GreenChannel]))/MagickSafePSNRRecipicol(10.0);
1421  if ((channel & BlueChannel) != 0)
1422  similarity[BlueChannel]=10.0*MagickSafeLog10(MagickSafeReciprocal(
1423  similarity[BlueChannel]))/MagickSafePSNRRecipicol(10.0);
1424  if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
1425  similarity[OpacityChannel]=10.0*MagickSafeLog10(MagickSafeReciprocal(
1426  similarity[OpacityChannel]))/MagickSafePSNRRecipicol(10.0);
1427  if (((channel & IndexChannel) != 0) && (image->colorspace == CMYKColorspace))
1428  similarity[BlackChannel]=10.0*MagickSafeLog10(MagickSafeReciprocal(
1429  similarity[BlackChannel]))/MagickSafePSNRRecipicol(10.0);
1430  similarity[CompositeChannels]=10.0*MagickSafeLog10(MagickSafeReciprocal(
1431  similarity[CompositeChannels]))/MagickSafePSNRRecipicol(10.0);
1432  return(status);
1433 }
1434 
1435 static MagickBooleanType GetPHASHSimilarity(const Image *image,
1436  const Image *reconstruct_image,const ChannelType channel,double *similarity,
1437  ExceptionInfo *exception)
1438 {
1439 #define PHASHNormalizationFactor 389.373723242
1440 
1441  ChannelPerceptualHash
1442  *image_phash,
1443  *reconstruct_phash;
1444 
1445  double
1446  error,
1447  difference;
1448 
1449  ssize_t
1450  i;
1451 
1452  /*
1453  Compute perceptual hash in the sRGB colorspace.
1454  */
1455  image_phash=GetImageChannelPerceptualHash(image,exception);
1456  if (image_phash == (ChannelPerceptualHash *) NULL)
1457  return(MagickFalse);
1458  reconstruct_phash=GetImageChannelPerceptualHash(reconstruct_image,exception);
1459  if (reconstruct_phash == (ChannelPerceptualHash *) NULL)
1460  {
1461  image_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(image_phash);
1462  return(MagickFalse);
1463  }
1464  for (i=0; i < MaximumNumberOfImageMoments; i++)
1465  {
1466  /*
1467  Compute sum of moment differences squared.
1468  */
1469  if ((channel & RedChannel) != 0)
1470  {
1471  error=reconstruct_phash[RedChannel].P[i]-image_phash[RedChannel].P[i];
1472  if (IsNaN(error) != 0)
1473  error=0.0;
1474  difference=error*error/PHASHNormalizationFactor;
1475  similarity[RedChannel]+=difference;
1476  similarity[CompositeChannels]+=difference;
1477  }
1478  if ((channel & GreenChannel) != 0)
1479  {
1480  error=reconstruct_phash[GreenChannel].P[i]-
1481  image_phash[GreenChannel].P[i];
1482  if (IsNaN(error) != 0)
1483  error=0.0;
1484  difference=error*error/PHASHNormalizationFactor;
1485  similarity[GreenChannel]+=difference;
1486  similarity[CompositeChannels]+=difference;
1487  }
1488  if ((channel & BlueChannel) != 0)
1489  {
1490  error=reconstruct_phash[BlueChannel].P[i]-image_phash[BlueChannel].P[i];
1491  if (IsNaN(error) != 0)
1492  error=0.0;
1493  difference=error*error/PHASHNormalizationFactor;
1494  similarity[BlueChannel]+=difference;
1495  similarity[CompositeChannels]+=difference;
1496  }
1497  if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse) &&
1498  (reconstruct_image->matte != MagickFalse))
1499  {
1500  error=reconstruct_phash[OpacityChannel].P[i]-
1501  image_phash[OpacityChannel].P[i];
1502  if (IsNaN(error) != 0)
1503  error=0.0;
1504  difference=error*error/PHASHNormalizationFactor;
1505  similarity[OpacityChannel]+=difference;
1506  similarity[CompositeChannels]+=difference;
1507  }
1508  if (((channel & IndexChannel) != 0) &&
1509  (image->colorspace == CMYKColorspace) &&
1510  (reconstruct_image->colorspace == CMYKColorspace))
1511  {
1512  error=reconstruct_phash[IndexChannel].P[i]-
1513  image_phash[IndexChannel].P[i];
1514  if (IsNaN(error) != 0)
1515  error=0.0;
1516  difference=error*error/PHASHNormalizationFactor;
1517  similarity[IndexChannel]+=difference;
1518  similarity[CompositeChannels]+=difference;
1519  }
1520  }
1521  /*
1522  Compute perceptual hash in the HCLP colorspace.
1523  */
1524  for (i=0; i < MaximumNumberOfImageMoments; i++)
1525  {
1526  /*
1527  Compute sum of moment differences squared.
1528  */
1529  if ((channel & RedChannel) != 0)
1530  {
1531  error=reconstruct_phash[RedChannel].Q[i]-image_phash[RedChannel].Q[i];
1532  if (IsNaN(error) != 0)
1533  error=0.0;
1534  difference=error*error/PHASHNormalizationFactor;
1535  similarity[RedChannel]+=difference;
1536  similarity[CompositeChannels]+=difference;
1537  }
1538  if ((channel & GreenChannel) != 0)
1539  {
1540  error=reconstruct_phash[GreenChannel].Q[i]-
1541  image_phash[GreenChannel].Q[i];
1542  if (IsNaN(error) != 0)
1543  error=0.0;
1544  difference=error*error/PHASHNormalizationFactor;
1545  similarity[GreenChannel]+=difference;
1546  similarity[CompositeChannels]+=difference;
1547  }
1548  if ((channel & BlueChannel) != 0)
1549  {
1550  error=reconstruct_phash[BlueChannel].Q[i]-image_phash[BlueChannel].Q[i];
1551  if (IsNaN(error) != 0)
1552  error=0.0;
1553  difference=error*error/PHASHNormalizationFactor;
1554  similarity[BlueChannel]+=difference;
1555  similarity[CompositeChannels]+=difference;
1556  }
1557  if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse) &&
1558  (reconstruct_image->matte != MagickFalse))
1559  {
1560  error=reconstruct_phash[OpacityChannel].Q[i]-
1561  image_phash[OpacityChannel].Q[i];
1562  if (IsNaN(error) != 0)
1563  error=0.0;
1564  difference=error*error/PHASHNormalizationFactor;
1565  similarity[OpacityChannel]+=difference;
1566  similarity[CompositeChannels]+=difference;
1567  }
1568  if (((channel & IndexChannel) != 0) &&
1569  (image->colorspace == CMYKColorspace) &&
1570  (reconstruct_image->colorspace == CMYKColorspace))
1571  {
1572  error=reconstruct_phash[IndexChannel].Q[i]-
1573  image_phash[IndexChannel].Q[i];
1574  if (IsNaN(error) != 0)
1575  error=0.0;
1576  difference=error*error/PHASHNormalizationFactor;
1577  similarity[IndexChannel]+=difference;
1578  similarity[CompositeChannels]+=difference;
1579  }
1580  }
1581  similarity[CompositeChannels]/=(double) GetNumberChannels(image,channel);
1582  /*
1583  Free resources.
1584  */
1585  reconstruct_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1586  reconstruct_phash);
1587  image_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(image_phash);
1588  return(MagickTrue);
1589 }
1590 
1591 static MagickBooleanType GetRMSESimilarity(const Image *image,
1592  const Image *reconstruct_image,const ChannelType channel,double *similarity,
1593  ExceptionInfo *exception)
1594 {
1595 #define RMSESquareRoot(x) sqrt((x) < 0.0 ? 0.0 : (x))
1596 
1597  MagickBooleanType
1598  status;
1599 
1600  status=GetMSESimilarity(image,reconstruct_image,channel,similarity,
1601  exception);
1602  if ((channel & RedChannel) != 0)
1603  similarity[RedChannel]=RMSESquareRoot(similarity[RedChannel]);
1604  if ((channel & GreenChannel) != 0)
1605  similarity[GreenChannel]=RMSESquareRoot(similarity[GreenChannel]);
1606  if ((channel & BlueChannel) != 0)
1607  similarity[BlueChannel]=RMSESquareRoot(similarity[BlueChannel]);
1608  if (((channel & OpacityChannel) != 0) &&
1609  (image->matte != MagickFalse))
1610  similarity[OpacityChannel]=RMSESquareRoot(similarity[OpacityChannel]);
1611  if (((channel & IndexChannel) != 0) &&
1612  (image->colorspace == CMYKColorspace))
1613  similarity[BlackChannel]=RMSESquareRoot(similarity[BlackChannel]);
1614  similarity[CompositeChannels]=RMSESquareRoot(similarity[CompositeChannels]);
1615  return(status);
1616 }
1617 
1618 MagickExport MagickBooleanType GetImageChannelDistortion(Image *image,
1619  const Image *reconstruct_image,const ChannelType channel,
1620  const MetricType metric,double *distortion,ExceptionInfo *exception)
1621 {
1622  double
1623  *channel_similarity;
1624 
1625  MagickBooleanType
1626  status;
1627 
1628  size_t
1629  length;
1630 
1631  assert(image != (Image *) NULL);
1632  assert(image->signature == MagickCoreSignature);
1633  assert(reconstruct_image != (const Image *) NULL);
1634  assert(reconstruct_image->signature == MagickCoreSignature);
1635  assert(distortion != (double *) NULL);
1636  if (IsEventLogging() != MagickFalse)
1637  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1638  *distortion=0.0;
1639  if (metric != PerceptualHashErrorMetric)
1640  if (ValidateImageMorphology(image,reconstruct_image) == MagickFalse)
1641  ThrowBinaryException(ImageError,"ImageMorphologyDiffers",image->filename);
1642  /*
1643  Get image distortion.
1644  */
1645  length=CompositeChannels+1UL;
1646  channel_similarity=(double *) AcquireQuantumMemory(length,
1647  sizeof(*channel_similarity));
1648  if (channel_similarity == (double *) NULL)
1649  ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1650  (void) memset(channel_similarity,0,length*sizeof(*channel_similarity));
1651  switch (metric)
1652  {
1653  case AbsoluteErrorMetric:
1654  {
1655  status=GetAESimilarity(image,reconstruct_image,channel,
1656  channel_similarity,exception);
1657  break;
1658  }
1659  case FuzzErrorMetric:
1660  {
1661  status=GetFUZZSimilarity(image,reconstruct_image,channel,
1662  channel_similarity,exception);
1663  break;
1664  }
1665  case MeanAbsoluteErrorMetric:
1666  {
1667  status=GetMAESimilarity(image,reconstruct_image,channel,
1668  channel_similarity,exception);
1669  break;
1670  }
1671  case MeanErrorPerPixelMetric:
1672  {
1673  status=GetMEPPSimilarity(image,reconstruct_image,channel,
1674  channel_similarity,exception);
1675  break;
1676  }
1677  case MeanSquaredErrorMetric:
1678  {
1679  status=GetMSESimilarity(image,reconstruct_image,channel,
1680  channel_similarity,exception);
1681  break;
1682  }
1683  case NormalizedCrossCorrelationErrorMetric:
1684  {
1685  status=GetNCCSimilarity(image,reconstruct_image,channel,
1686  channel_similarity,exception);
1687  break;
1688  }
1689  case PeakAbsoluteErrorMetric:
1690  {
1691  status=GetPASimilarity(image,reconstruct_image,channel,
1692  channel_similarity,exception);
1693  break;
1694  }
1695  case PeakSignalToNoiseRatioMetric:
1696  {
1697  status=GetPSNRSimilarity(image,reconstruct_image,channel,
1698  channel_similarity,exception);
1699  break;
1700  }
1701  case PerceptualHashErrorMetric:
1702  {
1703  status=GetPHASHSimilarity(image,reconstruct_image,channel,
1704  channel_similarity,exception);
1705  break;
1706  }
1707  case RootMeanSquaredErrorMetric:
1708  case UndefinedErrorMetric:
1709  default:
1710  {
1711  status=GetRMSESimilarity(image,reconstruct_image,channel,
1712  channel_similarity,exception);
1713  break;
1714  }
1715  }
1716  *distortion=channel_similarity[CompositeChannels];
1717  switch (metric)
1718  {
1719  case NormalizedCrossCorrelationErrorMetric:
1720  {
1721  *distortion=(1.0-(*distortion))/2.0;
1722  break;
1723  }
1724  default: break;
1725  }
1726  if (fabs(*distortion) < MagickEpsilon)
1727  *distortion=0.0;
1728  channel_similarity=(double *) RelinquishMagickMemory(channel_similarity);
1729  (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
1730  *distortion);
1731  return(status);
1732 }
1733 
1734 /*
1735 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1736 % %
1737 % %
1738 % %
1739 % G e t I m a g e C h a n n e l D i s t o r t i o n s %
1740 % %
1741 % %
1742 % %
1743 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1744 %
1745 % GetImageChannelDistortions() compares the image channels of an image to a
1746 % reconstructed image and returns the specified distortion metric for each
1747 % channel.
1748 %
1749 % The format of the GetImageChannelDistortions method is:
1750 %
1751 % double *GetImageChannelDistortions(const Image *image,
1752 % const Image *reconstruct_image,const MetricType metric,
1753 % ExceptionInfo *exception)
1754 %
1755 % A description of each parameter follows:
1756 %
1757 % o image: the image.
1758 %
1759 % o reconstruct_image: the reconstruct image.
1760 %
1761 % o metric: the metric.
1762 %
1763 % o exception: return any errors or warnings in this structure.
1764 %
1765 */
1766 MagickExport double *GetImageChannelDistortions(Image *image,
1767  const Image *reconstruct_image,const MetricType metric,
1768  ExceptionInfo *exception)
1769 {
1770  double
1771  *distortion,
1772  *similarity;
1773 
1774  MagickBooleanType
1775  status;
1776 
1777  size_t
1778  length;
1779 
1780  ssize_t
1781  i;
1782 
1783  assert(image != (Image *) NULL);
1784  assert(image->signature == MagickCoreSignature);
1785  assert(reconstruct_image != (const Image *) NULL);
1786  assert(reconstruct_image->signature == MagickCoreSignature);
1787  if (IsEventLogging() != MagickFalse)
1788  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1789  if (metric != PerceptualHashErrorMetric)
1790  if (ValidateImageMorphology(image,reconstruct_image) == MagickFalse)
1791  {
1792  (void) ThrowMagickException(&image->exception,GetMagickModule(),
1793  ImageError,"ImageMorphologyDiffers","`%s'",image->filename);
1794  return((double *) NULL);
1795  }
1796  /*
1797  Get image distortion.
1798  */
1799  length=CompositeChannels+1UL;
1800  similarity=(double *) AcquireQuantumMemory(length,
1801  sizeof(*similarity));
1802  if (similarity == (double *) NULL)
1803  ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1804  (void) memset(similarity,0,length*sizeof(*similarity));
1805  status=MagickTrue;
1806  switch (metric)
1807  {
1808  case AbsoluteErrorMetric:
1809  {
1810  status=GetAESimilarity(image,reconstruct_image,CompositeChannels,
1811  similarity,exception);
1812  break;
1813  }
1814  case FuzzErrorMetric:
1815  {
1816  status=GetFUZZSimilarity(image,reconstruct_image,CompositeChannels,
1817  similarity,exception);
1818  break;
1819  }
1820  case MeanAbsoluteErrorMetric:
1821  {
1822  status=GetMAESimilarity(image,reconstruct_image,CompositeChannels,
1823  similarity,exception);
1824  break;
1825  }
1826  case MeanErrorPerPixelMetric:
1827  {
1828  status=GetMEPPSimilarity(image,reconstruct_image,CompositeChannels,
1829  similarity,exception);
1830  break;
1831  }
1832  case MeanSquaredErrorMetric:
1833  {
1834  status=GetMSESimilarity(image,reconstruct_image,CompositeChannels,
1835  similarity,exception);
1836  break;
1837  }
1838  case NormalizedCrossCorrelationErrorMetric:
1839  {
1840  status=GetNCCSimilarity(image,reconstruct_image,CompositeChannels,
1841  similarity,exception);
1842  break;
1843  }
1844  case PeakAbsoluteErrorMetric:
1845  {
1846  status=GetPASimilarity(image,reconstruct_image,CompositeChannels,
1847  similarity,exception);
1848  break;
1849  }
1850  case PeakSignalToNoiseRatioMetric:
1851  {
1852  status=GetPSNRSimilarity(image,reconstruct_image,CompositeChannels,
1853  similarity,exception);
1854  break;
1855  }
1856  case PerceptualHashErrorMetric:
1857  {
1858  status=GetPHASHSimilarity(image,reconstruct_image,CompositeChannels,
1859  similarity,exception);
1860  break;
1861  }
1862  case RootMeanSquaredErrorMetric:
1863  case UndefinedErrorMetric:
1864  default:
1865  {
1866  status=GetRMSESimilarity(image,reconstruct_image,CompositeChannels,
1867  similarity,exception);
1868  break;
1869  }
1870  }
1871  if (status == MagickFalse)
1872  {
1873  similarity=(double *) RelinquishMagickMemory(similarity);
1874  return((double *) NULL);
1875  }
1876  distortion=similarity;
1877  switch (metric)
1878  {
1879  case NormalizedCrossCorrelationErrorMetric:
1880  {
1881  for (i=0; i <= (ssize_t) CompositeChannels; i++)
1882  distortion[i]=(1.0-distortion[i])/2.0;
1883  break;
1884  }
1885  default: break;
1886  }
1887  for (i=0; i <= (ssize_t) CompositeChannels; i++)
1888  if (fabs(distortion[i]) < MagickEpsilon)
1889  distortion[i]=0.0;
1890  return(distortion);
1891 }
1892 
1893 /*
1894 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1895 % %
1896 % %
1897 % %
1898 % I s I m a g e s E q u a l %
1899 % %
1900 % %
1901 % %
1902 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1903 %
1904 % IsImagesEqual() measures the difference between colors at each pixel
1905 % location of two images. A value other than 0 means the colors match
1906 % exactly. Otherwise an error measure is computed by summing over all
1907 % pixels in an image the distance squared in RGB space between each image
1908 % pixel and its corresponding pixel in the reconstruct image. The error
1909 % measure is assigned to these image members:
1910 %
1911 % o mean_error_per_pixel: The mean error for any single pixel in
1912 % the image.
1913 %
1914 % o normalized_mean_error: The normalized mean quantization error for
1915 % any single pixel in the image. This distance measure is normalized to
1916 % a range between 0 and 1. It is independent of the range of red, green,
1917 % and blue values in the image.
1918 %
1919 % o normalized_maximum_error: The normalized maximum quantization
1920 % error for any single pixel in the image. This distance measure is
1921 % normalized to a range between 0 and 1. It is independent of the range
1922 % of red, green, and blue values in your image.
1923 %
1924 % A small normalized mean square error, accessed as
1925 % image->normalized_mean_error, suggests the images are very similar in
1926 % spatial layout and color.
1927 %
1928 % The format of the IsImagesEqual method is:
1929 %
1930 % MagickBooleanType IsImagesEqual(Image *image,
1931 % const Image *reconstruct_image)
1932 %
1933 % A description of each parameter follows.
1934 %
1935 % o image: the image.
1936 %
1937 % o reconstruct_image: the reconstruct image.
1938 %
1939 */
1940 MagickExport MagickBooleanType IsImagesEqual(Image *image,
1941  const Image *reconstruct_image)
1942 {
1943  CacheView
1944  *image_view,
1945  *reconstruct_view;
1946 
1947  ExceptionInfo
1948  *exception;
1949 
1950  MagickBooleanType
1951  status;
1952 
1953  MagickRealType
1954  area,
1955  gamma,
1956  maximum_error,
1957  mean_error,
1958  mean_error_per_pixel;
1959 
1960  size_t
1961  columns,
1962  rows;
1963 
1964  ssize_t
1965  y;
1966 
1967  assert(image != (Image *) NULL);
1968  assert(image->signature == MagickCoreSignature);
1969  assert(reconstruct_image != (const Image *) NULL);
1970  assert(reconstruct_image->signature == MagickCoreSignature);
1971  exception=(&image->exception);
1972  if (ValidateImageMorphology(image,reconstruct_image) == MagickFalse)
1973  ThrowBinaryException(ImageError,"ImageMorphologyDiffers",image->filename);
1974  area=0.0;
1975  maximum_error=0.0;
1976  mean_error_per_pixel=0.0;
1977  mean_error=0.0;
1978  SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
1979  image_view=AcquireVirtualCacheView(image,exception);
1980  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1981  for (y=0; y < (ssize_t) rows; y++)
1982  {
1983  const IndexPacket
1984  *magick_restrict indexes,
1985  *magick_restrict reconstruct_indexes;
1986 
1987  const PixelPacket
1988  *magick_restrict p,
1989  *magick_restrict q;
1990 
1991  ssize_t
1992  x;
1993 
1994  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1995  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1996  if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1997  break;
1998  indexes=GetCacheViewVirtualIndexQueue(image_view);
1999  reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
2000  for (x=0; x < (ssize_t) columns; x++)
2001  {
2002  MagickRealType
2003  distance;
2004 
2005  distance=fabs((double) GetPixelRed(p)-(double) GetPixelRed(q));
2006  mean_error_per_pixel+=distance;
2007  mean_error+=distance*distance;
2008  if (distance > maximum_error)
2009  maximum_error=distance;
2010  area++;
2011  distance=fabs((double) GetPixelGreen(p)-(double) GetPixelGreen(q));
2012  mean_error_per_pixel+=distance;
2013  mean_error+=distance*distance;
2014  if (distance > maximum_error)
2015  maximum_error=distance;
2016  area++;
2017  distance=fabs((double) GetPixelBlue(p)-(double) GetPixelBlue(q));
2018  mean_error_per_pixel+=distance;
2019  mean_error+=distance*distance;
2020  if (distance > maximum_error)
2021  maximum_error=distance;
2022  area++;
2023  if (image->matte != MagickFalse)
2024  {
2025  distance=fabs((double) GetPixelOpacity(p)-(double)
2026  GetPixelOpacity(q));
2027  mean_error_per_pixel+=distance;
2028  mean_error+=distance*distance;
2029  if (distance > maximum_error)
2030  maximum_error=distance;
2031  area++;
2032  }
2033  if ((image->colorspace == CMYKColorspace) &&
2034  (reconstruct_image->colorspace == CMYKColorspace))
2035  {
2036  distance=fabs((double) GetPixelIndex(indexes+x)-(double)
2037  GetPixelIndex(reconstruct_indexes+x));
2038  mean_error_per_pixel+=distance;
2039  mean_error+=distance*distance;
2040  if (distance > maximum_error)
2041  maximum_error=distance;
2042  area++;
2043  }
2044  p++;
2045  q++;
2046  }
2047  }
2048  reconstruct_view=DestroyCacheView(reconstruct_view);
2049  image_view=DestroyCacheView(image_view);
2050  gamma=MagickSafeReciprocal(area);
2051  image->error.mean_error_per_pixel=gamma*mean_error_per_pixel;
2052  image->error.normalized_mean_error=gamma*QuantumScale*QuantumScale*mean_error;
2053  image->error.normalized_maximum_error=QuantumScale*maximum_error;
2054  status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
2055  return(status);
2056 }
2057 
2058 /*
2059 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2060 % %
2061 % %
2062 % %
2063 % S i m i l a r i t y I m a g e %
2064 % %
2065 % %
2066 % %
2067 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2068 %
2069 % SimilarityImage() compares the reference image of the image and returns the
2070 % best match offset. In addition, it returns a similarity image such that an
2071 % exact match location is completely white and if none of the pixels match,
2072 % black, otherwise some gray level in-between.
2073 %
2074 % The format of the SimilarityImageImage method is:
2075 %
2076 % Image *SimilarityImage(const Image *image,const Image *reference,
2077 % RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
2078 %
2079 % A description of each parameter follows:
2080 %
2081 % o image: the image.
2082 %
2083 % o reference: find an area of the image that closely resembles this image.
2084 %
2085 % o the best match offset of the reference image within the image.
2086 %
2087 % o similarity: the computed similarity between the images.
2088 %
2089 % o exception: return any errors or warnings in this structure.
2090 %
2091 */
2092 
2093 static double GetSimilarityMetric(const Image *image,
2094  const Image *reconstruct_image,const MetricType metric,
2095  const ssize_t x_offset,const ssize_t y_offset,ExceptionInfo *exception)
2096 {
2097  double
2098  *channel_similarity,
2099  similarity = 0.0;
2100 
2101  ExceptionInfo
2102  *sans_exception = AcquireExceptionInfo();
2103 
2104  Image
2105  *similarity_image;
2106 
2107  MagickBooleanType
2108  status = MagickTrue;
2109 
2110  RectangleInfo
2111  geometry;
2112 
2113  size_t
2114  length = CompositeChannels+1UL;
2115 
2116  SetGeometry(reconstruct_image,&geometry);
2117  geometry.x=x_offset;
2118  geometry.y=y_offset;
2119  similarity_image=CropImage(image,&geometry,sans_exception);
2120  sans_exception=DestroyExceptionInfo(sans_exception);
2121  if (similarity_image == (Image *) NULL)
2122  return(NAN);
2123  /*
2124  Get image distortion.
2125  */
2126  channel_similarity=(double *) AcquireQuantumMemory(length,
2127  sizeof(*channel_similarity));
2128  if (channel_similarity == (double *) NULL)
2129  ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
2130  (void) memset(channel_similarity,0,length*sizeof(*channel_similarity));
2131  switch (metric)
2132  {
2133  case AbsoluteErrorMetric:
2134  {
2135  status=GetAESimilarity(similarity_image,reconstruct_image,
2136  CompositeChannels,channel_similarity,exception);
2137  break;
2138  }
2139  case FuzzErrorMetric:
2140  {
2141  status=GetFUZZSimilarity(similarity_image,reconstruct_image,
2142  CompositeChannels,channel_similarity,exception);
2143  break;
2144  }
2145  case MeanAbsoluteErrorMetric:
2146  {
2147  status=GetMAESimilarity(similarity_image,reconstruct_image,
2148  CompositeChannels,channel_similarity,exception);
2149  break;
2150  }
2151  case MeanErrorPerPixelMetric:
2152  {
2153  status=GetMEPPSimilarity(similarity_image,reconstruct_image,
2154  CompositeChannels,channel_similarity,exception);
2155  break;
2156  }
2157  case MeanSquaredErrorMetric:
2158  {
2159  status=GetMSESimilarity(similarity_image,reconstruct_image,
2160  CompositeChannels,channel_similarity,exception);
2161  break;
2162  }
2163  case NormalizedCrossCorrelationErrorMetric:
2164  {
2165  status=GetNCCSimilarity(similarity_image,reconstruct_image,
2166  CompositeChannels,channel_similarity,exception);
2167  break;
2168  }
2169  case PeakAbsoluteErrorMetric:
2170  {
2171  status=GetPASimilarity(similarity_image,reconstruct_image,
2172  CompositeChannels,channel_similarity,exception);
2173  break;
2174  }
2175  case PeakSignalToNoiseRatioMetric:
2176  {
2177  status=GetPSNRSimilarity(similarity_image,reconstruct_image,
2178  CompositeChannels,channel_similarity,exception);
2179  break;
2180  }
2181  case PerceptualHashErrorMetric:
2182  {
2183  status=GetPHASHSimilarity(similarity_image,reconstruct_image,
2184  CompositeChannels,channel_similarity,exception);
2185  break;
2186  }
2187  case RootMeanSquaredErrorMetric:
2188  case UndefinedErrorMetric:
2189  default:
2190  {
2191  status=GetRMSESimilarity(similarity_image,reconstruct_image,
2192  CompositeChannels,channel_similarity,exception);
2193  break;
2194  }
2195  }
2196  similarity_image=DestroyImage(similarity_image);
2197  similarity=channel_similarity[CompositeChannels];
2198  channel_similarity=(double *) RelinquishMagickMemory(channel_similarity);
2199  if (status == MagickFalse)
2200  return(NAN);
2201  return(similarity);
2202 }
2203 
2204 MagickExport Image *SimilarityImage(Image *image,const Image *reference,
2205  RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
2206 {
2207  Image
2208  *similarity_image;
2209 
2210  similarity_image=SimilarityMetricImage(image,reference,
2211  RootMeanSquaredErrorMetric,offset,similarity_metric,exception);
2212  return(similarity_image);
2213 }
2214 
2215 MagickExport Image *SimilarityMetricImage(Image *image,const Image *reconstruct,
2216  const MetricType metric,RectangleInfo *offset,double *similarity_metric,
2217  ExceptionInfo *exception)
2218 {
2219 #define SimilarityImageTag "Similarity/Image"
2220 
2221  typedef struct
2222  {
2223  double
2224  similarity;
2225 
2226  ssize_t
2227  x,
2228  y;
2229  } SimilarityInfo;
2230 
2231  CacheView
2232  *similarity_view;
2233 
2234  const char
2235  *artifact;
2236 
2237  double
2238  similarity_threshold;
2239 
2240  Image
2241  *similarity_image = (Image *) NULL;
2242 
2243  MagickBooleanType
2244  status;
2245 
2246  MagickOffsetType
2247  progress;
2248 
2249  SimilarityInfo
2250  similarity_info = { 0 };
2251 
2252  ssize_t
2253  y;
2254 
2255  assert(image != (const Image *) NULL);
2256  assert(image->signature == MagickCoreSignature);
2257  assert(exception != (ExceptionInfo *) NULL);
2258  assert(exception->signature == MagickCoreSignature);
2259  assert(offset != (RectangleInfo *) NULL);
2260  if (IsEventLogging() != MagickFalse)
2261  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2262  SetGeometry(reconstruct,offset);
2263  *similarity_metric=0.0;
2264  offset->x=0;
2265  offset->y=0;
2266  if (ValidateImageMorphology(image,reconstruct) == MagickFalse)
2267  ThrowImageException(ImageError,"ImageMorphologyDiffers");
2268  if ((image->columns < reconstruct->columns) ||
2269  (image->rows < reconstruct->rows))
2270  {
2271  (void) ThrowMagickException(&image->exception,GetMagickModule(),
2272  OptionWarning,"GeometryDoesNotContainImage","`%s'",image->filename);
2273  return((Image *) NULL);
2274  }
2275  similarity_image=CloneImage(image,image->columns,image->rows,MagickTrue,
2276  exception);
2277  if (similarity_image == (Image *) NULL)
2278  return((Image *) NULL);
2279  similarity_image->depth=32;
2280  similarity_image->colorspace=GRAYColorspace;
2281  similarity_image->matte=MagickFalse;
2282  status=SetImageStorageClass(similarity_image,DirectClass);
2283  if (status == MagickFalse)
2284  {
2285  InheritException(exception,&similarity_image->exception);
2286  return(DestroyImage(similarity_image));
2287  }
2288  /*
2289  Measure similarity of reconstruction image against image.
2290  */
2291  similarity_threshold=DefaultSimilarityThreshold;
2292  artifact=GetImageArtifact(image,"compare:similarity-threshold");
2293  if (artifact != (const char *) NULL)
2294  similarity_threshold=StringToDouble(artifact,(char **) NULL);
2295  status=MagickTrue;
2296  similarity_info.similarity=GetSimilarityMetric(image,reconstruct,metric,
2297  similarity_info.x,similarity_info.y,exception);
2298  progress=0;
2299  similarity_view=AcquireVirtualCacheView(similarity_image,exception);
2300 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2301  #pragma omp parallel for schedule(static) shared(status,similarity_info) \
2302  magick_number_threads(image,reconstruct,similarity_image->rows << 2,1)
2303 #endif
2304  for (y=0; y < (ssize_t) similarity_image->rows; y++)
2305  {
2306  double
2307  similarity;
2308 
2309  MagickBooleanType
2310  threshold_trigger = MagickFalse;
2311 
2312  PixelPacket
2313  *magick_restrict q;
2314 
2315  SimilarityInfo
2316  channel_info = similarity_info;
2317 
2318  ssize_t
2319  x;
2320 
2321  if (status == MagickFalse)
2322  continue;
2323  if (threshold_trigger != MagickFalse)
2324  continue;
2325  q=QueueCacheViewAuthenticPixels(similarity_view,0,y,
2326  similarity_image->columns,1,exception);
2327  if (q == (PixelPacket *) NULL)
2328  {
2329  status=MagickFalse;
2330  continue;
2331  }
2332  for (x=0; x < (ssize_t) similarity_image->columns; x++)
2333  {
2334  MagickBooleanType
2335  update = MagickFalse;
2336 
2337  similarity=GetSimilarityMetric(image,reconstruct,metric,x,y,exception);
2338  switch (metric)
2339  {
2340  case NormalizedCrossCorrelationErrorMetric:
2341  case PeakSignalToNoiseRatioMetric:
2342  {
2343  if (similarity > channel_info.similarity)
2344  update=MagickTrue;
2345  break;
2346  }
2347  default:
2348  {
2349  if (similarity < channel_info.similarity)
2350  update=MagickTrue;
2351  break;
2352  }
2353  }
2354  if (update != MagickFalse)
2355  {
2356  channel_info.similarity=similarity;
2357  channel_info.x=x;
2358  channel_info.y=y;
2359  }
2360  switch (metric)
2361  {
2362  case NormalizedCrossCorrelationErrorMetric:
2363  case PeakSignalToNoiseRatioMetric:
2364  {
2365  SetPixelRed(q,ClampToQuantum((double) QuantumRange*similarity));
2366  break;
2367  }
2368  default:
2369  {
2370  SetPixelRed(q,ClampToQuantum((double) QuantumRange*(1.0-similarity)));
2371  break;
2372  }
2373  }
2374  SetPixelGreen(q,GetPixelRed(q));
2375  SetPixelBlue(q,GetPixelRed(q));
2376  q++;
2377  }
2378 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2379  #pragma omp critical (MagickCore_SimilarityMetricImage)
2380 #endif
2381  switch (metric)
2382  {
2383  case NormalizedCrossCorrelationErrorMetric:
2384  case PeakSignalToNoiseRatioMetric:
2385  {
2386  if (similarity_threshold != DefaultSimilarityThreshold)
2387  if (channel_info.similarity >= similarity_threshold)
2388  threshold_trigger=MagickTrue;
2389  if (channel_info.similarity >= similarity_info.similarity)
2390  similarity_info=channel_info;
2391  break;
2392  }
2393  default:
2394  {
2395  if (similarity_threshold != DefaultSimilarityThreshold)
2396  if (channel_info.similarity < similarity_threshold)
2397  threshold_trigger=MagickTrue;
2398  if (channel_info.similarity < similarity_info.similarity)
2399  similarity_info=channel_info;
2400  break;
2401  }
2402  }
2403  if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
2404  status=MagickFalse;
2405  if (image->progress_monitor != (MagickProgressMonitor) NULL)
2406  {
2407  MagickBooleanType
2408  proceed;
2409 
2410  progress++;
2411  proceed=SetImageProgress(image,SimilarityImageTag,progress,image->rows);
2412  if (proceed == MagickFalse)
2413  status=MagickFalse;
2414  }
2415  }
2416  similarity_view=DestroyCacheView(similarity_view);
2417  if (status == MagickFalse)
2418  similarity_image=DestroyImage(similarity_image);
2419  *similarity_metric=similarity_info.similarity;
2420  if (fabs(*similarity_metric) < MagickEpsilon)
2421  *similarity_metric=0.0;
2422  offset->x=similarity_info.x;
2423  offset->y=similarity_info.y;
2424  (void) FormatImageProperty((Image *) image,"similarity","%.*g",
2425  GetMagickPrecision(),*similarity_metric);
2426  (void) FormatImageProperty((Image *) image,"similarity.offset.x","%.*g",
2427  GetMagickPrecision(),(double) offset->x);
2428  (void) FormatImageProperty((Image *) image,"similarity.offset.y","%.*g",
2429  GetMagickPrecision(),(double) offset->y);
2430  return(similarity_image);
2431 }
_ExceptionInfo
Definition: exception.h:102
_CacheView
Definition: cache-view.c:65
_PixelPacket
Definition: pixel.h:131
_Image
Definition: image.h:133
_MagickPixelPacket
Definition: pixel.h:104
_RectangleInfo
Definition: geometry.h:123
_ChannelPerceptualHash
Definition: statistic.h:67
_ChannelStatistics
Definition: statistic.h:30