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 %    http://www.imagemagick.org/script/license.php                            %
  27 %                                                                             %
  28 %  Unless required by applicable law or agreed to in writing, software        %
  29 %  distributed under the License is distributed on an "AS IS" BASIS,          %
  30 %  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
  31 %  See the License for the specific language governing permissions and        %
  32 %  limitations under the License.                                             %
  33 %                                                                             %
  34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  35 %
  36 %
  37 %
  38 */
  39 \f
  40 /*
  41   Include declarations.
  42 */
  43 #include "MagickCore/studio.h"
  44 #include "MagickCore/artifact.h"
  45 #include "MagickCore/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   double
 136     fuzz;
 137
 138   const char
 139     *artifact;
 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   (void) SetImageMask(highlight_image,ReadPixelMask,(Image *) NULL,exception);
 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) SetImageAlphaChannel(highlight_image,OpaqueAlphaChannel,exception);
 207   (void) QueryColorCompliance("#f1001ecc",AllCompliance,&highlight,exception);
 208   artifact=GetImageArtifact(image,"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,"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,"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)
 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)
 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)
 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]=20.0*MagickLog10((double) 1.0/sqrt(distortion[i]));
1217   return(status);
1218 }
1219
1220 static MagickBooleanType GetPerceptualHashDistortion(const Image *image,
1221   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1222 {
1223   ChannelPerceptualHash
1224     *channel_phash,
1225     *reconstruct_phash;
1226
1227   const char
1228     *artifact;
1229
1230   MagickBooleanType
1231     normalize;
1232
1233   ssize_t
1234     channel;
1235
1236   /*
1237     Compute perceptual hash in the sRGB colorspace.
1238   */
1239   channel_phash=GetImagePerceptualHash(image,exception);
1240   if (channel_phash == (ChannelPerceptualHash *) NULL)
1241     return(MagickFalse);
1242   reconstruct_phash=GetImagePerceptualHash(reconstruct_image,exception);
1243   if (reconstruct_phash == (ChannelPerceptualHash *) NULL)
1244     {
1245       channel_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1246         channel_phash);
1247       return(MagickFalse);
1248     }
1249   artifact=GetImageArtifact(image,"phash:normalize");
1250   normalize=(artifact == (const char *) NULL) ||
1251     (IsStringTrue(artifact) == MagickFalse) ? MagickFalse : MagickTrue;
1252 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1253   #pragma omp parallel for schedule(static,4)
1254 #endif
1255   for (channel=0; channel < MaxPixelChannels; channel++)
1256   {
1257     double
1258       difference;
1259
1260     register ssize_t
1261       i;
1262
1263     difference=0.0;
1264     for (i=0; i < MaximumNumberOfImageMoments; i++)
1265     {
1266       double
1267         alpha,
1268         beta;
1269
1270       register ssize_t
1271         j;
1272
1273       for (j=0; j < (ssize_t) channel_phash[0].number_colorspaces; j++)
1274       {
1275         alpha=channel_phash[channel].phash[j][i];
1276         beta=reconstruct_phash[channel].phash[j][i];
1277         if (normalize == MagickFalse)
1278           difference+=(beta-alpha)*(beta-alpha);
1279         else
1280           difference=sqrt((beta-alpha)*(beta-alpha)/
1281             channel_phash[0].number_channels);
1282       }
1283     }
1284     distortion[channel]+=difference;
1285 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1286     #pragma omp critical (MagickCore_GetPerceptualHashDistortion)
1287 #endif
1288     distortion[CompositePixelChannel]+=difference;
1289   }
1290   /*
1291     Free resources.
1292   */
1293   reconstruct_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1294     reconstruct_phash);
1295   channel_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(channel_phash);
1296   return(MagickTrue);
1297 }
1298
1299 static MagickBooleanType GetRootMeanSquaredDistortion(const Image *image,
1300   const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1301 {
1302   MagickBooleanType
1303     status;
1304
1305   register ssize_t
1306     i;
1307
1308   status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1309   for (i=0; i <= MaxPixelChannels; i++)
1310     distortion[i]=sqrt(distortion[i]);
1311   return(status);
1312 }
1313
1314 MagickExport MagickBooleanType GetImageDistortion(Image *image,
1315   const Image *reconstruct_image,const MetricType metric,double *distortion,
1316   ExceptionInfo *exception)
1317 {
1318   double
1319     *channel_distortion;
1320
1321   MagickBooleanType
1322     status;
1323
1324   size_t
1325     length;
1326
1327   assert(image != (Image *) NULL);
1328   assert(image->signature == MagickCoreSignature);
1329   if (image->debug != MagickFalse)
1330     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1331   assert(reconstruct_image != (const Image *) NULL);
1332   assert(reconstruct_image->signature == MagickCoreSignature);
1333   assert(distortion != (double *) NULL);
1334   *distortion=0.0;
1335   if (image->debug != MagickFalse)
1336     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1337   /*
1338     Get image distortion.
1339   */
1340   length=MaxPixelChannels+1;
1341   channel_distortion=(double *) AcquireQuantumMemory(length,
1342     sizeof(*channel_distortion));
1343   if (channel_distortion == (double *) NULL)
1344     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1345   (void) ResetMagickMemory(channel_distortion,0,length*
1346     sizeof(*channel_distortion));
1347   switch (metric)
1348   {
1349     case AbsoluteErrorMetric:
1350     {
1351       status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1352         exception);
1353       break;
1354     }
1355     case FuzzErrorMetric:
1356     {
1357       status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1358         exception);
1359       break;
1360     }
1361     case MeanAbsoluteErrorMetric:
1362     {
1363       status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1364         channel_distortion,exception);
1365       break;
1366     }
1367     case MeanErrorPerPixelErrorMetric:
1368     {
1369       status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1370         exception);
1371       break;
1372     }
1373     case MeanSquaredErrorMetric:
1374     {
1375       status=GetMeanSquaredDistortion(image,reconstruct_image,
1376         channel_distortion,exception);
1377       break;
1378     }
1379     case NormalizedCrossCorrelationErrorMetric:
1380     default:
1381     {
1382       status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1383         channel_distortion,exception);
1384       break;
1385     }
1386     case PeakAbsoluteErrorMetric:
1387     {
1388       status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1389         channel_distortion,exception);
1390       break;
1391     }
1392     case PeakSignalToNoiseRatioErrorMetric:
1393     {
1394       status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1395         channel_distortion,exception);
1396       break;
1397     }
1398     case PerceptualHashErrorMetric:
1399     {
1400       status=GetPerceptualHashDistortion(image,reconstruct_image,
1401         channel_distortion,exception);
1402       break;
1403     }
1404     case RootMeanSquaredErrorMetric:
1405     {
1406       status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1407         channel_distortion,exception);
1408       break;
1409     }
1410   }
1411   *distortion=channel_distortion[CompositePixelChannel];
1412   channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1413   (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
1414     *distortion);
1415   return(status);
1416 }
1417 \f
1418 /*
1419 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1420 %                                                                             %
1421 %                                                                             %
1422 %                                                                             %
1423 %   G e t I m a g e D i s t o r t i o n s                                     %
1424 %                                                                             %
1425 %                                                                             %
1426 %                                                                             %
1427 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1428 %
1429 %  GetImageDistortions() compares the pixel channels of an image to a
1430 %  reconstructed image and returns the specified distortion metric for each
1431 %  channel.
1432 %
1433 %  The format of the GetImageDistortions method is:
1434 %
1435 %      double *GetImageDistortions(const Image *image,
1436 %        const Image *reconstruct_image,const MetricType metric,
1437 %        ExceptionInfo *exception)
1438 %
1439 %  A description of each parameter follows:
1440 %
1441 %    o image: the image.
1442 %
1443 %    o reconstruct_image: the reconstruct image.
1444 %
1445 %    o metric: the metric.
1446 %
1447 %    o exception: return any errors or warnings in this structure.
1448 %
1449 */
1450 MagickExport double *GetImageDistortions(Image *image,
1451   const Image *reconstruct_image,const MetricType metric,
1452   ExceptionInfo *exception)
1453 {
1454   double
1455     *channel_distortion;
1456
1457   MagickBooleanType
1458     status;
1459
1460   size_t
1461     length;
1462
1463   assert(image != (Image *) NULL);
1464   assert(image->signature == MagickCoreSignature);
1465   if (image->debug != MagickFalse)
1466     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1467   assert(reconstruct_image != (const Image *) NULL);
1468   assert(reconstruct_image->signature == MagickCoreSignature);
1469   if (image->debug != MagickFalse)
1470     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1471   /*
1472     Get image distortion.
1473   */
1474   length=MaxPixelChannels+1UL;
1475   channel_distortion=(double *) AcquireQuantumMemory(length,
1476     sizeof(*channel_distortion));
1477   if (channel_distortion == (double *) NULL)
1478     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1479   (void) ResetMagickMemory(channel_distortion,0,length*
1480     sizeof(*channel_distortion));
1481   status=MagickTrue;
1482   switch (metric)
1483   {
1484     case AbsoluteErrorMetric:
1485     {
1486       status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1487         exception);
1488       break;
1489     }
1490     case FuzzErrorMetric:
1491     {
1492       status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1493         exception);
1494       break;
1495     }
1496     case MeanAbsoluteErrorMetric:
1497     {
1498       status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1499         channel_distortion,exception);
1500       break;
1501     }
1502     case MeanErrorPerPixelErrorMetric:
1503     {
1504       status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1505         exception);
1506       break;
1507     }
1508     case MeanSquaredErrorMetric:
1509     {
1510       status=GetMeanSquaredDistortion(image,reconstruct_image,
1511         channel_distortion,exception);
1512       break;
1513     }
1514     case NormalizedCrossCorrelationErrorMetric:
1515     default:
1516     {
1517       status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1518         channel_distortion,exception);
1519       break;
1520     }
1521     case PeakAbsoluteErrorMetric:
1522     {
1523       status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1524         channel_distortion,exception);
1525       break;
1526     }
1527     case PeakSignalToNoiseRatioErrorMetric:
1528     {
1529       status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1530         channel_distortion,exception);
1531       break;
1532     }
1533     case PerceptualHashErrorMetric:
1534     {
1535       status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1536         channel_distortion,exception);
1537       break;
1538     }
1539     case RootMeanSquaredErrorMetric:
1540     {
1541       status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1542         channel_distortion,exception);
1543       break;
1544     }
1545   }
1546   if (status == MagickFalse)
1547     {
1548       channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1549       return((double *) NULL);
1550     }
1551   return(channel_distortion);
1552 }
1553 \f
1554 /*
1555 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1556 %                                                                             %
1557 %                                                                             %
1558 %                                                                             %
1559 %  I s I m a g e s E q u a l                                                  %
1560 %                                                                             %
1561 %                                                                             %
1562 %                                                                             %
1563 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1564 %
1565 %  IsImagesEqual() compare the pixels of two images and returns immediately
1566 %  if any pixel is not identical.
1567 %
1568 %  The format of the IsImagesEqual method is:
1569 %
1570 %      MagickBooleanType IsImagesEqual(const Image *image,
1571 %        const Image *reconstruct_image,ExceptionInfo *exception)
1572 %
1573 %  A description of each parameter follows.
1574 %
1575 %    o image: the image.
1576 %
1577 %    o reconstruct_image: the reconstruct image.
1578 %
1579 %    o exception: return any errors or warnings in this structure.
1580 %
1581 */
1582 MagickExport MagickBooleanType IsImagesEqual(const Image *image,
1583   const Image *reconstruct_image,ExceptionInfo *exception)
1584 {
1585   CacheView
1586     *image_view,
1587     *reconstruct_view;
1588
1589   size_t
1590     columns,
1591     rows;
1592
1593   ssize_t
1594     y;
1595
1596   assert(image != (Image *) NULL);
1597   assert(image->signature == MagickCoreSignature);
1598   assert(reconstruct_image != (const Image *) NULL);
1599   assert(reconstruct_image->signature == MagickCoreSignature);
1600   rows=MagickMax(image->rows,reconstruct_image->rows);
1601   columns=MagickMax(image->columns,reconstruct_image->columns);
1602   image_view=AcquireVirtualCacheView(image,exception);
1603   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1604   for (y=0; y < (ssize_t) rows; y++)
1605   {
1606     register const Quantum
1607       *magick_restrict p,
1608       *magick_restrict q;
1609
1610     register ssize_t
1611       x;
1612
1613     p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1614     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1615     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
1616       break;
1617     for (x=0; x < (ssize_t) columns; x++)
1618     {
1619       register ssize_t
1620         i;
1621
1622       if (GetPixelWriteMask(image,p) == 0)
1623         {
1624           p+=GetPixelChannels(image);
1625           q+=GetPixelChannels(reconstruct_image);
1626           continue;
1627         }
1628       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1629       {
1630         double
1631           distance;
1632
1633         PixelChannel channel=GetPixelChannelChannel(image,i);
1634         PixelTrait traits=GetPixelChannelTraits(image,channel);
1635         PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
1636           channel);
1637         if ((traits == UndefinedPixelTrait) ||
1638             (reconstruct_traits == UndefinedPixelTrait) ||
1639             ((reconstruct_traits & UpdatePixelTrait) == 0))
1640           continue;
1641         distance=fabs(p[i]-(double) GetPixelChannel(reconstruct_image,
1642           channel,q));
1643         if (distance >= MagickEpsilon)
1644           break;
1645       }
1646       if (i < (ssize_t) GetPixelChannels(image))
1647         break;
1648       p+=GetPixelChannels(image);
1649       q+=GetPixelChannels(reconstruct_image);
1650     }
1651     if (x < (ssize_t) columns)
1652       break;
1653   }
1654   reconstruct_view=DestroyCacheView(reconstruct_view);
1655   image_view=DestroyCacheView(image_view);
1656   return(y < (ssize_t) rows ? MagickFalse : MagickTrue);
1657 }
1658 \f
1659 /*
1660 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1661 %                                                                             %
1662 %                                                                             %
1663 %                                                                             %
1664 %  S e t I m a g e C o l o r M e t r i c                                      %
1665 %                                                                             %
1666 %                                                                             %
1667 %                                                                             %
1668 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1669 %
1670 %  SetImageColorMetric() measures the difference between colors at each pixel
1671 %  location of two images.  A value other than 0 means the colors match
1672 %  exactly.  Otherwise an error measure is computed by summing over all
1673 %  pixels in an image the distance squared in RGB space between each image
1674 %  pixel and its corresponding pixel in the reconstruct image.  The error
1675 %  measure is assigned to these image members:
1676 %
1677 %    o mean_error_per_pixel:  The mean error for any single pixel in
1678 %      the image.
1679 %
1680 %    o normalized_mean_error:  The normalized mean quantization error for
1681 %      any single pixel in the image.  This distance measure is normalized to
1682 %      a range between 0 and 1.  It is independent of the range of red, green,
1683 %      and blue values in the image.
1684 %
1685 %    o normalized_maximum_error:  The normalized maximum quantization
1686 %      error for any single pixel in the image.  This distance measure is
1687 %      normalized to a range between 0 and 1.  It is independent of the range
1688 %      of red, green, and blue values in your image.
1689 %
1690 %  A small normalized mean square error, accessed as
1691 %  image->normalized_mean_error, suggests the images are very similar in
1692 %  spatial layout and color.
1693 %
1694 %  The format of the SetImageColorMetric method is:
1695 %
1696 %      MagickBooleanType SetImageColorMetric(Image *image,
1697 %        const Image *reconstruct_image,ExceptionInfo *exception)
1698 %
1699 %  A description of each parameter follows.
1700 %
1701 %    o image: the image.
1702 %
1703 %    o reconstruct_image: the reconstruct image.
1704 %
1705 %    o exception: return any errors or warnings in this structure.
1706 %
1707 */
1708 MagickExport MagickBooleanType SetImageColorMetric(Image *image,
1709   const Image *reconstruct_image,ExceptionInfo *exception)
1710 {
1711   CacheView
1712     *image_view,
1713     *reconstruct_view;
1714
1715   double
1716     area,
1717     maximum_error,
1718     mean_error,
1719     mean_error_per_pixel;
1720
1721   MagickBooleanType
1722     status;
1723
1724   size_t
1725     columns,
1726     rows;
1727
1728   ssize_t
1729     y;
1730
1731   assert(image != (Image *) NULL);
1732   assert(image->signature == MagickCoreSignature);
1733   assert(reconstruct_image != (const Image *) NULL);
1734   assert(reconstruct_image->signature == MagickCoreSignature);
1735   area=0.0;
1736   maximum_error=0.0;
1737   mean_error_per_pixel=0.0;
1738   mean_error=0.0;
1739   rows=MagickMax(image->rows,reconstruct_image->rows);
1740   columns=MagickMax(image->columns,reconstruct_image->columns);
1741   image_view=AcquireVirtualCacheView(image,exception);
1742   reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1743   for (y=0; y < (ssize_t) rows; y++)
1744   {
1745     register const Quantum
1746       *magick_restrict p,
1747       *magick_restrict q;
1748
1749     register ssize_t
1750       x;
1751
1752     p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1753     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1754     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
1755       break;
1756     for (x=0; x < (ssize_t) columns; x++)
1757     {
1758       register ssize_t
1759         i;
1760
1761       if (GetPixelWriteMask(image,p) == 0)
1762         {
1763           p+=GetPixelChannels(image);
1764           q+=GetPixelChannels(reconstruct_image);
1765           continue;
1766         }
1767       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1768       {
1769         double
1770           distance;
1771
1772         PixelChannel channel=GetPixelChannelChannel(image,i);
1773         PixelTrait traits=GetPixelChannelTraits(image,channel);
1774         PixelTrait reconstruct_traits=GetPixelChannelTraits(reconstruct_image,
1775           channel);
1776         if ((traits == UndefinedPixelTrait) ||
1777             (reconstruct_traits == UndefinedPixelTrait) ||
1778             ((reconstruct_traits & UpdatePixelTrait) == 0))
1779           continue;
1780         distance=fabs(p[i]-(double) GetPixelChannel(reconstruct_image,
1781           channel,q));
1782         if (distance >= MagickEpsilon)
1783           {
1784             mean_error_per_pixel+=distance;
1785             mean_error+=distance*distance;
1786             if (distance > maximum_error)
1787               maximum_error=distance;
1788           }
1789         area++;
1790       }
1791       p+=GetPixelChannels(image);
1792       q+=GetPixelChannels(reconstruct_image);
1793     }
1794   }
1795   reconstruct_view=DestroyCacheView(reconstruct_view);
1796   image_view=DestroyCacheView(image_view);
1797   image->error.mean_error_per_pixel=(double) (mean_error_per_pixel/area);
1798   image->error.normalized_mean_error=(double) (QuantumScale*QuantumScale*
1799     mean_error/area);
1800   image->error.normalized_maximum_error=(double) (QuantumScale*maximum_error);
1801   status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
1802   return(status);
1803 }
1804 \f
1805 /*
1806 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1807 %                                                                             %
1808 %                                                                             %
1809 %                                                                             %
1810 %   S i m i l a r i t y I m a g e                                             %
1811 %                                                                             %
1812 %                                                                             %
1813 %                                                                             %
1814 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1815 %
1816 %  SimilarityImage() compares the reference image of the image and returns the
1817 %  best match offset.  In addition, it returns a similarity image such that an
1818 %  exact match location is completely white and if none of the pixels match,
1819 %  black, otherwise some gray level in-between.
1820 %
1821 %  The format of the SimilarityImageImage method is:
1822 %
1823 %      Image *SimilarityImage(const Image *image,const Image *reference,
1824 %        const MetricType metric,const double similarity_threshold,
1825 %        RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
1826 %
1827 %  A description of each parameter follows:
1828 %
1829 %    o image: the image.
1830 %
1831 %    o reference: find an area of the image that closely resembles this image.
1832 %
1833 %    o metric: the metric.
1834 %
1835 %    o similarity_threshold: minimum distortion for (sub)image match.
1836 %
1837 %    o offset: the best match offset of the reference image within the image.
1838 %
1839 %    o similarity: the computed similarity between the images.
1840 %
1841 %    o exception: return any errors or warnings in this structure.
1842 %
1843 */
1844
1845 static double GetSimilarityMetric(const Image *image,const Image *reference,
1846   const MetricType metric,const ssize_t x_offset,const ssize_t y_offset,
1847   ExceptionInfo *exception)
1848 {
1849   double
1850     distortion;
1851
1852   Image
1853     *similarity_image;
1854
1855   MagickBooleanType
1856     status;
1857
1858   RectangleInfo
1859     geometry;
1860
1861   SetGeometry(reference,&geometry);
1862   geometry.x=x_offset;
1863   geometry.y=y_offset;
1864   similarity_image=CropImage(image,&geometry,exception);
1865   if (similarity_image == (Image *) NULL)
1866     return(0.0);
1867   distortion=0.0;
1868   status=GetImageDistortion(similarity_image,reference,metric,&distortion,
1869     exception);
1870   similarity_image=DestroyImage(similarity_image);
1871   if (status == MagickFalse)
1872     return(0.0);
1873   return(distortion);
1874 }
1875
1876 MagickExport Image *SimilarityImage(const Image *image,const Image *reference,
1877   const MetricType metric,const double similarity_threshold,
1878   RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
1879 {
1880 #define SimilarityImageTag  "Similarity/Image"
1881
1882   CacheView
1883     *similarity_view;
1884
1885   Image
1886     *similarity_image;
1887
1888   MagickBooleanType
1889     status;
1890
1891   MagickOffsetType
1892     progress;
1893
1894   ssize_t
1895     y;
1896
1897   assert(image != (const Image *) NULL);
1898   assert(image->signature == MagickCoreSignature);
1899   if (image->debug != MagickFalse)
1900     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1901   assert(exception != (ExceptionInfo *) NULL);
1902   assert(exception->signature == MagickCoreSignature);
1903   assert(offset != (RectangleInfo *) NULL);
1904   SetGeometry(reference,offset);
1905   *similarity_metric=MagickMaximumValue;
1906   similarity_image=CloneImage(image,image->columns-reference->columns+1,
1907     image->rows-reference->rows+1,MagickTrue,exception);
1908   if (similarity_image == (Image *) NULL)
1909     return((Image *) NULL);
1910   status=SetImageStorageClass(similarity_image,DirectClass,exception);
1911   if (status == MagickFalse)
1912     {
1913       similarity_image=DestroyImage(similarity_image);
1914       return((Image *) NULL);
1915     }
1916   (void) SetImageAlphaChannel(similarity_image,DeactivateAlphaChannel,
1917     exception);
1918   /*
1919     Measure similarity of reference image against image.
1920   */
1921   status=MagickTrue;
1922   progress=0;
1923   similarity_view=AcquireAuthenticCacheView(similarity_image,exception);
1924 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1925   #pragma omp parallel for schedule(static,4) \
1926     shared(progress,status,similarity_metric) \
1927     magick_threads(image,image,image->rows,1)
1928 #endif
1929   for (y=0; y < (ssize_t) (image->rows-reference->rows+1); y++)
1930   {
1931     double
1932       similarity;
1933
1934     register Quantum
1935       *magick_restrict q;
1936
1937     register ssize_t
1938       x;
1939
1940     if (status == MagickFalse)
1941       continue;
1942 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1943     #pragma omp flush(similarity_metric)
1944 #endif
1945     if (*similarity_metric <= similarity_threshold)
1946       continue;
1947     q=GetCacheViewAuthenticPixels(similarity_view,0,y,similarity_image->columns,
1948       1,exception);
1949     if (q == (Quantum *) NULL)
1950       {
1951         status=MagickFalse;
1952         continue;
1953       }
1954     for (x=0; x < (ssize_t) (image->columns-reference->columns+1); x++)
1955     {
1956       register ssize_t
1957         i;
1958
1959 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1960       #pragma omp flush(similarity_metric)
1961 #endif
1962       if (*similarity_metric <= similarity_threshold)
1963         break;
1964       similarity=GetSimilarityMetric(image,reference,metric,x,y,exception);
1965 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1966       #pragma omp critical (MagickCore_SimilarityImage)
1967 #endif
1968       if ((metric == NormalizedCrossCorrelationErrorMetric) ||
1969           (metric == UndefinedErrorMetric))
1970         similarity=1.0-similarity;
1971       if (similarity < *similarity_metric)
1972         {
1973           offset->x=x;
1974           offset->y=y;
1975           *similarity_metric=similarity;
1976         }
1977       if (metric == PerceptualHashErrorMetric)
1978         similarity=MagickMin(0.01*similarity,1.0);
1979       if (GetPixelWriteMask(similarity_image,q) == 0)
1980         {
1981           SetPixelBackgoundColor(similarity_image,q);
1982           q+=GetPixelChannels(similarity_image);
1983           continue;
1984         }
1985       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1986       {
1987         PixelChannel channel=GetPixelChannelChannel(image,i);
1988         PixelTrait traits=GetPixelChannelTraits(image,channel);
1989         PixelTrait similarity_traits=GetPixelChannelTraits(similarity_image,
1990           channel);
1991         if ((traits == UndefinedPixelTrait) ||
1992             (similarity_traits == UndefinedPixelTrait) ||
1993             ((similarity_traits & UpdatePixelTrait) == 0))
1994           continue;
1995         SetPixelChannel(similarity_image,channel,ClampToQuantum(QuantumRange-
1996           QuantumRange*similarity),q);
1997       }
1998       q+=GetPixelChannels(similarity_image);
1999     }
2000     if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
2001       status=MagickFalse;
2002     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2003       {
2004         MagickBooleanType
2005           proceed;
2006
2007 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2008         #pragma omp critical (MagickCore_SimilarityImage)
2009 #endif
2010         proceed=SetImageProgress(image,SimilarityImageTag,progress++,
2011           image->rows);
2012         if (proceed == MagickFalse)
2013           status=MagickFalse;
2014       }
2015   }
2016   similarity_view=DestroyCacheView(similarity_view);
2017   if (status == MagickFalse)
2018     similarity_image=DestroyImage(similarity_image);
2019   return(similarity_image);
2020 }