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