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