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