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