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