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