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