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