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