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