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