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