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