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