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