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