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