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