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