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) && (GetRedPixelComponent(p) != q->red))
 266             difference=MagickTrue;
 267           if (((channel & GreenChannel) != 0) && (GetGreenPixelComponent(p) != q->green))
 268             difference=MagickTrue;
 269           if (((channel & BlueChannel) != 0) && (GetBluePixelComponent(p) != q->blue))
 270             difference=MagickTrue;
 271           if (((channel & OpacityChannel) != 0) &&
 272               (image->matte != MagickFalse) && (GetOpacityPixelComponent(p) != 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   MagickBooleanType
 359     status;
 360
 361   MagickPixelPacket
 362     zero;
 363
 364   ssize_t
 365     y;
 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   MagickBooleanType
 479     status;
 480
 481   register ssize_t
 482     i;
 483
 484   ssize_t
 485     y;
 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   MagickBooleanType
 595     status;
 596
 597   register ssize_t
 598     i;
 599
 600   ssize_t
 601     y;
 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   MagickBooleanType
 703     status;
 704
 705   MagickRealType
 706     alpha,
 707     area,
 708     beta,
 709     maximum_error,
 710     mean_error;
 711
 712   ssize_t
 713     y;
 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   MagickBooleanType
 832     status;
 833
 834   register ssize_t
 835     i;
 836
 837   ssize_t
 838     y;
 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   MagickBooleanType
1086     status;
1087
1088   ssize_t
1089     y;
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   status=MagickTrue;
1417   switch (metric)
1418   {
1419     case AbsoluteErrorMetric:
1420     {
1421       status=GetAbsoluteDistortion(image,reconstruct_image,AllChannels,
1422         channel_distortion,exception);
1423       break;
1424     }
1425     case FuzzErrorMetric:
1426     {
1427       status=GetFuzzDistortion(image,reconstruct_image,AllChannels,
1428         channel_distortion,exception);
1429       break;
1430     }
1431     case MeanAbsoluteErrorMetric:
1432     {
1433       status=GetMeanAbsoluteDistortion(image,reconstruct_image,AllChannels,
1434         channel_distortion,exception);
1435       break;
1436     }
1437     case MeanErrorPerPixelMetric:
1438     {
1439       status=GetMeanErrorPerPixel(image,reconstruct_image,AllChannels,
1440         channel_distortion,exception);
1441       break;
1442     }
1443     case MeanSquaredErrorMetric:
1444     {
1445       status=GetMeanSquaredDistortion(image,reconstruct_image,AllChannels,
1446         channel_distortion,exception);
1447       break;
1448     }
1449     case NormalizedCrossCorrelationErrorMetric:
1450     default:
1451     {
1452       status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1453         AllChannels,channel_distortion,exception);
1454       break;
1455     }
1456     case PeakAbsoluteErrorMetric:
1457     {
1458       status=GetPeakAbsoluteDistortion(image,reconstruct_image,AllChannels,
1459         channel_distortion,exception);
1460       break;
1461     }
1462     case PeakSignalToNoiseRatioMetric:
1463     {
1464       status=GetPeakSignalToNoiseRatio(image,reconstruct_image,AllChannels,
1465         channel_distortion,exception);
1466       break;
1467     }
1468     case RootMeanSquaredErrorMetric:
1469     {
1470       status=GetRootMeanSquaredDistortion(image,reconstruct_image,AllChannels,
1471         channel_distortion,exception);
1472       break;
1473     }
1474   }
1475   if (status == MagickFalse)
1476     {
1477       channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1478       return((double *) NULL);
1479     }
1480   return(channel_distortion);
1481 }
1482 \f
1483 /*
1484 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1485 %                                                                             %
1486 %                                                                             %
1487 %                                                                             %
1488 %  I s I m a g e s E q u a l                                                  %
1489 %                                                                             %
1490 %                                                                             %
1491 %                                                                             %
1492 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1493 %
1494 %  IsImagesEqual() measures the difference between colors at each pixel
1495 %  location of two images.  A value other than 0 means the colors match
1496 %  exactly.  Otherwise an error measure is computed by summing over all
1497 %  pixels in an image the distance squared in RGB space between each image
1498 %  pixel and its corresponding pixel in the reconstruct image.  The error
1499 %  measure is assigned to these image members:
1500 %
1501 %    o mean_error_per_pixel:  The mean error for any single pixel in
1502 %      the image.
1503 %
1504 %    o normalized_mean_error:  The normalized mean quantization error for
1505 %      any single pixel in the image.  This distance measure is normalized to
1506 %      a range between 0 and 1.  It is independent of the range of red, green,
1507 %      and blue values in the image.
1508 %
1509 %    o normalized_maximum_error:  The normalized maximum quantization
1510 %      error for any single pixel in the image.  This distance measure is
1511 %      normalized to a range between 0 and 1.  It is independent of the range
1512 %      of red, green, and blue values in your image.
1513 %
1514 %  A small normalized mean square error, accessed as
1515 %  image->normalized_mean_error, suggests the images are very similar in
1516 %  spatial layout and color.
1517 %
1518 %  The format of the IsImagesEqual method is:
1519 %
1520 %      MagickBooleanType IsImagesEqual(Image *image,
1521 %        const Image *reconstruct_image)
1522 %
1523 %  A description of each parameter follows.
1524 %
1525 %    o image: the image.
1526 %
1527 %    o reconstruct_image: the reconstruct image.
1528 %
1529 */
1530 MagickExport MagickBooleanType IsImagesEqual(Image *image,
1531   const Image *reconstruct_image)
1532 {
1533   CacheView
1534     *image_view,
1535     *reconstruct_view;
1536
1537   ExceptionInfo
1538     *exception;
1539
1540   MagickBooleanType
1541     status;
1542
1543   MagickRealType
1544     area,
1545     maximum_error,
1546     mean_error,
1547     mean_error_per_pixel;
1548
1549   ssize_t
1550     y;
1551
1552   assert(image != (Image *) NULL);
1553   assert(image->signature == MagickSignature);
1554   assert(reconstruct_image != (const Image *) NULL);
1555   assert(reconstruct_image->signature == MagickSignature);
1556   if ((reconstruct_image->columns != image->columns) ||
1557       (reconstruct_image->rows != image->rows))
1558     ThrowBinaryException(ImageError,"ImageSizeDiffers",image->filename);
1559   area=0.0;
1560   maximum_error=0.0;
1561   mean_error_per_pixel=0.0;
1562   mean_error=0.0;
1563   exception=(&image->exception);
1564   image_view=AcquireCacheView(image);
1565   reconstruct_view=AcquireCacheView(reconstruct_image);
1566   for (y=0; y < (ssize_t) image->rows; y++)
1567   {
1568     register const IndexPacket
1569       *restrict indexes,
1570       *restrict reconstruct_indexes;
1571
1572     register const PixelPacket
1573       *restrict p,
1574       *restrict q;
1575
1576     register ssize_t
1577       x;
1578
1579     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1580     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
1581       1,exception);
1582     if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1583       break;
1584     indexes=GetCacheViewVirtualIndexQueue(image_view);
1585     reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
1586     for (x=0; x < (ssize_t) image->columns; x++)
1587     {
1588       MagickRealType
1589         distance;
1590
1591       distance=fabs(p->red-(double) q->red);
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->green-(double) q->green);
1598       mean_error_per_pixel+=distance;
1599       mean_error+=distance*distance;
1600       if (distance > maximum_error)
1601         maximum_error=distance;
1602       area++;
1603       distance=fabs(p->blue-(double) q->blue);
1604       mean_error_per_pixel+=distance;
1605       mean_error+=distance*distance;
1606       if (distance > maximum_error)
1607         maximum_error=distance;
1608       area++;
1609       if (image->matte != MagickFalse)
1610         {
1611           distance=fabs(p->opacity-(double) q->opacity);
1612           mean_error_per_pixel+=distance;
1613           mean_error+=distance*distance;
1614           if (distance > maximum_error)
1615             maximum_error=distance;
1616           area++;
1617         }
1618       if ((image->colorspace == CMYKColorspace) &&
1619           (reconstruct_image->colorspace == CMYKColorspace))
1620         {
1621           distance=fabs(indexes[x]-(double) reconstruct_indexes[x]);
1622           mean_error_per_pixel+=distance;
1623           mean_error+=distance*distance;
1624           if (distance > maximum_error)
1625             maximum_error=distance;
1626           area++;
1627         }
1628       p++;
1629       q++;
1630     }
1631   }
1632   reconstruct_view=DestroyCacheView(reconstruct_view);
1633   image_view=DestroyCacheView(image_view);
1634   image->error.mean_error_per_pixel=(double) (mean_error_per_pixel/area);
1635   image->error.normalized_mean_error=(double) (QuantumScale*QuantumScale*
1636     mean_error/area);
1637   image->error.normalized_maximum_error=(double) (QuantumScale*maximum_error);
1638   status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
1639   return(status);
1640 }
1641 \f
1642 /*
1643 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1644 %                                                                             %
1645 %                                                                             %
1646 %                                                                             %
1647 %   S i m i l a r i t y I m a g e                                             %
1648 %                                                                             %
1649 %                                                                             %
1650 %                                                                             %
1651 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1652 %
1653 %  SimilarityImage() compares the reference image of the image and returns the
1654 %  best match offset.  In addition, it returns a similarity image such that an
1655 %  exact match location is completely white and if none of the pixels match,
1656 %  black, otherwise some gray level in-between.
1657 %
1658 %  The format of the SimilarityImageImage method is:
1659 %
1660 %      Image *SimilarityImage(const Image *image,const Image *reference,
1661 %        RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
1662 %
1663 %  A description of each parameter follows:
1664 %
1665 %    o image: the image.
1666 %
1667 %    o reference: find an area of the image that closely resembles this image.
1668 %
1669 %    o the best match offset of the reference image within the image.
1670 %
1671 %    o similarity: the computed similarity between the images.
1672 %
1673 %    o exception: return any errors or warnings in this structure.
1674 %
1675 */
1676
1677 static double GetNCCDistortion(const Image *image,
1678   const Image *reconstruct_image,
1679   const ChannelStatistics *reconstruct_statistics,ExceptionInfo *exception)
1680 {
1681 #define SimilarityImageTag  "Similarity/Image"
1682
1683   CacheView
1684     *image_view,
1685     *reconstruct_view;
1686
1687   ChannelStatistics
1688     *image_statistics;
1689
1690   double
1691     distortion;
1692
1693   MagickBooleanType
1694     status;
1695
1696   MagickRealType
1697     area,
1698     gamma;
1699
1700   ssize_t
1701     y;
1702
1703   unsigned long
1704     number_channels;
1705
1706   /*
1707     Normalize to account for variation due to lighting and exposure condition.
1708   */
1709   image_statistics=GetImageChannelStatistics(image,exception);
1710   status=MagickTrue;
1711   distortion=0.0;
1712   area=1.0/((MagickRealType) image->columns*image->rows);
1713   image_view=AcquireCacheView(image);
1714   reconstruct_view=AcquireCacheView(reconstruct_image);
1715   for (y=0; y < (ssize_t) image->rows; y++)
1716   {
1717     register const IndexPacket
1718       *restrict indexes,
1719       *restrict reconstruct_indexes;
1720
1721     register const PixelPacket
1722       *restrict p,
1723       *restrict q;
1724
1725     register ssize_t
1726       x;
1727
1728     if (status == MagickFalse)
1729       continue;
1730     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1731     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
1732       1,exception);
1733     if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1734       {
1735         status=MagickFalse;
1736         continue;
1737       }
1738     indexes=GetCacheViewVirtualIndexQueue(image_view);
1739     reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
1740     for (x=0; x < (ssize_t) image->columns; x++)
1741     {
1742       distortion+=area*QuantumScale*(p->red-
1743         image_statistics[RedChannel].mean)*(q->red-
1744         reconstruct_statistics[RedChannel].mean);
1745       distortion+=area*QuantumScale*(p->green-
1746         image_statistics[GreenChannel].mean)*(q->green-
1747         reconstruct_statistics[GreenChannel].mean);
1748       distortion+=area*QuantumScale*(p->blue-
1749         image_statistics[BlueChannel].mean)*(q->blue-
1750         reconstruct_statistics[BlueChannel].mean);
1751       if (image->matte != MagickFalse)
1752         distortion+=area*QuantumScale*(p->opacity-
1753           image_statistics[OpacityChannel].mean)*(q->opacity-
1754           reconstruct_statistics[OpacityChannel].mean);
1755       if ((image->colorspace == CMYKColorspace) &&
1756           (reconstruct_image->colorspace == CMYKColorspace))
1757         distortion+=area*QuantumScale*(indexes[x]-
1758           image_statistics[OpacityChannel].mean)*(reconstruct_indexes[x]-
1759           reconstruct_statistics[OpacityChannel].mean);
1760       p++;
1761       q++;
1762     }
1763   }
1764   reconstruct_view=DestroyCacheView(reconstruct_view);
1765   image_view=DestroyCacheView(image_view);
1766   /*
1767     Divide by the standard deviation.
1768   */
1769   gamma=image_statistics[AllChannels].standard_deviation*
1770     reconstruct_statistics[AllChannels].standard_deviation;
1771   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1772   distortion=QuantumRange*gamma*distortion;
1773   number_channels=3;
1774   if (image->matte != MagickFalse)
1775     number_channels++;
1776   if (image->colorspace == CMYKColorspace)
1777     number_channels++;
1778   distortion=sqrt(distortion/number_channels);
1779   /*
1780     Free resources.
1781   */
1782   image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1783     image_statistics);
1784   return(1.0-distortion);
1785 }
1786
1787 static double GetSimilarityMetric(const Image *image,const Image *reference,
1788   const ChannelStatistics *reference_statistics,const ssize_t x_offset,
1789   const ssize_t y_offset,ExceptionInfo *exception)
1790 {
1791   double
1792     distortion;
1793
1794   Image
1795     *similarity_image;
1796
1797   RectangleInfo
1798     geometry;
1799
1800   SetGeometry(reference,&geometry);
1801   geometry.x=x_offset;
1802   geometry.y=y_offset;
1803   similarity_image=CropImage(image,&geometry,exception);
1804   if (similarity_image == (Image *) NULL)
1805     return(0.0);
1806   distortion=GetNCCDistortion(reference,similarity_image,reference_statistics,
1807     exception);
1808   similarity_image=DestroyImage(similarity_image);
1809   return(distortion);
1810 }
1811
1812 MagickExport Image *SimilarityImage(Image *image,const Image *reference,
1813   RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
1814 {
1815 #define SimilarityImageTag  "Similarity/Image"
1816
1817   CacheView
1818     *similarity_view;
1819
1820   ChannelStatistics
1821     *reference_statistics;
1822
1823   Image
1824     *similarity_image;
1825
1826   MagickBooleanType
1827     status;
1828
1829   MagickOffsetType
1830     progress;
1831
1832   ssize_t
1833     y;
1834
1835   assert(image != (const Image *) NULL);
1836   assert(image->signature == MagickSignature);
1837   if (image->debug != MagickFalse)
1838     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1839   assert(exception != (ExceptionInfo *) NULL);
1840   assert(exception->signature == MagickSignature);
1841   assert(offset != (RectangleInfo *) NULL);
1842   SetGeometry(reference,offset);
1843   *similarity_metric=1.0;
1844   if ((reference->columns > image->columns) || (reference->rows > image->rows))
1845     ThrowImageException(ImageError,"ImageSizeDiffers");
1846   similarity_image=CloneImage(image,image->columns-reference->columns+1,
1847     image->rows-reference->rows+1,MagickTrue,exception);
1848   if (similarity_image == (Image *) NULL)
1849     return((Image *) NULL);
1850   if (SetImageStorageClass(similarity_image,DirectClass) == MagickFalse)
1851     {
1852       InheritException(exception,&similarity_image->exception);
1853       similarity_image=DestroyImage(similarity_image);
1854       return((Image *) NULL);
1855     }
1856   /*
1857     Measure similarity of reference image against image.
1858   */
1859   status=MagickTrue;
1860   progress=0;
1861   reference_statistics=GetImageChannelStatistics(reference,exception);
1862   similarity_view=AcquireCacheView(similarity_image);
1863 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1864   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1865 #endif
1866   for (y=0; y < (ssize_t) (image->rows-reference->rows+1); y++)
1867   {
1868     double
1869       similarity;
1870
1871     register ssize_t
1872       x;
1873
1874     register PixelPacket
1875       *restrict q;
1876
1877     if (status == MagickFalse)
1878       continue;
1879     q=GetCacheViewAuthenticPixels(similarity_view,0,y,similarity_image->columns,
1880       1,exception);
1881     if (q == (const PixelPacket *) NULL)
1882       {
1883         status=MagickFalse;
1884         continue;
1885       }
1886     for (x=0; x < (ssize_t) (image->columns-reference->columns+1); x++)
1887     {
1888       similarity=GetSimilarityMetric(image,reference,reference_statistics,x,y,
1889         exception);
1890 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1891   #pragma omp critical (MagickCore_SimilarityImage)
1892 #endif
1893       if (similarity < *similarity_metric)
1894         {
1895           *similarity_metric=similarity;
1896           offset->x=x;
1897           offset->y=y;
1898         }
1899       q->red=ClampToQuantum(QuantumRange-QuantumRange*similarity);
1900       q->green=q->red;
1901       q->blue=q->red;
1902       q++;
1903     }
1904     if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
1905       status=MagickFalse;
1906     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1907       {
1908         MagickBooleanType
1909           proceed;
1910
1911 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1912   #pragma omp critical (MagickCore_SimilarityImage)
1913 #endif
1914         proceed=SetImageProgress(image,SimilarityImageTag,progress++,
1915           image->rows);
1916         if (proceed == MagickFalse)
1917           status=MagickFalse;
1918       }
1919   }
1920   similarity_view=DestroyCacheView(similarity_view);
1921   reference_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1922     reference_statistics);
1923   return(similarity_image);
1924 }