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