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-2010 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   CacheView
 121     *highlight_view,
 122     *image_view,
 123     *reconstruct_view;
 124
 125   const char
 126     *artifact;
 127
 128   Image
 129     *difference_image,
 130     *highlight_image;
 131
 132   ssize_t
 133     y;
 134
 135   MagickBooleanType
 136     status;
 137
 138   MagickPixelPacket
 139     highlight,
 140     lowlight,
 141     zero;
 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(dynamic,4) shared(status)
 202 #endif
 203   for (y=0; y < (ssize_t) 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 ssize_t
 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 < (ssize_t) 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   CacheView
 353     *image_view,
 354     *reconstruct_view;
 355
 356   ssize_t
 357     y;
 358
 359   MagickBooleanType
 360     status;
 361
 362   MagickPixelPacket
 363     zero;
 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(dynamic,4) shared(status)
 374 #endif
 375   for (y=0; y < (ssize_t) 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 ssize_t
 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 < (ssize_t) 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 <= (ssize_t) 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 size_t GetNumberChannels(const Image *image,
 447   const ChannelType channel)
 448 {
 449   size_t
 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   CacheView
 473     *image_view,
 474     *reconstruct_view;
 475
 476   ssize_t
 477     y;
 478
 479   MagickBooleanType
 480     status;
 481
 482   register ssize_t
 483     i;
 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(dynamic,4) shared(status)
 490 #endif
 491   for (y=0; y < (ssize_t) 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 ssize_t
 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 < (ssize_t) 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 <= (ssize_t) 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 <= (ssize_t) 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   CacheView
 581     *image_view,
 582     *reconstruct_view;
 583
 584   ssize_t
 585     y;
 586
 587   MagickBooleanType
 588     status;
 589
 590   MagickRealType
 591     alpha,
 592     area,
 593     beta,
 594     maximum_error,
 595     mean_error;
 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 < (ssize_t) 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 ssize_t
 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 < (ssize_t) image->columns; x++)
 629     {
 630       MagickRealType
 631         distance;
 632
 633       if ((channel & OpacityChannel) != 0)
 634         {
 635           if (image->matte != MagickFalse)
 636             alpha=(MagickRealType) (QuantumScale*(GetAlphaPixelComponent(p)));
 637           if (reconstruct_image->matte != MagickFalse)
 638             beta=(MagickRealType) (QuantumScale*GetAlphaPixelComponent(q));
 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   CacheView
 710     *image_view,
 711     *reconstruct_view;
 712
 713   ssize_t
 714     y;
 715
 716   MagickBooleanType
 717     status;
 718
 719   register ssize_t
 720     i;
 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(dynamic,4) shared(status)
 727 #endif
 728   for (y=0; y < (ssize_t) 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 ssize_t
 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 < (ssize_t) 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 <= (ssize_t) 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 <= (ssize_t) 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 GetNormalizedCrossCorrelationError(const Image *image,
 815   const Image *reconstruct_image,const ChannelType channel,
 816   double *distortion,ExceptionInfo *exception)
 817 {
 818   CacheView
 819     *image_view,
 820     *reconstruct_view;
 821
 822   ssize_t
 823     y;
 824
 825   MagickBooleanType
 826     status;
 827
 828   register ssize_t
 829     i;
 830
 831   status=MagickTrue;
 832   image_view=AcquireCacheView(image);
 833   reconstruct_view=AcquireCacheView(reconstruct_image);
 834 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 835   #pragma omp parallel for schedule(dynamic,4) shared(status)
 836 #endif
 837   for (y=0; y < (ssize_t) image->rows; y++)
 838   {
 839     double
 840       channel_distortion[AllChannels+1];
 841
 842     register const IndexPacket
 843       *restrict indexes,
 844       *restrict reconstruct_indexes;
 845
 846     register const PixelPacket
 847       *restrict p,
 848       *restrict q;
 849
 850     register ssize_t
 851       i,
 852       x;
 853
 854     if (status == MagickFalse)
 855       continue;
 856     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
 857     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,
 858       reconstruct_image->columns,1,exception);
 859     if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
 860       {
 861         status=MagickFalse;
 862         continue;
 863       }
 864     indexes=GetCacheViewVirtualIndexQueue(image_view);
 865     reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
 866     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
 867     for (x=0; x < (ssize_t) image->columns; x++)
 868     {
 869       MagickRealType
 870         distance;
 871
 872       if ((channel & RedChannel) != 0)
 873         {
 874           distance=QuantumScale*fabs(p->red-(double) q->red);
 875           channel_distortion[RedChannel]+=distance;
 876           channel_distortion[AllChannels]+=distance;
 877         }
 878       if ((channel & GreenChannel) != 0)
 879         {
 880           distance=QuantumScale*fabs(p->green-(double) q->green);
 881           channel_distortion[GreenChannel]+=distance;
 882           channel_distortion[AllChannels]+=distance;
 883         }
 884       if ((channel & BlueChannel) != 0)
 885         {
 886           distance=QuantumScale*fabs(p->blue-(double) q->blue);
 887           channel_distortion[BlueChannel]+=distance;
 888           channel_distortion[AllChannels]+=distance;
 889         }
 890       if (((channel & OpacityChannel) != 0) &&
 891           (image->matte != MagickFalse))
 892         {
 893           distance=QuantumScale*fabs(p->opacity-(double) q->opacity);
 894           channel_distortion[OpacityChannel]+=distance;
 895           channel_distortion[AllChannels]+=distance;
 896         }
 897       if (((channel & IndexChannel) != 0) &&
 898           (image->colorspace == CMYKColorspace))
 899         {
 900           distance=QuantumScale*fabs(indexes[x]-(double)
 901             reconstruct_indexes[x]);
 902           channel_distortion[BlackChannel]+=distance;
 903           channel_distortion[AllChannels]+=distance;
 904         }
 905       p++;
 906       q++;
 907     }
 908 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 909   #pragma omp critical (MagickCore_GetNormalizedCrossCorrelationError)
 910 #endif
 911     for (i=0; i <= (ssize_t) AllChannels; i++)
 912       distortion[i]+=channel_distortion[i];
 913   }
 914   reconstruct_view=DestroyCacheView(reconstruct_view);
 915   image_view=DestroyCacheView(image_view);
 916   for (i=0; i <= (ssize_t) AllChannels; i++)
 917     distortion[i]/=((double) image->columns*image->rows);
 918   distortion[AllChannels]/=(double) GetNumberChannels(image,channel);
 919   return(status);
 920 }
 921
 922 static MagickBooleanType GetPeakAbsoluteError(const Image *image,
 923   const Image *reconstruct_image,const ChannelType channel,
 924   double *distortion,ExceptionInfo *exception)
 925 {
 926   CacheView
 927     *image_view,
 928     *reconstruct_view;
 929
 930   ssize_t
 931     y;
 932
 933   MagickBooleanType
 934     status;
 935
 936   status=MagickTrue;
 937   image_view=AcquireCacheView(image);
 938   reconstruct_view=AcquireCacheView(reconstruct_image);
 939 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 940   #pragma omp parallel for schedule(dynamic,4) shared(status)
 941 #endif
 942   for (y=0; y < (ssize_t) image->rows; y++)
 943   {
 944     double
 945       channel_distortion[AllChannels+1];
 946
 947     register const IndexPacket
 948       *restrict indexes,
 949       *restrict reconstruct_indexes;
 950
 951     register const PixelPacket
 952       *restrict p,
 953       *restrict q;
 954
 955     register ssize_t
 956       i,
 957       x;
 958
 959     if (status == MagickFalse)
 960       continue;
 961     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
 962     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,
 963       reconstruct_image->columns,1,exception);
 964     if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
 965       {
 966         status=MagickFalse;
 967         continue;
 968       }
 969     indexes=GetCacheViewVirtualIndexQueue(image_view);
 970     reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
 971     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
 972     for (x=0; x < (ssize_t) image->columns; x++)
 973     {
 974       MagickRealType
 975         distance;
 976
 977       if ((channel & RedChannel) != 0)
 978         {
 979           distance=QuantumScale*fabs(p->red-(double) q->red);
 980           if (distance > channel_distortion[RedChannel])
 981             channel_distortion[RedChannel]=distance;
 982           if (distance > channel_distortion[AllChannels])
 983             channel_distortion[AllChannels]=distance;
 984         }
 985       if ((channel & GreenChannel) != 0)
 986         {
 987           distance=QuantumScale*fabs(p->green-(double) q->green);
 988           if (distance > channel_distortion[GreenChannel])
 989             channel_distortion[GreenChannel]=distance;
 990           if (distance > channel_distortion[AllChannels])
 991             channel_distortion[AllChannels]=distance;
 992         }
 993       if ((channel & BlueChannel) != 0)
 994         {
 995           distance=QuantumScale*fabs(p->blue-(double) q->blue);
 996           if (distance > channel_distortion[BlueChannel])
 997             channel_distortion[BlueChannel]=distance;
 998           if (distance > channel_distortion[AllChannels])
 999             channel_distortion[AllChannels]=distance;
1000         }
1001       if (((channel & OpacityChannel) != 0) &&
1002           (image->matte != MagickFalse))
1003         {
1004           distance=QuantumScale*fabs(p->opacity-(double) q->opacity);
1005           if (distance > channel_distortion[OpacityChannel])
1006             channel_distortion[OpacityChannel]=distance;
1007           if (distance > channel_distortion[AllChannels])
1008             channel_distortion[AllChannels]=distance;
1009         }
1010       if (((channel & IndexChannel) != 0) &&
1011           (image->colorspace == CMYKColorspace) &&
1012           (reconstruct_image->colorspace == CMYKColorspace))
1013         {
1014           distance=QuantumScale*fabs(indexes[x]-(double)
1015             reconstruct_indexes[x]);
1016           if (distance > channel_distortion[BlackChannel])
1017             channel_distortion[BlackChannel]=distance;
1018           if (distance > channel_distortion[AllChannels])
1019             channel_distortion[AllChannels]=distance;
1020         }
1021       p++;
1022       q++;
1023     }
1024 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1025   #pragma omp critical (MagickCore_GetPeakAbsoluteError)
1026 #endif
1027     for (i=0; i <= (ssize_t) AllChannels; i++)
1028       if (channel_distortion[i] > distortion[i])
1029         distortion[i]=channel_distortion[i];
1030   }
1031   reconstruct_view=DestroyCacheView(reconstruct_view);
1032   image_view=DestroyCacheView(image_view);
1033   return(status);
1034 }
1035
1036 static MagickBooleanType GetPeakSignalToNoiseRatio(const Image *image,
1037   const Image *reconstruct_image,const ChannelType channel,
1038   double *distortion,ExceptionInfo *exception)
1039 {
1040   MagickBooleanType
1041     status;
1042
1043   status=GetMeanSquaredError(image,reconstruct_image,channel,distortion,
1044     exception);
1045   if ((channel & RedChannel) != 0)
1046     distortion[RedChannel]=20.0*log10((double) 1.0/sqrt(
1047       distortion[RedChannel]));
1048   if ((channel & GreenChannel) != 0)
1049     distortion[GreenChannel]=20.0*log10((double) 1.0/sqrt(
1050       distortion[GreenChannel]));
1051   if ((channel & BlueChannel) != 0)
1052     distortion[BlueChannel]=20.0*log10((double) 1.0/sqrt(
1053       distortion[BlueChannel]));
1054   if (((channel & OpacityChannel) != 0) &&
1055       (image->matte != MagickFalse))
1056     distortion[OpacityChannel]=20.0*log10((double) 1.0/sqrt(
1057       distortion[OpacityChannel]));
1058   if (((channel & IndexChannel) != 0) &&
1059       (image->colorspace == CMYKColorspace))
1060     distortion[BlackChannel]=20.0*log10((double) 1.0/sqrt(
1061       distortion[BlackChannel]));
1062   distortion[AllChannels]=20.0*log10((double) 1.0/sqrt(
1063     distortion[AllChannels]));
1064   return(status);
1065 }
1066
1067 static MagickBooleanType GetRootMeanSquaredError(const Image *image,
1068   const Image *reconstruct_image,const ChannelType channel,
1069   double *distortion,ExceptionInfo *exception)
1070 {
1071   MagickBooleanType
1072     status;
1073
1074   status=GetMeanSquaredError(image,reconstruct_image,channel,distortion,
1075     exception);
1076   if ((channel & RedChannel) != 0)
1077     distortion[RedChannel]=sqrt(distortion[RedChannel]);
1078   if ((channel & GreenChannel) != 0)
1079     distortion[GreenChannel]=sqrt(distortion[GreenChannel]);
1080   if ((channel & BlueChannel) != 0)
1081     distortion[BlueChannel]=sqrt(distortion[BlueChannel]);
1082   if (((channel & OpacityChannel) != 0) &&
1083       (image->matte != MagickFalse))
1084     distortion[OpacityChannel]=sqrt(distortion[OpacityChannel]);
1085   if (((channel & IndexChannel) != 0) &&
1086       (image->colorspace == CMYKColorspace))
1087     distortion[BlackChannel]=sqrt(distortion[BlackChannel]);
1088   distortion[AllChannels]=sqrt(distortion[AllChannels]);
1089   return(status);
1090 }
1091
1092 MagickExport MagickBooleanType GetImageChannelDistortion(Image *image,
1093   const Image *reconstruct_image,const ChannelType channel,
1094   const MetricType metric,double *distortion,ExceptionInfo *exception)
1095 {
1096   double
1097     *channel_distortion;
1098
1099   MagickBooleanType
1100     status;
1101
1102   size_t
1103     length;
1104
1105   assert(image != (Image *) NULL);
1106   assert(image->signature == MagickSignature);
1107   if (image->debug != MagickFalse)
1108     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1109   assert(reconstruct_image != (const Image *) NULL);
1110   assert(reconstruct_image->signature == MagickSignature);
1111   assert(distortion != (double *) NULL);
1112   *distortion=0.0;
1113   if (image->debug != MagickFalse)
1114     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1115   if ((reconstruct_image->columns != image->columns) ||
1116       (reconstruct_image->rows != image->rows))
1117     ThrowBinaryException(ImageError,"ImageSizeDiffers",image->filename);
1118   /*
1119     Get image distortion.
1120   */
1121   length=AllChannels+1UL;
1122   channel_distortion=(double *) AcquireQuantumMemory(length,
1123     sizeof(*channel_distortion));
1124   if (channel_distortion == (double *) NULL)
1125     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1126   (void) ResetMagickMemory(channel_distortion,0,length*
1127     sizeof(*channel_distortion));
1128   switch (metric)
1129   {
1130     case AbsoluteErrorMetric:
1131     {
1132       status=GetAbsoluteError(image,reconstruct_image,channel,
1133         channel_distortion,exception);
1134       break;
1135     }
1136     case MeanAbsoluteErrorMetric:
1137     {
1138       status=GetMeanAbsoluteError(image,reconstruct_image,channel,
1139         channel_distortion,exception);
1140       break;
1141     }
1142     case MeanErrorPerPixelMetric:
1143     {
1144       status=GetMeanErrorPerPixel(image,reconstruct_image,channel,
1145         channel_distortion,exception);
1146       break;
1147     }
1148     case MeanSquaredErrorMetric:
1149     {
1150       status=GetMeanSquaredError(image,reconstruct_image,channel,
1151         channel_distortion,exception);
1152       break;
1153     }
1154     case NormalizedCrossCorrelationErrorMetric:
1155     {
1156       status=GetNormalizedCrossCorrelationError(image,reconstruct_image,channel,
1157         channel_distortion,exception);
1158       break;
1159     }
1160     case PeakAbsoluteErrorMetric:
1161     default:
1162     {
1163       status=GetPeakAbsoluteError(image,reconstruct_image,channel,
1164         channel_distortion,exception);
1165       break;
1166     }
1167     case PeakSignalToNoiseRatioMetric:
1168     {
1169       status=GetPeakSignalToNoiseRatio(image,reconstruct_image,channel,
1170         channel_distortion,exception);
1171       break;
1172     }
1173     case RootMeanSquaredErrorMetric:
1174     {
1175       status=GetRootMeanSquaredError(image,reconstruct_image,channel,
1176         channel_distortion,exception);
1177       break;
1178     }
1179   }
1180   *distortion=channel_distortion[AllChannels];
1181   channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1182   return(status);
1183 }
1184 \f
1185 /*
1186 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1187 %                                                                             %
1188 %                                                                             %
1189 %                                                                             %
1190 %   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                       %
1191 %                                                                             %
1192 %                                                                             %
1193 %                                                                             %
1194 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1195 %
1196 %  GetImageChannelDistrortion() compares the image channels of an image to a
1197 %  reconstructed image and returns the specified distortion metric for each
1198 %  channel.
1199 %
1200 %  The format of the CompareImageChannels method is:
1201 %
1202 %      double *GetImageChannelDistortions(const Image *image,
1203 %        const Image *reconstruct_image,const MetricType metric,
1204 %        ExceptionInfo *exception)
1205 %
1206 %  A description of each parameter follows:
1207 %
1208 %    o image: the image.
1209 %
1210 %    o reconstruct_image: the reconstruct image.
1211 %
1212 %    o metric: the metric.
1213 %
1214 %    o exception: return any errors or warnings in this structure.
1215 %
1216 */
1217 MagickExport double *GetImageChannelDistortions(Image *image,
1218   const Image *reconstruct_image,const MetricType metric,
1219   ExceptionInfo *exception)
1220 {
1221   double
1222     *channel_distortion;
1223
1224   MagickBooleanType
1225     status;
1226
1227   size_t
1228     length;
1229
1230   assert(image != (Image *) NULL);
1231   assert(image->signature == MagickSignature);
1232   if (image->debug != MagickFalse)
1233     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1234   assert(reconstruct_image != (const Image *) NULL);
1235   assert(reconstruct_image->signature == MagickSignature);
1236   if (image->debug != MagickFalse)
1237     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1238   if ((reconstruct_image->columns != image->columns) ||
1239       (reconstruct_image->rows != image->rows))
1240     {
1241       (void) ThrowMagickException(&image->exception,GetMagickModule(),
1242         ImageError,"ImageSizeDiffers","`%s'",image->filename);
1243       return((double *) NULL);
1244     }
1245   /*
1246     Get image distortion.
1247   */
1248   length=AllChannels+1UL;
1249   channel_distortion=(double *) AcquireQuantumMemory(length,
1250     sizeof(*channel_distortion));
1251   if (channel_distortion == (double *) NULL)
1252     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1253   (void) ResetMagickMemory(channel_distortion,0,length*
1254     sizeof(*channel_distortion));
1255   switch (metric)
1256   {
1257     case AbsoluteErrorMetric:
1258     {
1259       status=GetAbsoluteError(image,reconstruct_image,AllChannels,
1260         channel_distortion,exception);
1261       break;
1262     }
1263     case MeanAbsoluteErrorMetric:
1264     {
1265       status=GetMeanAbsoluteError(image,reconstruct_image,AllChannels,
1266         channel_distortion,exception);
1267       break;
1268     }
1269     case MeanErrorPerPixelMetric:
1270     {
1271       status=GetMeanErrorPerPixel(image,reconstruct_image,AllChannels,
1272         channel_distortion,exception);
1273       break;
1274     }
1275     case MeanSquaredErrorMetric:
1276     {
1277       status=GetMeanSquaredError(image,reconstruct_image,AllChannels,
1278         channel_distortion,exception);
1279       break;
1280     }
1281     case NormalizedCrossCorrelationErrorMetric:
1282     {
1283       status=GetNormalizedCrossCorrelationError(image,reconstruct_image,
1284         AllChannels,channel_distortion,exception);
1285       break;
1286     }
1287     case PeakAbsoluteErrorMetric:
1288     default:
1289     {
1290       status=GetPeakAbsoluteError(image,reconstruct_image,AllChannels,
1291         channel_distortion,exception);
1292       break;
1293     }
1294     case PeakSignalToNoiseRatioMetric:
1295     {
1296       status=GetPeakSignalToNoiseRatio(image,reconstruct_image,AllChannels,
1297         channel_distortion,exception);
1298       break;
1299     }
1300     case RootMeanSquaredErrorMetric:
1301     {
1302       status=GetRootMeanSquaredError(image,reconstruct_image,AllChannels,
1303         channel_distortion,exception);
1304       break;
1305     }
1306   }
1307   return(channel_distortion);
1308 }
1309 \f
1310 /*
1311 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1312 %                                                                             %
1313 %                                                                             %
1314 %                                                                             %
1315 %  I s I m a g e s E q u a l                                                  %
1316 %                                                                             %
1317 %                                                                             %
1318 %                                                                             %
1319 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1320 %
1321 %  IsImagesEqual() measures the difference between colors at each pixel
1322 %  location of two images.  A value other than 0 means the colors match
1323 %  exactly.  Otherwise an error measure is computed by summing over all
1324 %  pixels in an image the distance squared in RGB space between each image
1325 %  pixel and its corresponding pixel in the reconstruct image.  The error
1326 %  measure is assigned to these image members:
1327 %
1328 %    o mean_error_per_pixel:  The mean error for any single pixel in
1329 %      the image.
1330 %
1331 %    o normalized_mean_error:  The normalized mean quantization error for
1332 %      any single pixel in the image.  This distance measure is normalized to
1333 %      a range between 0 and 1.  It is independent of the range of red, green,
1334 %      and blue values in the image.
1335 %
1336 %    o normalized_maximum_error:  The normalized maximum quantization
1337 %      error for any single pixel in the image.  This distance measure is
1338 %      normalized to a range between 0 and 1.  It is independent of the range
1339 %      of red, green, and blue values in your image.
1340 %
1341 %  A small normalized mean square error, accessed as
1342 %  image->normalized_mean_error, suggests the images are very similar in
1343 %  spatial layout and color.
1344 %
1345 %  The format of the IsImagesEqual method is:
1346 %
1347 %      MagickBooleanType IsImagesEqual(Image *image,
1348 %        const Image *reconstruct_image)
1349 %
1350 %  A description of each parameter follows.
1351 %
1352 %    o image: the image.
1353 %
1354 %    o reconstruct_image: the reconstruct image.
1355 %
1356 */
1357 MagickExport MagickBooleanType IsImagesEqual(Image *image,
1358   const Image *reconstruct_image)
1359 {
1360   CacheView
1361     *image_view,
1362     *reconstruct_view;
1363
1364   ExceptionInfo
1365     *exception;
1366
1367   ssize_t
1368     y;
1369
1370   MagickBooleanType
1371     status;
1372
1373   MagickRealType
1374     area,
1375     maximum_error,
1376     mean_error,
1377     mean_error_per_pixel;
1378
1379   assert(image != (Image *) NULL);
1380   assert(image->signature == MagickSignature);
1381   assert(reconstruct_image != (const Image *) NULL);
1382   assert(reconstruct_image->signature == MagickSignature);
1383   if ((reconstruct_image->columns != image->columns) ||
1384       (reconstruct_image->rows != image->rows))
1385     ThrowBinaryException(ImageError,"ImageSizeDiffers",image->filename);
1386   area=0.0;
1387   maximum_error=0.0;
1388   mean_error_per_pixel=0.0;
1389   mean_error=0.0;
1390   exception=(&image->exception);
1391   image_view=AcquireCacheView(image);
1392   reconstruct_view=AcquireCacheView(reconstruct_image);
1393   for (y=0; y < (ssize_t) image->rows; y++)
1394   {
1395     register const IndexPacket
1396       *restrict indexes,
1397       *restrict reconstruct_indexes;
1398
1399     register const PixelPacket
1400       *restrict p,
1401       *restrict q;
1402
1403     register ssize_t
1404       x;
1405
1406     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1407     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
1408       1,exception);
1409     if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1410       break;
1411     indexes=GetCacheViewVirtualIndexQueue(image_view);
1412     reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
1413     for (x=0; x < (ssize_t) image->columns; x++)
1414     {
1415       MagickRealType
1416         distance;
1417
1418       distance=fabs(p->red-(double) q->red);
1419       mean_error_per_pixel+=distance;
1420       mean_error+=distance*distance;
1421       if (distance > maximum_error)
1422         maximum_error=distance;
1423       area++;
1424       distance=fabs(p->green-(double) q->green);
1425       mean_error_per_pixel+=distance;
1426       mean_error+=distance*distance;
1427       if (distance > maximum_error)
1428         maximum_error=distance;
1429       area++;
1430       distance=fabs(p->blue-(double) q->blue);
1431       mean_error_per_pixel+=distance;
1432       mean_error+=distance*distance;
1433       if (distance > maximum_error)
1434         maximum_error=distance;
1435       area++;
1436       if (image->matte != MagickFalse)
1437         {
1438           distance=fabs(p->opacity-(double) q->opacity);
1439           mean_error_per_pixel+=distance;
1440           mean_error+=distance*distance;
1441           if (distance > maximum_error)
1442             maximum_error=distance;
1443           area++;
1444         }
1445       if ((image->colorspace == CMYKColorspace) &&
1446           (reconstruct_image->colorspace == CMYKColorspace))
1447         {
1448           distance=fabs(indexes[x]-(double) reconstruct_indexes[x]);
1449           mean_error_per_pixel+=distance;
1450           mean_error+=distance*distance;
1451           if (distance > maximum_error)
1452             maximum_error=distance;
1453           area++;
1454         }
1455       p++;
1456       q++;
1457     }
1458   }
1459   reconstruct_view=DestroyCacheView(reconstruct_view);
1460   image_view=DestroyCacheView(image_view);
1461   image->error.mean_error_per_pixel=(double) (mean_error_per_pixel/area);
1462   image->error.normalized_mean_error=(double) (QuantumScale*QuantumScale*
1463     mean_error/area);
1464   image->error.normalized_maximum_error=(double) (QuantumScale*maximum_error);
1465   status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
1466   return(status);
1467 }
1468 \f
1469 /*
1470 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1471 %                                                                             %
1472 %                                                                             %
1473 %                                                                             %
1474 %   S i m i l a r i t y I m a g e                                             %
1475 %                                                                             %
1476 %                                                                             %
1477 %                                                                             %
1478 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1479 %
1480 %  SimilarityImage() compares the reference image of the image and returns the
1481 %  best match offset.  In addition, it returns a similarity image such that an
1482 %  exact match location is completely white and if none of the pixels match,
1483 %  black, otherwise some gray level in-between.
1484 %
1485 %  The format of the SimilarityImageImage method is:
1486 %
1487 %      Image *SimilarityImage(const Image *image,const Image *reference,
1488 %        RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
1489 %
1490 %  A description of each parameter follows:
1491 %
1492 %    o image: the image.
1493 %
1494 %    o reference: find an area of the image that closely resembles this image.
1495 %
1496 %    o the best match offset of the reference image within the image.
1497 %
1498 %    o similarity: the computed similarity between the images.
1499 %
1500 %    o exception: return any errors or warnings in this structure.
1501 %
1502 */
1503
1504 static double GetSimilarityMetric(const Image *image,const Image *reference,
1505   const ssize_t x_offset,const ssize_t y_offset,ExceptionInfo *exception)
1506 {
1507   CacheView
1508     *image_view,
1509     *reference_view;
1510
1511   double
1512     similarity;
1513
1514   ssize_t
1515     y;
1516
1517   MagickBooleanType
1518     status;
1519
1520   /*
1521     Compute the similarity in pixels between two images.
1522   */
1523   status=MagickTrue;
1524   similarity=0.0;
1525   image_view=AcquireCacheView(image);
1526   reference_view=AcquireCacheView(reference);
1527   for (y=0; y < (ssize_t) reference->rows; y++)
1528   {
1529     register const IndexPacket
1530       *restrict indexes,
1531       *restrict reference_indexes;
1532
1533     register const PixelPacket
1534       *restrict p,
1535       *restrict q;
1536
1537     register ssize_t
1538       x;
1539
1540     if (status == MagickFalse)
1541       continue;
1542     p=GetCacheViewVirtualPixels(image_view,x_offset,y_offset+y,
1543       reference->columns,1,exception);
1544     q=GetCacheViewVirtualPixels(reference_view,0,y,reference->columns,1,
1545       exception);
1546     if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1547       {
1548         status=MagickFalse;
1549         continue;
1550       }
1551     indexes=GetCacheViewVirtualIndexQueue(image_view);
1552     reference_indexes=GetCacheViewVirtualIndexQueue(reference_view);
1553     for (x=0; x < (ssize_t) reference->columns; x++)
1554     {
1555       double
1556         thread_similarity;
1557
1558       MagickRealType
1559         distance;
1560
1561       thread_similarity=0.0;
1562       distance=QuantumScale*(p->red-(MagickRealType) q->red);
1563       thread_similarity+=distance*distance;
1564       distance=QuantumScale*(p->green-(MagickRealType) q->green);
1565       thread_similarity+=distance*distance;
1566       distance=QuantumScale*(p->blue-(MagickRealType) q->blue);
1567       thread_similarity+=distance*distance;
1568       if ((image->matte != MagickFalse) && (reference->matte != MagickFalse))
1569         {
1570           distance=QuantumScale*(p->opacity-(MagickRealType) q->opacity);
1571           thread_similarity+=distance*distance;
1572         }
1573       if ((image->colorspace == CMYKColorspace) &&
1574           (reference->colorspace == CMYKColorspace))
1575         {
1576           distance=QuantumScale*(indexes[x]-(MagickRealType)
1577             reference_indexes[x]);
1578           thread_similarity+=distance*distance;
1579         }
1580       similarity+=thread_similarity;
1581       p++;
1582       q++;
1583     }
1584   }
1585   reference_view=DestroyCacheView(reference_view);
1586   image_view=DestroyCacheView(image_view);
1587   if (status == MagickFalse)
1588     return(0.0);
1589   similarity/=((double) reference->columns*reference->rows);
1590   similarity/=(double) GetNumberChannels(reference,AllChannels);
1591   return(sqrt(similarity));
1592 }
1593
1594 MagickExport Image *SimilarityImage(Image *image,const Image *reference,
1595   RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
1596 {
1597 #define SimilarityImageTag  "Similarity/Image"
1598
1599   CacheView
1600     *similarity_view;
1601
1602   Image
1603     *similarity_image;
1604
1605   MagickBooleanType
1606     status;
1607
1608   MagickOffsetType
1609     progress;
1610
1611   ssize_t
1612     y;
1613
1614   assert(image != (const Image *) NULL);
1615   assert(image->signature == MagickSignature);
1616   if (image->debug != MagickFalse)
1617     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1618   assert(exception != (ExceptionInfo *) NULL);
1619   assert(exception->signature == MagickSignature);
1620   assert(offset != (RectangleInfo *) NULL);
1621   SetGeometry(reference,offset);
1622   *similarity_metric=1.0;
1623   if ((reference->columns > image->columns) || (reference->rows > image->rows))
1624     ThrowImageException(ImageError,"ImageSizeDiffers");
1625   similarity_image=CloneImage(image,image->columns-reference->columns+1,
1626     image->rows-reference->rows+1,MagickTrue,exception);
1627   if (similarity_image == (Image *) NULL)
1628     return((Image *) NULL);
1629   if (SetImageStorageClass(similarity_image,DirectClass) == MagickFalse)
1630     {
1631       InheritException(exception,&similarity_image->exception);
1632       similarity_image=DestroyImage(similarity_image);
1633       return((Image *) NULL);
1634     }
1635   /*
1636     Measure similarity of reference image against image.
1637   */
1638   status=MagickTrue;
1639   progress=0;
1640   similarity_view=AcquireCacheView(similarity_image);
1641 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1642   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1643 #endif
1644   for (y=0; y < (ssize_t) (image->rows-reference->rows+1); y++)
1645   {
1646     double
1647       similarity;
1648
1649     register ssize_t
1650       x;
1651
1652     register PixelPacket
1653       *restrict q;
1654
1655     if (status == MagickFalse)
1656       continue;
1657     q=GetCacheViewAuthenticPixels(similarity_view,0,y,
1658       similarity_image->columns,1,exception);
1659     if (q == (const PixelPacket *) NULL)
1660       {
1661         status=MagickFalse;
1662         continue;
1663       }
1664     for (x=0; x < (ssize_t) (image->columns-reference->columns+1); x++)
1665     {
1666       similarity=GetSimilarityMetric(image,reference,x,y,exception);
1667 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1668   #pragma omp critical (MagickCore_SimilarityImage)
1669 #endif
1670       if (similarity < *similarity_metric)
1671         {
1672           *similarity_metric=similarity;
1673           offset->x=x;
1674           offset->y=y;
1675         }
1676       q->red=ClampToQuantum(QuantumRange-QuantumRange*similarity);
1677       q->green=q->red;
1678       q->blue=q->red;
1679       q++;
1680     }
1681     if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
1682       status=MagickFalse;
1683     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1684       {
1685         MagickBooleanType
1686           proceed;
1687
1688 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1689   #pragma omp critical (MagickCore_SimilarityImage)
1690 #endif
1691         proceed=SetImageProgress(image,SimilarityImageTag,progress++,
1692           image->rows);
1693         if (proceed == MagickFalse)
1694           status=MagickFalse;
1695       }
1696   }
1697   similarity_view=DestroyCacheView(similarity_view);
1698   return(similarity_image);
1699 }