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