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