granicus.if.org Git - imagemagick/blob - MagickCore/compare.c

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