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

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