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 == MagickCoreSignature);
 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 == MagickCoreSignature);
 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       j,
 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 (j=0; j <= MaxPixelChannels; j++)
 463       distortion[j]+=channel_distortion[j];
 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     j;
 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       x;
 510
 511     if (status == MagickFalse)
 512       continue;
 513     p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
 514     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
 515     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
 516       {
 517         status=MagickFalse;
 518         continue;
 519       }
 520     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
 521     for (x=0; x < (ssize_t) columns; x++)
 522     {
 523       double
 524         Da,
 525         Sa;
 526
 527       register ssize_t
 528         i;
 529
 530       if (GetPixelReadMask(image,p) == 0)
 531         {
 532           p+=GetPixelChannels(image);
 533           q+=GetPixelChannels(reconstruct_image);
 534           continue;
 535         }
 536       Sa=QuantumScale*GetPixelAlpha(image,p);
 537       Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
 538       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 539       {
 540         double
 541           distance;
 542
 543         PixelChannel channel=GetPixelChannelChannel(image,i);
 544         PixelTrait traits=GetPixelChannelTraits(image,channel);
 545         PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
 546           channel);
 547         if ((traits == UndefinedPixelTrait) ||
 548             (reconstruct_traits == UndefinedPixelTrait) ||
 549             ((reconstruct_traits & UpdatePixelTrait) == 0))
 550           continue;
 551         distance=QuantumScale*(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
 552           channel,q));
 553         channel_distortion[i]+=distance*distance;
 554         channel_distortion[CompositePixelChannel]+=distance*distance;
 555       }
 556       p+=GetPixelChannels(image);
 557       q+=GetPixelChannels(reconstruct_image);
 558     }
 559 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 560     #pragma omp critical (MagickCore_GetFuzzDistortion)
 561 #endif
 562     for (j=0; j <= MaxPixelChannels; j++)
 563       distortion[j]+=channel_distortion[j];
 564   }
 565   reconstruct_view=DestroyCacheView(reconstruct_view);
 566   image_view=DestroyCacheView(image_view);
 567   for (j=0; j <= MaxPixelChannels; j++)
 568     distortion[j]/=((double) columns*rows);
 569   distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
 570   distortion[CompositePixelChannel]=sqrt(distortion[CompositePixelChannel]);
 571   return(status);
 572 }
 573
 574 static MagickBooleanType GetMeanAbsoluteDistortion(const Image *image,
 575   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
 576 {
 577   CacheView
 578     *image_view,
 579     *reconstruct_view;
 580
 581   MagickBooleanType
 582     status;
 583
 584   register ssize_t
 585     j;
 586
 587   size_t
 588     columns,
 589     rows;
 590
 591   ssize_t
 592     y;
 593
 594   status=MagickTrue;
 595   rows=MagickMax(image->rows,reconstruct_image->rows);
 596   columns=MagickMax(image->columns,reconstruct_image->columns);
 597   image_view=AcquireVirtualCacheView(image,exception);
 598   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
 599 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 600   #pragma omp parallel for schedule(static,4) shared(status) \
 601     magick_threads(image,image,rows,1)
 602 #endif
 603   for (y=0; y < (ssize_t) rows; y++)
 604   {
 605     double
 606       channel_distortion[MaxPixelChannels+1];
 607
 608     register const Quantum
 609       *restrict p,
 610       *restrict q;
 611
 612     register ssize_t
 613       x;
 614
 615     if (status == MagickFalse)
 616       continue;
 617     p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
 618     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
 619     if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
 620       {
 621         status=MagickFalse;
 622         continue;
 623       }
 624     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
 625     for (x=0; x < (ssize_t) columns; x++)
 626     {
 627       double
 628         Da,
 629         Sa;
 630
 631       register ssize_t
 632         i;
 633
 634       if (GetPixelReadMask(image,p) == 0)
 635         {
 636           p+=GetPixelChannels(image);
 637           q+=GetPixelChannels(reconstruct_image);
 638           continue;
 639         }
 640       Sa=QuantumScale*GetPixelAlpha(image,p);
 641       Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
 642       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 643       {
 644         double
 645           distance;
 646
 647         PixelChannel channel=GetPixelChannelChannel(image,i);
 648         PixelTrait traits=GetPixelChannelTraits(image,channel);
 649         PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
 650           channel);
 651         if ((traits == UndefinedPixelTrait) ||
 652             (reconstruct_traits == UndefinedPixelTrait) ||
 653             ((reconstruct_traits & UpdatePixelTrait) == 0))
 654           continue;
 655         distance=QuantumScale*fabs(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
 656           channel,q));
 657         channel_distortion[i]+=distance;
 658         channel_distortion[CompositePixelChannel]+=distance;
 659       }
 660       p+=GetPixelChannels(image);
 661       q+=GetPixelChannels(reconstruct_image);
 662     }
 663 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 664     #pragma omp critical (MagickCore_GetMeanAbsoluteError)
 665 #endif
 666     for (j=0; j <= MaxPixelChannels; j++)
 667       distortion[j]+=channel_distortion[j];
 668   }
 669   reconstruct_view=DestroyCacheView(reconstruct_view);
 670   image_view=DestroyCacheView(image_view);
 671   for (j=0; j <= MaxPixelChannels; j++)
 672     distortion[j]/=((double) columns*rows);
 673   distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
 674   return(status);
 675 }
 676
 677 static MagickBooleanType GetMeanErrorPerPixel(Image *image,
 678   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
 679 {
 680   CacheView
 681     *image_view,
 682     *reconstruct_view;
 683
 684   MagickBooleanType
 685     status;
 686
 687   double
 688     area,
 689     maximum_error,
 690     mean_error;
 691
 692   size_t
 693     columns,
 694     rows;
 695
 696   ssize_t
 697     y;
 698
 699   status=MagickTrue;
 700   area=0.0;
 701   maximum_error=0.0;
 702   mean_error=0.0;
 703   rows=MagickMax(image->rows,reconstruct_image->rows);
 704   columns=MagickMax(image->columns,reconstruct_image->columns);
 705   image_view=AcquireVirtualCacheView(image,exception);
 706   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
 707   for (y=0; y < (ssize_t) rows; y++)
 708   {
 709     register const Quantum
 710       *restrict p,
 711       *restrict q;
 712
 713     register ssize_t
 714       x;
 715
 716     p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
 717     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
 718     if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
 719       {
 720         status=MagickFalse;
 721         break;
 722       }
 723     for (x=0; x < (ssize_t) columns; x++)
 724     {
 725       double
 726         Da,
 727         Sa;
 728
 729       register ssize_t
 730         i;
 731
 732       if (GetPixelReadMask(image,p) == 0)
 733         {
 734           p+=GetPixelChannels(image);
 735           q+=GetPixelChannels(reconstruct_image);
 736           continue;
 737         }
 738       Sa=QuantumScale*GetPixelAlpha(image,p);
 739       Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
 740       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 741       {
 742         double
 743           distance;
 744
 745         PixelChannel channel=GetPixelChannelChannel(image,i);
 746         PixelTrait traits=GetPixelChannelTraits(image,channel);
 747         PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
 748           channel);
 749         if ((traits == UndefinedPixelTrait) ||
 750             (reconstruct_traits == UndefinedPixelTrait) ||
 751             ((reconstruct_traits & UpdatePixelTrait) == 0))
 752           continue;
 753         distance=fabs(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q));
 754         distortion[i]+=distance;
 755         distortion[CompositePixelChannel]+=distance;
 756         mean_error+=distance*distance;
 757         if (distance > maximum_error)
 758           maximum_error=distance;
 759         area++;
 760       }
 761       p+=GetPixelChannels(image);
 762       q+=GetPixelChannels(reconstruct_image);
 763     }
 764   }
 765   reconstruct_view=DestroyCacheView(reconstruct_view);
 766   image_view=DestroyCacheView(image_view);
 767   image->error.mean_error_per_pixel=distortion[CompositePixelChannel]/area;
 768   image->error.normalized_mean_error=QuantumScale*QuantumScale*mean_error/area;
 769   image->error.normalized_maximum_error=QuantumScale*maximum_error;
 770   return(status);
 771 }
 772
 773 static MagickBooleanType GetMeanSquaredDistortion(const Image *image,
 774   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
 775 {
 776   CacheView
 777     *image_view,
 778     *reconstruct_view;
 779
 780   MagickBooleanType
 781     status;
 782
 783   register ssize_t
 784     j;
 785
 786   size_t
 787     columns,
 788     rows;
 789
 790   ssize_t
 791     y;
 792
 793   status=MagickTrue;
 794   rows=MagickMax(image->rows,reconstruct_image->rows);
 795   columns=MagickMax(image->columns,reconstruct_image->columns);
 796   image_view=AcquireVirtualCacheView(image,exception);
 797   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
 798 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 799   #pragma omp parallel for schedule(static,4) shared(status) \
 800     magick_threads(image,image,rows,1)
 801 #endif
 802   for (y=0; y < (ssize_t) rows; y++)
 803   {
 804     double
 805       channel_distortion[MaxPixelChannels+1];
 806
 807     register const Quantum
 808       *restrict p,
 809       *restrict q;
 810
 811     register ssize_t
 812       x;
 813
 814     if (status == MagickFalse)
 815       continue;
 816     p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
 817     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
 818     if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
 819       {
 820         status=MagickFalse;
 821         continue;
 822       }
 823     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
 824     for (x=0; x < (ssize_t) columns; x++)
 825     {
 826       double
 827         Da,
 828         Sa;
 829
 830       register ssize_t
 831         i;
 832
 833       if (GetPixelReadMask(image,p) == 0)
 834         {
 835           p+=GetPixelChannels(image);
 836           q+=GetPixelChannels(reconstruct_image);
 837           continue;
 838         }
 839       Sa=QuantumScale*GetPixelAlpha(image,p);
 840       Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
 841       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 842       {
 843         double
 844           distance;
 845
 846         PixelChannel channel=GetPixelChannelChannel(image,i);
 847         PixelTrait traits=GetPixelChannelTraits(image,channel);
 848         PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
 849           channel);
 850         if ((traits == UndefinedPixelTrait) ||
 851             (reconstruct_traits == UndefinedPixelTrait) ||
 852             ((reconstruct_traits & UpdatePixelTrait) == 0))
 853           continue;
 854         distance=QuantumScale*(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
 855           channel,q));
 856         channel_distortion[i]+=distance*distance;
 857         channel_distortion[CompositePixelChannel]+=distance*distance;
 858       }
 859       p+=GetPixelChannels(image);
 860       q+=GetPixelChannels(reconstruct_image);
 861     }
 862 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 863     #pragma omp critical (MagickCore_GetMeanSquaredError)
 864 #endif
 865     for (j=0; j <= MaxPixelChannels; j++)
 866       distortion[j]+=channel_distortion[j];
 867   }
 868   reconstruct_view=DestroyCacheView(reconstruct_view);
 869   image_view=DestroyCacheView(image_view);
 870   for (j=0; j <= MaxPixelChannels; j++)
 871     distortion[j]/=((double) columns*rows);
 872   distortion[CompositePixelChannel]/=GetImageChannels(image);
 873   return(status);
 874 }
 875
 876 static MagickBooleanType GetNormalizedCrossCorrelationDistortion(
 877   const Image *image,const Image *reconstruct_image,double *distortion,
 878   ExceptionInfo *exception)
 879 {
 880 #define SimilarityImageTag  "Similarity/Image"
 881
 882   CacheView
 883     *image_view,
 884     *reconstruct_view;
 885
 886   ChannelStatistics
 887     *image_statistics,
 888     *reconstruct_statistics;
 889
 890   double
 891     area;
 892
 893   MagickBooleanType
 894     status;
 895
 896   MagickOffsetType
 897     progress;
 898
 899   register ssize_t
 900     i;
 901
 902   size_t
 903     columns,
 904     rows;
 905
 906   ssize_t
 907     y;
 908
 909   /*
 910     Normalize to account for variation due to lighting and exposure condition.
 911   */
 912   image_statistics=GetImageStatistics(image,exception);
 913   reconstruct_statistics=GetImageStatistics(reconstruct_image,exception);
 914   if ((image_statistics == (ChannelStatistics *) NULL) ||
 915       (reconstruct_statistics == (ChannelStatistics *) NULL))
 916     {
 917       if (image_statistics != (ChannelStatistics *) NULL)
 918         image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
 919           image_statistics);
 920       if (reconstruct_statistics != (ChannelStatistics *) NULL)
 921         reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
 922           reconstruct_statistics);
 923       return(MagickFalse);
 924     }
 925   status=MagickTrue;
 926   progress=0;
 927   for (i=0; i <= MaxPixelChannels; i++)
 928     distortion[i]=0.0;
 929   rows=MagickMax(image->rows,reconstruct_image->rows);
 930   columns=MagickMax(image->columns,reconstruct_image->columns);
 931   area=1.0/((double) columns*rows-1);
 932   image_view=AcquireVirtualCacheView(image,exception);
 933   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
 934   for (y=0; y < (ssize_t) rows; y++)
 935   {
 936     register const Quantum
 937       *restrict p,
 938       *restrict q;
 939
 940     register ssize_t
 941       x;
 942
 943     p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
 944     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
 945     if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
 946       {
 947         status=MagickFalse;
 948         break;
 949       }
 950     for (x=0; x < (ssize_t) columns; x++)
 951     {
 952       double
 953         Da,
 954         Sa;
 955
 956       if (GetPixelReadMask(image,p) == 0)
 957         {
 958           p+=GetPixelChannels(image);
 959           q+=GetPixelChannels(reconstruct_image);
 960           continue;
 961         }
 962       Sa=QuantumScale*(image->alpha_trait != UndefinedPixelTrait ?
 963         GetPixelAlpha(image,p) : OpaqueAlpha);
 964       Da=QuantumScale*(reconstruct_image->alpha_trait != UndefinedPixelTrait ?
 965         GetPixelAlpha(reconstruct_image,q) : OpaqueAlpha);
 966       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 967       {
 968         PixelChannel channel=GetPixelChannelChannel(image,i);
 969         PixelTrait traits=GetPixelChannelTraits(image,channel);
 970         PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
 971           channel);
 972         if ((traits == UndefinedPixelTrait) ||
 973             (reconstruct_traits == UndefinedPixelTrait) ||
 974             ((reconstruct_traits & UpdatePixelTrait) == 0))
 975           continue;
 976         if (channel == AlphaPixelChannel)
 977           {
 978             distortion[i]+=area*QuantumScale*(p[i]-
 979               image_statistics[channel].mean)*(GetPixelChannel(
 980               reconstruct_image,channel,q)-
 981               reconstruct_statistics[channel].mean);
 982           }
 983         else
 984           {
 985             distortion[i]+=area*QuantumScale*(Sa*p[i]-
 986               image_statistics[channel].mean)*(Da*GetPixelChannel(
 987               reconstruct_image,channel,q)-
 988               reconstruct_statistics[channel].mean);
 989           }
 990       }
 991       p+=GetPixelChannels(image);
 992       q+=GetPixelChannels(reconstruct_image);
 993     }
 994     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 995       {
 996         MagickBooleanType
 997           proceed;
 998
 999         proceed=SetImageProgress(image,SimilarityImageTag,progress++,rows);
1000         if (proceed == MagickFalse)
1001           {
1002             status=MagickFalse;
1003             break;
1004           }
1005       }
1006   }
1007   reconstruct_view=DestroyCacheView(reconstruct_view);
1008   image_view=DestroyCacheView(image_view);
1009   /*
1010     Divide by the standard deviation.
1011   */
1012   distortion[CompositePixelChannel]=0.0;
1013   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1014   {
1015     double
1016       gamma;
1017
1018     PixelChannel channel=GetPixelChannelChannel(image,i);
1019     gamma=image_statistics[channel].standard_deviation*
1020       reconstruct_statistics[channel].standard_deviation;
1021     gamma=PerceptibleReciprocal(gamma);
1022     distortion[i]=QuantumRange*gamma*distortion[i];
1023     distortion[CompositePixelChannel]+=distortion[i]*distortion[i];
1024   }
1025   distortion[CompositePixelChannel]=sqrt(distortion[CompositePixelChannel]/
1026     GetImageChannels(image));
1027   /*
1028     Free resources.
1029   */
1030   reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1031     reconstruct_statistics);
1032   image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1033     image_statistics);
1034   return(status);
1035 }
1036
1037 static MagickBooleanType GetPeakAbsoluteDistortion(const Image *image,
1038   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1039 {
1040   CacheView
1041     *image_view,
1042     *reconstruct_view;
1043
1044   MagickBooleanType
1045     status;
1046
1047   size_t
1048     columns,
1049     rows;
1050
1051   ssize_t
1052     y;
1053
1054   status=MagickTrue;
1055   rows=MagickMax(image->rows,reconstruct_image->rows);
1056   columns=MagickMax(image->columns,reconstruct_image->columns);
1057   image_view=AcquireVirtualCacheView(image,exception);
1058   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1059 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1060   #pragma omp parallel for schedule(static,4) shared(status) \
1061     magick_threads(image,image,rows,1)
1062 #endif
1063   for (y=0; y < (ssize_t) rows; y++)
1064   {
1065     double
1066       channel_distortion[MaxPixelChannels+1];
1067
1068     register const Quantum
1069       *restrict p,
1070       *restrict q;
1071
1072     register ssize_t
1073       j,
1074       x;
1075
1076     if (status == MagickFalse)
1077       continue;
1078     p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1079     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1080     if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1081       {
1082         status=MagickFalse;
1083         continue;
1084       }
1085     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
1086     for (x=0; x < (ssize_t) columns; x++)
1087     {
1088       double
1089         Da,
1090         Sa;
1091
1092       register ssize_t
1093         i;
1094
1095       if (GetPixelReadMask(image,p) == 0)
1096         {
1097           p+=GetPixelChannels(image);
1098           q+=GetPixelChannels(reconstruct_image);
1099           continue;
1100         }
1101       Sa=QuantumScale*GetPixelAlpha(image,p);
1102       Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
1103       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1104       {
1105         double
1106           distance;
1107
1108         PixelChannel channel=GetPixelChannelChannel(image,i);
1109         PixelTrait traits=GetPixelChannelTraits(image,channel);
1110         PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
1111           channel);
1112         if ((traits == UndefinedPixelTrait) ||
1113             (reconstruct_traits == UndefinedPixelTrait) ||
1114             ((reconstruct_traits & UpdatePixelTrait) == 0))
1115           continue;
1116         distance=QuantumScale*fabs(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
1117           channel,q));
1118         if (distance > channel_distortion[i])
1119           channel_distortion[i]=distance;
1120         if (distance > channel_distortion[CompositePixelChannel])
1121           channel_distortion[CompositePixelChannel]=distance;
1122       }
1123       p+=GetPixelChannels(image);
1124       q+=GetPixelChannels(reconstruct_image);
1125     }
1126 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1127     #pragma omp critical (MagickCore_GetPeakAbsoluteError)
1128 #endif
1129     for (j=0; j <= MaxPixelChannels; j++)
1130       if (channel_distortion[j] > distortion[j])
1131         distortion[j]=channel_distortion[j];
1132   }
1133   reconstruct_view=DestroyCacheView(reconstruct_view);
1134   image_view=DestroyCacheView(image_view);
1135   return(status);
1136 }
1137
1138 static inline double MagickLog10(const double x)
1139 {
1140 #define Log10Epsilon  (1.0e-11)
1141
1142  if (fabs(x) < Log10Epsilon)
1143    return(log10(Log10Epsilon));
1144  return(log10(fabs(x)));
1145 }
1146
1147 static MagickBooleanType GetPeakSignalToNoiseRatio(const Image *image,
1148   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1149 {
1150   MagickBooleanType
1151     status;
1152
1153   register ssize_t
1154     i;
1155
1156   status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1157   for (i=0; i <= MaxPixelChannels; i++)
1158     distortion[i]=20.0*MagickLog10((double) 1.0/sqrt(distortion[i]));
1159   return(status);
1160 }
1161
1162 static MagickBooleanType GetPerceptualHashDistortion(const Image *image,
1163   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1164 {
1165   ChannelPerceptualHash
1166     *image_phash,
1167     *reconstruct_phash;
1168
1169   ssize_t
1170     channel;
1171
1172   /*
1173     Compute perceptual hash in the sRGB colorspace.
1174   */
1175   image_phash=GetImagePerceptualHash(image,exception);
1176   if (image_phash == (ChannelPerceptualHash *) NULL)
1177     return(MagickFalse);
1178   reconstruct_phash=GetImagePerceptualHash(reconstruct_image,exception);
1179   if (reconstruct_phash == (ChannelPerceptualHash *) NULL)
1180     {
1181       image_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(image_phash);
1182       return(MagickFalse);
1183     }
1184 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1185   #pragma omp parallel for schedule(static,4)
1186 #endif
1187   for (channel=0; channel < MaxPixelChannels; channel++)
1188   {
1189     double
1190       difference;
1191
1192     register ssize_t
1193       i;
1194
1195     difference=0.0;
1196     for (i=0; i < MaximumNumberOfImageMoments; i++)
1197     {
1198       double
1199         alpha,
1200         beta;
1201
1202       alpha=image_phash[channel].srgb_hu_phash[i];
1203       beta=reconstruct_phash[channel].srgb_hu_phash[i];
1204       difference+=(beta-alpha)*(beta-alpha);
1205     }
1206     distortion[channel]+=difference;
1207 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1208     #pragma omp critical (MagickCore_GetPerceptualHashDistortion)
1209 #endif
1210     distortion[CompositePixelChannel]+=difference;
1211   }
1212   /*
1213     Compute perceptual hash in the HCLP colorspace.
1214   */
1215 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1216   #pragma omp parallel for schedule(static,4)
1217 #endif
1218   for (channel=0; channel < MaxPixelChannels; channel++)
1219   {
1220     double
1221       difference;
1222
1223     register ssize_t
1224       i;
1225
1226     difference=0.0;
1227     for (i=0; i < MaximumNumberOfImageMoments; i++)
1228     {
1229       double
1230         alpha,
1231         beta;
1232
1233       alpha=image_phash[channel].hclp_hu_phash[i];
1234       beta=reconstruct_phash[channel].hclp_hu_phash[i];
1235       difference+=(beta-alpha)*(beta-alpha);
1236     }
1237     distortion[channel]+=difference;
1238 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1239     #pragma omp critical (MagickCore_GetPerceptualHashDistortion)
1240 #endif
1241     distortion[CompositePixelChannel]+=difference;
1242   }
1243   /*
1244     Free resources.
1245   */
1246   reconstruct_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1247     reconstruct_phash);
1248   image_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(image_phash);
1249   return(MagickTrue);
1250 }
1251
1252 static MagickBooleanType GetRootMeanSquaredDistortion(const Image *image,
1253   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1254 {
1255   MagickBooleanType
1256     status;
1257
1258   register ssize_t
1259     i;
1260
1261   status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1262   for (i=0; i <= MaxPixelChannels; i++)
1263     distortion[i]=sqrt(distortion[i]);
1264   return(status);
1265 }
1266
1267 MagickExport MagickBooleanType GetImageDistortion(Image *image,
1268   const Image *reconstruct_image,const MetricType metric,double *distortion,
1269   ExceptionInfo *exception)
1270 {
1271   double
1272     *channel_distortion;
1273
1274   MagickBooleanType
1275     status;
1276
1277   size_t
1278     length;
1279
1280   assert(image != (Image *) NULL);
1281   assert(image->signature == MagickCoreSignature);
1282   if (image->debug != MagickFalse)
1283     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1284   assert(reconstruct_image != (const Image *) NULL);
1285   assert(reconstruct_image->signature == MagickCoreSignature);
1286   assert(distortion != (double *) NULL);
1287   *distortion=0.0;
1288   if (image->debug != MagickFalse)
1289     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1290   if (metric != PerceptualHashErrorMetric)
1291     if (ValidateImageMorphology(image,reconstruct_image) == MagickFalse)
1292       ThrowBinaryException(ImageError,"ImageMorphologyDiffers",image->filename);
1293   /*
1294     Get image distortion.
1295   */
1296   length=MaxPixelChannels+1;
1297   channel_distortion=(double *) AcquireQuantumMemory(length,
1298     sizeof(*channel_distortion));
1299   if (channel_distortion == (double *) NULL)
1300     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1301   (void) ResetMagickMemory(channel_distortion,0,length*
1302     sizeof(*channel_distortion));
1303   switch (metric)
1304   {
1305     case AbsoluteErrorMetric:
1306     {
1307       status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1308         exception);
1309       break;
1310     }
1311     case FuzzErrorMetric:
1312     {
1313       status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1314         exception);
1315       break;
1316     }
1317     case MeanAbsoluteErrorMetric:
1318     {
1319       status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1320         channel_distortion,exception);
1321       break;
1322     }
1323     case MeanErrorPerPixelErrorMetric:
1324     {
1325       status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1326         exception);
1327       break;
1328     }
1329     case MeanSquaredErrorMetric:
1330     {
1331       status=GetMeanSquaredDistortion(image,reconstruct_image,
1332         channel_distortion,exception);
1333       break;
1334     }
1335     case NormalizedCrossCorrelationErrorMetric:
1336     default:
1337     {
1338       status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1339         channel_distortion,exception);
1340       break;
1341     }
1342     case PeakAbsoluteErrorMetric:
1343     {
1344       status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1345         channel_distortion,exception);
1346       break;
1347     }
1348     case PeakSignalToNoiseRatioErrorMetric:
1349     {
1350       status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1351         channel_distortion,exception);
1352       break;
1353     }
1354     case PerceptualHashErrorMetric:
1355     {
1356       status=GetPerceptualHashDistortion(image,reconstruct_image,
1357         channel_distortion,exception);
1358       break;
1359     }
1360     case RootMeanSquaredErrorMetric:
1361     {
1362       status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1363         channel_distortion,exception);
1364       break;
1365     }
1366   }
1367   *distortion=channel_distortion[CompositePixelChannel];
1368   channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1369   (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
1370     *distortion);
1371   return(status);
1372 }
1373 \f
1374 /*
1375 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1376 %                                                                             %
1377 %                                                                             %
1378 %                                                                             %
1379 %   G e t I m a g e D i s t o r t i o n s                                     %
1380 %                                                                             %
1381 %                                                                             %
1382 %                                                                             %
1383 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1384 %
1385 %  GetImageDistortions() compares the pixel channels of an image to a
1386 %  reconstructed image and returns the specified distortion metric for each
1387 %  channel.
1388 %
1389 %  The format of the GetImageDistortions method is:
1390 %
1391 %      double *GetImageDistortions(const Image *image,
1392 %        const Image *reconstruct_image,const MetricType metric,
1393 %        ExceptionInfo *exception)
1394 %
1395 %  A description of each parameter follows:
1396 %
1397 %    o image: the image.
1398 %
1399 %    o reconstruct_image: the reconstruct image.
1400 %
1401 %    o metric: the metric.
1402 %
1403 %    o exception: return any errors or warnings in this structure.
1404 %
1405 */
1406 MagickExport double *GetImageDistortions(Image *image,
1407   const Image *reconstruct_image,const MetricType metric,
1408   ExceptionInfo *exception)
1409 {
1410   double
1411     *channel_distortion;
1412
1413   MagickBooleanType
1414     status;
1415
1416   size_t
1417     length;
1418
1419   assert(image != (Image *) NULL);
1420   assert(image->signature == MagickCoreSignature);
1421   if (image->debug != MagickFalse)
1422     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1423   assert(reconstruct_image != (const Image *) NULL);
1424   assert(reconstruct_image->signature == MagickCoreSignature);
1425   if (image->debug != MagickFalse)
1426     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1427   if (metric != PerceptualHashErrorMetric)
1428     if (ValidateImageMorphology(image,reconstruct_image) == MagickFalse)
1429       {
1430         (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1431           "ImageMorphologyDiffers","`%s'",image->filename);
1432         return((double *) NULL);
1433       }
1434   /*
1435     Get image distortion.
1436   */
1437   length=MaxPixelChannels+1UL;
1438   channel_distortion=(double *) AcquireQuantumMemory(length,
1439     sizeof(*channel_distortion));
1440   if (channel_distortion == (double *) NULL)
1441     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1442   (void) ResetMagickMemory(channel_distortion,0,length*
1443     sizeof(*channel_distortion));
1444   status=MagickTrue;
1445   switch (metric)
1446   {
1447     case AbsoluteErrorMetric:
1448     {
1449       status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1450         exception);
1451       break;
1452     }
1453     case FuzzErrorMetric:
1454     {
1455       status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1456         exception);
1457       break;
1458     }
1459     case MeanAbsoluteErrorMetric:
1460     {
1461       status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1462         channel_distortion,exception);
1463       break;
1464     }
1465     case MeanErrorPerPixelErrorMetric:
1466     {
1467       status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1468         exception);
1469       break;
1470     }
1471     case MeanSquaredErrorMetric:
1472     {
1473       status=GetMeanSquaredDistortion(image,reconstruct_image,
1474         channel_distortion,exception);
1475       break;
1476     }
1477     case NormalizedCrossCorrelationErrorMetric:
1478     default:
1479     {
1480       status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1481         channel_distortion,exception);
1482       break;
1483     }
1484     case PeakAbsoluteErrorMetric:
1485     {
1486       status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1487         channel_distortion,exception);
1488       break;
1489     }
1490     case PeakSignalToNoiseRatioErrorMetric:
1491     {
1492       status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1493         channel_distortion,exception);
1494       break;
1495     }
1496     case PerceptualHashErrorMetric:
1497     {
1498       status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1499         channel_distortion,exception);
1500       break;
1501     }
1502     case RootMeanSquaredErrorMetric:
1503     {
1504       status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1505         channel_distortion,exception);
1506       break;
1507     }
1508   }
1509   if (status == MagickFalse)
1510     {
1511       channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1512       return((double *) NULL);
1513     }
1514   return(channel_distortion);
1515 }
1516 \f
1517 /*
1518 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1519 %                                                                             %
1520 %                                                                             %
1521 %                                                                             %
1522 %  I s I m a g e s E q u a l                                                  %
1523 %                                                                             %
1524 %                                                                             %
1525 %                                                                             %
1526 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1527 %
1528 %  IsImagesEqual() measures the difference between colors at each pixel
1529 %  location of two images.  A value other than 0 means the colors match
1530 %  exactly.  Otherwise an error measure is computed by summing over all
1531 %  pixels in an image the distance squared in RGB space between each image
1532 %  pixel and its corresponding pixel in the reconstruct image.  The error
1533 %  measure is assigned to these image members:
1534 %
1535 %    o mean_error_per_pixel:  The mean error for any single pixel in
1536 %      the image.
1537 %
1538 %    o normalized_mean_error:  The normalized mean quantization error for
1539 %      any single pixel in the image.  This distance measure is normalized to
1540 %      a range between 0 and 1.  It is independent of the range of red, green,
1541 %      and blue values in the image.
1542 %
1543 %    o normalized_maximum_error:  The normalized maximum quantization
1544 %      error for any single pixel in the image.  This distance measure is
1545 %      normalized to a range between 0 and 1.  It is independent of the range
1546 %      of red, green, and blue values in your image.
1547 %
1548 %  A small normalized mean square error, accessed as
1549 %  image->normalized_mean_error, suggests the images are very similar in
1550 %  spatial layout and color.
1551 %
1552 %  The format of the IsImagesEqual method is:
1553 %
1554 %      MagickBooleanType IsImagesEqual(Image *image,
1555 %        const Image *reconstruct_image,ExceptionInfo *exception)
1556 %
1557 %  A description of each parameter follows.
1558 %
1559 %    o image: the image.
1560 %
1561 %    o reconstruct_image: the reconstruct image.
1562 %
1563 %    o exception: return any errors or warnings in this structure.
1564 %
1565 */
1566 MagickExport MagickBooleanType IsImagesEqual(Image *image,
1567   const Image *reconstruct_image,ExceptionInfo *exception)
1568 {
1569   CacheView
1570     *image_view,
1571     *reconstruct_view;
1572
1573   MagickBooleanType
1574     status;
1575
1576   double
1577     area,
1578     maximum_error,
1579     mean_error,
1580     mean_error_per_pixel;
1581
1582   size_t
1583     columns,
1584     rows;
1585
1586   ssize_t
1587     y;
1588
1589   assert(image != (Image *) NULL);
1590   assert(image->signature == MagickCoreSignature);
1591   assert(reconstruct_image != (const Image *) NULL);
1592   assert(reconstruct_image->signature == MagickCoreSignature);
1593   if (ValidateImageMorphology(image,reconstruct_image) == MagickFalse)
1594     ThrowBinaryException(ImageError,"ImageMorphologyDiffers",image->filename);
1595   area=0.0;
1596   maximum_error=0.0;
1597   mean_error_per_pixel=0.0;
1598   mean_error=0.0;
1599   rows=MagickMax(image->rows,reconstruct_image->rows);
1600   columns=MagickMax(image->columns,reconstruct_image->columns);
1601   image_view=AcquireVirtualCacheView(image,exception);
1602   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1603   for (y=0; y < (ssize_t) rows; y++)
1604   {
1605     register const Quantum
1606       *restrict p,
1607       *restrict q;
1608
1609     register ssize_t
1610       x;
1611
1612     p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1613     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1614     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
1615       break;
1616     for (x=0; x < (ssize_t) columns; x++)
1617     {
1618       register ssize_t
1619         i;
1620
1621       if (GetPixelReadMask(image,p) == 0)
1622         {
1623           p+=GetPixelChannels(image);
1624           q+=GetPixelChannels(reconstruct_image);
1625           continue;
1626         }
1627       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1628       {
1629         double
1630           distance;
1631
1632         PixelChannel channel=GetPixelChannelChannel(image,i);
1633         PixelTrait traits=GetPixelChannelTraits(image,channel);
1634         PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
1635           channel);
1636         if ((traits == UndefinedPixelTrait) ||
1637             (reconstruct_traits == UndefinedPixelTrait) ||
1638             ((reconstruct_traits & UpdatePixelTrait) == 0))
1639           continue;
1640         distance=fabs(p[i]-(double) GetPixelChannel(reconstruct_image,
1641           channel,q));
1642         if (distance >= MagickEpsilon)
1643           {
1644             mean_error_per_pixel+=distance;
1645             mean_error+=distance*distance;
1646             if (distance > maximum_error)
1647               maximum_error=distance;
1648           }
1649         area++;
1650       }
1651       p+=GetPixelChannels(image);
1652       q+=GetPixelChannels(reconstruct_image);
1653     }
1654   }
1655   reconstruct_view=DestroyCacheView(reconstruct_view);
1656   image_view=DestroyCacheView(image_view);
1657   image->error.mean_error_per_pixel=(double) (mean_error_per_pixel/area);
1658   image->error.normalized_mean_error=(double) (QuantumScale*QuantumScale*
1659     mean_error/area);
1660   image->error.normalized_maximum_error=(double) (QuantumScale*maximum_error);
1661   status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
1662   return(status);
1663 }
1664 \f
1665 /*
1666 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1667 %                                                                             %
1668 %                                                                             %
1669 %                                                                             %
1670 %   S i m i l a r i t y I m a g e                                             %
1671 %                                                                             %
1672 %                                                                             %
1673 %                                                                             %
1674 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1675 %
1676 %  SimilarityImage() compares the reference image of the image and returns the
1677 %  best match offset.  In addition, it returns a similarity image such that an
1678 %  exact match location is completely white and if none of the pixels match,
1679 %  black, otherwise some gray level in-between.
1680 %
1681 %  The format of the SimilarityImageImage method is:
1682 %
1683 %      Image *SimilarityImage(const Image *image,const Image *reference,
1684 %        const MetricType metric,const double similarity_threshold,
1685 %        RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
1686 %
1687 %  A description of each parameter follows:
1688 %
1689 %    o image: the image.
1690 %
1691 %    o reference: find an area of the image that closely resembles this image.
1692 %
1693 %    o metric: the metric.
1694 %
1695 %    o similarity_threshold: minimum distortion for (sub)image match.
1696 %
1697 %    o offset: the best match offset of the reference image within the image.
1698 %
1699 %    o similarity: the computed similarity between the images.
1700 %
1701 %    o exception: return any errors or warnings in this structure.
1702 %
1703 */
1704
1705 static double GetSimilarityMetric(const Image *image,const Image *reference,
1706   const MetricType metric,const ssize_t x_offset,const ssize_t y_offset,
1707   ExceptionInfo *exception)
1708 {
1709   double
1710     distortion;
1711
1712   Image
1713     *similarity_image;
1714
1715   MagickBooleanType
1716     status;
1717
1718   RectangleInfo
1719     geometry;
1720
1721   SetGeometry(reference,&geometry);
1722   geometry.x=x_offset;
1723   geometry.y=y_offset;
1724   similarity_image=CropImage(image,&geometry,exception);
1725   if (similarity_image == (Image *) NULL)
1726     return(0.0);
1727   distortion=0.0;
1728   status=GetImageDistortion(similarity_image,reference,metric,&distortion,
1729     exception);
1730   similarity_image=DestroyImage(similarity_image);
1731   if (status == MagickFalse)
1732     return(0.0);
1733   return(distortion);
1734 }
1735
1736 MagickExport Image *SimilarityImage(Image *image,const Image *reference,
1737   const MetricType metric,const double similarity_threshold,
1738   RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
1739 {
1740 #define SimilarityImageTag  "Similarity/Image"
1741
1742   CacheView
1743     *similarity_view;
1744
1745   Image
1746     *similarity_image;
1747
1748   MagickBooleanType
1749     status;
1750
1751   MagickOffsetType
1752     progress;
1753
1754   ssize_t
1755     y;
1756
1757   assert(image != (const Image *) NULL);
1758   assert(image->signature == MagickCoreSignature);
1759   if (image->debug != MagickFalse)
1760     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1761   assert(exception != (ExceptionInfo *) NULL);
1762   assert(exception->signature == MagickCoreSignature);
1763   assert(offset != (RectangleInfo *) NULL);
1764   SetGeometry(reference,offset);
1765   *similarity_metric=MagickMaximumValue;
1766   if (ValidateImageMorphology(image,reference) == MagickFalse)
1767     ThrowImageException(ImageError,"ImageMorphologyDiffers");
1768   similarity_image=CloneImage(image,image->columns-reference->columns+1,
1769     image->rows-reference->rows+1,MagickTrue,exception);
1770   if (similarity_image == (Image *) NULL)
1771     return((Image *) NULL);
1772   status=SetImageStorageClass(similarity_image,DirectClass,exception);
1773   if (status == MagickFalse)
1774     {
1775       similarity_image=DestroyImage(similarity_image);
1776       return((Image *) NULL);
1777     }
1778   (void) SetImageAlphaChannel(similarity_image,DeactivateAlphaChannel,
1779     exception);
1780   /*
1781     Measure similarity of reference image against image.
1782   */
1783   status=MagickTrue;
1784   progress=0;
1785   similarity_view=AcquireAuthenticCacheView(similarity_image,exception);
1786 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1787   #pragma omp parallel for schedule(static,4) \
1788     shared(progress,status,similarity_metric) \
1789     magick_threads(image,image,image->rows,1)
1790 #endif
1791   for (y=0; y < (ssize_t) (image->rows-reference->rows+1); y++)
1792   {
1793     double
1794       similarity;
1795
1796     register Quantum
1797       *restrict q;
1798
1799     register ssize_t
1800       x;
1801
1802     if (status == MagickFalse)
1803       continue;
1804 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1805     #pragma omp flush(similarity_metric)
1806 #endif
1807     if (*similarity_metric <= similarity_threshold)
1808       continue;
1809     q=GetCacheViewAuthenticPixels(similarity_view,0,y,similarity_image->columns,
1810       1,exception);
1811     if (q == (Quantum *) NULL)
1812       {
1813         status=MagickFalse;
1814         continue;
1815       }
1816     for (x=0; x < (ssize_t) (image->columns-reference->columns+1); x++)
1817     {
1818       register ssize_t
1819         i;
1820
1821 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1822       #pragma omp flush(similarity_metric)
1823 #endif
1824       if (*similarity_metric <= similarity_threshold)
1825         break;
1826       similarity=GetSimilarityMetric(image,reference,metric,x,y,exception);
1827 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1828       #pragma omp critical (MagickCore_SimilarityImage)
1829 #endif
1830       if ((metric == NormalizedCrossCorrelationErrorMetric) ||
1831           (metric == UndefinedErrorMetric))
1832         similarity=1.0-similarity;
1833       if (similarity < *similarity_metric)
1834         {
1835           offset->x=x;
1836           offset->y=y;
1837           *similarity_metric=similarity;
1838         }
1839       if (metric == PerceptualHashErrorMetric)
1840         similarity=MagickMin(0.01*similarity,1.0);
1841       if (GetPixelReadMask(similarity_image,q) == 0)
1842         {
1843           SetPixelBackgoundColor(similarity_image,q);
1844           q+=GetPixelChannels(similarity_image);
1845           continue;
1846         }
1847       for (i=0; i < (ssize_t) GetPixelChannels(similarity_image); i++)
1848       {
1849         PixelChannel channel=GetPixelChannelChannel(image,i);
1850         PixelTrait traits=GetPixelChannelTraits(image,channel);
1851         PixelTrait similarity_traits=GetPixelChannelTraits(similarity_image,
1852           channel);
1853         if ((traits == UndefinedPixelTrait) ||
1854             (similarity_traits == UndefinedPixelTrait) ||
1855             ((similarity_traits & UpdatePixelTrait) == 0))
1856           continue;
1857         SetPixelChannel(similarity_image,channel,ClampToQuantum(QuantumRange-
1858           QuantumRange*similarity),q);
1859       }
1860       q+=GetPixelChannels(similarity_image);
1861     }
1862     if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
1863       status=MagickFalse;
1864     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1865       {
1866         MagickBooleanType
1867           proceed;
1868
1869 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1870         #pragma omp critical (MagickCore_SimilarityImage)
1871 #endif
1872         proceed=SetImageProgress(image,SimilarityImageTag,progress++,
1873           image->rows);
1874         if (proceed == MagickFalse)
1875           status=MagickFalse;
1876       }
1877   }
1878   similarity_view=DestroyCacheView(similarity_view);
1879   if (status == MagickFalse)
1880     similarity_image=DestroyImage(similarity_image);
1881   return(similarity_image);
1882 }