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