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 #define Log10Epsilon  (1.0e-11)
1094
1095  if (fabs(x) < Log10Epsilon)
1096    return(log10(Log10Epsilon));
1097  return(log10(fabs(x)));
1098 }
1099
1100 static MagickBooleanType GetPeakSignalToNoiseRatio(const Image *image,
1101   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1102 {
1103   MagickBooleanType
1104     status;
1105
1106   register ssize_t
1107     i;
1108
1109   status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1110   for (i=0; i <= MaxPixelChannels; i++)
1111     distortion[i]=20.0*MagickLog10((double) 1.0/sqrt(distortion[i]));
1112   return(status);
1113 }
1114
1115 static Image *PerceptualImageHash(const Image *image,
1116   const ColorspaceType colorspace,ExceptionInfo *exception)
1117 {
1118   Image
1119     *phash_image;
1120
1121   MagickBooleanType
1122     status;
1123
1124   /*
1125     Transform colorspace then blur perceptual hash image.
1126   */
1127   phash_image=BlurImage(image,0.0,1.0,exception);
1128   if (phash_image == (Image *) NULL)
1129     return((Image *) NULL);
1130   phash_image->depth=8;
1131   status=TransformImageColorspace(phash_image,colorspace,exception);
1132   if (status == MagickFalse)
1133     phash_image=DestroyImage(phash_image);
1134   return(phash_image);
1135 }
1136
1137 static MagickBooleanType GetPerceptualHashDistortion(const Image *image,
1138   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1139 {
1140   ChannelMoments
1141     *image_moments,
1142     *reconstruct_moments;
1143
1144   Image
1145     *phash_image;
1146
1147   MagickBooleanType
1148     grayscale;
1149
1150   ssize_t
1151     channel;
1152
1153   /*
1154     Compute perceptual hash in the native image colorspace.
1155   */
1156   phash_image=PerceptualImageHash(image,image->colorspace,exception);
1157   if (phash_image == (Image *) NULL)
1158     return(MagickFalse);
1159   image_moments=GetImageMoments(phash_image,exception);
1160   phash_image=DestroyImage(phash_image);
1161   if (image_moments == (ChannelMoments *) NULL)
1162     return(MagickFalse);
1163   phash_image=PerceptualImageHash(reconstruct_image,
1164     reconstruct_image->colorspace,exception);
1165   if (phash_image == (Image *) NULL)
1166     {
1167       image_moments=(ChannelMoments *) RelinquishMagickMemory(image_moments);
1168       return(MagickFalse);
1169     }
1170   reconstruct_moments=GetImageMoments(phash_image,exception);
1171   phash_image=DestroyImage(phash_image);
1172   if (reconstruct_moments == (ChannelMoments *) NULL)
1173     {
1174       image_moments=(ChannelMoments *) RelinquishMagickMemory(image_moments);
1175       return(MagickFalse);
1176     }
1177   /*
1178     Compute sum of moment differences squared.
1179   */
1180 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1181   #pragma omp parallel for schedule(static,4)
1182 #endif
1183   for (channel=0; channel < MaxPixelChannels; channel++)
1184   {
1185     double
1186       difference;
1187
1188     register ssize_t
1189       i;
1190
1191     difference=0.0;
1192     for (i=0; i < 7; i++)
1193     {
1194       double
1195         alpha,
1196         beta;
1197
1198       alpha=MagickLog10(image_moments[channel].I[i]);
1199       beta=MagickLog10(reconstruct_moments[channel].I[i]);
1200       difference+=(beta-alpha)*(beta-alpha);
1201     }
1202     distortion[channel]+=difference;
1203 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1204     #pragma omp critical (MagickCore_GetPerceptualHashDistortion)
1205 #endif
1206     distortion[CompositePixelChannel]+=difference;
1207   }
1208   image_moments=(ChannelMoments *) RelinquishMagickMemory(image_moments);
1209   reconstruct_moments=(ChannelMoments *) RelinquishMagickMemory(
1210     reconstruct_moments);
1211   grayscale=(IsImageGray(image,exception) != MagickFalse) &&
1212     (IsImageGray(reconstruct_image,exception) != MagickFalse) ? MagickTrue :
1213     MagickFalse;
1214   if (grayscale != MagickFalse)
1215     return(MagickTrue);
1216   /*
1217     Compute perceptual hash in the HCLP colorspace.
1218   */
1219   phash_image=PerceptualImageHash(image,HCLpColorspace,exception);
1220   if (phash_image == (Image *) NULL)
1221     return(MagickFalse);
1222   image_moments=GetImageMoments(phash_image,exception);
1223   phash_image=DestroyImage(phash_image);
1224   if (image_moments == (ChannelMoments *) NULL)
1225     return(MagickFalse);
1226   phash_image=PerceptualImageHash(reconstruct_image,HCLpColorspace,exception);
1227   if (phash_image == (Image *) NULL)
1228     {
1229       image_moments=(ChannelMoments *) RelinquishMagickMemory(image_moments);
1230       return(MagickFalse);
1231     }
1232   reconstruct_moments=GetImageMoments(phash_image,exception);
1233   phash_image=DestroyImage(phash_image);
1234   if (reconstruct_moments == (ChannelMoments *) NULL)
1235     {
1236       image_moments=(ChannelMoments *) RelinquishMagickMemory(image_moments);
1237       return(MagickFalse);
1238     }
1239   /*
1240     Compute sum of moment differences squared.
1241   */
1242 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1243   #pragma omp parallel for schedule(static,4)
1244 #endif
1245   for (channel=0; channel < MaxPixelChannels; channel++)
1246   {
1247     double
1248       difference;
1249
1250     register ssize_t
1251       i;
1252
1253     difference=0.0;
1254     for (i=0; i < 7; i++)
1255     {
1256       double
1257         alpha,
1258         beta;
1259
1260       alpha=MagickLog10(image_moments[channel].I[i]);
1261       beta=MagickLog10(reconstruct_moments[channel].I[i]);
1262       difference+=(beta-alpha)*(beta-alpha);
1263     }
1264     distortion[channel]+=difference;
1265 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1266     #pragma omp critical (MagickCore_GetPerceptualHashDistortion)
1267 #endif
1268     distortion[CompositePixelChannel]+=difference;
1269   }
1270   image_moments=(ChannelMoments *) RelinquishMagickMemory(image_moments);
1271   reconstruct_moments=(ChannelMoments *) RelinquishMagickMemory(
1272     reconstruct_moments);
1273   return(MagickTrue);
1274 }
1275
1276
1277 static MagickBooleanType GetRootMeanSquaredDistortion(const Image *image,
1278   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1279 {
1280   MagickBooleanType
1281     status;
1282
1283   register ssize_t
1284     i;
1285
1286   status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1287   for (i=0; i <= MaxPixelChannels; i++)
1288     distortion[i]=sqrt(distortion[i]);
1289   return(status);
1290 }
1291
1292 MagickExport MagickBooleanType GetImageDistortion(Image *image,
1293   const Image *reconstruct_image,const MetricType metric,double *distortion,
1294   ExceptionInfo *exception)
1295 {
1296   double
1297     *channel_distortion;
1298
1299   MagickBooleanType
1300     status;
1301
1302   size_t
1303     length;
1304
1305   assert(image != (Image *) NULL);
1306   assert(image->signature == MagickSignature);
1307   if (image->debug != MagickFalse)
1308     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1309   assert(reconstruct_image != (const Image *) NULL);
1310   assert(reconstruct_image->signature == MagickSignature);
1311   assert(distortion != (double *) NULL);
1312   *distortion=0.0;
1313   if (image->debug != MagickFalse)
1314     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1315   if (metric != PerceptualHashErrorMetric)
1316     if ((reconstruct_image->columns != image->columns) ||
1317         (reconstruct_image->rows != image->rows))
1318       ThrowBinaryException(ImageError,"ImageSizeDiffers",image->filename);
1319   /*
1320     Get image distortion.
1321   */
1322   length=MaxPixelChannels+1;
1323   channel_distortion=(double *) AcquireQuantumMemory(length,
1324     sizeof(*channel_distortion));
1325   if (channel_distortion == (double *) NULL)
1326     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1327   (void) ResetMagickMemory(channel_distortion,0,length*
1328     sizeof(*channel_distortion));
1329   switch (metric)
1330   {
1331     case AbsoluteErrorMetric:
1332     {
1333       status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1334         exception);
1335       break;
1336     }
1337     case FuzzErrorMetric:
1338     {
1339       status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1340         exception);
1341       break;
1342     }
1343     case MeanAbsoluteErrorMetric:
1344     {
1345       status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1346         channel_distortion,exception);
1347       break;
1348     }
1349     case MeanErrorPerPixelErrorMetric:
1350     {
1351       status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1352         exception);
1353       break;
1354     }
1355     case MeanSquaredErrorMetric:
1356     {
1357       status=GetMeanSquaredDistortion(image,reconstruct_image,
1358         channel_distortion,exception);
1359       break;
1360     }
1361     case NormalizedCrossCorrelationErrorMetric:
1362     default:
1363     {
1364       status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1365         channel_distortion,exception);
1366       break;
1367     }
1368     case PeakAbsoluteErrorMetric:
1369     {
1370       status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1371         channel_distortion,exception);
1372       break;
1373     }
1374     case PeakSignalToNoiseRatioErrorMetric:
1375     {
1376       status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1377         channel_distortion,exception);
1378       break;
1379     }
1380     case PerceptualHashErrorMetric:
1381     {
1382       status=GetPerceptualHashDistortion(image,reconstruct_image,
1383         channel_distortion,exception);
1384       break;
1385     }
1386     case RootMeanSquaredErrorMetric:
1387     {
1388       status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1389         channel_distortion,exception);
1390       break;
1391     }
1392   }
1393   *distortion=channel_distortion[CompositePixelChannel];
1394   channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1395   (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
1396     *distortion);
1397   return(status);
1398 }
1399 \f
1400 /*
1401 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1402 %                                                                             %
1403 %                                                                             %
1404 %                                                                             %
1405 %   G e t I m a g e D i s t o r t i o n s                                     %
1406 %                                                                             %
1407 %                                                                             %
1408 %                                                                             %
1409 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1410 %
1411 %  GetImageDistortions() compares the pixel channels of an image to a
1412 %  reconstructed image and returns the specified distortion metric for each
1413 %  channel.
1414 %
1415 %  The format of the GetImageDistortions method is:
1416 %
1417 %      double *GetImageDistortions(const Image *image,
1418 %        const Image *reconstruct_image,const MetricType metric,
1419 %        ExceptionInfo *exception)
1420 %
1421 %  A description of each parameter follows:
1422 %
1423 %    o image: the image.
1424 %
1425 %    o reconstruct_image: the reconstruct image.
1426 %
1427 %    o metric: the metric.
1428 %
1429 %    o exception: return any errors or warnings in this structure.
1430 %
1431 */
1432 MagickExport double *GetImageDistortions(Image *image,
1433   const Image *reconstruct_image,const MetricType metric,
1434   ExceptionInfo *exception)
1435 {
1436   double
1437     *channel_distortion;
1438
1439   MagickBooleanType
1440     status;
1441
1442   size_t
1443     length;
1444
1445   assert(image != (Image *) NULL);
1446   assert(image->signature == MagickSignature);
1447   if (image->debug != MagickFalse)
1448     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1449   assert(reconstruct_image != (const Image *) NULL);
1450   assert(reconstruct_image->signature == MagickSignature);
1451   if (image->debug != MagickFalse)
1452     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1453   if (metric != PerceptualHashErrorMetric)
1454     if ((reconstruct_image->columns != image->columns) ||
1455         (reconstruct_image->rows != image->rows))
1456       {
1457         (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1458           "ImageSizeDiffers","`%s'",image->filename);
1459         return((double *) NULL);
1460       }
1461   /*
1462     Get image distortion.
1463   */
1464   length=MaxPixelChannels+1UL;
1465   channel_distortion=(double *) AcquireQuantumMemory(length,
1466     sizeof(*channel_distortion));
1467   if (channel_distortion == (double *) NULL)
1468     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1469   (void) ResetMagickMemory(channel_distortion,0,length*
1470     sizeof(*channel_distortion));
1471   status=MagickTrue;
1472   switch (metric)
1473   {
1474     case AbsoluteErrorMetric:
1475     {
1476       status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1477         exception);
1478       break;
1479     }
1480     case FuzzErrorMetric:
1481     {
1482       status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1483         exception);
1484       break;
1485     }
1486     case MeanAbsoluteErrorMetric:
1487     {
1488       status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1489         channel_distortion,exception);
1490       break;
1491     }
1492     case MeanErrorPerPixelErrorMetric:
1493     {
1494       status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1495         exception);
1496       break;
1497     }
1498     case MeanSquaredErrorMetric:
1499     {
1500       status=GetMeanSquaredDistortion(image,reconstruct_image,
1501         channel_distortion,exception);
1502       break;
1503     }
1504     case NormalizedCrossCorrelationErrorMetric:
1505     default:
1506     {
1507       status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1508         channel_distortion,exception);
1509       break;
1510     }
1511     case PeakAbsoluteErrorMetric:
1512     {
1513       status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1514         channel_distortion,exception);
1515       break;
1516     }
1517     case PeakSignalToNoiseRatioErrorMetric:
1518     {
1519       status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1520         channel_distortion,exception);
1521       break;
1522     }
1523     case PerceptualHashErrorMetric:
1524     {
1525       status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1526         channel_distortion,exception);
1527       break;
1528     }
1529     case RootMeanSquaredErrorMetric:
1530     {
1531       status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1532         channel_distortion,exception);
1533       break;
1534     }
1535   }
1536   if (status == MagickFalse)
1537     {
1538       channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1539       return((double *) NULL);
1540     }
1541   return(channel_distortion);
1542 }
1543 \f
1544 /*
1545 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1546 %                                                                             %
1547 %                                                                             %
1548 %                                                                             %
1549 %  I s I m a g e s E q u a l                                                  %
1550 %                                                                             %
1551 %                                                                             %
1552 %                                                                             %
1553 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1554 %
1555 %  IsImagesEqual() measures the difference between colors at each pixel
1556 %  location of two images.  A value other than 0 means the colors match
1557 %  exactly.  Otherwise an error measure is computed by summing over all
1558 %  pixels in an image the distance squared in RGB space between each image
1559 %  pixel and its corresponding pixel in the reconstruct image.  The error
1560 %  measure is assigned to these image members:
1561 %
1562 %    o mean_error_per_pixel:  The mean error for any single pixel in
1563 %      the image.
1564 %
1565 %    o normalized_mean_error:  The normalized mean quantization error for
1566 %      any single pixel in the image.  This distance measure is normalized to
1567 %      a range between 0 and 1.  It is independent of the range of red, green,
1568 %      and blue values in the image.
1569 %
1570 %    o normalized_maximum_error:  The normalized maximum quantization
1571 %      error for any single pixel in the image.  This distance measure is
1572 %      normalized to a range between 0 and 1.  It is independent of the range
1573 %      of red, green, and blue values in your image.
1574 %
1575 %  A small normalized mean square error, accessed as
1576 %  image->normalized_mean_error, suggests the images are very similar in
1577 %  spatial layout and color.
1578 %
1579 %  The format of the IsImagesEqual method is:
1580 %
1581 %      MagickBooleanType IsImagesEqual(Image *image,
1582 %        const Image *reconstruct_image,ExceptionInfo *exception)
1583 %
1584 %  A description of each parameter follows.
1585 %
1586 %    o image: the image.
1587 %
1588 %    o reconstruct_image: the reconstruct image.
1589 %
1590 %    o exception: return any errors or warnings in this structure.
1591 %
1592 */
1593 MagickExport MagickBooleanType IsImagesEqual(Image *image,
1594   const Image *reconstruct_image,ExceptionInfo *exception)
1595 {
1596   CacheView
1597     *image_view,
1598     *reconstruct_view;
1599
1600   MagickBooleanType
1601     status;
1602
1603   double
1604     area,
1605     maximum_error,
1606     mean_error,
1607     mean_error_per_pixel;
1608
1609   ssize_t
1610     y;
1611
1612   assert(image != (Image *) NULL);
1613   assert(image->signature == MagickSignature);
1614   assert(reconstruct_image != (const Image *) NULL);
1615   assert(reconstruct_image->signature == MagickSignature);
1616   if ((reconstruct_image->columns != image->columns) ||
1617       (reconstruct_image->rows != image->rows))
1618     ThrowBinaryException(ImageError,"ImageSizeDiffers",image->filename);
1619   area=0.0;
1620   maximum_error=0.0;
1621   mean_error_per_pixel=0.0;
1622   mean_error=0.0;
1623   image_view=AcquireVirtualCacheView(image,exception);
1624   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1625   for (y=0; y < (ssize_t) image->rows; y++)
1626   {
1627     register const Quantum
1628       *restrict p,
1629       *restrict q;
1630
1631     register ssize_t
1632       x;
1633
1634     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1635     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
1636       1,exception);
1637     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
1638       break;
1639     for (x=0; x < (ssize_t) image->columns; x++)
1640     {
1641       register ssize_t
1642         i;
1643
1644       if (GetPixelReadMask(image,p) == 0)
1645         {
1646           p+=GetPixelChannels(image);
1647           q+=GetPixelChannels(reconstruct_image);
1648           continue;
1649         }
1650       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1651       {
1652         double
1653           distance;
1654
1655         PixelChannel channel=GetPixelChannelChannel(image,i);
1656         PixelTrait traits=GetPixelChannelTraits(image,channel);
1657         PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
1658           channel);
1659         if ((traits == UndefinedPixelTrait) ||
1660             (reconstruct_traits == UndefinedPixelTrait) ||
1661             ((reconstruct_traits & UpdatePixelTrait) == 0))
1662           continue;
1663         distance=fabs(p[i]-(double) GetPixelChannel(reconstruct_image,
1664           channel,q));
1665         mean_error_per_pixel+=distance;
1666         mean_error+=distance*distance;
1667         if (distance > maximum_error)
1668           maximum_error=distance;
1669         area++;
1670       }
1671       p+=GetPixelChannels(image);
1672       q+=GetPixelChannels(reconstruct_image);
1673     }
1674   }
1675   reconstruct_view=DestroyCacheView(reconstruct_view);
1676   image_view=DestroyCacheView(image_view);
1677   image->error.mean_error_per_pixel=(double) (mean_error_per_pixel/area);
1678   image->error.normalized_mean_error=(double) (QuantumScale*QuantumScale*
1679     mean_error/area);
1680   image->error.normalized_maximum_error=(double) (QuantumScale*maximum_error);
1681   status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
1682   return(status);
1683 }
1684 \f
1685 /*
1686 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1687 %                                                                             %
1688 %                                                                             %
1689 %                                                                             %
1690 %   S i m i l a r i t y I m a g e                                             %
1691 %                                                                             %
1692 %                                                                             %
1693 %                                                                             %
1694 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1695 %
1696 %  SimilarityImage() compares the reference image of the image and returns the
1697 %  best match offset.  In addition, it returns a similarity image such that an
1698 %  exact match location is completely white and if none of the pixels match,
1699 %  black, otherwise some gray level in-between.
1700 %
1701 %  The format of the SimilarityImageImage method is:
1702 %
1703 %      Image *SimilarityImage(const Image *image,const Image *reference,
1704 %        const MetricType metric,const double similarity_threshold,
1705 %        RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
1706 %
1707 %  A description of each parameter follows:
1708 %
1709 %    o image: the image.
1710 %
1711 %    o reference: find an area of the image that closely resembles this image.
1712 %
1713 %    o metric: the metric.
1714 %
1715 %    o similarity_threshold: minimum distortion for (sub)image match.
1716 %
1717 %    o offset: the best match offset of the reference image within the image.
1718 %
1719 %    o similarity: the computed similarity between the images.
1720 %
1721 %    o exception: return any errors or warnings in this structure.
1722 %
1723 */
1724
1725 static double GetSimilarityMetric(const Image *image,const Image *reference,
1726   const MetricType metric,const ssize_t x_offset,const ssize_t y_offset,
1727   ExceptionInfo *exception)
1728 {
1729   double
1730     distortion;
1731
1732   Image
1733     *similarity_image;
1734
1735   MagickBooleanType
1736     status;
1737
1738   RectangleInfo
1739     geometry;
1740
1741   SetGeometry(reference,&geometry);
1742   geometry.x=x_offset;
1743   geometry.y=y_offset;
1744   similarity_image=CropImage(image,&geometry,exception);
1745   if (similarity_image == (Image *) NULL)
1746     return(0.0);
1747   distortion=0.0;
1748   status=GetImageDistortion(similarity_image,reference,metric,&distortion,
1749     exception);
1750   similarity_image=DestroyImage(similarity_image);
1751   if (status == MagickFalse)
1752     return(0.0);
1753   return(distortion);
1754 }
1755
1756 static inline double MagickMin(const double x,const double y)
1757 {
1758   if (x < y)
1759     return(x);
1760   return(y);
1761 }
1762
1763 MagickExport Image *SimilarityImage(Image *image,const Image *reference,
1764   const MetricType metric,const double similarity_threshold,
1765   RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
1766 {
1767 #define SimilarityImageTag  "Similarity/Image"
1768
1769   CacheView
1770     *similarity_view;
1771
1772   Image
1773     *similarity_image;
1774
1775   MagickBooleanType
1776     status;
1777
1778   MagickOffsetType
1779     progress;
1780
1781   ssize_t
1782     y;
1783
1784   assert(image != (const Image *) NULL);
1785   assert(image->signature == MagickSignature);
1786   if (image->debug != MagickFalse)
1787     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1788   assert(exception != (ExceptionInfo *) NULL);
1789   assert(exception->signature == MagickSignature);
1790   assert(offset != (RectangleInfo *) NULL);
1791   SetGeometry(reference,offset);
1792   *similarity_metric=MagickMaximumValue;
1793   if ((reference->columns > image->columns) || (reference->rows > image->rows))
1794     ThrowImageException(ImageError,"ImageSizeDiffers");
1795   similarity_image=CloneImage(image,image->columns-reference->columns+1,
1796     image->rows-reference->rows+1,MagickTrue,exception);
1797   if (similarity_image == (Image *) NULL)
1798     return((Image *) NULL);
1799   status=SetImageStorageClass(similarity_image,DirectClass,exception);
1800   if (status == MagickFalse)
1801     {
1802       similarity_image=DestroyImage(similarity_image);
1803       return((Image *) NULL);
1804     }
1805   (void) SetImageAlphaChannel(similarity_image,DeactivateAlphaChannel,
1806     exception);
1807   /*
1808     Measure similarity of reference image against image.
1809   */
1810   status=MagickTrue;
1811   progress=0;
1812   similarity_view=AcquireAuthenticCacheView(similarity_image,exception);
1813 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1814   #pragma omp parallel for schedule(static,4) \
1815     shared(progress,status,similarity_metric) \
1816     magick_threads(image,image,image->rows,1)
1817 #endif
1818   for (y=0; y < (ssize_t) (image->rows-reference->rows+1); y++)
1819   {
1820     double
1821       similarity;
1822
1823     register Quantum
1824       *restrict q;
1825
1826     register ssize_t
1827       x;
1828
1829     if (status == MagickFalse)
1830       continue;
1831 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1832     #pragma omp flush(similarity_metric)
1833 #endif
1834     if (*similarity_metric <= similarity_threshold)
1835       continue;
1836     q=GetCacheViewAuthenticPixels(similarity_view,0,y,similarity_image->columns,
1837       1,exception);
1838     if (q == (Quantum *) NULL)
1839       {
1840         status=MagickFalse;
1841         continue;
1842       }
1843     for (x=0; x < (ssize_t) (image->columns-reference->columns+1); x++)
1844     {
1845       register ssize_t
1846         i;
1847
1848 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1849       #pragma omp flush(similarity_metric)
1850 #endif
1851       if (*similarity_metric <= similarity_threshold)
1852         break;
1853       similarity=GetSimilarityMetric(image,reference,metric,x,y,exception);
1854 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1855       #pragma omp critical (MagickCore_SimilarityImage)
1856 #endif
1857       if (similarity < *similarity_metric)
1858         {
1859           offset->x=x;
1860           offset->y=y;
1861           *similarity_metric=similarity;
1862         }
1863       if (metric == PerceptualHashErrorMetric)
1864         similarity=MagickMin(0.01*similarity,1.0);
1865       if (GetPixelReadMask(similarity_image,q) == 0)
1866         {
1867           SetPixelBackgoundColor(similarity_image,q);
1868           q+=GetPixelChannels(similarity_image);
1869           continue;
1870         }
1871       for (i=0; i < (ssize_t) GetPixelChannels(similarity_image); i++)
1872       {
1873         PixelChannel channel=GetPixelChannelChannel(image,i);
1874         PixelTrait traits=GetPixelChannelTraits(image,channel);
1875         PixelTrait similarity_traits=GetPixelChannelTraits(similarity_image,
1876           channel);
1877         if ((traits == UndefinedPixelTrait) ||
1878             (similarity_traits == UndefinedPixelTrait) ||
1879             ((similarity_traits & UpdatePixelTrait) == 0))
1880           continue;
1881         SetPixelChannel(similarity_image,channel,ClampToQuantum(QuantumRange-
1882           QuantumRange*similarity),q);
1883       }
1884       q+=GetPixelChannels(similarity_image);
1885     }
1886     if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
1887       status=MagickFalse;
1888     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1889       {
1890         MagickBooleanType
1891           proceed;
1892
1893 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1894         #pragma omp critical (MagickCore_SimilarityImage)
1895 #endif
1896         proceed=SetImageProgress(image,SimilarityImageTag,progress++,
1897           image->rows);
1898         if (proceed == MagickFalse)
1899           status=MagickFalse;
1900       }
1901   }
1902   similarity_view=DestroyCacheView(similarity_view);
1903   if (status == MagickFalse)
1904     similarity_image=DestroyImage(similarity_image);
1905   return(similarity_image);
1906 }