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