granicus.if.org Git - imagemagick/blob - magick/compare.c

   1 /*
   2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   3 %                                                                             %
   4 %                                                                             %
   5 %                                                                             %
   6 %               CCCC   OOO   M   M  PPPP    AAA   RRRR    EEEEE               %
   7 %              C      O   O  MM MM  P   P  A   A  R   R   E                   %
   8 %              C      O   O  M M M  PPPP   AAAAA  RRRR    EEE                 %
   9 %              C      O   O  M   M  P      A   A  R R     E                   %
  10 %               CCCC   OOO   M   M  P      A   A  R  R    EEEEE               %
  11 %                                                                             %
  12 %                                                                             %
  13 %                      MagickCore Image Comparison Methods                    %
  14 %                                                                             %
  15 %                              Software Design                                %
  16 %                                John Cristy                                  %
  17 %                               December 2003                                 %
  18 %                                                                             %
  19 %                                                                             %
  20 %  Copyright 1999-2011 ImageMagick Studio LLC, a non-profit organization      %
  21 %  dedicated to making software imaging solutions freely available.           %
  22 %                                                                             %
  23 %  You may not use this file except in compliance with the License.  You may  %
  24 %  obtain a copy of the License at                                            %
  25 %                                                                             %
  26 %    http://www.imagemagick.org/script/license.php                            %
  27 %                                                                             %
  28 %  Unless required by applicable law or agreed to in writing, software        %
  29 %  distributed under the License is distributed on an "AS IS" BASIS,          %
  30 %  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
  31 %  See the License for the specific language governing permissions and        %
  32 %  limitations under the License.                                             %
  33 %                                                                             %
  34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  35 %
  36 %
  37 %
  38 */
  39 \f
  40 /*
  41   Include declarations.
  42 */
  43 #include "magick/studio.h"
  44 #include "magick/artifact.h"
  45 #include "magick/cache-view.h"
  46 #include "magick/client.h"
  47 #include "magick/color.h"
  48 #include "magick/color-private.h"
  49 #include "magick/colorspace.h"
  50 #include "magick/colorspace-private.h"
  51 #include "magick/compare.h"
  52 #include "magick/composite-private.h"
  53 #include "magick/constitute.h"
  54 #include "magick/exception-private.h"
  55 #include "magick/geometry.h"
  56 #include "magick/image-private.h"
  57 #include "magick/list.h"
  58 #include "magick/log.h"
  59 #include "magick/memory_.h"
  60 #include "magick/monitor.h"
  61 #include "magick/monitor-private.h"
  62 #include "magick/option.h"
  63 #include "magick/pixel-private.h"
  64 #include "magick/resource_.h"
  65 #include "magick/string_.h"
  66 #include "magick/statistic.h"
  67 #include "magick/transform.h"
  68 #include "magick/utility.h"
  69 #include "magick/version.h"
  70 \f
  71 /*
  72 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  73 %                                                                             %
  74 %                                                                             %
  75 %                                                                             %
  76 %   C o m p a r e I m a g e C h a n n e l s                                   %
  77 %                                                                             %
  78 %                                                                             %
  79 %                                                                             %
  80 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  81 %
  82 %  CompareImageChannels() compares one or more image channels of an image
  83 %  to a reconstructed image and returns the difference image.
  84 %
  85 %  The format of the CompareImageChannels method is:
  86 %
  87 %      Image *CompareImageChannels(const Image *image,
  88 %        const Image *reconstruct_image,const ChannelType channel,
  89 %        const MetricType metric,double *distortion,ExceptionInfo *exception)
  90 %
  91 %  A description of each parameter follows:
  92 %
  93 %    o image: the image.
  94 %
  95 %    o reconstruct_image: the reconstruct image.
  96 %
  97 %    o channel: the channel.
  98 %
  99 %    o metric: the metric.
 100 %
 101 %    o distortion: the computed distortion between the images.
 102 %
 103 %    o exception: return any errors or warnings in this structure.
 104 %
 105 */
 106
 107 MagickExport Image *CompareImages(Image *image,const Image *reconstruct_image,
 108   const MetricType metric,double *distortion,ExceptionInfo *exception)
 109 {
 110   Image
 111     *highlight_image;
 112
 113   highlight_image=CompareImageChannels(image,reconstruct_image,AllChannels,
 114     metric,distortion,exception);
 115   return(highlight_image);
 116 }
 117
 118 MagickExport Image *CompareImageChannels(Image *image,
 119   const Image *reconstruct_image,const ChannelType channel,
 120   const MetricType metric,double *distortion,ExceptionInfo *exception)
 121 {
 122   CacheView
 123     *highlight_view,
 124     *image_view,
 125     *reconstruct_view;
 126
 127   const char
 128     *artifact;
 129
 130   Image
 131     *difference_image,
 132     *highlight_image;
 133
 134   ssize_t
 135     y;
 136
 137   MagickBooleanType
 138     status;
 139
 140   MagickPixelPacket
 141     highlight,
 142     lowlight,
 143     zero;
 144
 145   assert(image != (Image *) NULL);
 146   assert(image->signature == MagickSignature);
 147   if (image->debug != MagickFalse)
 148     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 149   assert(reconstruct_image != (const Image *) NULL);
 150   assert(reconstruct_image->signature == MagickSignature);
 151   assert(distortion != (double *) NULL);
 152   *distortion=0.0;
 153   if (image->debug != MagickFalse)
 154     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 155   if ((reconstruct_image->columns != image->columns) ||
 156       (reconstruct_image->rows != image->rows))
 157     ThrowImageException(ImageError,"ImageSizeDiffers");
 158   status=GetImageChannelDistortion(image,reconstruct_image,channel,metric,
 159     distortion,exception);
 160   if (status == MagickFalse)
 161     return((Image *) NULL);
 162   difference_image=CloneImage(image,0,0,MagickTrue,exception);
 163   if (difference_image == (Image *) NULL)
 164     return((Image *) NULL);
 165   (void) SetImageAlphaChannel(difference_image,OpaqueAlphaChannel);
 166   highlight_image=CloneImage(image,image->columns,image->rows,MagickTrue,
 167     exception);
 168   if (highlight_image == (Image *) NULL)
 169     {
 170       difference_image=DestroyImage(difference_image);
 171       return((Image *) NULL);
 172     }
 173   if (SetImageStorageClass(highlight_image,DirectClass) == MagickFalse)
 174     {
 175       InheritException(exception,&highlight_image->exception);
 176       difference_image=DestroyImage(difference_image);
 177       highlight_image=DestroyImage(highlight_image);
 178       return((Image *) NULL);
 179     }
 180   (void) SetImageAlphaChannel(highlight_image,OpaqueAlphaChannel);
 181   (void) QueryMagickColor("#f1001ecc",&highlight,exception);
 182   artifact=GetImageArtifact(image,"highlight-color");
 183   if (artifact != (const char *) NULL)
 184     (void) QueryMagickColor(artifact,&highlight,exception);
 185   (void) QueryMagickColor("#ffffffcc",&lowlight,exception);
 186   artifact=GetImageArtifact(image,"lowlight-color");
 187   if (artifact != (const char *) NULL)
 188     (void) QueryMagickColor(artifact,&lowlight,exception);
 189   if (highlight_image->colorspace == CMYKColorspace)
 190     {
 191       ConvertRGBToCMYK(&highlight);
 192       ConvertRGBToCMYK(&lowlight);
 193     }
 194   /*
 195     Generate difference image.
 196   */
 197   status=MagickTrue;
 198   GetMagickPixelPacket(image,&zero);
 199   image_view=AcquireCacheView(image);
 200   reconstruct_view=AcquireCacheView(reconstruct_image);
 201   highlight_view=AcquireCacheView(highlight_image);
 202 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 203   #pragma omp parallel for schedule(dynamic,4) shared(status)
 204 #endif
 205   for (y=0; y < (ssize_t) image->rows; y++)
 206   {
 207     MagickBooleanType
 208       sync;
 209
 210     MagickPixelPacket
 211       pixel,
 212       reconstruct_pixel;
 213
 214     register const IndexPacket
 215       *restrict indexes,
 216       *restrict reconstruct_indexes;
 217
 218     register const PixelPacket
 219       *restrict p,
 220       *restrict q;
 221
 222     register IndexPacket
 223       *restrict highlight_indexes;
 224
 225     register ssize_t
 226       x;
 227
 228     register PixelPacket
 229       *restrict r;
 230
 231     if (status == MagickFalse)
 232       continue;
 233     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
 234     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
 235       1,exception);
 236     r=QueueCacheViewAuthenticPixels(highlight_view,0,y,highlight_image->columns,
 237       1,exception);
 238     if ((p == (const PixelPacket *) NULL) ||
 239         (q == (const PixelPacket *) NULL) || (r == (PixelPacket *) NULL))
 240       {
 241         status=MagickFalse;
 242         continue;
 243       }
 244     indexes=GetCacheViewVirtualIndexQueue(image_view);
 245     reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
 246     highlight_indexes=GetCacheViewAuthenticIndexQueue(highlight_view);
 247     pixel=zero;
 248     reconstruct_pixel=zero;
 249     for (x=0; x < (ssize_t) image->columns; x++)
 250     {
 251       MagickStatusType
 252         difference;
 253
 254       SetMagickPixelPacket(image,p,indexes+x,&pixel);
 255       SetMagickPixelPacket(reconstruct_image,q,reconstruct_indexes+x,
 256         &reconstruct_pixel);
 257       difference=MagickFalse;
 258       if (channel == AllChannels)
 259         {
 260           if (IsMagickColorSimilar(&pixel,&reconstruct_pixel) == MagickFalse)
 261             difference=MagickTrue;
 262         }
 263       else
 264         {
 265           if (((channel & RedChannel) != 0) && (p->red != q->red))
 266             difference=MagickTrue;
 267           if (((channel & GreenChannel) != 0) && (p->green != q->green))
 268             difference=MagickTrue;
 269           if (((channel & BlueChannel) != 0) && (p->blue != q->blue))
 270             difference=MagickTrue;
 271           if (((channel & OpacityChannel) != 0) &&
 272               (image->matte != MagickFalse) && (p->opacity != q->opacity))
 273             difference=MagickTrue;
 274           if ((((channel & IndexChannel) != 0) &&
 275                (image->colorspace == CMYKColorspace) &&
 276                (reconstruct_image->colorspace == CMYKColorspace)) &&
 277               (indexes[x] != reconstruct_indexes[x]))
 278             difference=MagickTrue;
 279         }
 280       if (difference != MagickFalse)
 281         SetPixelPacket(highlight_image,&highlight,r,highlight_indexes+x);
 282       else
 283         SetPixelPacket(highlight_image,&lowlight,r,highlight_indexes+x);
 284       p++;
 285       q++;
 286       r++;
 287     }
 288     sync=SyncCacheViewAuthenticPixels(highlight_view,exception);
 289     if (sync == MagickFalse)
 290       status=MagickFalse;
 291   }
 292   highlight_view=DestroyCacheView(highlight_view);
 293   reconstruct_view=DestroyCacheView(reconstruct_view);
 294   image_view=DestroyCacheView(image_view);
 295   (void) CompositeImage(difference_image,image->compose,highlight_image,0,0);
 296   highlight_image=DestroyImage(highlight_image);
 297   if (status == MagickFalse)
 298     difference_image=DestroyImage(difference_image);
 299   return(difference_image);
 300 }
 301 \f
 302 /*
 303 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 304 %                                                                             %
 305 %                                                                             %
 306 %                                                                             %
 307 %   G e t I m a g e C h a n n e l D i s t o r t i o n                         %
 308 %                                                                             %
 309 %                                                                             %
 310 %                                                                             %
 311 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 312 %
 313 %  GetImageChannelDistortion() compares one or more image channels of an image
 314 %  to a reconstructed image and returns the specified distortion metric.
 315 %
 316 %  The format of the CompareImageChannels method is:
 317 %
 318 %      MagickBooleanType GetImageChannelDistortion(const Image *image,
 319 %        const Image *reconstruct_image,const ChannelType channel,
 320 %        const MetricType metric,double *distortion,ExceptionInfo *exception)
 321 %
 322 %  A description of each parameter follows:
 323 %
 324 %    o image: the image.
 325 %
 326 %    o reconstruct_image: the reconstruct image.
 327 %
 328 %    o channel: the channel.
 329 %
 330 %    o metric: the metric.
 331 %
 332 %    o distortion: the computed distortion between the images.
 333 %
 334 %    o exception: return any errors or warnings in this structure.
 335 %
 336 */
 337
 338 MagickExport MagickBooleanType GetImageDistortion(Image *image,
 339   const Image *reconstruct_image,const MetricType metric,double *distortion,
 340   ExceptionInfo *exception)
 341 {
 342   MagickBooleanType
 343     status;
 344
 345   status=GetImageChannelDistortion(image,reconstruct_image,AllChannels,
 346     metric,distortion,exception);
 347   return(status);
 348 }
 349
 350 static MagickBooleanType GetAbsoluteDistortion(const Image *image,
 351   const Image *reconstruct_image,const ChannelType channel,double *distortion,
 352   ExceptionInfo *exception)
 353 {
 354   CacheView
 355     *image_view,
 356     *reconstruct_view;
 357
 358   ssize_t
 359     y;
 360
 361   MagickBooleanType
 362     status;
 363
 364   MagickPixelPacket
 365     zero;
 366
 367   /*
 368     Compute the absolute difference in pixels between two images.
 369   */
 370   status=MagickTrue;
 371   GetMagickPixelPacket(image,&zero);
 372   image_view=AcquireCacheView(image);
 373   reconstruct_view=AcquireCacheView(reconstruct_image);
 374 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 375   #pragma omp parallel for schedule(dynamic,4) shared(status)
 376 #endif
 377   for (y=0; y < (ssize_t) image->rows; y++)
 378   {
 379     double
 380       channel_distortion[AllChannels+1];
 381
 382     MagickPixelPacket
 383       pixel,
 384       reconstruct_pixel;
 385
 386     register const IndexPacket
 387       *restrict indexes,
 388       *restrict reconstruct_indexes;
 389
 390     register const PixelPacket
 391       *restrict p,
 392       *restrict q;
 393
 394     register ssize_t
 395       i,
 396       x;
 397
 398     if (status == MagickFalse)
 399       continue;
 400     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
 401     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
 402       1,exception);
 403     if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
 404       {
 405         status=MagickFalse;
 406         continue;
 407       }
 408     indexes=GetCacheViewVirtualIndexQueue(image_view);
 409     reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
 410     pixel=zero;
 411     reconstruct_pixel=pixel;
 412     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
 413     for (x=0; x < (ssize_t) image->columns; x++)
 414     {
 415       SetMagickPixelPacket(image,p,indexes+x,&pixel);
 416       SetMagickPixelPacket(reconstruct_image,q,reconstruct_indexes+x,
 417         &reconstruct_pixel);
 418       if (IsMagickColorSimilar(&pixel,&reconstruct_pixel) == MagickFalse)
 419         {
 420           if ((channel & RedChannel) != 0)
 421             channel_distortion[RedChannel]++;
 422           if ((channel & GreenChannel) != 0)
 423             channel_distortion[GreenChannel]++;
 424           if ((channel & BlueChannel) != 0)
 425             channel_distortion[BlueChannel]++;
 426           if (((channel & OpacityChannel) != 0) &&
 427               (image->matte != MagickFalse))
 428             channel_distortion[OpacityChannel]++;
 429           if (((channel & IndexChannel) != 0) &&
 430               (image->colorspace == CMYKColorspace))
 431             channel_distortion[BlackChannel]++;
 432           channel_distortion[AllChannels]++;
 433         }
 434       p++;
 435       q++;
 436     }
 437 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 438   #pragma omp critical (MagickCore_GetAbsoluteError)
 439 #endif
 440     for (i=0; i <= (ssize_t) AllChannels; i++)
 441       distortion[i]+=channel_distortion[i];
 442   }
 443   reconstruct_view=DestroyCacheView(reconstruct_view);
 444   image_view=DestroyCacheView(image_view);
 445   return(status);
 446 }
 447
 448 static size_t GetNumberChannels(const Image *image,
 449   const ChannelType channel)
 450 {
 451   size_t
 452     channels;
 453
 454   channels=0;
 455   if ((channel & RedChannel) != 0)
 456     channels++;
 457   if ((channel & GreenChannel) != 0)
 458     channels++;
 459   if ((channel & BlueChannel) != 0)
 460     channels++;
 461   if (((channel & OpacityChannel) != 0) &&
 462        (image->matte != MagickFalse))
 463     channels++;
 464   if (((channel & IndexChannel) != 0) &&
 465       (image->colorspace == CMYKColorspace))
 466     channels++;
 467   return(channels);
 468 }
 469
 470 static MagickBooleanType 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   MagickBooleanType
 831     status;
 832
 833   MagickOffsetType
 834     progress;
 835
 836   MagickRealType
 837     area;
 838
 839   register ssize_t
 840     i;
 841
 842   ssize_t
 843     y;
 844
 845   /*
 846     Normalize to account for variation due to lighting and exposure condition.
 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       gamma;
 930
 931     gamma=image_statistics[i].standard_deviation*
 932       reconstruct_statistics[i].standard_deviation;
 933     gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
 934     distortion[i]=QuantumRange*gamma*distortion[i];
 935   }
 936   distortion[AllChannels]=0.0;
 937   if ((channel & RedChannel) != 0)
 938     distortion[AllChannels]+=distortion[RedChannel]*distortion[RedChannel];
 939   if ((channel & GreenChannel) != 0)
 940     distortion[AllChannels]+=distortion[GreenChannel]*distortion[GreenChannel];
 941   if ((channel & BlueChannel) != 0)
 942     distortion[AllChannels]+=distortion[BlueChannel]*distortion[BlueChannel];
 943   if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
 944     distortion[AllChannels]+=distortion[OpacityChannel]*
 945       distortion[OpacityChannel];
 946   if (((channel & IndexChannel) != 0) &&
 947       (image->colorspace == CMYKColorspace))
 948     distortion[AllChannels]+=distortion[BlackChannel]*distortion[BlackChannel];
 949   distortion[AllChannels]=sqrt(distortion[AllChannels]/GetNumberChannels(image,
 950     channel));
 951   /*
 952     Free resources.
 953   */
 954   reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
 955     reconstruct_statistics);
 956   image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
 957     image_statistics);
 958   return(status);
 959 }
 960
 961 static MagickBooleanType GetPeakAbsoluteDistortion(const Image *image,
 962   const Image *reconstruct_image,const ChannelType channel,
 963   double *distortion,ExceptionInfo *exception)
 964 {
 965   CacheView
 966     *image_view,
 967     *reconstruct_view;
 968
 969   ssize_t
 970     y;
 971
 972   MagickBooleanType
 973     status;
 974
 975   status=MagickTrue;
 976   image_view=AcquireCacheView(image);
 977   reconstruct_view=AcquireCacheView(reconstruct_image);
 978 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 979   #pragma omp parallel for schedule(dynamic,4) shared(status)
 980 #endif
 981   for (y=0; y < (ssize_t) image->rows; y++)
 982   {
 983     double
 984       channel_distortion[AllChannels+1];
 985
 986     register const IndexPacket
 987       *restrict indexes,
 988       *restrict reconstruct_indexes;
 989
 990     register const PixelPacket
 991       *restrict p,
 992       *restrict q;
 993
 994     register ssize_t
 995       i,
 996       x;
 997
 998     if (status == MagickFalse)
 999       continue;
1000     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1001     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,
1002       reconstruct_image->columns,1,exception);
1003     if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1004       {
1005         status=MagickFalse;
1006         continue;
1007       }
1008     indexes=GetCacheViewVirtualIndexQueue(image_view);
1009     reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
1010     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
1011     for (x=0; x < (ssize_t) image->columns; x++)
1012     {
1013       MagickRealType
1014         distance;
1015
1016       if ((channel & RedChannel) != 0)
1017         {
1018           distance=QuantumScale*fabs(p->red-(double) q->red);
1019           if (distance > channel_distortion[RedChannel])
1020             channel_distortion[RedChannel]=distance;
1021           if (distance > channel_distortion[AllChannels])
1022             channel_distortion[AllChannels]=distance;
1023         }
1024       if ((channel & GreenChannel) != 0)
1025         {
1026           distance=QuantumScale*fabs(p->green-(double) q->green);
1027           if (distance > channel_distortion[GreenChannel])
1028             channel_distortion[GreenChannel]=distance;
1029           if (distance > channel_distortion[AllChannels])
1030             channel_distortion[AllChannels]=distance;
1031         }
1032       if ((channel & BlueChannel) != 0)
1033         {
1034           distance=QuantumScale*fabs(p->blue-(double) q->blue);
1035           if (distance > channel_distortion[BlueChannel])
1036             channel_distortion[BlueChannel]=distance;
1037           if (distance > channel_distortion[AllChannels])
1038             channel_distortion[AllChannels]=distance;
1039         }
1040       if (((channel & OpacityChannel) != 0) &&
1041           (image->matte != MagickFalse))
1042         {
1043           distance=QuantumScale*fabs(p->opacity-(double) q->opacity);
1044           if (distance > channel_distortion[OpacityChannel])
1045             channel_distortion[OpacityChannel]=distance;
1046           if (distance > channel_distortion[AllChannels])
1047             channel_distortion[AllChannels]=distance;
1048         }
1049       if (((channel & IndexChannel) != 0) &&
1050           (image->colorspace == CMYKColorspace) &&
1051           (reconstruct_image->colorspace == CMYKColorspace))
1052         {
1053           distance=QuantumScale*fabs(indexes[x]-(double)
1054             reconstruct_indexes[x]);
1055           if (distance > channel_distortion[BlackChannel])
1056             channel_distortion[BlackChannel]=distance;
1057           if (distance > channel_distortion[AllChannels])
1058             channel_distortion[AllChannels]=distance;
1059         }
1060       p++;
1061       q++;
1062     }
1063 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1064   #pragma omp critical (MagickCore_GetPeakAbsoluteError)
1065 #endif
1066     for (i=0; i <= (ssize_t) AllChannels; i++)
1067       if (channel_distortion[i] > distortion[i])
1068         distortion[i]=channel_distortion[i];
1069   }
1070   reconstruct_view=DestroyCacheView(reconstruct_view);
1071   image_view=DestroyCacheView(image_view);
1072   return(status);
1073 }
1074
1075 static MagickBooleanType GetPeakSignalToNoiseRatio(const Image *image,
1076   const Image *reconstruct_image,const ChannelType channel,
1077   double *distortion,ExceptionInfo *exception)
1078 {
1079   MagickBooleanType
1080     status;
1081
1082   status=GetMeanSquaredDistortion(image,reconstruct_image,channel,distortion,
1083     exception);
1084   if ((channel & RedChannel) != 0)
1085     distortion[RedChannel]=20.0*log10((double) 1.0/sqrt(
1086       distortion[RedChannel]));
1087   if ((channel & GreenChannel) != 0)
1088     distortion[GreenChannel]=20.0*log10((double) 1.0/sqrt(
1089       distortion[GreenChannel]));
1090   if ((channel & BlueChannel) != 0)
1091     distortion[BlueChannel]=20.0*log10((double) 1.0/sqrt(
1092       distortion[BlueChannel]));
1093   if (((channel & OpacityChannel) != 0) &&
1094       (image->matte != MagickFalse))
1095     distortion[OpacityChannel]=20.0*log10((double) 1.0/sqrt(
1096       distortion[OpacityChannel]));
1097   if (((channel & IndexChannel) != 0) &&
1098       (image->colorspace == CMYKColorspace))
1099     distortion[BlackChannel]=20.0*log10((double) 1.0/sqrt(
1100       distortion[BlackChannel]));
1101   distortion[AllChannels]=20.0*log10((double) 1.0/sqrt(
1102     distortion[AllChannels]));
1103   return(status);
1104 }
1105
1106 static MagickBooleanType GetRootMeanSquaredDistortion(const Image *image,
1107   const Image *reconstruct_image,const ChannelType channel,
1108   double *distortion,ExceptionInfo *exception)
1109 {
1110   MagickBooleanType
1111     status;
1112
1113   status=GetMeanSquaredDistortion(image,reconstruct_image,channel,distortion,
1114     exception);
1115   if ((channel & RedChannel) != 0)
1116     distortion[RedChannel]=sqrt(distortion[RedChannel]);
1117   if ((channel & GreenChannel) != 0)
1118     distortion[GreenChannel]=sqrt(distortion[GreenChannel]);
1119   if ((channel & BlueChannel) != 0)
1120     distortion[BlueChannel]=sqrt(distortion[BlueChannel]);
1121   if (((channel & OpacityChannel) != 0) &&
1122       (image->matte != MagickFalse))
1123     distortion[OpacityChannel]=sqrt(distortion[OpacityChannel]);
1124   if (((channel & IndexChannel) != 0) &&
1125       (image->colorspace == CMYKColorspace))
1126     distortion[BlackChannel]=sqrt(distortion[BlackChannel]);
1127   distortion[AllChannels]=sqrt(distortion[AllChannels]);
1128   return(status);
1129 }
1130
1131 MagickExport MagickBooleanType GetImageChannelDistortion(Image *image,
1132   const Image *reconstruct_image,const ChannelType channel,
1133   const MetricType metric,double *distortion,ExceptionInfo *exception)
1134 {
1135   double
1136     *channel_distortion;
1137
1138   MagickBooleanType
1139     status;
1140
1141   size_t
1142     length;
1143
1144   assert(image != (Image *) NULL);
1145   assert(image->signature == MagickSignature);
1146   if (image->debug != MagickFalse)
1147     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1148   assert(reconstruct_image != (const Image *) NULL);
1149   assert(reconstruct_image->signature == MagickSignature);
1150   assert(distortion != (double *) NULL);
1151   *distortion=0.0;
1152   if (image->debug != MagickFalse)
1153     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1154   if ((reconstruct_image->columns != image->columns) ||
1155       (reconstruct_image->rows != image->rows))
1156     ThrowBinaryException(ImageError,"ImageSizeDiffers",image->filename);
1157   /*
1158     Get image distortion.
1159   */
1160   length=AllChannels+1UL;
1161   channel_distortion=(double *) AcquireQuantumMemory(length,
1162     sizeof(*channel_distortion));
1163   if (channel_distortion == (double *) NULL)
1164     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1165   (void) ResetMagickMemory(channel_distortion,0,length*
1166     sizeof(*channel_distortion));
1167   switch (metric)
1168   {
1169     case AbsoluteErrorMetric:
1170     {
1171       status=GetAbsoluteDistortion(image,reconstruct_image,channel,
1172         channel_distortion,exception);
1173       break;
1174     }
1175     case MeanAbsoluteErrorMetric:
1176     {
1177       status=GetMeanAbsoluteDistortion(image,reconstruct_image,channel,
1178         channel_distortion,exception);
1179       break;
1180     }
1181     case MeanErrorPerPixelMetric:
1182     {
1183       status=GetMeanErrorPerPixel(image,reconstruct_image,channel,
1184         channel_distortion,exception);
1185       break;
1186     }
1187     case MeanSquaredErrorMetric:
1188     {
1189       status=GetMeanSquaredDistortion(image,reconstruct_image,channel,
1190         channel_distortion,exception);
1191       break;
1192     }
1193     case NormalizedCrossCorrelationErrorMetric:
1194     default:
1195     {
1196       status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1197         channel,channel_distortion,exception);
1198       break;
1199     }
1200     case PeakAbsoluteErrorMetric:
1201     {
1202       status=GetPeakAbsoluteDistortion(image,reconstruct_image,channel,
1203         channel_distortion,exception);
1204       break;
1205     }
1206     case PeakSignalToNoiseRatioMetric:
1207     {
1208       status=GetPeakSignalToNoiseRatio(image,reconstruct_image,channel,
1209         channel_distortion,exception);
1210       break;
1211     }
1212     case RootMeanSquaredErrorMetric:
1213     {
1214       status=GetRootMeanSquaredDistortion(image,reconstruct_image,channel,
1215         channel_distortion,exception);
1216       break;
1217     }
1218   }
1219   *distortion=channel_distortion[AllChannels];
1220   channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1221   return(status);
1222 }
1223 \f
1224 /*
1225 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1226 %                                                                             %
1227 %                                                                             %
1228 %                                                                             %
1229 %   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                       %
1230 %                                                                             %
1231 %                                                                             %
1232 %                                                                             %
1233 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1234 %
1235 %  GetImageChannelDistrortion() compares the image channels of an image to a
1236 %  reconstructed image and returns the specified distortion metric for each
1237 %  channel.
1238 %
1239 %  The format of the CompareImageChannels method is:
1240 %
1241 %      double *GetImageChannelDistortions(const Image *image,
1242 %        const Image *reconstruct_image,const MetricType metric,
1243 %        ExceptionInfo *exception)
1244 %
1245 %  A description of each parameter follows:
1246 %
1247 %    o image: the image.
1248 %
1249 %    o reconstruct_image: the reconstruct image.
1250 %
1251 %    o metric: the metric.
1252 %
1253 %    o exception: return any errors or warnings in this structure.
1254 %
1255 */
1256 MagickExport double *GetImageChannelDistortions(Image *image,
1257   const Image *reconstruct_image,const MetricType metric,
1258   ExceptionInfo *exception)
1259 {
1260   double
1261     *channel_distortion;
1262
1263   MagickBooleanType
1264     status;
1265
1266   size_t
1267     length;
1268
1269   assert(image != (Image *) NULL);
1270   assert(image->signature == MagickSignature);
1271   if (image->debug != MagickFalse)
1272     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1273   assert(reconstruct_image != (const Image *) NULL);
1274   assert(reconstruct_image->signature == MagickSignature);
1275   if (image->debug != MagickFalse)
1276     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1277   if ((reconstruct_image->columns != image->columns) ||
1278       (reconstruct_image->rows != image->rows))
1279     {
1280       (void) ThrowMagickException(&image->exception,GetMagickModule(),
1281         ImageError,"ImageSizeDiffers","`%s'",image->filename);
1282       return((double *) NULL);
1283     }
1284   /*
1285     Get image distortion.
1286   */
1287   length=AllChannels+1UL;
1288   channel_distortion=(double *) AcquireQuantumMemory(length,
1289     sizeof(*channel_distortion));
1290   if (channel_distortion == (double *) NULL)
1291     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1292   (void) ResetMagickMemory(channel_distortion,0,length*
1293     sizeof(*channel_distortion));
1294   switch (metric)
1295   {
1296     case AbsoluteErrorMetric:
1297     {
1298       status=GetAbsoluteDistortion(image,reconstruct_image,AllChannels,
1299         channel_distortion,exception);
1300       break;
1301     }
1302     case MeanAbsoluteErrorMetric:
1303     {
1304       status=GetMeanAbsoluteDistortion(image,reconstruct_image,AllChannels,
1305         channel_distortion,exception);
1306       break;
1307     }
1308     case MeanErrorPerPixelMetric:
1309     {
1310       status=GetMeanErrorPerPixel(image,reconstruct_image,AllChannels,
1311         channel_distortion,exception);
1312       break;
1313     }
1314     case MeanSquaredErrorMetric:
1315     {
1316       status=GetMeanSquaredDistortion(image,reconstruct_image,AllChannels,
1317         channel_distortion,exception);
1318       break;
1319     }
1320     case NormalizedCrossCorrelationErrorMetric:
1321     default:
1322     {
1323       status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1324         AllChannels,channel_distortion,exception);
1325       break;
1326     }
1327     case PeakAbsoluteErrorMetric:
1328     {
1329       status=GetPeakAbsoluteDistortion(image,reconstruct_image,AllChannels,
1330         channel_distortion,exception);
1331       break;
1332     }
1333     case PeakSignalToNoiseRatioMetric:
1334     {
1335       status=GetPeakSignalToNoiseRatio(image,reconstruct_image,AllChannels,
1336         channel_distortion,exception);
1337       break;
1338     }
1339     case RootMeanSquaredErrorMetric:
1340     {
1341       status=GetRootMeanSquaredDistortion(image,reconstruct_image,AllChannels,
1342         channel_distortion,exception);
1343       break;
1344     }
1345   }
1346   return(channel_distortion);
1347 }
1348 \f
1349 /*
1350 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1351 %                                                                             %
1352 %                                                                             %
1353 %                                                                             %
1354 %  I s I m a g e s E q u a l                                                  %
1355 %                                                                             %
1356 %                                                                             %
1357 %                                                                             %
1358 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1359 %
1360 %  IsImagesEqual() measures the difference between colors at each pixel
1361 %  location of two images.  A value other than 0 means the colors match
1362 %  exactly.  Otherwise an error measure is computed by summing over all
1363 %  pixels in an image the distance squared in RGB space between each image
1364 %  pixel and its corresponding pixel in the reconstruct image.  The error
1365 %  measure is assigned to these image members:
1366 %
1367 %    o mean_error_per_pixel:  The mean error for any single pixel in
1368 %      the image.
1369 %
1370 %    o normalized_mean_error:  The normalized mean quantization error for
1371 %      any single pixel in the image.  This distance measure is normalized to
1372 %      a range between 0 and 1.  It is independent of the range of red, green,
1373 %      and blue values in the image.
1374 %
1375 %    o normalized_maximum_error:  The normalized maximum quantization
1376 %      error for any single pixel in the image.  This distance measure is
1377 %      normalized to a range between 0 and 1.  It is independent of the range
1378 %      of red, green, and blue values in your image.
1379 %
1380 %  A small normalized mean square error, accessed as
1381 %  image->normalized_mean_error, suggests the images are very similar in
1382 %  spatial layout and color.
1383 %
1384 %  The format of the IsImagesEqual method is:
1385 %
1386 %      MagickBooleanType IsImagesEqual(Image *image,
1387 %        const Image *reconstruct_image)
1388 %
1389 %  A description of each parameter follows.
1390 %
1391 %    o image: the image.
1392 %
1393 %    o reconstruct_image: the reconstruct image.
1394 %
1395 */
1396 MagickExport MagickBooleanType IsImagesEqual(Image *image,
1397   const Image *reconstruct_image)
1398 {
1399   CacheView
1400     *image_view,
1401     *reconstruct_view;
1402
1403   ExceptionInfo
1404     *exception;
1405
1406   ssize_t
1407     y;
1408
1409   MagickBooleanType
1410     status;
1411
1412   MagickRealType
1413     area,
1414     maximum_error,
1415     mean_error,
1416     mean_error_per_pixel;
1417
1418   assert(image != (Image *) NULL);
1419   assert(image->signature == MagickSignature);
1420   assert(reconstruct_image != (const Image *) NULL);
1421   assert(reconstruct_image->signature == MagickSignature);
1422   if ((reconstruct_image->columns != image->columns) ||
1423       (reconstruct_image->rows != image->rows))
1424     ThrowBinaryException(ImageError,"ImageSizeDiffers",image->filename);
1425   area=0.0;
1426   maximum_error=0.0;
1427   mean_error_per_pixel=0.0;
1428   mean_error=0.0;
1429   exception=(&image->exception);
1430   image_view=AcquireCacheView(image);
1431   reconstruct_view=AcquireCacheView(reconstruct_image);
1432   for (y=0; y < (ssize_t) image->rows; y++)
1433   {
1434     register const IndexPacket
1435       *restrict indexes,
1436       *restrict reconstruct_indexes;
1437
1438     register const PixelPacket
1439       *restrict p,
1440       *restrict q;
1441
1442     register ssize_t
1443       x;
1444
1445     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1446     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
1447       1,exception);
1448     if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1449       break;
1450     indexes=GetCacheViewVirtualIndexQueue(image_view);
1451     reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
1452     for (x=0; x < (ssize_t) image->columns; x++)
1453     {
1454       MagickRealType
1455         distance;
1456
1457       distance=fabs(p->red-(double) q->red);
1458       mean_error_per_pixel+=distance;
1459       mean_error+=distance*distance;
1460       if (distance > maximum_error)
1461         maximum_error=distance;
1462       area++;
1463       distance=fabs(p->green-(double) q->green);
1464       mean_error_per_pixel+=distance;
1465       mean_error+=distance*distance;
1466       if (distance > maximum_error)
1467         maximum_error=distance;
1468       area++;
1469       distance=fabs(p->blue-(double) q->blue);
1470       mean_error_per_pixel+=distance;
1471       mean_error+=distance*distance;
1472       if (distance > maximum_error)
1473         maximum_error=distance;
1474       area++;
1475       if (image->matte != MagickFalse)
1476         {
1477           distance=fabs(p->opacity-(double) q->opacity);
1478           mean_error_per_pixel+=distance;
1479           mean_error+=distance*distance;
1480           if (distance > maximum_error)
1481             maximum_error=distance;
1482           area++;
1483         }
1484       if ((image->colorspace == CMYKColorspace) &&
1485           (reconstruct_image->colorspace == CMYKColorspace))
1486         {
1487           distance=fabs(indexes[x]-(double) reconstruct_indexes[x]);
1488           mean_error_per_pixel+=distance;
1489           mean_error+=distance*distance;
1490           if (distance > maximum_error)
1491             maximum_error=distance;
1492           area++;
1493         }
1494       p++;
1495       q++;
1496     }
1497   }
1498   reconstruct_view=DestroyCacheView(reconstruct_view);
1499   image_view=DestroyCacheView(image_view);
1500   image->error.mean_error_per_pixel=(double) (mean_error_per_pixel/area);
1501   image->error.normalized_mean_error=(double) (QuantumScale*QuantumScale*
1502     mean_error/area);
1503   image->error.normalized_maximum_error=(double) (QuantumScale*maximum_error);
1504   status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
1505   return(status);
1506 }
1507 \f
1508 /*
1509 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1510 %                                                                             %
1511 %                                                                             %
1512 %                                                                             %
1513 %   S i m i l a r i t y I m a g e                                             %
1514 %                                                                             %
1515 %                                                                             %
1516 %                                                                             %
1517 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1518 %
1519 %  SimilarityImage() compares the reference image of the image and returns the
1520 %  best match offset.  In addition, it returns a similarity image such that an
1521 %  exact match location is completely white and if none of the pixels match,
1522 %  black, otherwise some gray level in-between.
1523 %
1524 %  The format of the SimilarityImageImage method is:
1525 %
1526 %      Image *SimilarityImage(const Image *image,const Image *reference,
1527 %        RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
1528 %
1529 %  A description of each parameter follows:
1530 %
1531 %    o image: the image.
1532 %
1533 %    o reference: find an area of the image that closely resembles this image.
1534 %
1535 %    o the best match offset of the reference image within the image.
1536 %
1537 %    o similarity: the computed similarity between the images.
1538 %
1539 %    o exception: return any errors or warnings in this structure.
1540 %
1541 */
1542
1543 static double GetNCCDistortion(const Image *image,
1544   const Image *reconstruct_image,
1545   const ChannelStatistics *reconstruct_statistics,ExceptionInfo *exception)
1546 {
1547 #define SimilarityImageTag  "Similarity/Image"
1548
1549   CacheView
1550     *image_view,
1551     *reconstruct_view;
1552
1553   ChannelStatistics
1554     *image_statistics;
1555
1556   double
1557     distortion;
1558
1559   MagickBooleanType
1560     status;
1561
1562   MagickRealType
1563     area,
1564     gamma;
1565
1566   ssize_t
1567     y;
1568
1569   unsigned long
1570     number_channels;
1571
1572   /*
1573     Normalize to account for variation due to lighting and exposure condition.
1574   */
1575   image_statistics=GetImageChannelStatistics(image,exception);
1576   status=MagickTrue;
1577   distortion=0.0;
1578   area=1.0/((MagickRealType) image->columns*image->rows);
1579   image_view=AcquireCacheView(image);
1580   reconstruct_view=AcquireCacheView(reconstruct_image);
1581   for (y=0; y < (ssize_t) image->rows; y++)
1582   {
1583     register const IndexPacket
1584       *restrict indexes,
1585       *restrict reconstruct_indexes;
1586
1587     register const PixelPacket
1588       *restrict p,
1589       *restrict q;
1590
1591     register ssize_t
1592       x;
1593
1594     if (status == MagickFalse)
1595       continue;
1596     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1597     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
1598       1,exception);
1599     if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1600       {
1601         status=MagickFalse;
1602         continue;
1603       }
1604     indexes=GetCacheViewVirtualIndexQueue(image_view);
1605     reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
1606     for (x=0; x < (ssize_t) image->columns; x++)
1607     {
1608       distortion+=area*QuantumScale*(p->red-
1609         image_statistics[RedChannel].mean)*(q->red-
1610         reconstruct_statistics[RedChannel].mean);
1611       distortion+=area*QuantumScale*(p->green-
1612         image_statistics[GreenChannel].mean)*(q->green-
1613         reconstruct_statistics[GreenChannel].mean);
1614       distortion+=area*QuantumScale*(p->blue-
1615         image_statistics[BlueChannel].mean)*(q->blue-
1616         reconstruct_statistics[BlueChannel].mean);
1617       if (image->matte != MagickFalse)
1618         distortion+=area*QuantumScale*(p->opacity-
1619           image_statistics[OpacityChannel].mean)*(q->opacity-
1620           reconstruct_statistics[OpacityChannel].mean);
1621       if ((image->colorspace == CMYKColorspace) &&
1622           (reconstruct_image->colorspace == CMYKColorspace))
1623         distortion+=area*QuantumScale*(indexes[x]-
1624           image_statistics[OpacityChannel].mean)*(reconstruct_indexes[x]-
1625           reconstruct_statistics[OpacityChannel].mean);
1626       p++;
1627       q++;
1628     }
1629   }
1630   reconstruct_view=DestroyCacheView(reconstruct_view);
1631   image_view=DestroyCacheView(image_view);
1632   /*
1633     Divide by the standard deviation.
1634   */
1635   gamma=image_statistics[AllChannels].standard_deviation*
1636     reconstruct_statistics[AllChannels].standard_deviation;
1637   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1638   distortion=QuantumRange*gamma*distortion;
1639   number_channels=3;
1640   if (image->matte != MagickFalse)
1641     number_channels++;
1642   if (image->colorspace == CMYKColorspace)
1643     number_channels++;
1644   distortion=sqrt(distortion/number_channels);
1645   /*
1646     Free resources.
1647   */
1648   image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1649     image_statistics);
1650   return(1.0-distortion);
1651 }
1652
1653 static double GetSimilarityMetric(const Image *image,const Image *reference,
1654   const ChannelStatistics *reference_statistics,const ssize_t x_offset,
1655   const ssize_t y_offset,ExceptionInfo *exception)
1656 {
1657   double
1658     distortion;
1659
1660   Image
1661     *similarity_image;
1662
1663   RectangleInfo
1664     geometry;
1665
1666   SetGeometry(reference,&geometry);
1667   geometry.x=x_offset;
1668   geometry.y=y_offset;
1669   similarity_image=CropImage(image,&geometry,exception);
1670   if (similarity_image == (Image *) NULL)
1671     return(0.0);
1672   distortion=GetNCCDistortion(reference,similarity_image,reference_statistics,
1673     exception);
1674   similarity_image=DestroyImage(similarity_image);
1675   return(distortion);
1676 }
1677
1678 MagickExport Image *SimilarityImage(Image *image,const Image *reference,
1679   RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
1680 {
1681 #define SimilarityImageTag  "Similarity/Image"
1682
1683   CacheView
1684     *similarity_view;
1685
1686   ChannelStatistics
1687     *reference_statistics;
1688
1689   Image
1690     *similarity_image;
1691
1692   MagickBooleanType
1693     status;
1694
1695   MagickOffsetType
1696     progress;
1697
1698   ssize_t
1699     y;
1700
1701   assert(image != (const Image *) NULL);
1702   assert(image->signature == MagickSignature);
1703   if (image->debug != MagickFalse)
1704     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1705   assert(exception != (ExceptionInfo *) NULL);
1706   assert(exception->signature == MagickSignature);
1707   assert(offset != (RectangleInfo *) NULL);
1708   SetGeometry(reference,offset);
1709   *similarity_metric=1.0;
1710   if ((reference->columns > image->columns) || (reference->rows > image->rows))
1711     ThrowImageException(ImageError,"ImageSizeDiffers");
1712   similarity_image=CloneImage(image,image->columns-reference->columns+1,
1713     image->rows-reference->rows+1,MagickTrue,exception);
1714   if (similarity_image == (Image *) NULL)
1715     return((Image *) NULL);
1716   if (SetImageStorageClass(similarity_image,DirectClass) == MagickFalse)
1717     {
1718       InheritException(exception,&similarity_image->exception);
1719       similarity_image=DestroyImage(similarity_image);
1720       return((Image *) NULL);
1721     }
1722   /*
1723     Measure similarity of reference image against image.
1724   */
1725   status=MagickTrue;
1726   progress=0;
1727   reference_statistics=GetImageChannelStatistics(reference,exception);
1728   similarity_view=AcquireCacheView(similarity_image);
1729 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1730   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1731 #endif
1732   for (y=0; y < (ssize_t) (image->rows-reference->rows+1); y++)
1733   {
1734     double
1735       similarity;
1736
1737     register ssize_t
1738       x;
1739
1740     register PixelPacket
1741       *restrict q;
1742
1743     if (status == MagickFalse)
1744       continue;
1745     q=GetCacheViewAuthenticPixels(similarity_view,0,y,similarity_image->columns,
1746       1,exception);
1747     if (q == (const PixelPacket *) NULL)
1748       {
1749         status=MagickFalse;
1750         continue;
1751       }
1752     for (x=0; x < (ssize_t) (image->columns-reference->columns+1); x++)
1753     {
1754       similarity=GetSimilarityMetric(image,reference,reference_statistics,x,y,
1755         exception);
1756 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1757   #pragma omp critical (MagickCore_SimilarityImage)
1758 #endif
1759       if (similarity < *similarity_metric)
1760         {
1761           *similarity_metric=similarity;
1762           offset->x=x;
1763           offset->y=y;
1764         }
1765       q->red=ClampToQuantum(QuantumRange-QuantumRange*similarity);
1766       q->green=q->red;
1767       q->blue=q->red;
1768       q++;
1769     }
1770     if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
1771       status=MagickFalse;
1772     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1773       {
1774         MagickBooleanType
1775           proceed;
1776
1777 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1778   #pragma omp critical (MagickCore_SimilarityImage)
1779 #endif
1780         proceed=SetImageProgress(image,SimilarityImageTag,progress++,
1781           image->rows);
1782         if (proceed == MagickFalse)
1783           status=MagickFalse;
1784       }
1785   }
1786   similarity_view=DestroyCacheView(similarity_view);
1787   reference_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1788     reference_statistics);
1789   return(similarity_image);
1790 }