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