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

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