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 %                                  Cristy                                     %
  17 %                               December 2003                                 %
  18 %                                                                             %
  19 %                                                                             %
  20 %  Copyright 1999-2019 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 %    https://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/attribute.h"
  46 #include "MagickCore/cache-view.h"
  47 #include "MagickCore/channel.h"
  48 #include "MagickCore/client.h"
  49 #include "MagickCore/color.h"
  50 #include "MagickCore/color-private.h"
  51 #include "MagickCore/colorspace.h"
  52 #include "MagickCore/colorspace-private.h"
  53 #include "MagickCore/compare.h"
  54 #include "MagickCore/composite-private.h"
  55 #include "MagickCore/constitute.h"
  56 #include "MagickCore/exception-private.h"
  57 #include "MagickCore/geometry.h"
  58 #include "MagickCore/image-private.h"
  59 #include "MagickCore/list.h"
  60 #include "MagickCore/log.h"
  61 #include "MagickCore/memory_.h"
  62 #include "MagickCore/monitor.h"
  63 #include "MagickCore/monitor-private.h"
  64 #include "MagickCore/option.h"
  65 #include "MagickCore/pixel-accessor.h"
  66 #include "MagickCore/property.h"
  67 #include "MagickCore/resource_.h"
  68 #include "MagickCore/string_.h"
  69 #include "MagickCore/statistic.h"
  70 #include "MagickCore/string-private.h"
  71 #include "MagickCore/thread-private.h"
  72 #include "MagickCore/transform.h"
  73 #include "MagickCore/utility.h"
  74 #include "MagickCore/version.h"
  75 \f
  76 /*
  77 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  78 %                                                                             %
  79 %                                                                             %
  80 %                                                                             %
  81 %   C o m p a r e I m a g e                                                   %
  82 %                                                                             %
  83 %                                                                             %
  84 %                                                                             %
  85 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  86 %
  87 %  CompareImages() compares one or more pixel channels of an image to a
  88 %  reconstructed image and returns the difference image.
  89 %
  90 %  The format of the CompareImages method is:
  91 %
  92 %      Image *CompareImages(const Image *image,const Image *reconstruct_image,
  93 %        const MetricType metric,double *distortion,ExceptionInfo *exception)
  94 %
  95 %  A description of each parameter follows:
  96 %
  97 %    o image: the image.
  98 %
  99 %    o reconstruct_image: the reconstruct image.
 100 %
 101 %    o metric: the metric.
 102 %
 103 %    o distortion: the computed distortion between the images.
 104 %
 105 %    o exception: return any errors or warnings in this structure.
 106 %
 107 */
 108
 109 static size_t GetImageChannels(const Image *image)
 110 {
 111   register ssize_t
 112     i;
 113
 114   size_t
 115     channels;
 116
 117   channels=0;
 118   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 119   {
 120     PixelChannel channel = GetPixelChannelChannel(image,i);
 121     PixelTrait traits = GetPixelChannelTraits(image,channel);
 122     if ((traits & UpdatePixelTrait) != 0)
 123       channels++;
 124   }
 125   return(channels == 0 ? (size_t) 1 : channels);
 126 }
 127
 128 MagickExport Image *CompareImages(Image *image,const Image *reconstruct_image,
 129   const MetricType metric,double *distortion,ExceptionInfo *exception)
 130 {
 131   CacheView
 132     *highlight_view,
 133     *image_view,
 134     *reconstruct_view;
 135
 136   const char
 137     *artifact;
 138
 139   double
 140     fuzz;
 141
 142   Image
 143     *clone_image,
 144     *difference_image,
 145     *highlight_image;
 146
 147   MagickBooleanType
 148     status;
 149
 150   PixelInfo
 151     highlight,
 152     lowlight,
 153     masklight;
 154
 155   RectangleInfo
 156     geometry;
 157
 158   size_t
 159     columns,
 160     rows;
 161
 162   ssize_t
 163     y;
 164
 165   assert(image != (Image *) NULL);
 166   assert(image->signature == MagickCoreSignature);
 167   if (image->debug != MagickFalse)
 168     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 169   assert(reconstruct_image != (const Image *) NULL);
 170   assert(reconstruct_image->signature == MagickCoreSignature);
 171   assert(distortion != (double *) NULL);
 172   *distortion=0.0;
 173   if (image->debug != MagickFalse)
 174     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 175   status=GetImageDistortion(image,reconstruct_image,metric,distortion,
 176     exception);
 177   if (status == MagickFalse)
 178     return((Image *) NULL);
 179   columns=MagickMax(image->columns,reconstruct_image->columns);
 180   rows=MagickMax(image->rows,reconstruct_image->rows);
 181   SetGeometry(image,&geometry);
 182   geometry.width=columns;
 183   geometry.height=rows;
 184   clone_image=CloneImage(image,0,0,MagickTrue,exception);
 185   if (clone_image == (Image *) NULL)
 186     return((Image *) NULL);
 187   (void) SetImageMask(clone_image,ReadPixelMask,(Image *) NULL,exception);
 188   difference_image=ExtentImage(clone_image,&geometry,exception);
 189   clone_image=DestroyImage(clone_image);
 190   if (difference_image == (Image *) NULL)
 191     return((Image *) NULL);
 192   (void) SetImageAlphaChannel(difference_image,OpaqueAlphaChannel,exception);
 193   highlight_image=CloneImage(image,columns,rows,MagickTrue,exception);
 194   if (highlight_image == (Image *) NULL)
 195     {
 196       difference_image=DestroyImage(difference_image);
 197       return((Image *) NULL);
 198     }
 199   status=SetImageStorageClass(highlight_image,DirectClass,exception);
 200   if (status == MagickFalse)
 201     {
 202       difference_image=DestroyImage(difference_image);
 203       highlight_image=DestroyImage(highlight_image);
 204       return((Image *) NULL);
 205     }
 206   (void) SetImageMask(highlight_image,ReadPixelMask,(Image *) NULL,exception);
 207   (void) SetImageAlphaChannel(highlight_image,OpaqueAlphaChannel,exception);
 208   (void) QueryColorCompliance("#f1001ecc",AllCompliance,&highlight,exception);
 209   artifact=GetImageArtifact(image,"compare:highlight-color");
 210   if (artifact != (const char *) NULL)
 211     (void) QueryColorCompliance(artifact,AllCompliance,&highlight,exception);
 212   (void) QueryColorCompliance("#ffffffcc",AllCompliance,&lowlight,exception);
 213   artifact=GetImageArtifact(image,"compare:lowlight-color");
 214   if (artifact != (const char *) NULL)
 215     (void) QueryColorCompliance(artifact,AllCompliance,&lowlight,exception);
 216   (void) QueryColorCompliance("#888888cc",AllCompliance,&masklight,exception);
 217   artifact=GetImageArtifact(image,"compare:masklight-color");
 218   if (artifact != (const char *) NULL)
 219     (void) QueryColorCompliance(artifact,AllCompliance,&masklight,exception);
 220   /*
 221     Generate difference image.
 222   */
 223   status=MagickTrue;
 224   fuzz=GetFuzzyColorDistance(image,reconstruct_image);
 225   image_view=AcquireVirtualCacheView(image,exception);
 226   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
 227   highlight_view=AcquireAuthenticCacheView(highlight_image,exception);
 228 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 229   #pragma omp parallel for schedule(static) shared(status) \
 230     magick_number_threads(image,highlight_image,rows,1)
 231 #endif
 232   for (y=0; y < (ssize_t) rows; y++)
 233   {
 234     MagickBooleanType
 235       sync;
 236
 237     register const Quantum
 238       *magick_restrict p,
 239       *magick_restrict q;
 240
 241     register Quantum
 242       *magick_restrict r;
 243
 244     register ssize_t
 245       x;
 246
 247     if (status == MagickFalse)
 248       continue;
 249     p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
 250     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
 251     r=QueueCacheViewAuthenticPixels(highlight_view,0,y,columns,1,exception);
 252     if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL) ||
 253         (r == (Quantum *) NULL))
 254       {
 255         status=MagickFalse;
 256         continue;
 257       }
 258     for (x=0; x < (ssize_t) columns; x++)
 259     {
 260       double
 261         Da,
 262         Sa;
 263
 264       MagickStatusType
 265         difference;
 266
 267       register ssize_t
 268         i;
 269
 270       if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
 271           (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
 272         {
 273           SetPixelViaPixelInfo(highlight_image,&masklight,r);
 274           p+=GetPixelChannels(image);
 275           q+=GetPixelChannels(reconstruct_image);
 276           r+=GetPixelChannels(highlight_image);
 277           continue;
 278         }
 279       difference=MagickFalse;
 280       Sa=QuantumScale*GetPixelAlpha(image,p);
 281       Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
 282       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 283       {
 284         double
 285           distance;
 286
 287         PixelChannel channel = GetPixelChannelChannel(image,i);
 288         PixelTrait traits = GetPixelChannelTraits(image,channel);
 289         PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
 290           channel);
 291         if ((traits == UndefinedPixelTrait) ||
 292             (reconstruct_traits == UndefinedPixelTrait) ||
 293             ((reconstruct_traits & UpdatePixelTrait) == 0))
 294           continue;
 295         if (channel == AlphaPixelChannel)
 296           distance=(double) p[i]-GetPixelChannel(reconstruct_image,channel,q);
 297         else
 298           distance=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
 299         if ((distance*distance) > fuzz)
 300           {
 301             difference=MagickTrue;
 302             break;
 303           }
 304       }
 305       if (difference == MagickFalse)
 306         SetPixelViaPixelInfo(highlight_image,&lowlight,r);
 307       else
 308         SetPixelViaPixelInfo(highlight_image,&highlight,r);
 309       p+=GetPixelChannels(image);
 310       q+=GetPixelChannels(reconstruct_image);
 311       r+=GetPixelChannels(highlight_image);
 312     }
 313     sync=SyncCacheViewAuthenticPixels(highlight_view,exception);
 314     if (sync == MagickFalse)
 315       status=MagickFalse;
 316   }
 317   highlight_view=DestroyCacheView(highlight_view);
 318   reconstruct_view=DestroyCacheView(reconstruct_view);
 319   image_view=DestroyCacheView(image_view);
 320   (void) CompositeImage(difference_image,highlight_image,image->compose,
 321     MagickTrue,0,0,exception);
 322   highlight_image=DestroyImage(highlight_image);
 323   if (status == MagickFalse)
 324     difference_image=DestroyImage(difference_image);
 325   return(difference_image);
 326 }
 327 \f
 328 /*
 329 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 330 %                                                                             %
 331 %                                                                             %
 332 %                                                                             %
 333 %   G e t I m a g e D i s t o r t i o n                                       %
 334 %                                                                             %
 335 %                                                                             %
 336 %                                                                             %
 337 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 338 %
 339 %  GetImageDistortion() compares one or more pixel channels of an image to a
 340 %  reconstructed image and returns the specified distortion metric.
 341 %
 342 %  The format of the GetImageDistortion method is:
 343 %
 344 %      MagickBooleanType GetImageDistortion(const Image *image,
 345 %        const Image *reconstruct_image,const MetricType metric,
 346 %        double *distortion,ExceptionInfo *exception)
 347 %
 348 %  A description of each parameter follows:
 349 %
 350 %    o image: the image.
 351 %
 352 %    o reconstruct_image: the reconstruct image.
 353 %
 354 %    o metric: the metric.
 355 %
 356 %    o distortion: the computed distortion between the images.
 357 %
 358 %    o exception: return any errors or warnings in this structure.
 359 %
 360 */
 361
 362 static MagickBooleanType GetAbsoluteDistortion(const Image *image,
 363   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
 364 {
 365   CacheView
 366     *image_view,
 367     *reconstruct_view;
 368
 369   double
 370     fuzz;
 371
 372   MagickBooleanType
 373     status;
 374
 375   size_t
 376     columns,
 377     rows;
 378
 379   ssize_t
 380     y;
 381
 382   /*
 383     Compute the absolute difference in pixels between two images.
 384   */
 385   status=MagickTrue;
 386   fuzz=(double) MagickMin(GetPixelChannels(image),
 387     GetPixelChannels(reconstruct_image))*
 388     GetFuzzyColorDistance(image,reconstruct_image);
 389   rows=MagickMax(image->rows,reconstruct_image->rows);
 390   columns=MagickMax(image->columns,reconstruct_image->columns);
 391   image_view=AcquireVirtualCacheView(image,exception);
 392   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
 393 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 394   #pragma omp parallel for schedule(static) shared(status) \
 395     magick_number_threads(image,image,rows,1)
 396 #endif
 397   for (y=0; y < (ssize_t) rows; y++)
 398   {
 399     double
 400       channel_distortion[MaxPixelChannels+1];
 401
 402     register const Quantum
 403       *magick_restrict p,
 404       *magick_restrict q;
 405
 406     register ssize_t
 407       j,
 408       x;
 409
 410     if (status == MagickFalse)
 411       continue;
 412     p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
 413     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
 414     if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
 415       {
 416         status=MagickFalse;
 417         continue;
 418       }
 419     (void) memset(channel_distortion,0,sizeof(channel_distortion));
 420     for (x=0; x < (ssize_t) columns; x++)
 421     {
 422       double
 423         Da,
 424         distance,
 425         Sa;
 426
 427       MagickBooleanType
 428         difference;
 429
 430       register ssize_t
 431         i;
 432
 433       if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
 434           (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
 435         {
 436           p+=GetPixelChannels(image);
 437           q+=GetPixelChannels(reconstruct_image);
 438           continue;
 439         }
 440       difference=MagickFalse;
 441       distance=0.0;
 442       Sa=QuantumScale*GetPixelAlpha(image,p);
 443       Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
 444       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 445       {
 446         double
 447           pixel;
 448
 449         PixelChannel channel = GetPixelChannelChannel(image,i);
 450         PixelTrait traits = GetPixelChannelTraits(image,channel);
 451         PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
 452           channel);
 453         if ((traits == UndefinedPixelTrait) ||
 454             (reconstruct_traits == UndefinedPixelTrait) ||
 455             ((reconstruct_traits & UpdatePixelTrait) == 0))
 456           continue;
 457         if (channel == AlphaPixelChannel)
 458           pixel=(double) p[i]-GetPixelChannel(reconstruct_image,channel,q);
 459         else
 460           pixel=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
 461         distance+=pixel*pixel;
 462         if (distance > fuzz)
 463           {
 464             channel_distortion[i]++;
 465             difference=MagickTrue;
 466           }
 467       }
 468       if (difference != MagickFalse)
 469         channel_distortion[CompositePixelChannel]++;
 470       p+=GetPixelChannels(image);
 471       q+=GetPixelChannels(reconstruct_image);
 472     }
 473 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 474     #pragma omp critical (MagickCore_GetAbsoluteDistortion)
 475 #endif
 476     for (j=0; j <= MaxPixelChannels; j++)
 477       distortion[j]+=channel_distortion[j];
 478   }
 479   reconstruct_view=DestroyCacheView(reconstruct_view);
 480   image_view=DestroyCacheView(image_view);
 481   return(status);
 482 }
 483
 484 static MagickBooleanType GetFuzzDistortion(const Image *image,
 485   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
 486 {
 487   CacheView
 488     *image_view,
 489     *reconstruct_view;
 490
 491   double
 492     area;
 493
 494   MagickBooleanType
 495     status;
 496
 497   register ssize_t
 498     j;
 499
 500   size_t
 501     columns,
 502     rows;
 503
 504   ssize_t
 505     y;
 506
 507   status=MagickTrue;
 508   rows=MagickMax(image->rows,reconstruct_image->rows);
 509   columns=MagickMax(image->columns,reconstruct_image->columns);
 510   area=0.0;
 511   image_view=AcquireVirtualCacheView(image,exception);
 512   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
 513 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 514   #pragma omp parallel for schedule(static) shared(status) \
 515     magick_number_threads(image,image,rows,1) reduction(+:area)
 516 #endif
 517   for (y=0; y < (ssize_t) rows; y++)
 518   {
 519     double
 520       channel_distortion[MaxPixelChannels+1];
 521
 522     register const Quantum
 523       *magick_restrict p,
 524       *magick_restrict q;
 525
 526     register ssize_t
 527       x;
 528
 529     if (status == MagickFalse)
 530       continue;
 531     p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
 532     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
 533     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
 534       {
 535         status=MagickFalse;
 536         continue;
 537       }
 538     (void) memset(channel_distortion,0,sizeof(channel_distortion));
 539     for (x=0; x < (ssize_t) columns; x++)
 540     {
 541       double
 542         Da,
 543         Sa;
 544
 545       register ssize_t
 546         i;
 547
 548       if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
 549           (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
 550         {
 551           p+=GetPixelChannels(image);
 552           q+=GetPixelChannels(reconstruct_image);
 553           continue;
 554         }
 555       Sa=QuantumScale*GetPixelAlpha(image,p);
 556       Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
 557       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 558       {
 559         double
 560           distance;
 561
 562         PixelChannel channel = GetPixelChannelChannel(image,i);
 563         PixelTrait traits = GetPixelChannelTraits(image,channel);
 564         PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
 565           channel);
 566         if ((traits == UndefinedPixelTrait) ||
 567             (reconstruct_traits == UndefinedPixelTrait) ||
 568             ((reconstruct_traits & UpdatePixelTrait) == 0))
 569           continue;
 570         if (channel == AlphaPixelChannel)
 571           distance=QuantumScale*(p[i]-GetPixelChannel(reconstruct_image,
 572             channel,q));
 573         else
 574           distance=QuantumScale*(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
 575             channel,q));
 576         channel_distortion[i]+=distance*distance;
 577         channel_distortion[CompositePixelChannel]+=distance*distance;
 578       }
 579       area++;
 580       p+=GetPixelChannels(image);
 581       q+=GetPixelChannels(reconstruct_image);
 582     }
 583 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 584     #pragma omp critical (MagickCore_GetFuzzDistortion)
 585 #endif
 586     for (j=0; j <= MaxPixelChannels; j++)
 587       distortion[j]+=channel_distortion[j];
 588   }
 589   reconstruct_view=DestroyCacheView(reconstruct_view);
 590   image_view=DestroyCacheView(image_view);
 591   area=PerceptibleReciprocal(area);
 592   for (j=0; j <= MaxPixelChannels; j++)
 593     distortion[j]*=area;
 594   distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
 595   distortion[CompositePixelChannel]=sqrt(distortion[CompositePixelChannel]);
 596   return(status);
 597 }
 598
 599 static MagickBooleanType GetMeanAbsoluteDistortion(const Image *image,
 600   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
 601 {
 602   CacheView
 603     *image_view,
 604     *reconstruct_view;
 605
 606   double
 607     area;
 608
 609   MagickBooleanType
 610     status;
 611
 612   register ssize_t
 613     j;
 614
 615   size_t
 616     columns,
 617     rows;
 618
 619   ssize_t
 620     y;
 621
 622   status=MagickTrue;
 623   rows=MagickMax(image->rows,reconstruct_image->rows);
 624   columns=MagickMax(image->columns,reconstruct_image->columns);
 625   area=0.0;
 626   image_view=AcquireVirtualCacheView(image,exception);
 627   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
 628 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 629   #pragma omp parallel for schedule(static) shared(status) \
 630     magick_number_threads(image,image,rows,1) reduction(+:area)
 631 #endif
 632   for (y=0; y < (ssize_t) rows; y++)
 633   {
 634     double
 635       channel_distortion[MaxPixelChannels+1];
 636
 637     register const Quantum
 638       *magick_restrict p,
 639       *magick_restrict q;
 640
 641     register ssize_t
 642       x;
 643
 644     if (status == MagickFalse)
 645       continue;
 646     p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
 647     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
 648     if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
 649       {
 650         status=MagickFalse;
 651         continue;
 652       }
 653     (void) memset(channel_distortion,0,sizeof(channel_distortion));
 654     for (x=0; x < (ssize_t) columns; x++)
 655     {
 656       double
 657         Da,
 658         Sa;
 659
 660       register ssize_t
 661         i;
 662
 663       if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
 664           (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
 665         {
 666           p+=GetPixelChannels(image);
 667           q+=GetPixelChannels(reconstruct_image);
 668           continue;
 669         }
 670       Sa=QuantumScale*GetPixelAlpha(image,p);
 671       Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
 672       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 673       {
 674         double
 675           distance;
 676
 677         PixelChannel channel = GetPixelChannelChannel(image,i);
 678         PixelTrait traits = GetPixelChannelTraits(image,channel);
 679         PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
 680           channel);
 681         if ((traits == UndefinedPixelTrait) ||
 682             (reconstruct_traits == UndefinedPixelTrait) ||
 683             ((reconstruct_traits & UpdatePixelTrait) == 0))
 684           continue;
 685         if (channel == AlphaPixelChannel)
 686           distance=QuantumScale*fabs((double) p[i]-
 687             GetPixelChannel(reconstruct_image,channel,q));
 688         else
 689           distance=QuantumScale*fabs(Sa*p[i]-Da*
 690             GetPixelChannel(reconstruct_image,channel,q));
 691         channel_distortion[i]+=distance;
 692         channel_distortion[CompositePixelChannel]+=distance;
 693       }
 694       area++;
 695       p+=GetPixelChannels(image);
 696       q+=GetPixelChannels(reconstruct_image);
 697     }
 698 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 699     #pragma omp critical (MagickCore_GetMeanAbsoluteError)
 700 #endif
 701     for (j=0; j <= MaxPixelChannels; j++)
 702       distortion[j]+=channel_distortion[j];
 703   }
 704   reconstruct_view=DestroyCacheView(reconstruct_view);
 705   image_view=DestroyCacheView(image_view);
 706   area=PerceptibleReciprocal(area);
 707   for (j=0; j <= MaxPixelChannels; j++)
 708     distortion[j]*=area;
 709   distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
 710   return(status);
 711 }
 712
 713 static MagickBooleanType GetMeanErrorPerPixel(Image *image,
 714   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
 715 {
 716   CacheView
 717     *image_view,
 718     *reconstruct_view;
 719
 720   MagickBooleanType
 721     status;
 722
 723   double
 724     area,
 725     maximum_error,
 726     mean_error;
 727
 728   size_t
 729     columns,
 730     rows;
 731
 732   ssize_t
 733     y;
 734
 735   status=MagickTrue;
 736   area=0.0;
 737   maximum_error=0.0;
 738   mean_error=0.0;
 739   rows=MagickMax(image->rows,reconstruct_image->rows);
 740   columns=MagickMax(image->columns,reconstruct_image->columns);
 741   image_view=AcquireVirtualCacheView(image,exception);
 742   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
 743   for (y=0; y < (ssize_t) rows; y++)
 744   {
 745     register const Quantum
 746       *magick_restrict p,
 747       *magick_restrict q;
 748
 749     register ssize_t
 750       x;
 751
 752     p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
 753     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
 754     if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
 755       {
 756         status=MagickFalse;
 757         break;
 758       }
 759     for (x=0; x < (ssize_t) columns; x++)
 760     {
 761       double
 762         Da,
 763         Sa;
 764
 765       register ssize_t
 766         i;
 767
 768       if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
 769           (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
 770         {
 771           p+=GetPixelChannels(image);
 772           q+=GetPixelChannels(reconstruct_image);
 773           continue;
 774         }
 775       Sa=QuantumScale*GetPixelAlpha(image,p);
 776       Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
 777       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 778       {
 779         double
 780           distance;
 781
 782         PixelChannel channel = GetPixelChannelChannel(image,i);
 783         PixelTrait traits = GetPixelChannelTraits(image,channel);
 784         PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
 785           channel);
 786         if ((traits == UndefinedPixelTrait) ||
 787             (reconstruct_traits == UndefinedPixelTrait) ||
 788             ((reconstruct_traits & UpdatePixelTrait) == 0))
 789           continue;
 790         if (channel == AlphaPixelChannel)
 791           distance=fabs((double) p[i]-
 792             GetPixelChannel(reconstruct_image,channel,q));
 793         else
 794           distance=fabs(Sa*p[i]-Da*
 795             GetPixelChannel(reconstruct_image,channel,q));
 796         distortion[i]+=distance;
 797         distortion[CompositePixelChannel]+=distance;
 798         mean_error+=distance*distance;
 799         if (distance > maximum_error)
 800           maximum_error=distance;
 801         area++;
 802       }
 803       p+=GetPixelChannels(image);
 804       q+=GetPixelChannels(reconstruct_image);
 805     }
 806   }
 807   reconstruct_view=DestroyCacheView(reconstruct_view);
 808   image_view=DestroyCacheView(image_view);
 809   image->error.mean_error_per_pixel=distortion[CompositePixelChannel]/area;
 810   image->error.normalized_mean_error=QuantumScale*QuantumScale*mean_error/area;
 811   image->error.normalized_maximum_error=QuantumScale*maximum_error;
 812   return(status);
 813 }
 814
 815 static MagickBooleanType GetMeanSquaredDistortion(const Image *image,
 816   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
 817 {
 818   CacheView
 819     *image_view,
 820     *reconstruct_view;
 821
 822   double
 823     area;
 824
 825   MagickBooleanType
 826     status;
 827
 828   register ssize_t
 829     j;
 830
 831   size_t
 832     columns,
 833     rows;
 834
 835   ssize_t
 836     y;
 837
 838   status=MagickTrue;
 839   rows=MagickMax(image->rows,reconstruct_image->rows);
 840   columns=MagickMax(image->columns,reconstruct_image->columns);
 841   area=0.0;
 842   image_view=AcquireVirtualCacheView(image,exception);
 843   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
 844 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 845   #pragma omp parallel for schedule(static) shared(status) \
 846     magick_number_threads(image,image,rows,1) reduction(+:area)
 847 #endif
 848   for (y=0; y < (ssize_t) rows; y++)
 849   {
 850     double
 851       channel_distortion[MaxPixelChannels+1];
 852
 853     register const Quantum
 854       *magick_restrict p,
 855       *magick_restrict q;
 856
 857     register ssize_t
 858       x;
 859
 860     if (status == MagickFalse)
 861       continue;
 862     p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
 863     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
 864     if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
 865       {
 866         status=MagickFalse;
 867         continue;
 868       }
 869     (void) memset(channel_distortion,0,sizeof(channel_distortion));
 870     for (x=0; x < (ssize_t) columns; x++)
 871     {
 872       double
 873         Da,
 874         Sa;
 875
 876       register ssize_t
 877         i;
 878
 879       if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
 880           (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
 881         {
 882           p+=GetPixelChannels(image);
 883           q+=GetPixelChannels(reconstruct_image);
 884           continue;
 885         }
 886       Sa=QuantumScale*GetPixelAlpha(image,p);
 887       Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
 888       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 889       {
 890         double
 891           distance;
 892
 893         PixelChannel channel = GetPixelChannelChannel(image,i);
 894         PixelTrait traits = GetPixelChannelTraits(image,channel);
 895         PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
 896           channel);
 897         if ((traits == UndefinedPixelTrait) ||
 898             (reconstruct_traits == UndefinedPixelTrait) ||
 899             ((reconstruct_traits & UpdatePixelTrait) == 0))
 900           continue;
 901         if (channel == AlphaPixelChannel)
 902           distance=QuantumScale*(p[i]-GetPixelChannel(reconstruct_image,
 903             channel,q));
 904         else
 905           distance=QuantumScale*(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
 906             channel,q));
 907         channel_distortion[i]+=distance*distance;
 908         channel_distortion[CompositePixelChannel]+=distance*distance;
 909       }
 910       area++;
 911       p+=GetPixelChannels(image);
 912       q+=GetPixelChannels(reconstruct_image);
 913     }
 914 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 915     #pragma omp critical (MagickCore_GetMeanSquaredError)
 916 #endif
 917     for (j=0; j <= MaxPixelChannels; j++)
 918       distortion[j]+=channel_distortion[j];
 919   }
 920   reconstruct_view=DestroyCacheView(reconstruct_view);
 921   image_view=DestroyCacheView(image_view);
 922   area=PerceptibleReciprocal(area);
 923   for (j=0; j <= MaxPixelChannels; j++)
 924     distortion[j]*=area;
 925   distortion[CompositePixelChannel]/=GetImageChannels(image);
 926   return(status);
 927 }
 928
 929 static MagickBooleanType GetNormalizedCrossCorrelationDistortion(
 930   const Image *image,const Image *reconstruct_image,double *distortion,
 931   ExceptionInfo *exception)
 932 {
 933 #define SimilarityImageTag  "Similarity/Image"
 934
 935   CacheView
 936     *image_view,
 937     *reconstruct_view;
 938
 939   ChannelStatistics
 940     *image_statistics,
 941     *reconstruct_statistics;
 942
 943   double
 944     area;
 945
 946   MagickBooleanType
 947     status;
 948
 949   MagickOffsetType
 950     progress;
 951
 952   register ssize_t
 953     i;
 954
 955   size_t
 956     columns,
 957     rows;
 958
 959   ssize_t
 960     y;
 961
 962   /*
 963     Normalize to account for variation due to lighting and exposure condition.
 964   */
 965   image_statistics=GetImageStatistics(image,exception);
 966   reconstruct_statistics=GetImageStatistics(reconstruct_image,exception);
 967   if ((image_statistics == (ChannelStatistics *) NULL) ||
 968       (reconstruct_statistics == (ChannelStatistics *) NULL))
 969     {
 970       if (image_statistics != (ChannelStatistics *) NULL)
 971         image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
 972           image_statistics);
 973       if (reconstruct_statistics != (ChannelStatistics *) NULL)
 974         reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
 975           reconstruct_statistics);
 976       return(MagickFalse);
 977     }
 978   status=MagickTrue;
 979   progress=0;
 980   for (i=0; i <= MaxPixelChannels; i++)
 981     distortion[i]=0.0;
 982   rows=MagickMax(image->rows,reconstruct_image->rows);
 983   columns=MagickMax(image->columns,reconstruct_image->columns);
 984   area=0.0;
 985   image_view=AcquireVirtualCacheView(image,exception);
 986   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
 987   for (y=0; y < (ssize_t) rows; y++)
 988   {
 989     register const Quantum
 990       *magick_restrict p,
 991       *magick_restrict q;
 992
 993     register ssize_t
 994       x;
 995
 996     p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
 997     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
 998     if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
 999       {
1000         status=MagickFalse;
1001         break;
1002       }
1003     for (x=0; x < (ssize_t) columns; x++)
1004     {
1005       if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1006           (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1007         {
1008           p+=GetPixelChannels(image);
1009           q+=GetPixelChannels(reconstruct_image);
1010           continue;
1011         }
1012       area++;
1013       p+=GetPixelChannels(image);
1014       q+=GetPixelChannels(reconstruct_image);
1015     }
1016   }
1017   area=PerceptibleReciprocal(area);
1018   for (y=0; y < (ssize_t) rows; y++)
1019   {
1020     register const Quantum
1021       *magick_restrict p,
1022       *magick_restrict q;
1023
1024     register ssize_t
1025       x;
1026
1027     p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1028     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1029     if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1030       {
1031         status=MagickFalse;
1032         break;
1033       }
1034     for (x=0; x < (ssize_t) columns; x++)
1035     {
1036       double
1037         Da,
1038         Sa;
1039
1040       if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1041           (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1042         {
1043           p+=GetPixelChannels(image);
1044           q+=GetPixelChannels(reconstruct_image);
1045           continue;
1046         }
1047       Sa=QuantumScale*GetPixelAlpha(image,p);
1048       Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
1049       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1050       {
1051         PixelChannel channel = GetPixelChannelChannel(image,i);
1052         PixelTrait traits = GetPixelChannelTraits(image,channel);
1053         PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1054           channel);
1055         if ((traits == UndefinedPixelTrait) ||
1056             (reconstruct_traits == UndefinedPixelTrait) ||
1057             ((reconstruct_traits & UpdatePixelTrait) == 0))
1058           continue;
1059         if (channel == AlphaPixelChannel)
1060           {
1061             distortion[i]+=area*QuantumScale*(p[i]-
1062               image_statistics[channel].mean)*(GetPixelChannel(
1063               reconstruct_image,channel,q)-
1064               reconstruct_statistics[channel].mean);
1065           }
1066         else
1067           {
1068             distortion[i]+=area*QuantumScale*(Sa*p[i]-
1069               image_statistics[channel].mean)*(Da*GetPixelChannel(
1070               reconstruct_image,channel,q)-
1071               reconstruct_statistics[channel].mean);
1072           }
1073       }
1074       p+=GetPixelChannels(image);
1075       q+=GetPixelChannels(reconstruct_image);
1076     }
1077     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1078       {
1079         MagickBooleanType
1080           proceed;
1081
1082 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1083         #pragma omp atomic
1084 #endif
1085         progress++;
1086         proceed=SetImageProgress(image,SimilarityImageTag,progress,rows);
1087         if (proceed == MagickFalse)
1088           {
1089             status=MagickFalse;
1090             break;
1091           }
1092       }
1093   }
1094   reconstruct_view=DestroyCacheView(reconstruct_view);
1095   image_view=DestroyCacheView(image_view);
1096   /*
1097     Divide by the standard deviation.
1098   */
1099   distortion[CompositePixelChannel]=0.0;
1100   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1101   {
1102     double
1103       gamma;
1104
1105     PixelChannel channel = GetPixelChannelChannel(image,i);
1106     gamma=image_statistics[channel].standard_deviation*
1107       reconstruct_statistics[channel].standard_deviation;
1108     gamma=PerceptibleReciprocal(gamma);
1109     distortion[i]=QuantumRange*gamma*distortion[i];
1110     distortion[CompositePixelChannel]+=distortion[i]*distortion[i];
1111   }
1112   distortion[CompositePixelChannel]=sqrt(distortion[CompositePixelChannel]/
1113     GetImageChannels(image));
1114   /*
1115     Free resources.
1116   */
1117   reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1118     reconstruct_statistics);
1119   image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1120     image_statistics);
1121   return(status);
1122 }
1123
1124 static MagickBooleanType GetPeakAbsoluteDistortion(const Image *image,
1125   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1126 {
1127   CacheView
1128     *image_view,
1129     *reconstruct_view;
1130
1131   MagickBooleanType
1132     status;
1133
1134   size_t
1135     columns,
1136     rows;
1137
1138   ssize_t
1139     y;
1140
1141   status=MagickTrue;
1142   rows=MagickMax(image->rows,reconstruct_image->rows);
1143   columns=MagickMax(image->columns,reconstruct_image->columns);
1144   image_view=AcquireVirtualCacheView(image,exception);
1145   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1146 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1147   #pragma omp parallel for schedule(static) shared(status) \
1148     magick_number_threads(image,image,rows,1)
1149 #endif
1150   for (y=0; y < (ssize_t) rows; y++)
1151   {
1152     double
1153       channel_distortion[MaxPixelChannels+1];
1154
1155     register const Quantum
1156       *magick_restrict p,
1157       *magick_restrict q;
1158
1159     register ssize_t
1160       j,
1161       x;
1162
1163     if (status == MagickFalse)
1164       continue;
1165     p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1166     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1167     if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1168       {
1169         status=MagickFalse;
1170         continue;
1171       }
1172     (void) memset(channel_distortion,0,sizeof(channel_distortion));
1173     for (x=0; x < (ssize_t) columns; x++)
1174     {
1175       double
1176         Da,
1177         Sa;
1178
1179       register ssize_t
1180         i;
1181
1182       if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1183           (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1184         {
1185           p+=GetPixelChannels(image);
1186           q+=GetPixelChannels(reconstruct_image);
1187           continue;
1188         }
1189       Sa=QuantumScale*GetPixelAlpha(image,p);
1190       Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
1191       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1192       {
1193         double
1194           distance;
1195
1196         PixelChannel channel = GetPixelChannelChannel(image,i);
1197         PixelTrait traits = GetPixelChannelTraits(image,channel);
1198         PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1199           channel);
1200         if ((traits == UndefinedPixelTrait) ||
1201             (reconstruct_traits == UndefinedPixelTrait) ||
1202             ((reconstruct_traits & UpdatePixelTrait) == 0))
1203           continue;
1204         if (channel == AlphaPixelChannel)
1205           distance=QuantumScale*fabs((double) p[i]-
1206             GetPixelChannel(reconstruct_image,channel,q));
1207         else
1208           distance=QuantumScale*fabs(Sa*p[i]-Da*
1209             GetPixelChannel(reconstruct_image,channel,q));
1210         if (distance > channel_distortion[i])
1211           channel_distortion[i]=distance;
1212         if (distance > channel_distortion[CompositePixelChannel])
1213           channel_distortion[CompositePixelChannel]=distance;
1214       }
1215       p+=GetPixelChannels(image);
1216       q+=GetPixelChannels(reconstruct_image);
1217     }
1218 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1219     #pragma omp critical (MagickCore_GetPeakAbsoluteError)
1220 #endif
1221     for (j=0; j <= MaxPixelChannels; j++)
1222       if (channel_distortion[j] > distortion[j])
1223         distortion[j]=channel_distortion[j];
1224   }
1225   reconstruct_view=DestroyCacheView(reconstruct_view);
1226   image_view=DestroyCacheView(image_view);
1227   return(status);
1228 }
1229
1230 static inline double MagickLog10(const double x)
1231 {
1232 #define Log10Epsilon  (1.0e-11)
1233
1234  if (fabs(x) < Log10Epsilon)
1235    return(log10(Log10Epsilon));
1236  return(log10(fabs(x)));
1237 }
1238
1239 static MagickBooleanType GetPeakSignalToNoiseRatio(const Image *image,
1240   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1241 {
1242   MagickBooleanType
1243     status;
1244
1245   register ssize_t
1246     i;
1247
1248   status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1249   for (i=0; i <= MaxPixelChannels; i++)
1250     if (fabs(distortion[i]) < MagickEpsilon)
1251       distortion[i]=INFINITY;
1252     else
1253       distortion[i]=10.0*MagickLog10(1.0)-10.0*MagickLog10(distortion[i]);
1254   return(status);
1255 }
1256
1257 static MagickBooleanType GetPerceptualHashDistortion(const Image *image,
1258   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1259 {
1260   ChannelPerceptualHash
1261     *channel_phash,
1262     *reconstruct_phash;
1263
1264   const char
1265     *artifact;
1266
1267   MagickBooleanType
1268     normalize;
1269
1270   ssize_t
1271     channel;
1272
1273   /*
1274     Compute perceptual hash in the sRGB colorspace.
1275   */
1276   channel_phash=GetImagePerceptualHash(image,exception);
1277   if (channel_phash == (ChannelPerceptualHash *) NULL)
1278     return(MagickFalse);
1279   reconstruct_phash=GetImagePerceptualHash(reconstruct_image,exception);
1280   if (reconstruct_phash == (ChannelPerceptualHash *) NULL)
1281     {
1282       channel_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1283         channel_phash);
1284       return(MagickFalse);
1285     }
1286   artifact=GetImageArtifact(image,"phash:normalize");
1287   normalize=(artifact == (const char *) NULL) ||
1288     (IsStringTrue(artifact) == MagickFalse) ? MagickFalse : MagickTrue;
1289 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1290   #pragma omp parallel for schedule(static)
1291 #endif
1292   for (channel=0; channel < MaxPixelChannels; channel++)
1293   {
1294     double
1295       difference;
1296
1297     register ssize_t
1298       i;
1299
1300     difference=0.0;
1301     for (i=0; i < MaximumNumberOfImageMoments; i++)
1302     {
1303       double
1304         alpha,
1305         beta;
1306
1307       register ssize_t
1308         j;
1309
1310       for (j=0; j < (ssize_t) channel_phash[0].number_colorspaces; j++)
1311       {
1312         alpha=channel_phash[channel].phash[j][i];
1313         beta=reconstruct_phash[channel].phash[j][i];
1314         if (normalize == MagickFalse)
1315           difference+=(beta-alpha)*(beta-alpha);
1316         else
1317           difference=sqrt((beta-alpha)*(beta-alpha)/
1318             channel_phash[0].number_channels);
1319       }
1320     }
1321     distortion[channel]+=difference;
1322 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1323     #pragma omp critical (MagickCore_GetPerceptualHashDistortion)
1324 #endif
1325     distortion[CompositePixelChannel]+=difference;
1326   }
1327   /*
1328     Free resources.
1329   */
1330   reconstruct_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1331     reconstruct_phash);
1332   channel_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(channel_phash);
1333   return(MagickTrue);
1334 }
1335
1336 static MagickBooleanType GetRootMeanSquaredDistortion(const Image *image,
1337   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1338 {
1339   MagickBooleanType
1340     status;
1341
1342   register ssize_t
1343     i;
1344
1345   status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1346   for (i=0; i <= MaxPixelChannels; i++)
1347     distortion[i]=sqrt(distortion[i]);
1348   return(status);
1349 }
1350
1351 static MagickBooleanType GetStructuralSimilarityDistortion(const Image *image,
1352   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1353 {
1354 #define SSIMRadius  5.0
1355 #define SSIMSigma  1.5
1356 #define SSIMBlocksize  8
1357 #define SSIMK1  0.01
1358 #define SSIMK2  0.03
1359 #define SSIML  1.0
1360
1361   CacheView
1362     *image_view,
1363     *reconstruct_view;
1364
1365   char
1366     geometry[MagickPathExtent];
1367
1368   const char
1369     *artifact;
1370
1371   double
1372     c1,
1373     c2,
1374     radius,
1375     sigma;
1376
1377   KernelInfo
1378     *kernel_info;
1379
1380   MagickBooleanType
1381     status;
1382
1383   register ssize_t
1384     i;
1385
1386   size_t
1387     columns,
1388     rows;
1389
1390   ssize_t
1391     y;
1392
1393   /*
1394     Compute structural similarity index @
1395     https://en.wikipedia.org/wiki/Structural_similarity.
1396   */
1397   radius=SSIMRadius;
1398   artifact=GetImageArtifact(image,"compare:ssim-radius");
1399   if (artifact != (const char *) NULL)
1400     radius=StringToDouble(artifact,(char **) NULL);
1401   sigma=SSIMSigma;
1402   artifact=GetImageArtifact(image,"compare:ssim-sigma");
1403   if (artifact != (const char *) NULL)
1404     sigma=StringToDouble(artifact,(char **) NULL);
1405   (void) FormatLocaleString(geometry,MagickPathExtent,"gaussian:%.20gx%.20g",
1406     radius,sigma);
1407   kernel_info=AcquireKernelInfo(geometry,exception);
1408   if (kernel_info == (KernelInfo *) NULL)
1409     ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
1410       image->filename);
1411   c1=pow(SSIMK1*SSIML,2.0);
1412   artifact=GetImageArtifact(image,"compare:ssim-k1");
1413   if (artifact != (const char *) NULL)
1414     c1=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
1415   c2=pow(SSIMK2*SSIML,2.0);
1416   artifact=GetImageArtifact(image,"compare:ssim-k2");
1417   if (artifact != (const char *) NULL)
1418     c2=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
1419   status=MagickTrue;
1420   rows=MagickMax(image->rows,reconstruct_image->rows);
1421   columns=MagickMax(image->columns,reconstruct_image->columns);
1422   image_view=AcquireVirtualCacheView(image,exception);
1423   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1424 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1425   #pragma omp parallel for schedule(static) shared(status) \
1426     magick_number_threads(image,reconstruct_image,rows,1)
1427 #endif
1428   for (y=0; y < (ssize_t) rows; y++)
1429   {
1430     double
1431       channel_distortion[MaxPixelChannels+1];
1432
1433     register const Quantum
1434       *magick_restrict p,
1435       *magick_restrict q;
1436
1437     register ssize_t
1438       i,
1439       x;
1440
1441     if (status == MagickFalse)
1442       continue;
1443     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) kernel_info->width/2L),y-
1444       ((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
1445       kernel_info->height,exception);
1446     q=GetCacheViewVirtualPixels(reconstruct_view,-((ssize_t) kernel_info->width/
1447       2L),y-((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
1448       kernel_info->height,exception);
1449     if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1450       {
1451         status=MagickFalse;
1452         continue;
1453       }
1454     (void) memset(channel_distortion,0,sizeof(channel_distortion));
1455     for (x=0; x < (ssize_t) columns; x++)
1456     {
1457       double
1458         x_pixel_mu[MaxPixelChannels+1],
1459         x_pixel_sigma_squared[MaxPixelChannels+1],
1460         xy_sigma[MaxPixelChannels+1],
1461         y_pixel_mu[MaxPixelChannels+1],
1462         y_pixel_sigma_squared[MaxPixelChannels+1];
1463
1464       register const Quantum
1465         *magick_restrict reference,
1466         *magick_restrict target;
1467
1468       register MagickRealType
1469         *k;
1470
1471       ssize_t
1472         v;
1473
1474       (void) memset(x_pixel_mu,0,sizeof(x_pixel_mu));
1475       (void) memset(x_pixel_sigma_squared,0,sizeof(x_pixel_sigma_squared));
1476       (void) memset(xy_sigma,0,sizeof(xy_sigma));
1477       (void) memset(x_pixel_sigma_squared,0,sizeof(y_pixel_sigma_squared));
1478       (void) memset(y_pixel_mu,0,sizeof(y_pixel_mu));
1479       (void) memset(y_pixel_sigma_squared,0,sizeof(y_pixel_sigma_squared));
1480       k=kernel_info->values;
1481       reference=p;
1482       target=q;
1483       for (v=0; v < (ssize_t) kernel_info->height; v++)
1484       {
1485         register ssize_t
1486           u;
1487
1488         for (u=0; u < (ssize_t) kernel_info->width; u++)
1489         {
1490           for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1491           {
1492             double
1493               x_pixel,
1494               y_pixel;
1495
1496             PixelChannel channel = GetPixelChannelChannel(image,i);
1497             PixelTrait traits = GetPixelChannelTraits(image,channel);
1498             PixelTrait reconstruct_traits = GetPixelChannelTraits(
1499               reconstruct_image,channel);
1500             if ((traits == UndefinedPixelTrait) ||
1501                 (reconstruct_traits == UndefinedPixelTrait) ||
1502                 ((reconstruct_traits & UpdatePixelTrait) == 0))
1503               continue;
1504             x_pixel=QuantumScale*reference[i];
1505             x_pixel_mu[i]+=(*k)*x_pixel;
1506             x_pixel_sigma_squared[i]+=(*k)*x_pixel*x_pixel;
1507             y_pixel=QuantumScale*
1508               GetPixelChannel(reconstruct_image,channel,target);
1509             y_pixel_mu[i]+=(*k)*y_pixel;
1510             y_pixel_sigma_squared[i]+=(*k)*y_pixel*y_pixel;
1511             xy_sigma[i]+=(*k)*x_pixel*y_pixel;
1512           }
1513           k++;
1514           reference+=GetPixelChannels(image);
1515           target+=GetPixelChannels(reconstruct_image);
1516         }
1517         reference+=GetPixelChannels(image)*columns;
1518         target+=GetPixelChannels(reconstruct_image)*columns;
1519       }
1520       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1521       {
1522         double
1523           ssim,
1524           x_pixel_mu_squared,
1525           x_pixel_sigmas_squared,
1526           xy_mu,
1527           xy_sigmas,
1528           y_pixel_mu_squared,
1529           y_pixel_sigmas_squared;
1530
1531         PixelChannel channel = GetPixelChannelChannel(image,i);
1532         PixelTrait traits = GetPixelChannelTraits(image,channel);
1533         PixelTrait reconstruct_traits = GetPixelChannelTraits(
1534           reconstruct_image,channel);
1535         if ((traits == UndefinedPixelTrait) ||
1536             (reconstruct_traits == UndefinedPixelTrait) ||
1537             ((reconstruct_traits & UpdatePixelTrait) == 0))
1538           continue;
1539         x_pixel_mu_squared=x_pixel_mu[i]*x_pixel_mu[i];
1540         y_pixel_mu_squared=y_pixel_mu[i]*y_pixel_mu[i];
1541         xy_mu=x_pixel_mu[i]*y_pixel_mu[i];
1542         xy_sigmas=xy_sigma[i]-xy_mu;
1543         x_pixel_sigmas_squared=x_pixel_sigma_squared[i]-x_pixel_mu_squared;
1544         y_pixel_sigmas_squared=y_pixel_sigma_squared[i]-y_pixel_mu_squared;
1545         ssim=((2.0*xy_mu+c1)*(2.0*xy_sigmas+c2))/
1546           ((x_pixel_mu_squared+y_pixel_mu_squared+c1)*
1547            (x_pixel_sigmas_squared+y_pixel_sigmas_squared+c2));
1548         channel_distortion[i]+=ssim;
1549         channel_distortion[CompositePixelChannel]+=ssim;
1550       }
1551       p+=GetPixelChannels(image);
1552       q+=GetPixelChannels(reconstruct_image);
1553     }
1554 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1555     #pragma omp critical (MagickCore_GetStructuralSimilarityDistortion)
1556 #endif
1557     for (i=0; i <= MaxPixelChannels; i++)
1558       distortion[i]+=channel_distortion[i];
1559   }
1560   image_view=DestroyCacheView(image_view);
1561   reconstruct_view=DestroyCacheView(reconstruct_view);
1562   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1563   {
1564     PixelChannel channel = GetPixelChannelChannel(image,i);
1565     PixelTrait traits = GetPixelChannelTraits(image,channel);
1566     if ((traits == UndefinedPixelTrait) || ((traits & UpdatePixelTrait) == 0))
1567       continue;
1568     distortion[i]/=((double) columns*rows);
1569   }
1570   distortion[CompositePixelChannel]/=((double) columns*rows);
1571   distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
1572   kernel_info=DestroyKernelInfo(kernel_info);
1573   return(status);
1574 }
1575
1576 static MagickBooleanType GetStructuralDisimilarityDistortion(const Image *image,
1577   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1578 {
1579   MagickBooleanType
1580     status;
1581
1582   register ssize_t
1583     i;
1584
1585   status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1586     distortion,exception);
1587   for (i=0; i <= MaxPixelChannels; i++)
1588     distortion[i]=(1.0-(distortion[i]))/2.0;
1589   return(status);
1590 }
1591
1592 MagickExport MagickBooleanType GetImageDistortion(Image *image,
1593   const Image *reconstruct_image,const MetricType metric,double *distortion,
1594   ExceptionInfo *exception)
1595 {
1596   double
1597     *channel_distortion;
1598
1599   MagickBooleanType
1600     status;
1601
1602   size_t
1603     length;
1604
1605   assert(image != (Image *) NULL);
1606   assert(image->signature == MagickCoreSignature);
1607   if (image->debug != MagickFalse)
1608     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1609   assert(reconstruct_image != (const Image *) NULL);
1610   assert(reconstruct_image->signature == MagickCoreSignature);
1611   assert(distortion != (double *) NULL);
1612   *distortion=0.0;
1613   if (image->debug != MagickFalse)
1614     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1615   /*
1616     Get image distortion.
1617   */
1618   length=MaxPixelChannels+1;
1619   channel_distortion=(double *) AcquireQuantumMemory(length,
1620     sizeof(*channel_distortion));
1621   if (channel_distortion == (double *) NULL)
1622     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1623   (void) memset(channel_distortion,0,length*
1624     sizeof(*channel_distortion));
1625   switch (metric)
1626   {
1627     case AbsoluteErrorMetric:
1628     {
1629       status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1630         exception);
1631       break;
1632     }
1633     case FuzzErrorMetric:
1634     {
1635       status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1636         exception);
1637       break;
1638     }
1639     case MeanAbsoluteErrorMetric:
1640     {
1641       status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1642         channel_distortion,exception);
1643       break;
1644     }
1645     case MeanErrorPerPixelErrorMetric:
1646     {
1647       status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1648         exception);
1649       break;
1650     }
1651     case MeanSquaredErrorMetric:
1652     {
1653       status=GetMeanSquaredDistortion(image,reconstruct_image,
1654         channel_distortion,exception);
1655       break;
1656     }
1657     case NormalizedCrossCorrelationErrorMetric:
1658     default:
1659     {
1660       status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1661         channel_distortion,exception);
1662       break;
1663     }
1664     case PeakAbsoluteErrorMetric:
1665     {
1666       status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1667         channel_distortion,exception);
1668       break;
1669     }
1670     case PeakSignalToNoiseRatioErrorMetric:
1671     {
1672       status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1673         channel_distortion,exception);
1674       break;
1675     }
1676     case PerceptualHashErrorMetric:
1677     {
1678       status=GetPerceptualHashDistortion(image,reconstruct_image,
1679         channel_distortion,exception);
1680       break;
1681     }
1682     case RootMeanSquaredErrorMetric:
1683     {
1684       status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1685         channel_distortion,exception);
1686       break;
1687     }
1688     case StructuralSimilarityErrorMetric:
1689     {
1690       status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1691         channel_distortion,exception);
1692       break;
1693     }
1694     case StructuralDissimilarityErrorMetric:
1695     {
1696       status=GetStructuralDisimilarityDistortion(image,reconstruct_image,
1697         channel_distortion,exception);
1698       break;
1699     }
1700   }
1701   *distortion=channel_distortion[CompositePixelChannel];
1702   channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1703   (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
1704     *distortion);
1705   return(status);
1706 }
1707 \f
1708 /*
1709 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1710 %                                                                             %
1711 %                                                                             %
1712 %                                                                             %
1713 %   G e t I m a g e D i s t o r t i o n s                                     %
1714 %                                                                             %
1715 %                                                                             %
1716 %                                                                             %
1717 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1718 %
1719 %  GetImageDistortions() compares the pixel channels of an image to a
1720 %  reconstructed image and returns the specified distortion metric for each
1721 %  channel.
1722 %
1723 %  The format of the GetImageDistortions method is:
1724 %
1725 %      double *GetImageDistortions(const Image *image,
1726 %        const Image *reconstruct_image,const MetricType metric,
1727 %        ExceptionInfo *exception)
1728 %
1729 %  A description of each parameter follows:
1730 %
1731 %    o image: the image.
1732 %
1733 %    o reconstruct_image: the reconstruct image.
1734 %
1735 %    o metric: the metric.
1736 %
1737 %    o exception: return any errors or warnings in this structure.
1738 %
1739 */
1740 MagickExport double *GetImageDistortions(Image *image,
1741   const Image *reconstruct_image,const MetricType metric,
1742   ExceptionInfo *exception)
1743 {
1744   double
1745     *channel_distortion;
1746
1747   MagickBooleanType
1748     status;
1749
1750   size_t
1751     length;
1752
1753   assert(image != (Image *) NULL);
1754   assert(image->signature == MagickCoreSignature);
1755   if (image->debug != MagickFalse)
1756     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1757   assert(reconstruct_image != (const Image *) NULL);
1758   assert(reconstruct_image->signature == MagickCoreSignature);
1759   if (image->debug != MagickFalse)
1760     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1761   /*
1762     Get image distortion.
1763   */
1764   length=MaxPixelChannels+1UL;
1765   channel_distortion=(double *) AcquireQuantumMemory(length,
1766     sizeof(*channel_distortion));
1767   if (channel_distortion == (double *) NULL)
1768     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1769   (void) memset(channel_distortion,0,length*
1770     sizeof(*channel_distortion));
1771   status=MagickTrue;
1772   switch (metric)
1773   {
1774     case AbsoluteErrorMetric:
1775     {
1776       status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1777         exception);
1778       break;
1779     }
1780     case FuzzErrorMetric:
1781     {
1782       status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1783         exception);
1784       break;
1785     }
1786     case MeanAbsoluteErrorMetric:
1787     {
1788       status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1789         channel_distortion,exception);
1790       break;
1791     }
1792     case MeanErrorPerPixelErrorMetric:
1793     {
1794       status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1795         exception);
1796       break;
1797     }
1798     case MeanSquaredErrorMetric:
1799     {
1800       status=GetMeanSquaredDistortion(image,reconstruct_image,
1801         channel_distortion,exception);
1802       break;
1803     }
1804     case NormalizedCrossCorrelationErrorMetric:
1805     default:
1806     {
1807       status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1808         channel_distortion,exception);
1809       break;
1810     }
1811     case PeakAbsoluteErrorMetric:
1812     {
1813       status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1814         channel_distortion,exception);
1815       break;
1816     }
1817     case PeakSignalToNoiseRatioErrorMetric:
1818     {
1819       status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1820         channel_distortion,exception);
1821       break;
1822     }
1823     case PerceptualHashErrorMetric:
1824     {
1825       status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1826         channel_distortion,exception);
1827       break;
1828     }
1829     case RootMeanSquaredErrorMetric:
1830     {
1831       status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1832         channel_distortion,exception);
1833       break;
1834     }
1835     case StructuralSimilarityErrorMetric:
1836     {
1837       status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1838         channel_distortion,exception);
1839       break;
1840     }
1841     case StructuralDissimilarityErrorMetric:
1842     {
1843       status=GetStructuralDisimilarityDistortion(image,reconstruct_image,
1844         channel_distortion,exception);
1845       break;
1846     }
1847   }
1848   if (status == MagickFalse)
1849     {
1850       channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1851       return((double *) NULL);
1852     }
1853   return(channel_distortion);
1854 }
1855 \f
1856 /*
1857 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1858 %                                                                             %
1859 %                                                                             %
1860 %                                                                             %
1861 %  I s I m a g e s E q u a l                                                  %
1862 %                                                                             %
1863 %                                                                             %
1864 %                                                                             %
1865 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1866 %
1867 %  IsImagesEqual() compare the pixels of two images and returns immediately
1868 %  if any pixel is not identical.
1869 %
1870 %  The format of the IsImagesEqual method is:
1871 %
1872 %      MagickBooleanType IsImagesEqual(const Image *image,
1873 %        const Image *reconstruct_image,ExceptionInfo *exception)
1874 %
1875 %  A description of each parameter follows.
1876 %
1877 %    o image: the image.
1878 %
1879 %    o reconstruct_image: the reconstruct image.
1880 %
1881 %    o exception: return any errors or warnings in this structure.
1882 %
1883 */
1884 MagickExport MagickBooleanType IsImagesEqual(const Image *image,
1885   const Image *reconstruct_image,ExceptionInfo *exception)
1886 {
1887   CacheView
1888     *image_view,
1889     *reconstruct_view;
1890
1891   size_t
1892     columns,
1893     rows;
1894
1895   ssize_t
1896     y;
1897
1898   assert(image != (Image *) NULL);
1899   assert(image->signature == MagickCoreSignature);
1900   assert(reconstruct_image != (const Image *) NULL);
1901   assert(reconstruct_image->signature == MagickCoreSignature);
1902   rows=MagickMax(image->rows,reconstruct_image->rows);
1903   columns=MagickMax(image->columns,reconstruct_image->columns);
1904   image_view=AcquireVirtualCacheView(image,exception);
1905   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1906   for (y=0; y < (ssize_t) rows; y++)
1907   {
1908     register const Quantum
1909       *magick_restrict p,
1910       *magick_restrict q;
1911
1912     register ssize_t
1913       x;
1914
1915     p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1916     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1917     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
1918       break;
1919     for (x=0; x < (ssize_t) columns; x++)
1920     {
1921       register ssize_t
1922         i;
1923
1924       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1925       {
1926         double
1927           distance;
1928
1929         PixelChannel channel = GetPixelChannelChannel(image,i);
1930         PixelTrait traits = GetPixelChannelTraits(image,channel);
1931         PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1932           channel);
1933         if ((traits == UndefinedPixelTrait) ||
1934             (reconstruct_traits == UndefinedPixelTrait) ||
1935             ((reconstruct_traits & UpdatePixelTrait) == 0))
1936           continue;
1937         distance=fabs(p[i]-(double) GetPixelChannel(reconstruct_image,
1938           channel,q));
1939         if (distance >= MagickEpsilon)
1940           break;
1941       }
1942       if (i < (ssize_t) GetPixelChannels(image))
1943         break;
1944       p+=GetPixelChannels(image);
1945       q+=GetPixelChannels(reconstruct_image);
1946     }
1947     if (x < (ssize_t) columns)
1948       break;
1949   }
1950   reconstruct_view=DestroyCacheView(reconstruct_view);
1951   image_view=DestroyCacheView(image_view);
1952   return(y < (ssize_t) rows ? MagickFalse : MagickTrue);
1953 }
1954 \f
1955 /*
1956 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1957 %                                                                             %
1958 %                                                                             %
1959 %                                                                             %
1960 %  S e t I m a g e C o l o r M e t r i c                                      %
1961 %                                                                             %
1962 %                                                                             %
1963 %                                                                             %
1964 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1965 %
1966 %  SetImageColorMetric() measures the difference between colors at each pixel
1967 %  location of two images.  A value other than 0 means the colors match
1968 %  exactly.  Otherwise an error measure is computed by summing over all
1969 %  pixels in an image the distance squared in RGB space between each image
1970 %  pixel and its corresponding pixel in the reconstruct image.  The error
1971 %  measure is assigned to these image members:
1972 %
1973 %    o mean_error_per_pixel:  The mean error for any single pixel in
1974 %      the image.
1975 %
1976 %    o normalized_mean_error:  The normalized mean quantization error for
1977 %      any single pixel in the image.  This distance measure is normalized to
1978 %      a range between 0 and 1.  It is independent of the range of red, green,
1979 %      and blue values in the image.
1980 %
1981 %    o normalized_maximum_error:  The normalized maximum quantization
1982 %      error for any single pixel in the image.  This distance measure is
1983 %      normalized to a range between 0 and 1.  It is independent of the range
1984 %      of red, green, and blue values in your image.
1985 %
1986 %  A small normalized mean square error, accessed as
1987 %  image->normalized_mean_error, suggests the images are very similar in
1988 %  spatial layout and color.
1989 %
1990 %  The format of the SetImageColorMetric method is:
1991 %
1992 %      MagickBooleanType SetImageColorMetric(Image *image,
1993 %        const Image *reconstruct_image,ExceptionInfo *exception)
1994 %
1995 %  A description of each parameter follows.
1996 %
1997 %    o image: the image.
1998 %
1999 %    o reconstruct_image: the reconstruct image.
2000 %
2001 %    o exception: return any errors or warnings in this structure.
2002 %
2003 */
2004 MagickExport MagickBooleanType SetImageColorMetric(Image *image,
2005   const Image *reconstruct_image,ExceptionInfo *exception)
2006 {
2007   CacheView
2008     *image_view,
2009     *reconstruct_view;
2010
2011   double
2012     area,
2013     maximum_error,
2014     mean_error,
2015     mean_error_per_pixel;
2016
2017   MagickBooleanType
2018     status;
2019
2020   size_t
2021     columns,
2022     rows;
2023
2024   ssize_t
2025     y;
2026
2027   assert(image != (Image *) NULL);
2028   assert(image->signature == MagickCoreSignature);
2029   assert(reconstruct_image != (const Image *) NULL);
2030   assert(reconstruct_image->signature == MagickCoreSignature);
2031   area=0.0;
2032   maximum_error=0.0;
2033   mean_error_per_pixel=0.0;
2034   mean_error=0.0;
2035   rows=MagickMax(image->rows,reconstruct_image->rows);
2036   columns=MagickMax(image->columns,reconstruct_image->columns);
2037   image_view=AcquireVirtualCacheView(image,exception);
2038   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
2039   for (y=0; y < (ssize_t) rows; y++)
2040   {
2041     register const Quantum
2042       *magick_restrict p,
2043       *magick_restrict q;
2044
2045     register ssize_t
2046       x;
2047
2048     p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
2049     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
2050     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2051       break;
2052     for (x=0; x < (ssize_t) columns; x++)
2053     {
2054       register ssize_t
2055         i;
2056
2057       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2058       {
2059         double
2060           distance;
2061
2062         PixelChannel channel = GetPixelChannelChannel(image,i);
2063         PixelTrait traits = GetPixelChannelTraits(image,channel);
2064         PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2065           channel);
2066         if ((traits == UndefinedPixelTrait) ||
2067             (reconstruct_traits == UndefinedPixelTrait) ||
2068             ((reconstruct_traits & UpdatePixelTrait) == 0))
2069           continue;
2070         distance=fabs(p[i]-(double) GetPixelChannel(reconstruct_image,
2071           channel,q));
2072         if (distance >= MagickEpsilon)
2073           {
2074             mean_error_per_pixel+=distance;
2075             mean_error+=distance*distance;
2076             if (distance > maximum_error)
2077               maximum_error=distance;
2078           }
2079         area++;
2080       }
2081       p+=GetPixelChannels(image);
2082       q+=GetPixelChannels(reconstruct_image);
2083     }
2084   }
2085   reconstruct_view=DestroyCacheView(reconstruct_view);
2086   image_view=DestroyCacheView(image_view);
2087   image->error.mean_error_per_pixel=(double) (mean_error_per_pixel/area);
2088   image->error.normalized_mean_error=(double) (QuantumScale*QuantumScale*
2089     mean_error/area);
2090   image->error.normalized_maximum_error=(double) (QuantumScale*maximum_error);
2091   status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
2092   return(status);
2093 }
2094 \f
2095 /*
2096 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2097 %                                                                             %
2098 %                                                                             %
2099 %                                                                             %
2100 %   S i m i l a r i t y I m a g e                                             %
2101 %                                                                             %
2102 %                                                                             %
2103 %                                                                             %
2104 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2105 %
2106 %  SimilarityImage() compares the reference image of the image and returns the
2107 %  best match offset.  In addition, it returns a similarity image such that an
2108 %  exact match location is completely white and if none of the pixels match,
2109 %  black, otherwise some gray level in-between.
2110 %
2111 %  The format of the SimilarityImageImage method is:
2112 %
2113 %      Image *SimilarityImage(const Image *image,const Image *reference,
2114 %        const MetricType metric,const double similarity_threshold,
2115 %        RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
2116 %
2117 %  A description of each parameter follows:
2118 %
2119 %    o image: the image.
2120 %
2121 %    o reference: find an area of the image that closely resembles this image.
2122 %
2123 %    o metric: the metric.
2124 %
2125 %    o similarity_threshold: minimum distortion for (sub)image match.
2126 %
2127 %    o offset: the best match offset of the reference image within the image.
2128 %
2129 %    o similarity: the computed similarity between the images.
2130 %
2131 %    o exception: return any errors or warnings in this structure.
2132 %
2133 */
2134
2135 static double GetSimilarityMetric(const Image *image,const Image *reference,
2136   const MetricType metric,const ssize_t x_offset,const ssize_t y_offset,
2137   ExceptionInfo *exception)
2138 {
2139   double
2140     distortion;
2141
2142   Image
2143     *similarity_image;
2144
2145   MagickBooleanType
2146     status;
2147
2148   RectangleInfo
2149     geometry;
2150
2151   SetGeometry(reference,&geometry);
2152   geometry.x=x_offset;
2153   geometry.y=y_offset;
2154   similarity_image=CropImage(image,&geometry,exception);
2155   if (similarity_image == (Image *) NULL)
2156     return(0.0);
2157   distortion=0.0;
2158   status=GetImageDistortion(similarity_image,reference,metric,&distortion,
2159     exception);
2160   similarity_image=DestroyImage(similarity_image);
2161   if (status == MagickFalse)
2162     return(0.0);
2163   return(distortion);
2164 }
2165
2166 MagickExport Image *SimilarityImage(const Image *image,const Image *reference,
2167   const MetricType metric,const double similarity_threshold,
2168   RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
2169 {
2170 #define SimilarityImageTag  "Similarity/Image"
2171
2172   CacheView
2173     *similarity_view;
2174
2175   Image
2176     *similarity_image;
2177
2178   MagickBooleanType
2179     status;
2180
2181   MagickOffsetType
2182     progress;
2183
2184   ssize_t
2185     y;
2186
2187   assert(image != (const Image *) NULL);
2188   assert(image->signature == MagickCoreSignature);
2189   if (image->debug != MagickFalse)
2190     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2191   assert(exception != (ExceptionInfo *) NULL);
2192   assert(exception->signature == MagickCoreSignature);
2193   assert(offset != (RectangleInfo *) NULL);
2194   SetGeometry(reference,offset);
2195   *similarity_metric=MagickMaximumValue;
2196   similarity_image=CloneImage(image,image->columns-reference->columns+1,
2197     image->rows-reference->rows+1,MagickTrue,exception);
2198   if (similarity_image == (Image *) NULL)
2199     return((Image *) NULL);
2200   status=SetImageStorageClass(similarity_image,DirectClass,exception);
2201   if (status == MagickFalse)
2202     {
2203       similarity_image=DestroyImage(similarity_image);
2204       return((Image *) NULL);
2205     }
2206   (void) SetImageAlphaChannel(similarity_image,DeactivateAlphaChannel,
2207     exception);
2208   /*
2209     Measure similarity of reference image against image.
2210   */
2211   status=MagickTrue;
2212   progress=0;
2213   similarity_view=AcquireAuthenticCacheView(similarity_image,exception);
2214 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2215   #pragma omp parallel for schedule(static) \
2216     shared(progress,status,similarity_metric) \
2217     magick_number_threads(image,image,image->rows-reference->rows+1,1)
2218 #endif
2219   for (y=0; y < (ssize_t) (image->rows-reference->rows+1); y++)
2220   {
2221     double
2222       similarity;
2223
2224     register Quantum
2225       *magick_restrict q;
2226
2227     register ssize_t
2228       x;
2229
2230     if (status == MagickFalse)
2231       continue;
2232 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2233     #pragma omp flush(similarity_metric)
2234 #endif
2235     if (*similarity_metric <= similarity_threshold)
2236       continue;
2237     q=GetCacheViewAuthenticPixels(similarity_view,0,y,similarity_image->columns,
2238       1,exception);
2239     if (q == (Quantum *) NULL)
2240       {
2241         status=MagickFalse;
2242         continue;
2243       }
2244     for (x=0; x < (ssize_t) (image->columns-reference->columns+1); x++)
2245     {
2246       register ssize_t
2247         i;
2248
2249 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2250       #pragma omp flush(similarity_metric)
2251 #endif
2252       if (*similarity_metric <= similarity_threshold)
2253         break;
2254       similarity=GetSimilarityMetric(image,reference,metric,x,y,exception);
2255 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2256       #pragma omp critical (MagickCore_SimilarityImage)
2257 #endif
2258       if ((metric == NormalizedCrossCorrelationErrorMetric) ||
2259           (metric == UndefinedErrorMetric))
2260         similarity=1.0-similarity;
2261       if (similarity < *similarity_metric)
2262         {
2263           offset->x=x;
2264           offset->y=y;
2265           *similarity_metric=similarity;
2266         }
2267       if (metric == PerceptualHashErrorMetric)
2268         similarity=MagickMin(0.01*similarity,1.0);
2269       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2270       {
2271         PixelChannel channel = GetPixelChannelChannel(image,i);
2272         PixelTrait traits = GetPixelChannelTraits(image,channel);
2273         PixelTrait similarity_traits=GetPixelChannelTraits(similarity_image,
2274           channel);
2275         if ((traits == UndefinedPixelTrait) ||
2276             (similarity_traits == UndefinedPixelTrait) ||
2277             ((similarity_traits & UpdatePixelTrait) == 0))
2278           continue;
2279         SetPixelChannel(similarity_image,channel,ClampToQuantum(QuantumRange-
2280           QuantumRange*similarity),q);
2281       }
2282       q+=GetPixelChannels(similarity_image);
2283     }
2284     if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
2285       status=MagickFalse;
2286     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2287       {
2288         MagickBooleanType
2289           proceed;
2290
2291 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2292         #pragma omp atomic
2293 #endif
2294         progress++;
2295         proceed=SetImageProgress(image,SimilarityImageTag,progress,image->rows);
2296         if (proceed == MagickFalse)
2297           status=MagickFalse;
2298       }
2299   }
2300   similarity_view=DestroyCacheView(similarity_view);
2301   if (status == MagickFalse)
2302     similarity_image=DestroyImage(similarity_image);
2303   return(similarity_image);
2304 }