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