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