MagickCore  6.9.12-54
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  if (fabs(gamma) >= MagickEpsilon)
1231  {
1232  gamma=PerceptibleReciprocal(gamma);
1233  distortion[i]=QuantumRange*gamma*distortion[i];
1234  }
1235  }
1236  distortion[CompositeChannels]=0.0;
1237  if ((channel & RedChannel) != 0)
1238  distortion[CompositeChannels]+=distortion[RedChannel]*
1239  distortion[RedChannel];
1240  if ((channel & GreenChannel) != 0)
1241  distortion[CompositeChannels]+=distortion[GreenChannel]*
1242  distortion[GreenChannel];
1243  if ((channel & BlueChannel) != 0)
1244  distortion[CompositeChannels]+=distortion[BlueChannel]*
1245  distortion[BlueChannel];
1246  if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
1247  distortion[CompositeChannels]+=distortion[OpacityChannel]*
1248  distortion[OpacityChannel];
1249  if (((channel & IndexChannel) != 0) && (image->colorspace == CMYKColorspace))
1250  distortion[CompositeChannels]+=distortion[BlackChannel]*
1251  distortion[BlackChannel];
1252  distortion[CompositeChannels]=sqrt(distortion[CompositeChannels]/
1253  GetNumberChannels(image,channel));
1254  /*
1255  Free resources.
1256  */
1257  reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1258  reconstruct_statistics);
1259  image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1260  image_statistics);
1261  return(status);
1262 }
1263 
1264 static MagickBooleanType GetPeakAbsoluteDistortion(const Image *image,
1265  const Image *reconstruct_image,const ChannelType channel,
1266  double *distortion,ExceptionInfo *exception)
1267 {
1268  CacheView
1269  *image_view,
1270  *reconstruct_view;
1271 
1272  MagickBooleanType
1273  status;
1274 
1275  size_t
1276  columns,
1277  rows;
1278 
1279  ssize_t
1280  y;
1281 
1282  status=MagickTrue;
1283  rows=MagickMax(image->rows,reconstruct_image->rows);
1284  columns=MagickMax(image->columns,reconstruct_image->columns);
1285  image_view=AcquireVirtualCacheView(image,exception);
1286  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1287 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1288  #pragma omp parallel for schedule(static) shared(status) \
1289  magick_number_threads(image,image,rows,1)
1290 #endif
1291  for (y=0; y < (ssize_t) rows; y++)
1292  {
1293  double
1294  channel_distortion[CompositeChannels+1];
1295 
1296  const IndexPacket
1297  *magick_restrict indexes,
1298  *magick_restrict reconstruct_indexes;
1299 
1300  const PixelPacket
1301  *magick_restrict p,
1302  *magick_restrict q;
1303 
1304  ssize_t
1305  i,
1306  x;
1307 
1308  if (status == MagickFalse)
1309  continue;
1310  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1311  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1312  if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1313  {
1314  status=MagickFalse;
1315  continue;
1316  }
1317  indexes=GetCacheViewVirtualIndexQueue(image_view);
1318  reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
1319  (void) memset(channel_distortion,0,sizeof(channel_distortion));
1320  for (x=0; x < (ssize_t) columns; x++)
1321  {
1322  MagickRealType
1323  distance,
1324  Da,
1325  Sa;
1326 
1327  Sa=QuantumScale*(image->matte != MagickFalse ? GetPixelAlpha(p) :
1328  (QuantumRange-OpaqueOpacity));
1329  Da=QuantumScale*(reconstruct_image->matte != MagickFalse ?
1330  GetPixelAlpha(q) : (QuantumRange-OpaqueOpacity));
1331  if ((channel & RedChannel) != 0)
1332  {
1333  distance=QuantumScale*fabs((double) (Sa*GetPixelRed(p)-Da*
1334  GetPixelRed(q)));
1335  if (distance > channel_distortion[RedChannel])
1336  channel_distortion[RedChannel]=distance;
1337  if (distance > channel_distortion[CompositeChannels])
1338  channel_distortion[CompositeChannels]=distance;
1339  }
1340  if ((channel & GreenChannel) != 0)
1341  {
1342  distance=QuantumScale*fabs((double) (Sa*GetPixelGreen(p)-Da*
1343  GetPixelGreen(q)));
1344  if (distance > channel_distortion[GreenChannel])
1345  channel_distortion[GreenChannel]=distance;
1346  if (distance > channel_distortion[CompositeChannels])
1347  channel_distortion[CompositeChannels]=distance;
1348  }
1349  if ((channel & BlueChannel) != 0)
1350  {
1351  distance=QuantumScale*fabs((double) (Sa*GetPixelBlue(p)-Da*
1352  GetPixelBlue(q)));
1353  if (distance > channel_distortion[BlueChannel])
1354  channel_distortion[BlueChannel]=distance;
1355  if (distance > channel_distortion[CompositeChannels])
1356  channel_distortion[CompositeChannels]=distance;
1357  }
1358  if (((channel & OpacityChannel) != 0) &&
1359  (image->matte != MagickFalse))
1360  {
1361  distance=QuantumScale*fabs((double) (GetPixelOpacity(p)-(double)
1362  GetPixelOpacity(q)));
1363  if (distance > channel_distortion[OpacityChannel])
1364  channel_distortion[OpacityChannel]=distance;
1365  if (distance > channel_distortion[CompositeChannels])
1366  channel_distortion[CompositeChannels]=distance;
1367  }
1368  if (((channel & IndexChannel) != 0) &&
1369  (image->colorspace == CMYKColorspace) &&
1370  (reconstruct_image->colorspace == CMYKColorspace))
1371  {
1372  distance=QuantumScale*fabs((double) (Sa*GetPixelIndex(indexes+x)-Da*
1373  GetPixelIndex(reconstruct_indexes+x)));
1374  if (distance > channel_distortion[BlackChannel])
1375  channel_distortion[BlackChannel]=distance;
1376  if (distance > channel_distortion[CompositeChannels])
1377  channel_distortion[CompositeChannels]=distance;
1378  }
1379  p++;
1380  q++;
1381  }
1382 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1383  #pragma omp critical (MagickCore_GetPeakAbsoluteError)
1384 #endif
1385  for (i=0; i <= (ssize_t) CompositeChannels; i++)
1386  if (channel_distortion[i] > distortion[i])
1387  distortion[i]=channel_distortion[i];
1388  }
1389  reconstruct_view=DestroyCacheView(reconstruct_view);
1390  image_view=DestroyCacheView(image_view);
1391  return(status);
1392 }
1393 
1394 static inline double MagickLog10(const double x)
1395 {
1396 #define Log10Epsilon (1.0e-11)
1397 
1398  if (fabs(x) < Log10Epsilon)
1399  return(log10(Log10Epsilon));
1400  return(log10(fabs(x)));
1401 }
1402 
1403 static MagickBooleanType GetPeakSignalToNoiseRatio(const Image *image,
1404  const Image *reconstruct_image,const ChannelType channel,
1405  double *distortion,ExceptionInfo *exception)
1406 {
1407  MagickBooleanType
1408  status;
1409 
1410  status=GetMeanSquaredDistortion(image,reconstruct_image,channel,distortion,
1411  exception);
1412  if ((channel & RedChannel) != 0)
1413  {
1414  if (fabs(distortion[RedChannel]) < MagickEpsilon)
1415  distortion[RedChannel]=INFINITY;
1416  else
1417  distortion[RedChannel]=10.0*MagickLog10(1.0)-10.0*
1418  MagickLog10(distortion[RedChannel]);
1419  }
1420  if ((channel & GreenChannel) != 0)
1421  {
1422  if (fabs(distortion[GreenChannel]) < MagickEpsilon)
1423  distortion[GreenChannel]=INFINITY;
1424  else
1425  distortion[GreenChannel]=10.0*MagickLog10(1.0)-10.0*
1426  MagickLog10(distortion[GreenChannel]);
1427  }
1428  if ((channel & BlueChannel) != 0)
1429  {
1430  if (fabs(distortion[BlueChannel]) < MagickEpsilon)
1431  distortion[BlueChannel]=INFINITY;
1432  else
1433  distortion[BlueChannel]=10.0*MagickLog10(1.0)-10.0*
1434  MagickLog10(distortion[BlueChannel]);
1435  }
1436  if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
1437  {
1438  if (fabs(distortion[OpacityChannel]) < MagickEpsilon)
1439  distortion[OpacityChannel]=INFINITY;
1440  else
1441  distortion[OpacityChannel]=10.0*MagickLog10(1.0)-10.0*
1442  MagickLog10(distortion[OpacityChannel]);
1443  }
1444  if (((channel & IndexChannel) != 0) && (image->colorspace == CMYKColorspace))
1445  {
1446  if (fabs(distortion[BlackChannel]) < MagickEpsilon)
1447  distortion[BlackChannel]=INFINITY;
1448  else
1449  distortion[BlackChannel]=10.0*MagickLog10(1.0)-10.0*
1450  MagickLog10(distortion[BlackChannel]);
1451  }
1452  if (fabs(distortion[CompositeChannels]) < MagickEpsilon)
1453  distortion[CompositeChannels]=INFINITY;
1454  else
1455  distortion[CompositeChannels]=10.0*MagickLog10(1.0)-10.0*
1456  MagickLog10(distortion[CompositeChannels]);
1457  return(status);
1458 }
1459 
1460 static MagickBooleanType GetPerceptualHashDistortion(const Image *image,
1461  const Image *reconstruct_image,const ChannelType channel,double *distortion,
1462  ExceptionInfo *exception)
1463 {
1465  *image_phash,
1466  *reconstruct_phash;
1467 
1468  double
1469  difference;
1470 
1471  ssize_t
1472  i;
1473 
1474  /*
1475  Compute perceptual hash in the sRGB colorspace.
1476  */
1477  image_phash=GetImageChannelPerceptualHash(image,exception);
1478  if (image_phash == (ChannelPerceptualHash *) NULL)
1479  return(MagickFalse);
1480  reconstruct_phash=GetImageChannelPerceptualHash(reconstruct_image,exception);
1481  if (reconstruct_phash == (ChannelPerceptualHash *) NULL)
1482  {
1483  image_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(image_phash);
1484  return(MagickFalse);
1485  }
1486  for (i=0; i < MaximumNumberOfImageMoments; i++)
1487  {
1488  /*
1489  Compute sum of moment differences squared.
1490  */
1491  if ((channel & RedChannel) != 0)
1492  {
1493  difference=reconstruct_phash[RedChannel].P[i]-
1494  image_phash[RedChannel].P[i];
1495  distortion[RedChannel]+=difference*difference;
1496  distortion[CompositeChannels]+=difference*difference;
1497  }
1498  if ((channel & GreenChannel) != 0)
1499  {
1500  difference=reconstruct_phash[GreenChannel].P[i]-
1501  image_phash[GreenChannel].P[i];
1502  distortion[GreenChannel]+=difference*difference;
1503  distortion[CompositeChannels]+=difference*difference;
1504  }
1505  if ((channel & BlueChannel) != 0)
1506  {
1507  difference=reconstruct_phash[BlueChannel].P[i]-
1508  image_phash[BlueChannel].P[i];
1509  distortion[BlueChannel]+=difference*difference;
1510  distortion[CompositeChannels]+=difference*difference;
1511  }
1512  if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse) &&
1513  (reconstruct_image->matte != MagickFalse))
1514  {
1515  difference=reconstruct_phash[OpacityChannel].P[i]-
1516  image_phash[OpacityChannel].P[i];
1517  distortion[OpacityChannel]+=difference*difference;
1518  distortion[CompositeChannels]+=difference*difference;
1519  }
1520  if (((channel & IndexChannel) != 0) &&
1521  (image->colorspace == CMYKColorspace) &&
1522  (reconstruct_image->colorspace == CMYKColorspace))
1523  {
1524  difference=reconstruct_phash[IndexChannel].P[i]-
1525  image_phash[IndexChannel].P[i];
1526  distortion[IndexChannel]+=difference*difference;
1527  distortion[CompositeChannels]+=difference*difference;
1528  }
1529  }
1530  /*
1531  Compute perceptual hash in the HCLP colorspace.
1532  */
1533  for (i=0; i < MaximumNumberOfImageMoments; i++)
1534  {
1535  /*
1536  Compute sum of moment differences squared.
1537  */
1538  if ((channel & RedChannel) != 0)
1539  {
1540  difference=reconstruct_phash[RedChannel].Q[i]-
1541  image_phash[RedChannel].Q[i];
1542  distortion[RedChannel]+=difference*difference;
1543  distortion[CompositeChannels]+=difference*difference;
1544  }
1545  if ((channel & GreenChannel) != 0)
1546  {
1547  difference=reconstruct_phash[GreenChannel].Q[i]-
1548  image_phash[GreenChannel].Q[i];
1549  distortion[GreenChannel]+=difference*difference;
1550  distortion[CompositeChannels]+=difference*difference;
1551  }
1552  if ((channel & BlueChannel) != 0)
1553  {
1554  difference=reconstruct_phash[BlueChannel].Q[i]-
1555  image_phash[BlueChannel].Q[i];
1556  distortion[BlueChannel]+=difference*difference;
1557  distortion[CompositeChannels]+=difference*difference;
1558  }
1559  if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse) &&
1560  (reconstruct_image->matte != MagickFalse))
1561  {
1562  difference=reconstruct_phash[OpacityChannel].Q[i]-
1563  image_phash[OpacityChannel].Q[i];
1564  distortion[OpacityChannel]+=difference*difference;
1565  distortion[CompositeChannels]+=difference*difference;
1566  }
1567  if (((channel & IndexChannel) != 0) &&
1568  (image->colorspace == CMYKColorspace) &&
1569  (reconstruct_image->colorspace == CMYKColorspace))
1570  {
1571  difference=reconstruct_phash[IndexChannel].Q[i]-
1572  image_phash[IndexChannel].Q[i];
1573  distortion[IndexChannel]+=difference*difference;
1574  distortion[CompositeChannels]+=difference*difference;
1575  }
1576  }
1577  /*
1578  Free resources.
1579  */
1580  reconstruct_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1581  reconstruct_phash);
1582  image_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(image_phash);
1583  return(MagickTrue);
1584 }
1585 
1586 static MagickBooleanType GetRootMeanSquaredDistortion(const Image *image,
1587  const Image *reconstruct_image,const ChannelType channel,double *distortion,
1588  ExceptionInfo *exception)
1589 {
1590  MagickBooleanType
1591  status;
1592 
1593  status=GetMeanSquaredDistortion(image,reconstruct_image,channel,distortion,
1594  exception);
1595  if ((channel & RedChannel) != 0)
1596  distortion[RedChannel]=sqrt(distortion[RedChannel]);
1597  if ((channel & GreenChannel) != 0)
1598  distortion[GreenChannel]=sqrt(distortion[GreenChannel]);
1599  if ((channel & BlueChannel) != 0)
1600  distortion[BlueChannel]=sqrt(distortion[BlueChannel]);
1601  if (((channel & OpacityChannel) != 0) &&
1602  (image->matte != MagickFalse))
1603  distortion[OpacityChannel]=sqrt(distortion[OpacityChannel]);
1604  if (((channel & IndexChannel) != 0) &&
1605  (image->colorspace == CMYKColorspace))
1606  distortion[BlackChannel]=sqrt(distortion[BlackChannel]);
1607  distortion[CompositeChannels]=sqrt(distortion[CompositeChannels]);
1608  return(status);
1609 }
1610 
1611 MagickExport MagickBooleanType GetImageChannelDistortion(Image *image,
1612  const Image *reconstruct_image,const ChannelType channel,
1613  const MetricType metric,double *distortion,ExceptionInfo *exception)
1614 {
1615  double
1616  *channel_distortion;
1617 
1618  MagickBooleanType
1619  status;
1620 
1621  size_t
1622  length;
1623 
1624  assert(image != (Image *) NULL);
1625  assert(image->signature == MagickCoreSignature);
1626  assert(reconstruct_image != (const Image *) NULL);
1627  assert(reconstruct_image->signature == MagickCoreSignature);
1628  assert(distortion != (double *) NULL);
1629  if (IsEventLogging() != MagickFalse)
1630  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1631  *distortion=0.0;
1632  if (metric != PerceptualHashErrorMetric)
1633  if (ValidateImageMorphology(image,reconstruct_image) == MagickFalse)
1634  ThrowBinaryException(ImageError,"ImageMorphologyDiffers",image->filename);
1635  /*
1636  Get image distortion.
1637  */
1638  length=CompositeChannels+1UL;
1639  channel_distortion=(double *) AcquireQuantumMemory(length,
1640  sizeof(*channel_distortion));
1641  if (channel_distortion == (double *) NULL)
1642  ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1643  (void) memset(channel_distortion,0,length*sizeof(*channel_distortion));
1644  switch (metric)
1645  {
1646  case AbsoluteErrorMetric:
1647  {
1648  status=GetAbsoluteDistortion(image,reconstruct_image,channel,
1649  channel_distortion,exception);
1650  break;
1651  }
1652  case FuzzErrorMetric:
1653  {
1654  status=GetFuzzDistortion(image,reconstruct_image,channel,
1655  channel_distortion,exception);
1656  break;
1657  }
1658  case MeanAbsoluteErrorMetric:
1659  {
1660  status=GetMeanAbsoluteDistortion(image,reconstruct_image,channel,
1661  channel_distortion,exception);
1662  break;
1663  }
1664  case MeanErrorPerPixelMetric:
1665  {
1666  status=GetMeanErrorPerPixel(image,reconstruct_image,channel,
1667  channel_distortion,exception);
1668  break;
1669  }
1670  case MeanSquaredErrorMetric:
1671  {
1672  status=GetMeanSquaredDistortion(image,reconstruct_image,channel,
1673  channel_distortion,exception);
1674  break;
1675  }
1676  case NormalizedCrossCorrelationErrorMetric:
1677  default:
1678  {
1679  status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1680  channel,channel_distortion,exception);
1681  break;
1682  }
1683  case PeakAbsoluteErrorMetric:
1684  {
1685  status=GetPeakAbsoluteDistortion(image,reconstruct_image,channel,
1686  channel_distortion,exception);
1687  break;
1688  }
1689  case PeakSignalToNoiseRatioMetric:
1690  {
1691  status=GetPeakSignalToNoiseRatio(image,reconstruct_image,channel,
1692  channel_distortion,exception);
1693  break;
1694  }
1695  case PerceptualHashErrorMetric:
1696  {
1697  status=GetPerceptualHashDistortion(image,reconstruct_image,channel,
1698  channel_distortion,exception);
1699  break;
1700  }
1701  case RootMeanSquaredErrorMetric:
1702  {
1703  status=GetRootMeanSquaredDistortion(image,reconstruct_image,channel,
1704  channel_distortion,exception);
1705  break;
1706  }
1707  }
1708  *distortion=channel_distortion[CompositeChannels];
1709  channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1710  (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
1711  *distortion);
1712  return(status);
1713 }
1714 
1715 /*
1716 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1717 % %
1718 % %
1719 % %
1720 % 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 %
1721 % %
1722 % %
1723 % %
1724 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1725 %
1726 % GetImageChannelDistortions() compares the image channels of an image to a
1727 % reconstructed image and returns the specified distortion metric for each
1728 % channel.
1729 %
1730 % The format of the GetImageChannelDistortions method is:
1731 %
1732 % double *GetImageChannelDistortions(const Image *image,
1733 % const Image *reconstruct_image,const MetricType metric,
1734 % ExceptionInfo *exception)
1735 %
1736 % A description of each parameter follows:
1737 %
1738 % o image: the image.
1739 %
1740 % o reconstruct_image: the reconstruct image.
1741 %
1742 % o metric: the metric.
1743 %
1744 % o exception: return any errors or warnings in this structure.
1745 %
1746 */
1747 MagickExport double *GetImageChannelDistortions(Image *image,
1748  const Image *reconstruct_image,const MetricType metric,
1749  ExceptionInfo *exception)
1750 {
1751  double
1752  *channel_distortion;
1753 
1754  MagickBooleanType
1755  status;
1756 
1757  size_t
1758  length;
1759 
1760  assert(image != (Image *) NULL);
1761  assert(image->signature == MagickCoreSignature);
1762  assert(reconstruct_image != (const Image *) NULL);
1763  assert(reconstruct_image->signature == MagickCoreSignature);
1764  if (IsEventLogging() != MagickFalse)
1765  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1766  if (metric != PerceptualHashErrorMetric)
1767  if (ValidateImageMorphology(image,reconstruct_image) == MagickFalse)
1768  {
1769  (void) ThrowMagickException(&image->exception,GetMagickModule(),
1770  ImageError,"ImageMorphologyDiffers","`%s'",image->filename);
1771  return((double *) NULL);
1772  }
1773  /*
1774  Get image distortion.
1775  */
1776  length=CompositeChannels+1UL;
1777  channel_distortion=(double *) AcquireQuantumMemory(length,
1778  sizeof(*channel_distortion));
1779  if (channel_distortion == (double *) NULL)
1780  ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1781  (void) memset(channel_distortion,0,length*
1782  sizeof(*channel_distortion));
1783  status=MagickTrue;
1784  switch (metric)
1785  {
1786  case AbsoluteErrorMetric:
1787  {
1788  status=GetAbsoluteDistortion(image,reconstruct_image,CompositeChannels,
1789  channel_distortion,exception);
1790  break;
1791  }
1792  case FuzzErrorMetric:
1793  {
1794  status=GetFuzzDistortion(image,reconstruct_image,CompositeChannels,
1795  channel_distortion,exception);
1796  break;
1797  }
1798  case MeanAbsoluteErrorMetric:
1799  {
1800  status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1801  CompositeChannels,channel_distortion,exception);
1802  break;
1803  }
1804  case MeanErrorPerPixelMetric:
1805  {
1806  status=GetMeanErrorPerPixel(image,reconstruct_image,CompositeChannels,
1807  channel_distortion,exception);
1808  break;
1809  }
1810  case MeanSquaredErrorMetric:
1811  {
1812  status=GetMeanSquaredDistortion(image,reconstruct_image,CompositeChannels,
1813  channel_distortion,exception);
1814  break;
1815  }
1816  case NormalizedCrossCorrelationErrorMetric:
1817  default:
1818  {
1819  status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1820  CompositeChannels,channel_distortion,exception);
1821  break;
1822  }
1823  case PeakAbsoluteErrorMetric:
1824  {
1825  status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1826  CompositeChannels,channel_distortion,exception);
1827  break;
1828  }
1829  case PeakSignalToNoiseRatioMetric:
1830  {
1831  status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1832  CompositeChannels,channel_distortion,exception);
1833  break;
1834  }
1835  case PerceptualHashErrorMetric:
1836  {
1837  status=GetPerceptualHashDistortion(image,reconstruct_image,
1838  CompositeChannels,channel_distortion,exception);
1839  break;
1840  }
1841  case RootMeanSquaredErrorMetric:
1842  {
1843  status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1844  CompositeChannels,channel_distortion,exception);
1845  break;
1846  }
1847  }
1848  if (status == MagickFalse)
1849  {
1850  channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1851  return((double *) NULL);
1852  }
1853  return(channel_distortion);
1854 }
1855 
1856 /*
1857 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1858 % %
1859 % %
1860 % %
1861 % I s I m a g e s E q u a l %
1862 % %
1863 % %
1864 % %
1865 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1866 %
1867 % IsImagesEqual() measures the difference between colors at each pixel
1868 % location of two images. A value other than 0 means the colors match
1869 % exactly. Otherwise an error measure is computed by summing over all
1870 % pixels in an image the distance squared in RGB space between each image
1871 % pixel and its corresponding pixel in the reconstruct image. The error
1872 % measure is assigned to these image members:
1873 %
1874 % o mean_error_per_pixel: The mean error for any single pixel in
1875 % the image.
1876 %
1877 % o normalized_mean_error: The normalized mean quantization error for
1878 % any single pixel in the image. This distance measure is normalized to
1879 % a range between 0 and 1. It is independent of the range of red, green,
1880 % and blue values in the image.
1881 %
1882 % o normalized_maximum_error: The normalized maximum quantization
1883 % error for any single pixel in the image. This distance measure is
1884 % normalized to a range between 0 and 1. It is independent of the range
1885 % of red, green, and blue values in your image.
1886 %
1887 % A small normalized mean square error, accessed as
1888 % image->normalized_mean_error, suggests the images are very similar in
1889 % spatial layout and color.
1890 %
1891 % The format of the IsImagesEqual method is:
1892 %
1893 % MagickBooleanType IsImagesEqual(Image *image,
1894 % const Image *reconstruct_image)
1895 %
1896 % A description of each parameter follows.
1897 %
1898 % o image: the image.
1899 %
1900 % o reconstruct_image: the reconstruct image.
1901 %
1902 */
1903 MagickExport MagickBooleanType IsImagesEqual(Image *image,
1904  const Image *reconstruct_image)
1905 {
1906  CacheView
1907  *image_view,
1908  *reconstruct_view;
1909 
1911  *exception;
1912 
1913  MagickBooleanType
1914  status;
1915 
1916  MagickRealType
1917  area,
1918  gamma,
1919  maximum_error,
1920  mean_error,
1921  mean_error_per_pixel;
1922 
1923  size_t
1924  columns,
1925  rows;
1926 
1927  ssize_t
1928  y;
1929 
1930  assert(image != (Image *) NULL);
1931  assert(image->signature == MagickCoreSignature);
1932  assert(reconstruct_image != (const Image *) NULL);
1933  assert(reconstruct_image->signature == MagickCoreSignature);
1934  exception=(&image->exception);
1935  if (ValidateImageMorphology(image,reconstruct_image) == MagickFalse)
1936  ThrowBinaryException(ImageError,"ImageMorphologyDiffers",image->filename);
1937  area=0.0;
1938  maximum_error=0.0;
1939  mean_error_per_pixel=0.0;
1940  mean_error=0.0;
1941  rows=MagickMax(image->rows,reconstruct_image->rows);
1942  columns=MagickMax(image->columns,reconstruct_image->columns);
1943  image_view=AcquireVirtualCacheView(image,exception);
1944  reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1945  for (y=0; y < (ssize_t) rows; y++)
1946  {
1947  const IndexPacket
1948  *magick_restrict indexes,
1949  *magick_restrict reconstruct_indexes;
1950 
1951  const PixelPacket
1952  *magick_restrict p,
1953  *magick_restrict q;
1954 
1955  ssize_t
1956  x;
1957 
1958  p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1959  q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1960  if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1961  break;
1962  indexes=GetCacheViewVirtualIndexQueue(image_view);
1963  reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
1964  for (x=0; x < (ssize_t) columns; x++)
1965  {
1966  MagickRealType
1967  distance;
1968 
1969  distance=fabs((double) (GetPixelRed(p)-(double) GetPixelRed(q)));
1970  mean_error_per_pixel+=distance;
1971  mean_error+=distance*distance;
1972  if (distance > maximum_error)
1973  maximum_error=distance;
1974  area++;
1975  distance=fabs((double) (GetPixelGreen(p)-(double) GetPixelGreen(q)));
1976  mean_error_per_pixel+=distance;
1977  mean_error+=distance*distance;
1978  if (distance > maximum_error)
1979  maximum_error=distance;
1980  area++;
1981  distance=fabs((double) (GetPixelBlue(p)-(double) GetPixelBlue(q)));
1982  mean_error_per_pixel+=distance;
1983  mean_error+=distance*distance;
1984  if (distance > maximum_error)
1985  maximum_error=distance;
1986  area++;
1987  if (image->matte != MagickFalse)
1988  {
1989  distance=fabs((double) (GetPixelOpacity(p)-(double)
1990  GetPixelOpacity(q)));
1991  mean_error_per_pixel+=distance;
1992  mean_error+=distance*distance;
1993  if (distance > maximum_error)
1994  maximum_error=distance;
1995  area++;
1996  }
1997  if ((image->colorspace == CMYKColorspace) &&
1998  (reconstruct_image->colorspace == CMYKColorspace))
1999  {
2000  distance=fabs((double) (GetPixelIndex(indexes+x)-(double)
2001  GetPixelIndex(reconstruct_indexes+x)));
2002  mean_error_per_pixel+=distance;
2003  mean_error+=distance*distance;
2004  if (distance > maximum_error)
2005  maximum_error=distance;
2006  area++;
2007  }
2008  p++;
2009  q++;
2010  }
2011  }
2012  reconstruct_view=DestroyCacheView(reconstruct_view);
2013  image_view=DestroyCacheView(image_view);
2014  gamma=PerceptibleReciprocal(area);
2015  image->error.mean_error_per_pixel=gamma*mean_error_per_pixel;
2016  image->error.normalized_mean_error=gamma*QuantumScale*QuantumScale*mean_error;
2017  image->error.normalized_maximum_error=QuantumScale*maximum_error;
2018  status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
2019  return(status);
2020 }
2021 
2022 /*
2023 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2024 % %
2025 % %
2026 % %
2027 % S i m i l a r i t y I m a g e %
2028 % %
2029 % %
2030 % %
2031 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2032 %
2033 % SimilarityImage() compares the reference image of the image and returns the
2034 % best match offset. In addition, it returns a similarity image such that an
2035 % exact match location is completely white and if none of the pixels match,
2036 % black, otherwise some gray level in-between.
2037 %
2038 % The format of the SimilarityImageImage method is:
2039 %
2040 % Image *SimilarityImage(const Image *image,const Image *reference,
2041 % RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
2042 %
2043 % A description of each parameter follows:
2044 %
2045 % o image: the image.
2046 %
2047 % o reference: find an area of the image that closely resembles this image.
2048 %
2049 % o the best match offset of the reference image within the image.
2050 %
2051 % o similarity: the computed similarity between the images.
2052 %
2053 % o exception: return any errors or warnings in this structure.
2054 %
2055 */
2056 
2057 static double GetSimilarityMetric(const Image *image,const Image *reference,
2058  const MetricType metric,const ssize_t x_offset,const ssize_t y_offset,
2059  ExceptionInfo *exception)
2060 {
2061  double
2062  distortion;
2063 
2064  Image
2065  *similarity_image;
2066 
2067  MagickBooleanType
2068  status;
2069 
2071  geometry;
2072 
2073  SetGeometry(reference,&geometry);
2074  geometry.x=x_offset;
2075  geometry.y=y_offset;
2076  similarity_image=CropImage(image,&geometry,exception);
2077  if (similarity_image == (Image *) NULL)
2078  return(0.0);
2079  distortion=0.0;
2080  status=GetImageDistortion(similarity_image,reference,metric,&distortion,
2081  exception);
2082  (void) status;
2083  similarity_image=DestroyImage(similarity_image);
2084  return(distortion);
2085 }
2086 
2087 MagickExport Image *SimilarityImage(Image *image,const Image *reference,
2088  RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
2089 {
2090  Image
2091  *similarity_image;
2092 
2093  similarity_image=SimilarityMetricImage(image,reference,
2094  RootMeanSquaredErrorMetric,offset,similarity_metric,exception);
2095  return(similarity_image);
2096 }
2097 
2098 MagickExport Image *SimilarityMetricImage(Image *image,const Image *reference,
2099  const MetricType metric,RectangleInfo *offset,double *similarity_metric,
2100  ExceptionInfo *exception)
2101 {
2102 #define SimilarityImageTag "Similarity/Image"
2103 
2104  CacheView
2105  *similarity_view;
2106 
2107  const char
2108  *artifact;
2109 
2110  double
2111  similarity_threshold;
2112 
2113  Image
2114  *similarity_image;
2115 
2116  MagickBooleanType
2117  status;
2118 
2119  MagickOffsetType
2120  progress;
2121 
2122  ssize_t
2123  y;
2124 
2125  assert(image != (const Image *) NULL);
2126  assert(image->signature == MagickCoreSignature);
2127  assert(exception != (ExceptionInfo *) NULL);
2128  assert(exception->signature == MagickCoreSignature);
2129  assert(offset != (RectangleInfo *) NULL);
2130  if (IsEventLogging() != MagickFalse)
2131  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2132  SetGeometry(reference,offset);
2133  *similarity_metric=MagickMaximumValue;
2134  if (ValidateImageMorphology(image,reference) == MagickFalse)
2135  ThrowImageException(ImageError,"ImageMorphologyDiffers");
2136  similarity_image=CloneImage(image,image->columns-reference->columns+1,
2137  image->rows-reference->rows+1,MagickTrue,exception);
2138  if (similarity_image == (Image *) NULL)
2139  return((Image *) NULL);
2140  if (SetImageStorageClass(similarity_image,DirectClass) == MagickFalse)
2141  {
2142  InheritException(exception,&similarity_image->exception);
2143  similarity_image=DestroyImage(similarity_image);
2144  return((Image *) NULL);
2145  }
2146  (void) SetImageAlphaChannel(similarity_image,DeactivateAlphaChannel);
2147  /*
2148  Measure similarity of reference image against image.
2149  */
2150  similarity_threshold=(-1.0);
2151  artifact=GetImageArtifact(image,"compare:similarity-threshold");
2152  if (artifact != (const char *) NULL)
2153  similarity_threshold=StringToDouble(artifact,(char **) NULL);
2154  status=MagickTrue;
2155  progress=0;
2156  similarity_view=AcquireVirtualCacheView(similarity_image,exception);
2157 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2158  #pragma omp parallel for schedule(static) \
2159  shared(progress,status,similarity_metric) \
2160  magick_number_threads(image,image,image->rows-reference->rows+1,1)
2161 #endif
2162  for (y=0; y < (ssize_t) (image->rows-reference->rows+1); y++)
2163  {
2164  double
2165  similarity;
2166 
2167  ssize_t
2168  x;
2169 
2170  PixelPacket
2171  *magick_restrict q;
2172 
2173  if (status == MagickFalse)
2174  continue;
2175 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2176  #pragma omp flush(similarity_metric)
2177 #endif
2178  if (*similarity_metric <= similarity_threshold)
2179  continue;
2180  q=GetCacheViewAuthenticPixels(similarity_view,0,y,similarity_image->columns,
2181  1,exception);
2182  if (q == (const PixelPacket *) NULL)
2183  {
2184  status=MagickFalse;
2185  continue;
2186  }
2187  for (x=0; x < (ssize_t) (image->columns-reference->columns+1); x++)
2188  {
2189 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2190  #pragma omp flush(similarity_metric)
2191 #endif
2192  if (*similarity_metric <= similarity_threshold)
2193  break;
2194  similarity=GetSimilarityMetric(image,reference,metric,x,y,exception);
2195 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2196  #pragma omp critical (MagickCore_SimilarityImage)
2197 #endif
2198  if ((metric == NormalizedCrossCorrelationErrorMetric) ||
2199  (metric == UndefinedErrorMetric))
2200  similarity=1.0-similarity;
2201  if (similarity < *similarity_metric)
2202  {
2203  *similarity_metric=similarity;
2204  offset->x=x;
2205  offset->y=y;
2206  }
2207  if (metric == PerceptualHashErrorMetric)
2208  similarity=MagickMin(0.01*similarity,1.0);
2209  SetPixelRed(q,ClampToQuantum(QuantumRange-QuantumRange*similarity));
2210  SetPixelGreen(q,GetPixelRed(q));
2211  SetPixelBlue(q,GetPixelRed(q));
2212  q++;
2213  }
2214  if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
2215  status=MagickFalse;
2216  if (image->progress_monitor != (MagickProgressMonitor) NULL)
2217  {
2218  MagickBooleanType
2219  proceed;
2220 
2221 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2222  #pragma omp atomic
2223 #endif
2224  progress++;
2225  proceed=SetImageProgress(image,SimilarityImageTag,progress,image->rows);
2226  if (proceed == MagickFalse)
2227  status=MagickFalse;
2228  }
2229  }
2230  similarity_view=DestroyCacheView(similarity_view);
2231  if (status == MagickFalse)
2232  similarity_image=DestroyImage(similarity_image);
2233  return(similarity_image);
2234 }
Definition: image.h:152