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