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

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