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