1 /* 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % % 4 % % 5 % % 6 % SSSSS TTTTT AAA TTTTT IIIII SSSSS TTTTT IIIII CCCC % 7 % SS T A A T I SS T I C % 8 % SSS T AAAAA T I SSS T I C % 9 % SS T A A T I SS T I C % 10 % SSSSS T A A T IIIII SSSSS T IIIII CCCC % 11 % % 12 % % 13 % MagickCore Image Statistical Methods % 14 % % 15 % Software Design % 16 % Cristy % 17 % July 1992 % 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/property.h" 46 #include "MagickCore/accelerate-private.h" 47 #include "MagickCore/animate.h" 48 #include "MagickCore/blob.h" 49 #include "MagickCore/blob-private.h" 50 #include "MagickCore/cache.h" 51 #include "MagickCore/cache-private.h" 52 #include "MagickCore/cache-view.h" 53 #include "MagickCore/client.h" 54 #include "MagickCore/color.h" 55 #include "MagickCore/color-private.h" 56 #include "MagickCore/colorspace.h" 57 #include "MagickCore/colorspace-private.h" 58 #include "MagickCore/composite.h" 59 #include "MagickCore/composite-private.h" 60 #include "MagickCore/compress.h" 61 #include "MagickCore/constitute.h" 62 #include "MagickCore/display.h" 63 #include "MagickCore/draw.h" 64 #include "MagickCore/enhance.h" 65 #include "MagickCore/exception.h" 66 #include "MagickCore/exception-private.h" 67 #include "MagickCore/gem.h" 68 #include "MagickCore/gem-private.h" 69 #include "MagickCore/geometry.h" 70 #include "MagickCore/list.h" 71 #include "MagickCore/image-private.h" 72 #include "MagickCore/magic.h" 73 #include "MagickCore/magick.h" 74 #include "MagickCore/memory_.h" 75 #include "MagickCore/module.h" 76 #include "MagickCore/monitor.h" 77 #include "MagickCore/monitor-private.h" 78 #include "MagickCore/option.h" 79 #include "MagickCore/paint.h" 80 #include "MagickCore/pixel-accessor.h" 81 #include "MagickCore/profile.h" 82 #include "MagickCore/quantize.h" 83 #include "MagickCore/quantum-private.h" 84 #include "MagickCore/random_.h" 85 #include "MagickCore/random-private.h" 86 #include "MagickCore/resource_.h" 87 #include "MagickCore/segment.h" 88 #include "MagickCore/semaphore.h" 89 #include "MagickCore/signature-private.h" 90 #include "MagickCore/statistic.h" 91 #include "MagickCore/string_.h" 92 #include "MagickCore/thread-private.h" 93 #include "MagickCore/timer.h" 94 #include "MagickCore/utility.h" 95 #include "MagickCore/version.h" 96 97 /* 99 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 100 % % 101 % % 102 % % 103 % E v a l u a t e I m a g e % 104 % % 105 % % 106 % % 107 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 108 % 109 % EvaluateImage() applies a value to the image with an arithmetic, relational, 110 % or logical operator to an image. Use these operations to lighten or darken 111 % an image, to increase or decrease contrast in an image, or to produce the 112 % "negative" of an image. 113 % 114 % The format of the EvaluateImage method is: 115 % 116 % MagickBooleanType EvaluateImage(Image *image, 117 % const MagickEvaluateOperator op,const double value, 118 % ExceptionInfo *exception) 119 % MagickBooleanType EvaluateImages(Image *images, 120 % const MagickEvaluateOperator op,const double value, 121 % ExceptionInfo *exception) 122 % 123 % A description of each parameter follows: 124 % 125 % o image: the image. 126 % 127 % o op: A channel op. 128 % 129 % o value: A value value. 130 % 131 % o exception: return any errors or warnings in this structure. 132 % 133 */ 134 135 typedef struct _PixelChannels 136 { 137 double 138 channel[CompositePixelChannel]; 139 } PixelChannels; 140 141 static PixelChannels **DestroyPixelThreadSet(PixelChannels **pixels) 142 { 143 register ssize_t 144 i; 145 146 assert(pixels != (PixelChannels **) NULL); 147 for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++) 148 if (pixels[i] != (PixelChannels *) NULL) 149 pixels[i]=(PixelChannels *) RelinquishMagickMemory(pixels[i]); 150 pixels=(PixelChannels **) RelinquishMagickMemory(pixels); 151 return(pixels); 152 } 153 154 static PixelChannels **AcquirePixelThreadSet(const Image *image) 155 { 156 PixelChannels 157 **pixels; 158 159 register ssize_t 160 i; 161 162 size_t 163 number_threads; 164 165 number_threads=(size_t) GetMagickResourceLimit(ThreadResource); 166 pixels=(PixelChannels **) AcquireQuantumMemory(number_threads, 167 sizeof(*pixels)); 168 if (pixels == (PixelChannels **) NULL) 169 return((PixelChannels **) NULL); 170 (void) ResetMagickMemory(pixels,0,number_threads*sizeof(*pixels)); 171 for (i=0; i < (ssize_t) number_threads; i++) 172 { 173 register ssize_t 174 j; 175 176 pixels[i]=(PixelChannels *) AcquireQuantumMemory(image->columns, 177 sizeof(**pixels)); 178 if (pixels[i] == (PixelChannels *) NULL) 179 return(DestroyPixelThreadSet(pixels)); 180 for (j=0; j < (ssize_t) image->columns; j++) 181 { 182 register ssize_t 183 k; 184 185 for (k=0; k < MaxPixelChannels; k++) 186 pixels[i][j].channel[k]=0.0; 187 } 188 } 189 return(pixels); 190 } 191 192 static inline double EvaluateMax(const double x,const double y) 193 { 194 if (x > y) 195 return(x); 196 return(y); 197 } 198 199 #if defined(__cplusplus) || defined(c_plusplus) 200 extern "C" { 201 #endif 202 203 static int IntensityCompare(const void *x,const void *y) 204 { 205 const PixelChannels 206 *color_1, 207 *color_2; 208 209 double 210 distance; 211 212 register ssize_t 213 i; 214 215 color_1=(const PixelChannels *) x; 216 color_2=(const PixelChannels *) y; 217 distance=0.0; 218 for (i=0; i < MaxPixelChannels; i++) 219 distance+=color_1->channel[i]-(double) color_2->channel[i]; 220 return(distance < 0 ? -1 : distance > 0 ? 1 : 0); 221 } 222 223 #if defined(__cplusplus) || defined(c_plusplus) 224 } 225 #endif 226 227 static double ApplyEvaluateOperator(RandomInfo *random_info,const Quantum pixel, 228 const MagickEvaluateOperator op,const double value) 229 { 230 double 231 result; 232 233 result=0.0; 234 switch (op) 235 { 236 case UndefinedEvaluateOperator: 237 break; 238 case AbsEvaluateOperator: 239 { 240 result=(double) fabs((double) (pixel+value)); 241 break; 242 } 243 case AddEvaluateOperator: 244 { 245 result=(double) (pixel+value); 246 break; 247 } 248 case AddModulusEvaluateOperator: 249 { 250 /* 251 This returns a 'floored modulus' of the addition which is a positive 252 result. It differs from % or fmod() that returns a 'truncated modulus' 253 result, where floor() is replaced by trunc() and could return a 254 negative result (which is clipped). 255 */ 256 result=pixel+value; 257 result-=(QuantumRange+1.0)*floor((double) result/(QuantumRange+1.0)); 258 break; 259 } 260 case AndEvaluateOperator: 261 { 262 result=(double) ((size_t) pixel & (size_t) (value+0.5)); 263 break; 264 } 265 case CosineEvaluateOperator: 266 { 267 result=(double) (QuantumRange*(0.5*cos((double) (2.0*MagickPI* 268 QuantumScale*pixel*value))+0.5)); 269 break; 270 } 271 case DivideEvaluateOperator: 272 { 273 result=pixel/(value == 0.0 ? 1.0 : value); 274 break; 275 } 276 case ExponentialEvaluateOperator: 277 { 278 result=(double) (QuantumRange*exp((double) (value*QuantumScale*pixel))); 279 break; 280 } 281 case GaussianNoiseEvaluateOperator: 282 { 283 result=(double) GenerateDifferentialNoise(random_info,pixel, 284 GaussianNoise,value); 285 break; 286 } 287 case ImpulseNoiseEvaluateOperator: 288 { 289 result=(double) GenerateDifferentialNoise(random_info,pixel,ImpulseNoise, 290 value); 291 break; 292 } 293 case LaplacianNoiseEvaluateOperator: 294 { 295 result=(double) GenerateDifferentialNoise(random_info,pixel, 296 LaplacianNoise,value); 297 break; 298 } 299 case LeftShiftEvaluateOperator: 300 { 301 result=(double) ((size_t) pixel << (size_t) (value+0.5)); 302 break; 303 } 304 case LogEvaluateOperator: 305 { 306 if ((QuantumScale*pixel) >= MagickEpsilon) 307 result=(double) (QuantumRange*log((double) (QuantumScale*value*pixel+ 308 1.0))/log((double) (value+1.0))); 309 break; 310 } 311 case MaxEvaluateOperator: 312 { 313 result=(double) EvaluateMax((double) pixel,value); 314 break; 315 } 316 case MeanEvaluateOperator: 317 { 318 result=(double) (pixel+value); 319 break; 320 } 321 case MedianEvaluateOperator: 322 { 323 result=(double) (pixel+value); 324 break; 325 } 326 case MinEvaluateOperator: 327 { 328 result=(double) MagickMin((double) pixel,value); 329 break; 330 } 331 case MultiplicativeNoiseEvaluateOperator: 332 { 333 result=(double) GenerateDifferentialNoise(random_info,pixel, 334 MultiplicativeGaussianNoise,value); 335 break; 336 } 337 case MultiplyEvaluateOperator: 338 { 339 result=(double) (value*pixel); 340 break; 341 } 342 case OrEvaluateOperator: 343 { 344 result=(double) ((size_t) pixel | (size_t) (value+0.5)); 345 break; 346 } 347 case PoissonNoiseEvaluateOperator: 348 { 349 result=(double) GenerateDifferentialNoise(random_info,pixel,PoissonNoise, 350 value); 351 break; 352 } 353 case PowEvaluateOperator: 354 { 355 result=(double) (QuantumRange*pow((double) (QuantumScale*pixel),(double) 356 value)); 357 break; 358 } 359 case RightShiftEvaluateOperator: 360 { 361 result=(double) ((size_t) pixel >> (size_t) (value+0.5)); 362 break; 363 } 364 case RootMeanSquareEvaluateOperator: 365 { 366 result=(double) (pixel*pixel+value); 367 break; 368 } 369 case SetEvaluateOperator: 370 { 371 result=value; 372 break; 373 } 374 case SineEvaluateOperator: 375 { 376 result=(double) (QuantumRange*(0.5*sin((double) (2.0*MagickPI* 377 QuantumScale*pixel*value))+0.5)); 378 break; 379 } 380 case SubtractEvaluateOperator: 381 { 382 result=(double) (pixel-value); 383 break; 384 } 385 case SumEvaluateOperator: 386 { 387 result=(double) (pixel+value); 388 break; 389 } 390 case ThresholdEvaluateOperator: 391 { 392 result=(double) (((double) pixel <= value) ? 0 : QuantumRange); 393 break; 394 } 395 case ThresholdBlackEvaluateOperator: 396 { 397 result=(double) (((double) pixel <= value) ? 0 : pixel); 398 break; 399 } 400 case ThresholdWhiteEvaluateOperator: 401 { 402 result=(double) (((double) pixel > value) ? QuantumRange : pixel); 403 break; 404 } 405 case UniformNoiseEvaluateOperator: 406 { 407 result=(double) GenerateDifferentialNoise(random_info,pixel,UniformNoise, 408 value); 409 break; 410 } 411 case XorEvaluateOperator: 412 { 413 result=(double) ((size_t) pixel ^ (size_t) (value+0.5)); 414 break; 415 } 416 } 417 return(result); 418 } 419 420 MagickExport Image *EvaluateImages(const Image *images, 421 const MagickEvaluateOperator op,ExceptionInfo *exception) 422 { 423 #define EvaluateImageTag "Evaluate/Image" 424 425 CacheView 426 *evaluate_view; 427 428 Image 429 *image; 430 431 MagickBooleanType 432 status; 433 434 MagickOffsetType 435 progress; 436 437 PixelChannels 438 **magick_restrict evaluate_pixels; 439 440 RandomInfo 441 **magick_restrict random_info; 442 443 size_t 444 number_images; 445 446 ssize_t 447 y; 448 449 #if defined(MAGICKCORE_OPENMP_SUPPORT) 450 unsigned long 451 key; 452 #endif 453 454 assert(images != (Image *) NULL); 455 assert(images->signature == MagickCoreSignature); 456 if (images->debug != MagickFalse) 457 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename); 458 assert(exception != (ExceptionInfo *) NULL); 459 assert(exception->signature == MagickCoreSignature); 460 image=CloneImage(images,images->columns,images->rows,MagickTrue, 461 exception); 462 if (image == (Image *) NULL) 463 return((Image *) NULL); 464 if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse) 465 { 466 image=DestroyImage(image); 467 return((Image *) NULL); 468 } 469 number_images=GetImageListLength(images); 470 evaluate_pixels=AcquirePixelThreadSet(images); 471 if (evaluate_pixels == (PixelChannels **) NULL) 472 { 473 image=DestroyImage(image); 474 (void) ThrowMagickException(exception,GetMagickModule(), 475 ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename); 476 return((Image *) NULL); 477 } 478 /* 479 Evaluate image pixels. 480 */ 481 status=MagickTrue; 482 progress=0; 483 random_info=AcquireRandomInfoThreadSet(); 484 evaluate_view=AcquireAuthenticCacheView(image,exception); 485 if (op == MedianEvaluateOperator) 486 { 487 #if defined(MAGICKCORE_OPENMP_SUPPORT) 488 key=GetRandomSecretKey(random_info[0]); 489 #pragma omp parallel for schedule(static,4) shared(progress,status) \ 490 magick_threads(image,images,image->rows,key == ~0UL) 491 #endif 492 for (y=0; y < (ssize_t) image->rows; y++) 493 { 494 CacheView 495 *image_view; 496 497 const Image 498 *next; 499 500 const int 501 id = GetOpenMPThreadId(); 502 503 register PixelChannels 504 *evaluate_pixel; 505 506 register Quantum 507 *magick_restrict q; 508 509 register ssize_t 510 x; 511 512 if (status == MagickFalse) 513 continue; 514 q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,image->columns,1, 515 exception); 516 if (q == (Quantum *) NULL) 517 { 518 status=MagickFalse; 519 continue; 520 } 521 evaluate_pixel=evaluate_pixels[id]; 522 for (x=0; x < (ssize_t) image->columns; x++) 523 { 524 register ssize_t 525 j, 526 k; 527 528 for (j=0; j < (ssize_t) number_images; j++) 529 for (k=0; k < MaxPixelChannels; k++) 530 evaluate_pixel[j].channel[k]=0.0; 531 next=images; 532 for (j=0; j < (ssize_t) number_images; j++) 533 { 534 register const Quantum 535 *p; 536 537 register ssize_t 538 i; 539 540 image_view=AcquireVirtualCacheView(next,exception); 541 p=GetCacheViewVirtualPixels(image_view,x,y,1,1,exception); 542 if (p == (const Quantum *) NULL) 543 { 544 image_view=DestroyCacheView(image_view); 545 break; 546 } 547 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 548 { 549 PixelChannel channel=GetPixelChannelChannel(image,i); 550 PixelTrait evaluate_traits=GetPixelChannelTraits(image,channel); 551 PixelTrait traits=GetPixelChannelTraits(next,channel); 552 if ((traits == UndefinedPixelTrait) || 553 (evaluate_traits == UndefinedPixelTrait)) 554 continue; 555 if ((evaluate_traits & UpdatePixelTrait) == 0) 556 continue; 557 evaluate_pixel[j].channel[i]=ApplyEvaluateOperator( 558 random_info[id],GetPixelChannel(image,channel,p),op, 559 evaluate_pixel[j].channel[i]); 560 } 561 image_view=DestroyCacheView(image_view); 562 next=GetNextImageInList(next); 563 } 564 qsort((void *) evaluate_pixel,number_images,sizeof(*evaluate_pixel), 565 IntensityCompare); 566 for (k=0; k < (ssize_t) GetPixelChannels(image); k++) 567 q[k]=ClampToQuantum(evaluate_pixel[j/2].channel[k]); 568 q+=GetPixelChannels(image); 569 } 570 if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse) 571 status=MagickFalse; 572 if (images->progress_monitor != (MagickProgressMonitor) NULL) 573 { 574 MagickBooleanType 575 proceed; 576 577 #if defined(MAGICKCORE_OPENMP_SUPPORT) 578 #pragma omp critical (MagickCore_EvaluateImages) 579 #endif 580 proceed=SetImageProgress(images,EvaluateImageTag,progress++, 581 image->rows); 582 if (proceed == MagickFalse) 583 status=MagickFalse; 584 } 585 } 586 } 587 else 588 { 589 #if defined(MAGICKCORE_OPENMP_SUPPORT) 590 key=GetRandomSecretKey(random_info[0]); 591 #pragma omp parallel for schedule(static,4) shared(progress,status) \ 592 magick_threads(image,images,image->rows,key == ~0UL) 593 #endif 594 for (y=0; y < (ssize_t) image->rows; y++) 595 { 596 CacheView 597 *image_view; 598 599 const Image 600 *next; 601 602 const int 603 id = GetOpenMPThreadId(); 604 605 register ssize_t 606 i, 607 x; 608 609 register PixelChannels 610 *evaluate_pixel; 611 612 register Quantum 613 *magick_restrict q; 614 615 ssize_t 616 j; 617 618 if (status == MagickFalse) 619 continue; 620 q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,image->columns,1, 621 exception); 622 if (q == (Quantum *) NULL) 623 { 624 status=MagickFalse; 625 continue; 626 } 627 evaluate_pixel=evaluate_pixels[id]; 628 for (j=0; j < (ssize_t) image->columns; j++) 629 for (i=0; i < MaxPixelChannels; i++) 630 evaluate_pixel[j].channel[i]=0.0; 631 next=images; 632 for (j=0; j < (ssize_t) number_images; j++) 633 { 634 register const Quantum 635 *p; 636 637 image_view=AcquireVirtualCacheView(next,exception); 638 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1, 639 exception); 640 if (p == (const Quantum *) NULL) 641 { 642 image_view=DestroyCacheView(image_view); 643 break; 644 } 645 for (x=0; x < (ssize_t) image->columns; x++) 646 { 647 register ssize_t 648 i; 649 650 if (GetPixelReadMask(next,p) == 0) 651 { 652 p+=GetPixelChannels(next); 653 continue; 654 } 655 for (i=0; i < (ssize_t) GetPixelChannels(next); i++) 656 { 657 PixelChannel channel=GetPixelChannelChannel(image,i); 658 PixelTrait traits=GetPixelChannelTraits(next,channel); 659 PixelTrait evaluate_traits=GetPixelChannelTraits(image,channel); 660 if ((traits == UndefinedPixelTrait) || 661 (evaluate_traits == UndefinedPixelTrait)) 662 continue; 663 if ((traits & UpdatePixelTrait) == 0) 664 continue; 665 evaluate_pixel[x].channel[i]=ApplyEvaluateOperator( 666 random_info[id],GetPixelChannel(image,channel,p),j == 0 ? 667 AddEvaluateOperator : op,evaluate_pixel[x].channel[i]); 668 } 669 p+=GetPixelChannels(next); 670 } 671 image_view=DestroyCacheView(image_view); 672 next=GetNextImageInList(next); 673 } 674 for (x=0; x < (ssize_t) image->columns; x++) 675 { 676 register ssize_t 677 i; 678 679 switch (op) 680 { 681 case MeanEvaluateOperator: 682 { 683 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 684 evaluate_pixel[x].channel[i]/=(double) number_images; 685 break; 686 } 687 case MultiplyEvaluateOperator: 688 { 689 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 690 { 691 register ssize_t 692 j; 693 694 for (j=0; j < (ssize_t) (number_images-1); j++) 695 evaluate_pixel[x].channel[i]*=QuantumScale; 696 } 697 break; 698 } 699 case RootMeanSquareEvaluateOperator: 700 { 701 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 702 evaluate_pixel[x].channel[i]=sqrt(evaluate_pixel[x].channel[i]/ 703 number_images); 704 break; 705 } 706 default: 707 break; 708 } 709 } 710 for (x=0; x < (ssize_t) image->columns; x++) 711 { 712 register ssize_t 713 i; 714 715 if (GetPixelReadMask(image,q) == 0) 716 { 717 q+=GetPixelChannels(image); 718 continue; 719 } 720 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 721 { 722 PixelChannel channel=GetPixelChannelChannel(image,i); 723 PixelTrait traits=GetPixelChannelTraits(image,channel); 724 if (traits == UndefinedPixelTrait) 725 continue; 726 if ((traits & UpdatePixelTrait) == 0) 727 continue; 728 q[i]=ClampToQuantum(evaluate_pixel[x].channel[i]); 729 } 730 q+=GetPixelChannels(image); 731 } 732 if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse) 733 status=MagickFalse; 734 if (images->progress_monitor != (MagickProgressMonitor) NULL) 735 { 736 MagickBooleanType 737 proceed; 738 739 #if defined(MAGICKCORE_OPENMP_SUPPORT) 740 #pragma omp critical (MagickCore_EvaluateImages) 741 #endif 742 proceed=SetImageProgress(images,EvaluateImageTag,progress++, 743 image->rows); 744 if (proceed == MagickFalse) 745 status=MagickFalse; 746 } 747 } 748 } 749 evaluate_view=DestroyCacheView(evaluate_view); 750 evaluate_pixels=DestroyPixelThreadSet(evaluate_pixels); 751 random_info=DestroyRandomInfoThreadSet(random_info); 752 if (status == MagickFalse) 753 image=DestroyImage(image); 754 return(image); 755 } 756 757 MagickExport MagickBooleanType EvaluateImage(Image *image, 758 const MagickEvaluateOperator op,const double value,ExceptionInfo *exception) 759 { 760 CacheView 761 *image_view; 762 763 MagickBooleanType 764 status; 765 766 MagickOffsetType 767 progress; 768 769 RandomInfo 770 **magick_restrict random_info; 771 772 ssize_t 773 y; 774 775 #if defined(MAGICKCORE_OPENMP_SUPPORT) 776 unsigned long 777 key; 778 #endif 779 780 assert(image != (Image *) NULL); 781 assert(image->signature == MagickCoreSignature); 782 if (image->debug != MagickFalse) 783 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 784 assert(exception != (ExceptionInfo *) NULL); 785 assert(exception->signature == MagickCoreSignature); 786 if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse) 787 return(MagickFalse); 788 status=MagickTrue; 789 progress=0; 790 random_info=AcquireRandomInfoThreadSet(); 791 image_view=AcquireAuthenticCacheView(image,exception); 792 #if defined(MAGICKCORE_OPENMP_SUPPORT) 793 key=GetRandomSecretKey(random_info[0]); 794 #pragma omp parallel for schedule(static,4) shared(progress,status) \ 795 magick_threads(image,image,image->rows,key == ~0UL) 796 #endif 797 for (y=0; y < (ssize_t) image->rows; y++) 798 { 799 const int 800 id = GetOpenMPThreadId(); 801 802 register Quantum 803 *magick_restrict q; 804 805 register ssize_t 806 x; 807 808 if (status == MagickFalse) 809 continue; 810 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception); 811 if (q == (Quantum *) NULL) 812 { 813 status=MagickFalse; 814 continue; 815 } 816 for (x=0; x < (ssize_t) image->columns; x++) 817 { 818 double 819 result; 820 821 register ssize_t 822 i; 823 824 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 825 { 826 PixelChannel channel=GetPixelChannelChannel(image,i); 827 PixelTrait traits=GetPixelChannelTraits(image,channel); 828 if (traits == UndefinedPixelTrait) 829 continue; 830 if (((traits & CopyPixelTrait) != 0) || 831 (GetPixelReadMask(image,q) == 0)) 832 continue; 833 result=ApplyEvaluateOperator(random_info[id],q[i],op,value); 834 if (op == MeanEvaluateOperator) 835 result/=2.0; 836 q[i]=ClampToQuantum(result); 837 } 838 q+=GetPixelChannels(image); 839 } 840 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse) 841 status=MagickFalse; 842 if (image->progress_monitor != (MagickProgressMonitor) NULL) 843 { 844 MagickBooleanType 845 proceed; 846 847 #if defined(MAGICKCORE_OPENMP_SUPPORT) 848 #pragma omp critical (MagickCore_EvaluateImage) 849 #endif 850 proceed=SetImageProgress(image,EvaluateImageTag,progress++,image->rows); 851 if (proceed == MagickFalse) 852 status=MagickFalse; 853 } 854 } 855 image_view=DestroyCacheView(image_view); 856 random_info=DestroyRandomInfoThreadSet(random_info); 857 return(status); 858 } 859 860 /* 862 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 863 % % 864 % % 865 % % 866 % F u n c t i o n I m a g e % 867 % % 868 % % 869 % % 870 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 871 % 872 % FunctionImage() applies a value to the image with an arithmetic, relational, 873 % or logical operator to an image. Use these operations to lighten or darken 874 % an image, to increase or decrease contrast in an image, or to produce the 875 % "negative" of an image. 876 % 877 % The format of the FunctionImage method is: 878 % 879 % MagickBooleanType FunctionImage(Image *image, 880 % const MagickFunction function,const ssize_t number_parameters, 881 % const double *parameters,ExceptionInfo *exception) 882 % 883 % A description of each parameter follows: 884 % 885 % o image: the image. 886 % 887 % o function: A channel function. 888 % 889 % o parameters: one or more parameters. 890 % 891 % o exception: return any errors or warnings in this structure. 892 % 893 */ 894 895 static Quantum ApplyFunction(Quantum pixel,const MagickFunction function, 896 const size_t number_parameters,const double *parameters, 897 ExceptionInfo *exception) 898 { 899 double 900 result; 901 902 register ssize_t 903 i; 904 905 (void) exception; 906 result=0.0; 907 switch (function) 908 { 909 case PolynomialFunction: 910 { 911 /* 912 Polynomial: polynomial constants, highest to lowest order (e.g. c0*x^3+ 913 c1*x^2+c2*x+c3). 914 */ 915 result=0.0; 916 for (i=0; i < (ssize_t) number_parameters; i++) 917 result=result*QuantumScale*pixel+parameters[i]; 918 result*=QuantumRange; 919 break; 920 } 921 case SinusoidFunction: 922 { 923 double 924 amplitude, 925 bias, 926 frequency, 927 phase; 928 929 /* 930 Sinusoid: frequency, phase, amplitude, bias. 931 */ 932 frequency=(number_parameters >= 1) ? parameters[0] : 1.0; 933 phase=(number_parameters >= 2) ? parameters[1] : 0.0; 934 amplitude=(number_parameters >= 3) ? parameters[2] : 0.5; 935 bias=(number_parameters >= 4) ? parameters[3] : 0.5; 936 result=(double) (QuantumRange*(amplitude*sin((double) (2.0* 937 MagickPI*(frequency*QuantumScale*pixel+phase/360.0)))+bias)); 938 break; 939 } 940 case ArcsinFunction: 941 { 942 double 943 bias, 944 center, 945 range, 946 width; 947 948 /* 949 Arcsin (peged at range limits for invalid results): width, center, 950 range, and bias. 951 */ 952 width=(number_parameters >= 1) ? parameters[0] : 1.0; 953 center=(number_parameters >= 2) ? parameters[1] : 0.5; 954 range=(number_parameters >= 3) ? parameters[2] : 1.0; 955 bias=(number_parameters >= 4) ? parameters[3] : 0.5; 956 result=2.0/width*(QuantumScale*pixel-center); 957 if ( result <= -1.0 ) 958 result=bias-range/2.0; 959 else 960 if (result >= 1.0) 961 result=bias+range/2.0; 962 else 963 result=(double) (range/MagickPI*asin((double) result)+bias); 964 result*=QuantumRange; 965 break; 966 } 967 case ArctanFunction: 968 { 969 double 970 center, 971 bias, 972 range, 973 slope; 974 975 /* 976 Arctan: slope, center, range, and bias. 977 */ 978 slope=(number_parameters >= 1) ? parameters[0] : 1.0; 979 center=(number_parameters >= 2) ? parameters[1] : 0.5; 980 range=(number_parameters >= 3) ? parameters[2] : 1.0; 981 bias=(number_parameters >= 4) ? parameters[3] : 0.5; 982 result=(double) (MagickPI*slope*(QuantumScale*pixel-center)); 983 result=(double) (QuantumRange*(range/MagickPI*atan((double) 984 result)+bias)); 985 break; 986 } 987 case UndefinedFunction: 988 break; 989 } 990 return(ClampToQuantum(result)); 991 } 992 993 MagickExport MagickBooleanType FunctionImage(Image *image, 994 const MagickFunction function,const size_t number_parameters, 995 const double *parameters,ExceptionInfo *exception) 996 { 997 #define FunctionImageTag "Function/Image " 998 999 CacheView 1000 *image_view; 1001 1002 MagickBooleanType 1003 status; 1004 1005 MagickOffsetType 1006 progress; 1007 1008 ssize_t 1009 y; 1010 1011 assert(image != (Image *) NULL); 1012 assert(image->signature == MagickCoreSignature); 1013 if (image->debug != MagickFalse) 1014 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 1015 assert(exception != (ExceptionInfo *) NULL); 1016 assert(exception->signature == MagickCoreSignature); 1017 #if defined(MAGICKCORE_OPENCL_SUPPORT) 1018 if (AccelerateFunctionImage(image,function,number_parameters,parameters, 1019 exception) != MagickFalse) 1020 return(MagickTrue); 1021 #endif 1022 if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse) 1023 return(MagickFalse); 1024 status=MagickTrue; 1025 progress=0; 1026 image_view=AcquireAuthenticCacheView(image,exception); 1027 #if defined(MAGICKCORE_OPENMP_SUPPORT) 1028 #pragma omp parallel for schedule(static,4) shared(progress,status) \ 1029 magick_threads(image,image,image->rows,1) 1030 #endif 1031 for (y=0; y < (ssize_t) image->rows; y++) 1032 { 1033 register Quantum 1034 *magick_restrict q; 1035 1036 register ssize_t 1037 x; 1038 1039 if (status == MagickFalse) 1040 continue; 1041 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception); 1042 if (q == (Quantum *) NULL) 1043 { 1044 status=MagickFalse; 1045 continue; 1046 } 1047 for (x=0; x < (ssize_t) image->columns; x++) 1048 { 1049 register ssize_t 1050 i; 1051 1052 if (GetPixelReadMask(image,q) == 0) 1053 { 1054 q+=GetPixelChannels(image); 1055 continue; 1056 } 1057 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 1058 { 1059 PixelChannel channel=GetPixelChannelChannel(image,i); 1060 PixelTrait traits=GetPixelChannelTraits(image,channel); 1061 if (traits == UndefinedPixelTrait) 1062 continue; 1063 if ((traits & UpdatePixelTrait) == 0) 1064 continue; 1065 q[i]=ApplyFunction(q[i],function,number_parameters,parameters, 1066 exception); 1067 } 1068 q+=GetPixelChannels(image); 1069 } 1070 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse) 1071 status=MagickFalse; 1072 if (image->progress_monitor != (MagickProgressMonitor) NULL) 1073 { 1074 MagickBooleanType 1075 proceed; 1076 1077 #if defined(MAGICKCORE_OPENMP_SUPPORT) 1078 #pragma omp critical (MagickCore_FunctionImage) 1079 #endif 1080 proceed=SetImageProgress(image,FunctionImageTag,progress++,image->rows); 1081 if (proceed == MagickFalse) 1082 status=MagickFalse; 1083 } 1084 } 1085 image_view=DestroyCacheView(image_view); 1086 return(status); 1087 } 1088 1089 /* 1091 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1092 % % 1093 % % 1094 % % 1095 % G e t I m a g e E n t r o p y % 1096 % % 1097 % % 1098 % % 1099 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1100 % 1101 % GetImageEntropy() returns the entropy of one or more image channels. 1102 % 1103 % The format of the GetImageEntropy method is: 1104 % 1105 % MagickBooleanType GetImageEntropy(const Image *image,double *entropy, 1106 % ExceptionInfo *exception) 1107 % 1108 % A description of each parameter follows: 1109 % 1110 % o image: the image. 1111 % 1112 % o entropy: the average entropy of the selected channels. 1113 % 1114 % o exception: return any errors or warnings in this structure. 1115 % 1116 */ 1117 MagickExport MagickBooleanType GetImageEntropy(const Image *image, 1118 double *entropy,ExceptionInfo *exception) 1119 { 1120 double 1121 area; 1122 1123 ChannelStatistics 1124 *channel_statistics; 1125 1126 register ssize_t 1127 i; 1128 1129 assert(image != (Image *) NULL); 1130 assert(image->signature == MagickCoreSignature); 1131 if (image->debug != MagickFalse) 1132 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 1133 channel_statistics=GetImageStatistics(image,exception); 1134 if (channel_statistics == (ChannelStatistics *) NULL) 1135 return(MagickFalse); 1136 area=0.0; 1137 channel_statistics[CompositePixelChannel].entropy=0.0; 1138 channel_statistics[CompositePixelChannel].standard_deviation=0.0; 1139 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 1140 { 1141 PixelChannel channel=GetPixelChannelChannel(image,i); 1142 PixelTrait traits=GetPixelChannelTraits(image,channel); 1143 if (traits == UndefinedPixelTrait) 1144 continue; 1145 if ((traits & UpdatePixelTrait) == 0) 1146 continue; 1147 channel_statistics[CompositePixelChannel].entropy+= 1148 channel_statistics[i].entropy; 1149 area++; 1150 } 1151 if (area > MagickEpsilon) 1152 { 1153 channel_statistics[CompositePixelChannel].entropy/=area; 1154 channel_statistics[CompositePixelChannel].standard_deviation= 1155 sqrt(channel_statistics[CompositePixelChannel].standard_deviation/area); 1156 } 1157 *entropy=channel_statistics[CompositePixelChannel].entropy; 1158 channel_statistics=(ChannelStatistics *) RelinquishMagickMemory( 1159 channel_statistics); 1160 return(MagickTrue); 1161 } 1162 1163 /* 1165 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1166 % % 1167 % % 1168 % % 1169 % G e t I m a g e E x t r e m a % 1170 % % 1171 % % 1172 % % 1173 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1174 % 1175 % GetImageExtrema() returns the extrema of one or more image channels. 1176 % 1177 % The format of the GetImageExtrema method is: 1178 % 1179 % MagickBooleanType GetImageExtrema(const Image *image,size_t *minima, 1180 % size_t *maxima,ExceptionInfo *exception) 1181 % 1182 % A description of each parameter follows: 1183 % 1184 % o image: the image. 1185 % 1186 % o minima: the minimum value in the channel. 1187 % 1188 % o maxima: the maximum value in the channel. 1189 % 1190 % o exception: return any errors or warnings in this structure. 1191 % 1192 */ 1193 MagickExport MagickBooleanType GetImageExtrema(const Image *image, 1194 size_t *minima,size_t *maxima,ExceptionInfo *exception) 1195 { 1196 double 1197 max, 1198 min; 1199 1200 MagickBooleanType 1201 status; 1202 1203 assert(image != (Image *) NULL); 1204 assert(image->signature == MagickCoreSignature); 1205 if (image->debug != MagickFalse) 1206 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 1207 status=GetImageRange(image,&min,&max,exception); 1208 *minima=(size_t) ceil(min-0.5); 1209 *maxima=(size_t) floor(max+0.5); 1210 return(status); 1211 } 1212 1213 /* 1215 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1216 % % 1217 % % 1218 % % 1219 % G e t I m a g e K u r t o s i s % 1220 % % 1221 % % 1222 % % 1223 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1224 % 1225 % GetImageKurtosis() returns the kurtosis and skewness of one or more image 1226 % channels. 1227 % 1228 % The format of the GetImageKurtosis method is: 1229 % 1230 % MagickBooleanType GetImageKurtosis(const Image *image,double *kurtosis, 1231 % double *skewness,ExceptionInfo *exception) 1232 % 1233 % A description of each parameter follows: 1234 % 1235 % o image: the image. 1236 % 1237 % o kurtosis: the kurtosis of the channel. 1238 % 1239 % o skewness: the skewness of the channel. 1240 % 1241 % o exception: return any errors or warnings in this structure. 1242 % 1243 */ 1244 MagickExport MagickBooleanType GetImageKurtosis(const Image *image, 1245 double *kurtosis,double *skewness,ExceptionInfo *exception) 1246 { 1247 CacheView 1248 *image_view; 1249 1250 double 1251 area, 1252 mean, 1253 standard_deviation, 1254 sum_squares, 1255 sum_cubes, 1256 sum_fourth_power; 1257 1258 MagickBooleanType 1259 status; 1260 1261 ssize_t 1262 y; 1263 1264 assert(image != (Image *) NULL); 1265 assert(image->signature == MagickCoreSignature); 1266 if (image->debug != MagickFalse) 1267 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 1268 status=MagickTrue; 1269 *kurtosis=0.0; 1270 *skewness=0.0; 1271 area=0.0; 1272 mean=0.0; 1273 standard_deviation=0.0; 1274 sum_squares=0.0; 1275 sum_cubes=0.0; 1276 sum_fourth_power=0.0; 1277 image_view=AcquireVirtualCacheView(image,exception); 1278 #if defined(MAGICKCORE_OPENMP_SUPPORT) 1279 #pragma omp parallel for schedule(static,4) shared(status) \ 1280 magick_threads(image,image,image->rows,1) 1281 #endif 1282 for (y=0; y < (ssize_t) image->rows; y++) 1283 { 1284 register const Quantum 1285 *magick_restrict p; 1286 1287 register ssize_t 1288 x; 1289 1290 if (status == MagickFalse) 1291 continue; 1292 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception); 1293 if (p == (const Quantum *) NULL) 1294 { 1295 status=MagickFalse; 1296 continue; 1297 } 1298 for (x=0; x < (ssize_t) image->columns; x++) 1299 { 1300 register ssize_t 1301 i; 1302 1303 if (GetPixelReadMask(image,p) == 0) 1304 { 1305 p+=GetPixelChannels(image); 1306 continue; 1307 } 1308 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 1309 { 1310 PixelChannel channel=GetPixelChannelChannel(image,i); 1311 PixelTrait traits=GetPixelChannelTraits(image,channel); 1312 if (traits == UndefinedPixelTrait) 1313 continue; 1314 if ((traits & UpdatePixelTrait) == 0) 1315 continue; 1316 #if defined(MAGICKCORE_OPENMP_SUPPORT) 1317 #pragma omp critical (MagickCore_GetImageKurtosis) 1318 #endif 1319 { 1320 mean+=p[i]; 1321 sum_squares+=(double) p[i]*p[i]; 1322 sum_cubes+=(double) p[i]*p[i]*p[i]; 1323 sum_fourth_power+=(double) p[i]*p[i]*p[i]*p[i]; 1324 area++; 1325 } 1326 } 1327 p+=GetPixelChannels(image); 1328 } 1329 } 1330 image_view=DestroyCacheView(image_view); 1331 if (area != 0.0) 1332 { 1333 mean/=area; 1334 sum_squares/=area; 1335 sum_cubes/=area; 1336 sum_fourth_power/=area; 1337 } 1338 standard_deviation=sqrt(sum_squares-(mean*mean)); 1339 if (standard_deviation != 0.0) 1340 { 1341 *kurtosis=sum_fourth_power-4.0*mean*sum_cubes+6.0*mean*mean*sum_squares- 1342 3.0*mean*mean*mean*mean; 1343 *kurtosis/=standard_deviation*standard_deviation*standard_deviation* 1344 standard_deviation; 1345 *kurtosis-=3.0; 1346 *skewness=sum_cubes-3.0*mean*sum_squares+2.0*mean*mean*mean; 1347 *skewness/=standard_deviation*standard_deviation*standard_deviation; 1348 } 1349 return(status); 1350 } 1351 1352 /* 1354 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1355 % % 1356 % % 1357 % % 1358 % G e t I m a g e M e a n % 1359 % % 1360 % % 1361 % % 1362 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1363 % 1364 % GetImageMean() returns the mean and standard deviation of one or more image 1365 % channels. 1366 % 1367 % The format of the GetImageMean method is: 1368 % 1369 % MagickBooleanType GetImageMean(const Image *image,double *mean, 1370 % double *standard_deviation,ExceptionInfo *exception) 1371 % 1372 % A description of each parameter follows: 1373 % 1374 % o image: the image. 1375 % 1376 % o mean: the average value in the channel. 1377 % 1378 % o standard_deviation: the standard deviation of the channel. 1379 % 1380 % o exception: return any errors or warnings in this structure. 1381 % 1382 */ 1383 MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean, 1384 double *standard_deviation,ExceptionInfo *exception) 1385 { 1386 double 1387 area; 1388 1389 ChannelStatistics 1390 *channel_statistics; 1391 1392 register ssize_t 1393 i; 1394 1395 assert(image != (Image *) NULL); 1396 assert(image->signature == MagickCoreSignature); 1397 if (image->debug != MagickFalse) 1398 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 1399 channel_statistics=GetImageStatistics(image,exception); 1400 if (channel_statistics == (ChannelStatistics *) NULL) 1401 return(MagickFalse); 1402 area=0.0; 1403 channel_statistics[CompositePixelChannel].mean=0.0; 1404 channel_statistics[CompositePixelChannel].standard_deviation=0.0; 1405 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 1406 { 1407 PixelChannel channel=GetPixelChannelChannel(image,i); 1408 PixelTrait traits=GetPixelChannelTraits(image,channel); 1409 if (traits == UndefinedPixelTrait) 1410 continue; 1411 if ((traits & UpdatePixelTrait) == 0) 1412 continue; 1413 channel_statistics[CompositePixelChannel].mean+=channel_statistics[i].mean; 1414 channel_statistics[CompositePixelChannel].standard_deviation+= 1415 channel_statistics[i].variance-channel_statistics[i].mean* 1416 channel_statistics[i].mean; 1417 area++; 1418 } 1419 if (area > MagickEpsilon) 1420 { 1421 channel_statistics[CompositePixelChannel].mean/=area; 1422 channel_statistics[CompositePixelChannel].standard_deviation= 1423 sqrt(channel_statistics[CompositePixelChannel].standard_deviation/area); 1424 } 1425 *mean=channel_statistics[CompositePixelChannel].mean; 1426 *standard_deviation= 1427 channel_statistics[CompositePixelChannel].standard_deviation; 1428 channel_statistics=(ChannelStatistics *) RelinquishMagickMemory( 1429 channel_statistics); 1430 return(MagickTrue); 1431 } 1432 1433 /* 1435 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1436 % % 1437 % % 1438 % % 1439 % G e t I m a g e M o m e n t s % 1440 % % 1441 % % 1442 % % 1443 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1444 % 1445 % GetImageMoments() returns the normalized moments of one or more image 1446 % channels. 1447 % 1448 % The format of the GetImageMoments method is: 1449 % 1450 % ChannelMoments *GetImageMoments(const Image *image, 1451 % ExceptionInfo *exception) 1452 % 1453 % A description of each parameter follows: 1454 % 1455 % o image: the image. 1456 % 1457 % o exception: return any errors or warnings in this structure. 1458 % 1459 */ 1460 1461 static size_t GetImageChannels(const Image *image) 1462 { 1463 register ssize_t 1464 i; 1465 1466 size_t 1467 channels; 1468 1469 channels=0; 1470 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 1471 { 1472 PixelChannel channel=GetPixelChannelChannel(image,i); 1473 PixelTrait traits=GetPixelChannelTraits(image,channel); 1474 if ((traits & UpdatePixelTrait) != 0) 1475 channels++; 1476 } 1477 return((size_t) (channels == 0 ? 1 : channels)); 1478 } 1479 1480 MagickExport ChannelMoments *GetImageMoments(const Image *image, 1481 ExceptionInfo *exception) 1482 { 1483 #define MaxNumberImageMoments 8 1484 1485 CacheView 1486 *image_view; 1487 1488 ChannelMoments 1489 *channel_moments; 1490 1491 double 1492 M00[MaxPixelChannels+1], 1493 M01[MaxPixelChannels+1], 1494 M02[MaxPixelChannels+1], 1495 M03[MaxPixelChannels+1], 1496 M10[MaxPixelChannels+1], 1497 M11[MaxPixelChannels+1], 1498 M12[MaxPixelChannels+1], 1499 M20[MaxPixelChannels+1], 1500 M21[MaxPixelChannels+1], 1501 M22[MaxPixelChannels+1], 1502 M30[MaxPixelChannels+1]; 1503 1504 PointInfo 1505 centroid[MaxPixelChannels+1]; 1506 1507 ssize_t 1508 channel, 1509 y; 1510 1511 assert(image != (Image *) NULL); 1512 assert(image->signature == MagickCoreSignature); 1513 if (image->debug != MagickFalse) 1514 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 1515 channel_moments=(ChannelMoments *) AcquireQuantumMemory(MaxPixelChannels+1, 1516 sizeof(*channel_moments)); 1517 if (channel_moments == (ChannelMoments *) NULL) 1518 return(channel_moments); 1519 (void) ResetMagickMemory(channel_moments,0,(MaxPixelChannels+1)* 1520 sizeof(*channel_moments)); 1521 (void) ResetMagickMemory(centroid,0,sizeof(centroid)); 1522 (void) ResetMagickMemory(M00,0,sizeof(M00)); 1523 (void) ResetMagickMemory(M01,0,sizeof(M01)); 1524 (void) ResetMagickMemory(M02,0,sizeof(M02)); 1525 (void) ResetMagickMemory(M03,0,sizeof(M03)); 1526 (void) ResetMagickMemory(M10,0,sizeof(M10)); 1527 (void) ResetMagickMemory(M11,0,sizeof(M11)); 1528 (void) ResetMagickMemory(M12,0,sizeof(M12)); 1529 (void) ResetMagickMemory(M20,0,sizeof(M20)); 1530 (void) ResetMagickMemory(M21,0,sizeof(M21)); 1531 (void) ResetMagickMemory(M22,0,sizeof(M22)); 1532 (void) ResetMagickMemory(M30,0,sizeof(M30)); 1533 image_view=AcquireVirtualCacheView(image,exception); 1534 for (y=0; y < (ssize_t) image->rows; y++) 1535 { 1536 register const Quantum 1537 *magick_restrict p; 1538 1539 register ssize_t 1540 x; 1541 1542 /* 1543 Compute center of mass (centroid). 1544 */ 1545 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception); 1546 if (p == (const Quantum *) NULL) 1547 break; 1548 for (x=0; x < (ssize_t) image->columns; x++) 1549 { 1550 register ssize_t 1551 i; 1552 1553 if (GetPixelReadMask(image,p) == 0) 1554 { 1555 p+=GetPixelChannels(image); 1556 continue; 1557 } 1558 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 1559 { 1560 PixelChannel channel=GetPixelChannelChannel(image,i); 1561 PixelTrait traits=GetPixelChannelTraits(image,channel); 1562 if (traits == UndefinedPixelTrait) 1563 continue; 1564 if ((traits & UpdatePixelTrait) == 0) 1565 continue; 1566 M00[channel]+=QuantumScale*p[i]; 1567 M00[MaxPixelChannels]+=QuantumScale*p[i]; 1568 M10[channel]+=x*QuantumScale*p[i]; 1569 M10[MaxPixelChannels]+=x*QuantumScale*p[i]; 1570 M01[channel]+=y*QuantumScale*p[i]; 1571 M01[MaxPixelChannels]+=y*QuantumScale*p[i]; 1572 } 1573 p+=GetPixelChannels(image); 1574 } 1575 } 1576 for (channel=0; channel <= MaxPixelChannels; channel++) 1577 { 1578 /* 1579 Compute center of mass (centroid). 1580 */ 1581 if (M00[channel] < MagickEpsilon) 1582 { 1583 M00[channel]+=MagickEpsilon; 1584 centroid[channel].x=(double) image->columns/2.0; 1585 centroid[channel].y=(double) image->rows/2.0; 1586 continue; 1587 } 1588 M00[channel]+=MagickEpsilon; 1589 centroid[channel].x=M10[channel]/M00[channel]; 1590 centroid[channel].y=M01[channel]/M00[channel]; 1591 } 1592 for (y=0; y < (ssize_t) image->rows; y++) 1593 { 1594 register const Quantum 1595 *magick_restrict p; 1596 1597 register ssize_t 1598 x; 1599 1600 /* 1601 Compute the image moments. 1602 */ 1603 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception); 1604 if (p == (const Quantum *) NULL) 1605 break; 1606 for (x=0; x < (ssize_t) image->columns; x++) 1607 { 1608 register ssize_t 1609 i; 1610 1611 if (GetPixelReadMask(image,p) == 0) 1612 { 1613 p+=GetPixelChannels(image); 1614 continue; 1615 } 1616 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 1617 { 1618 PixelChannel channel=GetPixelChannelChannel(image,i); 1619 PixelTrait traits=GetPixelChannelTraits(image,channel); 1620 if (traits == UndefinedPixelTrait) 1621 continue; 1622 if ((traits & UpdatePixelTrait) == 0) 1623 continue; 1624 M11[channel]+=(x-centroid[channel].x)*(y-centroid[channel].y)* 1625 QuantumScale*p[i]; 1626 M11[MaxPixelChannels]+=(x-centroid[channel].x)*(y-centroid[channel].y)* 1627 QuantumScale*p[i]; 1628 M20[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)* 1629 QuantumScale*p[i]; 1630 M20[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)* 1631 QuantumScale*p[i]; 1632 M02[channel]+=(y-centroid[channel].y)*(y-centroid[channel].y)* 1633 QuantumScale*p[i]; 1634 M02[MaxPixelChannels]+=(y-centroid[channel].y)*(y-centroid[channel].y)* 1635 QuantumScale*p[i]; 1636 M21[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)* 1637 (y-centroid[channel].y)*QuantumScale*p[i]; 1638 M21[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)* 1639 (y-centroid[channel].y)*QuantumScale*p[i]; 1640 M12[channel]+=(x-centroid[channel].x)*(y-centroid[channel].y)* 1641 (y-centroid[channel].y)*QuantumScale*p[i]; 1642 M12[MaxPixelChannels]+=(x-centroid[channel].x)*(y-centroid[channel].y)* 1643 (y-centroid[channel].y)*QuantumScale*p[i]; 1644 M22[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)* 1645 (y-centroid[channel].y)*(y-centroid[channel].y)*QuantumScale*p[i]; 1646 M22[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)* 1647 (y-centroid[channel].y)*(y-centroid[channel].y)*QuantumScale*p[i]; 1648 M30[channel]+=(x-centroid[channel].x)*(x-centroid[channel].x)* 1649 (x-centroid[channel].x)*QuantumScale*p[i]; 1650 M30[MaxPixelChannels]+=(x-centroid[channel].x)*(x-centroid[channel].x)* 1651 (x-centroid[channel].x)*QuantumScale*p[i]; 1652 M03[channel]+=(y-centroid[channel].y)*(y-centroid[channel].y)* 1653 (y-centroid[channel].y)*QuantumScale*p[i]; 1654 M03[MaxPixelChannels]+=(y-centroid[channel].y)*(y-centroid[channel].y)* 1655 (y-centroid[channel].y)*QuantumScale*p[i]; 1656 } 1657 p+=GetPixelChannels(image); 1658 } 1659 } 1660 M00[MaxPixelChannels]/=GetImageChannels(image); 1661 M01[MaxPixelChannels]/=GetImageChannels(image); 1662 M02[MaxPixelChannels]/=GetImageChannels(image); 1663 M03[MaxPixelChannels]/=GetImageChannels(image); 1664 M10[MaxPixelChannels]/=GetImageChannels(image); 1665 M11[MaxPixelChannels]/=GetImageChannels(image); 1666 M12[MaxPixelChannels]/=GetImageChannels(image); 1667 M20[MaxPixelChannels]/=GetImageChannels(image); 1668 M21[MaxPixelChannels]/=GetImageChannels(image); 1669 M22[MaxPixelChannels]/=GetImageChannels(image); 1670 M30[MaxPixelChannels]/=GetImageChannels(image); 1671 for (channel=0; channel <= MaxPixelChannels; channel++) 1672 { 1673 /* 1674 Compute elliptical angle, major and minor axes, eccentricity, & intensity. 1675 */ 1676 channel_moments[channel].centroid=centroid[channel]; 1677 channel_moments[channel].ellipse_axis.x=sqrt((2.0/M00[channel])* 1678 ((M20[channel]+M02[channel])+sqrt(4.0*M11[channel]*M11[channel]+ 1679 (M20[channel]-M02[channel])*(M20[channel]-M02[channel])))); 1680 channel_moments[channel].ellipse_axis.y=sqrt((2.0/M00[channel])* 1681 ((M20[channel]+M02[channel])-sqrt(4.0*M11[channel]*M11[channel]+ 1682 (M20[channel]-M02[channel])*(M20[channel]-M02[channel])))); 1683 channel_moments[channel].ellipse_angle=RadiansToDegrees(0.5*atan(2.0* 1684 M11[channel]/(M20[channel]-M02[channel]+MagickEpsilon))); 1685 channel_moments[channel].ellipse_eccentricity=sqrt(1.0-( 1686 channel_moments[channel].ellipse_axis.y/ 1687 (channel_moments[channel].ellipse_axis.x+MagickEpsilon))); 1688 channel_moments[channel].ellipse_intensity=M00[channel]/ 1689 (MagickPI*channel_moments[channel].ellipse_axis.x* 1690 channel_moments[channel].ellipse_axis.y+MagickEpsilon); 1691 } 1692 for (channel=0; channel <= MaxPixelChannels; channel++) 1693 { 1694 /* 1695 Normalize image moments. 1696 */ 1697 M10[channel]=0.0; 1698 M01[channel]=0.0; 1699 M11[channel]/=pow(M00[channel],1.0+(1.0+1.0)/2.0); 1700 M20[channel]/=pow(M00[channel],1.0+(2.0+0.0)/2.0); 1701 M02[channel]/=pow(M00[channel],1.0+(0.0+2.0)/2.0); 1702 M21[channel]/=pow(M00[channel],1.0+(2.0+1.0)/2.0); 1703 M12[channel]/=pow(M00[channel],1.0+(1.0+2.0)/2.0); 1704 M22[channel]/=pow(M00[channel],1.0+(2.0+2.0)/2.0); 1705 M30[channel]/=pow(M00[channel],1.0+(3.0+0.0)/2.0); 1706 M03[channel]/=pow(M00[channel],1.0+(0.0+3.0)/2.0); 1707 M00[channel]=1.0; 1708 } 1709 image_view=DestroyCacheView(image_view); 1710 for (channel=0; channel <= MaxPixelChannels; channel++) 1711 { 1712 /* 1713 Compute Hu invariant moments. 1714 */ 1715 channel_moments[channel].invariant[0]=M20[channel]+M02[channel]; 1716 channel_moments[channel].invariant[1]=(M20[channel]-M02[channel])* 1717 (M20[channel]-M02[channel])+4.0*M11[channel]*M11[channel]; 1718 channel_moments[channel].invariant[2]=(M30[channel]-3.0*M12[channel])* 1719 (M30[channel]-3.0*M12[channel])+(3.0*M21[channel]-M03[channel])* 1720 (3.0*M21[channel]-M03[channel]); 1721 channel_moments[channel].invariant[3]=(M30[channel]+M12[channel])* 1722 (M30[channel]+M12[channel])+(M21[channel]+M03[channel])* 1723 (M21[channel]+M03[channel]); 1724 channel_moments[channel].invariant[4]=(M30[channel]-3.0*M12[channel])* 1725 (M30[channel]+M12[channel])*((M30[channel]+M12[channel])* 1726 (M30[channel]+M12[channel])-3.0*(M21[channel]+M03[channel])* 1727 (M21[channel]+M03[channel]))+(3.0*M21[channel]-M03[channel])* 1728 (M21[channel]+M03[channel])*(3.0*(M30[channel]+M12[channel])* 1729 (M30[channel]+M12[channel])-(M21[channel]+M03[channel])* 1730 (M21[channel]+M03[channel])); 1731 channel_moments[channel].invariant[5]=(M20[channel]-M02[channel])* 1732 ((M30[channel]+M12[channel])*(M30[channel]+M12[channel])- 1733 (M21[channel]+M03[channel])*(M21[channel]+M03[channel]))+ 1734 4.0*M11[channel]*(M30[channel]+M12[channel])*(M21[channel]+M03[channel]); 1735 channel_moments[channel].invariant[6]=(3.0*M21[channel]-M03[channel])* 1736 (M30[channel]+M12[channel])*((M30[channel]+M12[channel])* 1737 (M30[channel]+M12[channel])-3.0*(M21[channel]+M03[channel])* 1738 (M21[channel]+M03[channel]))-(M30[channel]-3*M12[channel])* 1739 (M21[channel]+M03[channel])*(3.0*(M30[channel]+M12[channel])* 1740 (M30[channel]+M12[channel])-(M21[channel]+M03[channel])* 1741 (M21[channel]+M03[channel])); 1742 channel_moments[channel].invariant[7]=M11[channel]*((M30[channel]+ 1743 M12[channel])*(M30[channel]+M12[channel])-(M03[channel]+M21[channel])* 1744 (M03[channel]+M21[channel]))-(M20[channel]-M02[channel])* 1745 (M30[channel]+M12[channel])*(M03[channel]+M21[channel]); 1746 } 1747 if (y < (ssize_t) image->rows) 1748 channel_moments=(ChannelMoments *) RelinquishMagickMemory(channel_moments); 1749 return(channel_moments); 1750 } 1751 1752 /* 1754 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1755 % % 1756 % % 1757 % % 1758 % G e t I m a g e C h a n n e l P e r c e p t u a l H a s h % 1759 % % 1760 % % 1761 % % 1762 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1763 % 1764 % GetImagePerceptualHash() returns the perceptual hash of one or more 1765 % image channels. 1766 % 1767 % The format of the GetImagePerceptualHash method is: 1768 % 1769 % ChannelPerceptualHash *GetImagePerceptualHash(const Image *image, 1770 % ExceptionInfo *exception) 1771 % 1772 % A description of each parameter follows: 1773 % 1774 % o image: the image. 1775 % 1776 % o exception: return any errors or warnings in this structure. 1777 % 1778 */ 1779 1780 static inline double MagickLog10(const double x) 1781 { 1782 #define Log10Epsilon (1.0e-11) 1783 1784 if (fabs(x) < Log10Epsilon) 1785 return(log10(Log10Epsilon)); 1786 return(log10(fabs(x))); 1787 } 1788 1789 MagickExport ChannelPerceptualHash *GetImagePerceptualHash( 1790 const Image *image,ExceptionInfo *exception) 1791 { 1792 ChannelMoments 1793 *moments; 1794 1795 ChannelPerceptualHash 1796 *perceptual_hash; 1797 1798 Image 1799 *hash_image; 1800 1801 MagickBooleanType 1802 status; 1803 1804 register ssize_t 1805 i; 1806 1807 ssize_t 1808 channel; 1809 1810 /* 1811 Blur then transform to sRGB colorspace. 1812 */ 1813 hash_image=BlurImage(image,0.0,1.0,exception); 1814 if (hash_image == (Image *) NULL) 1815 return((ChannelPerceptualHash *) NULL); 1816 hash_image->depth=8; 1817 status=TransformImageColorspace(hash_image,sRGBColorspace,exception); 1818 if (status == MagickFalse) 1819 return((ChannelPerceptualHash *) NULL); 1820 moments=GetImageMoments(hash_image,exception); 1821 hash_image=DestroyImage(hash_image); 1822 if (moments == (ChannelMoments *) NULL) 1823 return((ChannelPerceptualHash *) NULL); 1824 perceptual_hash=(ChannelPerceptualHash *) AcquireQuantumMemory( 1825 MaxPixelChannels+1UL,sizeof(*perceptual_hash)); 1826 if (perceptual_hash == (ChannelPerceptualHash *) NULL) 1827 return((ChannelPerceptualHash *) NULL); 1828 for (channel=0; channel <= MaxPixelChannels; channel++) 1829 for (i=0; i < MaximumNumberOfImageMoments; i++) 1830 perceptual_hash[channel].srgb_hu_phash[i]= 1831 (-MagickLog10(moments[channel].invariant[i])); 1832 moments=(ChannelMoments *) RelinquishMagickMemory(moments); 1833 /* 1834 Blur then transform to HCLp colorspace. 1835 */ 1836 hash_image=BlurImage(image,0.0,1.0,exception); 1837 if (hash_image == (Image *) NULL) 1838 { 1839 perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory( 1840 perceptual_hash); 1841 return((ChannelPerceptualHash *) NULL); 1842 } 1843 hash_image->depth=8; 1844 status=TransformImageColorspace(hash_image,HCLpColorspace,exception); 1845 if (status == MagickFalse) 1846 { 1847 perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory( 1848 perceptual_hash); 1849 return((ChannelPerceptualHash *) NULL); 1850 } 1851 moments=GetImageMoments(hash_image,exception); 1852 hash_image=DestroyImage(hash_image); 1853 if (moments == (ChannelMoments *) NULL) 1854 { 1855 perceptual_hash=(ChannelPerceptualHash *) RelinquishMagickMemory( 1856 perceptual_hash); 1857 return((ChannelPerceptualHash *) NULL); 1858 } 1859 for (channel=0; channel <= MaxPixelChannels; channel++) 1860 for (i=0; i < MaximumNumberOfImageMoments; i++) 1861 perceptual_hash[channel].hclp_hu_phash[i]= 1862 (-MagickLog10(moments[channel].invariant[i])); 1863 moments=(ChannelMoments *) RelinquishMagickMemory(moments); 1864 return(perceptual_hash); 1865 } 1866 1867 /* 1869 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1870 % % 1871 % % 1872 % % 1873 % G e t I m a g e R a n g e % 1874 % % 1875 % % 1876 % % 1877 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1878 % 1879 % GetImageRange() returns the range of one or more image channels. 1880 % 1881 % The format of the GetImageRange method is: 1882 % 1883 % MagickBooleanType GetImageRange(const Image *image,double *minima, 1884 % double *maxima,ExceptionInfo *exception) 1885 % 1886 % A description of each parameter follows: 1887 % 1888 % o image: the image. 1889 % 1890 % o minima: the minimum value in the channel. 1891 % 1892 % o maxima: the maximum value in the channel. 1893 % 1894 % o exception: return any errors or warnings in this structure. 1895 % 1896 */ 1897 MagickExport MagickBooleanType GetImageRange(const Image *image,double *minima, 1898 double *maxima,ExceptionInfo *exception) 1899 { 1900 CacheView 1901 *image_view; 1902 1903 MagickBooleanType 1904 initialize, 1905 status; 1906 1907 ssize_t 1908 y; 1909 1910 assert(image != (Image *) NULL); 1911 assert(image->signature == MagickCoreSignature); 1912 if (image->debug != MagickFalse) 1913 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 1914 status=MagickTrue; 1915 initialize=MagickTrue; 1916 *maxima=0.0; 1917 *minima=0.0; 1918 image_view=AcquireVirtualCacheView(image,exception); 1919 #if defined(MAGICKCORE_OPENMP_SUPPORT) 1920 #pragma omp parallel for schedule(static,4) shared(status,initialize) \ 1921 magick_threads(image,image,image->rows,1) 1922 #endif 1923 for (y=0; y < (ssize_t) image->rows; y++) 1924 { 1925 double 1926 row_maxima = 0.0, 1927 row_minima = 0.0; 1928 1929 MagickBooleanType 1930 row_initialize; 1931 1932 register const Quantum 1933 *magick_restrict p; 1934 1935 register ssize_t 1936 x; 1937 1938 if (status == MagickFalse) 1939 continue; 1940 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception); 1941 if (p == (const Quantum *) NULL) 1942 { 1943 status=MagickFalse; 1944 continue; 1945 } 1946 row_initialize=MagickTrue; 1947 for (x=0; x < (ssize_t) image->columns; x++) 1948 { 1949 register ssize_t 1950 i; 1951 1952 if (GetPixelReadMask(image,p) == 0) 1953 { 1954 p+=GetPixelChannels(image); 1955 continue; 1956 } 1957 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 1958 { 1959 PixelChannel channel=GetPixelChannelChannel(image,i); 1960 PixelTrait traits=GetPixelChannelTraits(image,channel); 1961 if (traits == UndefinedPixelTrait) 1962 continue; 1963 if ((traits & UpdatePixelTrait) == 0) 1964 continue; 1965 if (row_initialize != MagickFalse) 1966 { 1967 row_minima=(double) p[i]; 1968 row_maxima=(double) p[i]; 1969 row_initialize=MagickFalse; 1970 } 1971 else 1972 { 1973 if ((double) p[i] < row_minima) 1974 row_minima=(double) p[i]; 1975 if ((double) p[i] > row_maxima) 1976 row_maxima=(double) p[i]; 1977 } 1978 } 1979 p+=GetPixelChannels(image); 1980 } 1981 #if defined(MAGICKCORE_OPENMP_SUPPORT) 1982 #pragma omp critical (MagickCore_GetImageRange) 1983 #endif 1984 { 1985 if (initialize != MagickFalse) 1986 { 1987 *minima=row_minima; 1988 *maxima=row_maxima; 1989 initialize=MagickFalse; 1990 } 1991 else 1992 { 1993 if (row_minima < *minima) 1994 *minima=row_minima; 1995 if (row_maxima > *maxima) 1996 *maxima=row_maxima; 1997 } 1998 } 1999 } 2000 image_view=DestroyCacheView(image_view); 2001 return(status); 2002 } 2003 2004 /* 2006 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2007 % % 2008 % % 2009 % % 2010 % G e t I m a g e S t a t i s t i c s % 2011 % % 2012 % % 2013 % % 2014 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2015 % 2016 % GetImageStatistics() returns statistics for each channel in the image. The 2017 % statistics include the channel depth, its minima, maxima, mean, standard 2018 % deviation, kurtosis and skewness. You can access the red channel mean, for 2019 % example, like this: 2020 % 2021 % channel_statistics=GetImageStatistics(image,exception); 2022 % red_mean=channel_statistics[RedPixelChannel].mean; 2023 % 2024 % Use MagickRelinquishMemory() to free the statistics buffer. 2025 % 2026 % The format of the GetImageStatistics method is: 2027 % 2028 % ChannelStatistics *GetImageStatistics(const Image *image, 2029 % ExceptionInfo *exception) 2030 % 2031 % A description of each parameter follows: 2032 % 2033 % o image: the image. 2034 % 2035 % o exception: return any errors or warnings in this structure. 2036 % 2037 */ 2038 MagickExport ChannelStatistics *GetImageStatistics(const Image *image, 2039 ExceptionInfo *exception) 2040 { 2041 ChannelStatistics 2042 *channel_statistics; 2043 2044 double 2045 *histogram; 2046 2047 MagickStatusType 2048 status; 2049 2050 QuantumAny 2051 range; 2052 2053 register ssize_t 2054 i; 2055 2056 size_t 2057 channels, 2058 depth; 2059 2060 ssize_t 2061 y; 2062 2063 assert(image != (Image *) NULL); 2064 assert(image->signature == MagickCoreSignature); 2065 if (image->debug != MagickFalse) 2066 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 2067 histogram=(double *) AcquireQuantumMemory(MaxMap+1UL,MaxPixelChannels* 2068 sizeof(*histogram)); 2069 channel_statistics=(ChannelStatistics *) AcquireQuantumMemory( 2070 MaxPixelChannels+1,sizeof(*channel_statistics)); 2071 if ((channel_statistics == (ChannelStatistics *) NULL) || 2072 (histogram == (double *) NULL)) 2073 { 2074 if (histogram != (double *) NULL) 2075 histogram=(double *) RelinquishMagickMemory(histogram); 2076 if (channel_statistics != (ChannelStatistics *) NULL) 2077 channel_statistics=(ChannelStatistics *) RelinquishMagickMemory( 2078 channel_statistics); 2079 return(channel_statistics); 2080 } 2081 (void) ResetMagickMemory(channel_statistics,0,(MaxPixelChannels+1)* 2082 sizeof(*channel_statistics)); 2083 for (i=0; i <= (ssize_t) MaxPixelChannels; i++) 2084 { 2085 channel_statistics[i].depth=1; 2086 channel_statistics[i].maxima=(-MagickMaximumValue); 2087 channel_statistics[i].minima=MagickMaximumValue; 2088 } 2089 (void) ResetMagickMemory(histogram,0,(MaxMap+1)*MaxPixelChannels* 2090 sizeof(*histogram)); 2091 for (y=0; y < (ssize_t) image->rows; y++) 2092 { 2093 register const Quantum 2094 *magick_restrict p; 2095 2096 register ssize_t 2097 x; 2098 2099 p=GetVirtualPixels(image,0,y,image->columns,1,exception); 2100 if (p == (const Quantum *) NULL) 2101 break; 2102 for (x=0; x < (ssize_t) image->columns; x++) 2103 { 2104 register ssize_t 2105 i; 2106 2107 if (GetPixelReadMask(image,p) == 0) 2108 { 2109 p+=GetPixelChannels(image); 2110 continue; 2111 } 2112 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 2113 { 2114 PixelChannel channel=GetPixelChannelChannel(image,i); 2115 PixelTrait traits=GetPixelChannelTraits(image,channel); 2116 if (traits == UndefinedPixelTrait) 2117 continue; 2118 if (channel_statistics[channel].depth != MAGICKCORE_QUANTUM_DEPTH) 2119 { 2120 depth=channel_statistics[channel].depth; 2121 range=GetQuantumRange(depth); 2122 status=p[i] != ScaleAnyToQuantum(ScaleQuantumToAny(p[i],range), 2123 range) ? MagickTrue : MagickFalse; 2124 if (status != MagickFalse) 2125 { 2126 channel_statistics[channel].depth++; 2127 i--; 2128 continue; 2129 } 2130 } 2131 if ((double) p[i] < channel_statistics[channel].minima) 2132 channel_statistics[channel].minima=(double) p[i]; 2133 if ((double) p[i] > channel_statistics[channel].maxima) 2134 channel_statistics[channel].maxima=(double) p[i]; 2135 channel_statistics[channel].sum+=p[i]; 2136 channel_statistics[channel].sum_squared+=(double) p[i]*p[i]; 2137 channel_statistics[channel].sum_cubed+=(double) p[i]*p[i]*p[i]; 2138 channel_statistics[channel].sum_fourth_power+=(double) p[i]*p[i]*p[i]* 2139 p[i]; 2140 channel_statistics[channel].area++; 2141 histogram[GetPixelChannels(image)*ScaleQuantumToMap( 2142 ClampToQuantum((double) p[i]))+i]++; 2143 } 2144 p+=GetPixelChannels(image); 2145 } 2146 } 2147 for (i=0; i < (ssize_t) MaxPixelChannels; i++) 2148 { 2149 double 2150 area, 2151 number_bins; 2152 2153 register ssize_t 2154 j; 2155 2156 area=PerceptibleReciprocal(channel_statistics[i].area); 2157 channel_statistics[i].sum*=area; 2158 channel_statistics[i].sum_squared*=area; 2159 channel_statistics[i].sum_cubed*=area; 2160 channel_statistics[i].sum_fourth_power*=area; 2161 channel_statistics[i].mean=channel_statistics[i].sum; 2162 channel_statistics[i].variance=channel_statistics[i].sum_squared; 2163 channel_statistics[i].standard_deviation=sqrt( 2164 channel_statistics[i].variance-(channel_statistics[i].mean* 2165 channel_statistics[i].mean)); 2166 number_bins=0.0; 2167 for (j=0; j < (ssize_t) (MaxMap+1U); j++) 2168 if (histogram[GetPixelChannels(image)*j+i] > 0.0) 2169 number_bins++; 2170 for (j=0; j < (ssize_t) (MaxMap+1U); j++) 2171 { 2172 double 2173 count; 2174 2175 count=histogram[GetPixelChannels(image)*j+i]*area; 2176 channel_statistics[i].entropy+=-count*MagickLog10(count)/ 2177 MagickLog10(number_bins); 2178 } 2179 } 2180 for (i=0; i < (ssize_t) MaxPixelChannels; i++) 2181 { 2182 PixelTrait traits=GetPixelChannelTraits(image,(PixelChannel) i); 2183 if ((traits & UpdatePixelTrait) == 0) 2184 continue; 2185 channel_statistics[CompositePixelChannel].area+=channel_statistics[i].area; 2186 channel_statistics[CompositePixelChannel].minima=MagickMin( 2187 channel_statistics[CompositePixelChannel].minima, 2188 channel_statistics[i].minima); 2189 channel_statistics[CompositePixelChannel].maxima=EvaluateMax( 2190 channel_statistics[CompositePixelChannel].maxima, 2191 channel_statistics[i].maxima); 2192 channel_statistics[CompositePixelChannel].sum+=channel_statistics[i].sum; 2193 channel_statistics[CompositePixelChannel].sum_squared+= 2194 channel_statistics[i].sum_squared; 2195 channel_statistics[CompositePixelChannel].sum_cubed+= 2196 channel_statistics[i].sum_cubed; 2197 channel_statistics[CompositePixelChannel].sum_fourth_power+= 2198 channel_statistics[i].sum_fourth_power; 2199 channel_statistics[CompositePixelChannel].mean+=channel_statistics[i].mean; 2200 channel_statistics[CompositePixelChannel].variance+= 2201 channel_statistics[i].variance-channel_statistics[i].mean* 2202 channel_statistics[i].mean; 2203 channel_statistics[CompositePixelChannel].standard_deviation+= 2204 channel_statistics[i].variance-channel_statistics[i].mean* 2205 channel_statistics[i].mean; 2206 if (channel_statistics[i].entropy > MagickEpsilon) 2207 channel_statistics[CompositePixelChannel].entropy+= 2208 channel_statistics[i].entropy; 2209 } 2210 channels=GetImageChannels(image); 2211 channel_statistics[CompositePixelChannel].area/=channels; 2212 channel_statistics[CompositePixelChannel].sum/=channels; 2213 channel_statistics[CompositePixelChannel].sum_squared/=channels; 2214 channel_statistics[CompositePixelChannel].sum_cubed/=channels; 2215 channel_statistics[CompositePixelChannel].sum_fourth_power/=channels; 2216 channel_statistics[CompositePixelChannel].mean/=channels; 2217 channel_statistics[CompositePixelChannel].variance/=channels; 2218 channel_statistics[CompositePixelChannel].standard_deviation= 2219 sqrt(channel_statistics[CompositePixelChannel].standard_deviation/channels); 2220 channel_statistics[CompositePixelChannel].kurtosis/=channels; 2221 channel_statistics[CompositePixelChannel].skewness/=channels; 2222 channel_statistics[CompositePixelChannel].entropy/=channels; 2223 for (i=0; i <= (ssize_t) MaxPixelChannels; i++) 2224 { 2225 double 2226 standard_deviation; 2227 2228 if (channel_statistics[i].standard_deviation == 0.0) 2229 continue; 2230 standard_deviation=PerceptibleReciprocal( 2231 channel_statistics[i].standard_deviation); 2232 channel_statistics[i].skewness=(channel_statistics[i].sum_cubed-3.0* 2233 channel_statistics[i].mean*channel_statistics[i].sum_squared+2.0* 2234 channel_statistics[i].mean*channel_statistics[i].mean* 2235 channel_statistics[i].mean)*(standard_deviation*standard_deviation* 2236 standard_deviation); 2237 channel_statistics[i].kurtosis=(channel_statistics[i].sum_fourth_power-4.0* 2238 channel_statistics[i].mean*channel_statistics[i].sum_cubed+6.0* 2239 channel_statistics[i].mean*channel_statistics[i].mean* 2240 channel_statistics[i].sum_squared-3.0*channel_statistics[i].mean* 2241 channel_statistics[i].mean*1.0*channel_statistics[i].mean* 2242 channel_statistics[i].mean)*(standard_deviation*standard_deviation* 2243 standard_deviation*standard_deviation)-3.0; 2244 } 2245 histogram=(double *) RelinquishMagickMemory(histogram); 2246 if (y < (ssize_t) image->rows) 2247 channel_statistics=(ChannelStatistics *) RelinquishMagickMemory( 2248 channel_statistics); 2249 return(channel_statistics); 2250 } 2251 2252 /* 2254 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2255 % % 2256 % % 2257 % % 2258 % P o l y n o m i a l I m a g e % 2259 % % 2260 % % 2261 % % 2262 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2263 % 2264 % PolynomialImage() returns a new image where each pixel is the sum of the 2265 % pixels in the image sequence after applying its corresponding terms 2266 % (coefficient and degree pairs). 2267 % 2268 % The format of the PolynomialImage method is: 2269 % 2270 % Image *PolynomialImage(const Image *images,const size_t number_terms, 2271 % const double *terms,ExceptionInfo *exception) 2272 % 2273 % A description of each parameter follows: 2274 % 2275 % o images: the image sequence. 2276 % 2277 % o number_terms: the number of terms in the list. The actual list length 2278 % is 2 x number_terms + 1 (the constant). 2279 % 2280 % o terms: the list of polynomial coefficients and degree pairs and a 2281 % constant. 2282 % 2283 % o exception: return any errors or warnings in this structure. 2284 % 2285 */ 2286 2287 MagickExport Image *PolynomialImage(const Image *images, 2288 const size_t number_terms,const double *terms,ExceptionInfo *exception) 2289 { 2290 #define PolynomialImageTag "Polynomial/Image" 2291 2292 CacheView 2293 *polynomial_view; 2294 2295 Image 2296 *image; 2297 2298 MagickBooleanType 2299 status; 2300 2301 MagickOffsetType 2302 progress; 2303 2304 PixelChannels 2305 **magick_restrict polynomial_pixels; 2306 2307 size_t 2308 number_images; 2309 2310 ssize_t 2311 y; 2312 2313 assert(images != (Image *) NULL); 2314 assert(images->signature == MagickCoreSignature); 2315 if (images->debug != MagickFalse) 2316 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename); 2317 assert(exception != (ExceptionInfo *) NULL); 2318 assert(exception->signature == MagickCoreSignature); 2319 image=CloneImage(images,images->columns,images->rows,MagickTrue, 2320 exception); 2321 if (image == (Image *) NULL) 2322 return((Image *) NULL); 2323 if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse) 2324 { 2325 image=DestroyImage(image); 2326 return((Image *) NULL); 2327 } 2328 number_images=GetImageListLength(images); 2329 polynomial_pixels=AcquirePixelThreadSet(images); 2330 if (polynomial_pixels == (PixelChannels **) NULL) 2331 { 2332 image=DestroyImage(image); 2333 (void) ThrowMagickException(exception,GetMagickModule(), 2334 ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename); 2335 return((Image *) NULL); 2336 } 2337 /* 2338 Polynomial image pixels. 2339 */ 2340 status=MagickTrue; 2341 progress=0; 2342 polynomial_view=AcquireAuthenticCacheView(image,exception); 2343 #if defined(MAGICKCORE_OPENMP_SUPPORT) 2344 #pragma omp parallel for schedule(static,4) shared(progress,status) \ 2345 magick_threads(image,image,image->rows,1) 2346 #endif 2347 for (y=0; y < (ssize_t) image->rows; y++) 2348 { 2349 CacheView 2350 *image_view; 2351 2352 const Image 2353 *next; 2354 2355 const int 2356 id = GetOpenMPThreadId(); 2357 2358 register ssize_t 2359 i, 2360 x; 2361 2362 register PixelChannels 2363 *polynomial_pixel; 2364 2365 register Quantum 2366 *magick_restrict q; 2367 2368 ssize_t 2369 j; 2370 2371 if (status == MagickFalse) 2372 continue; 2373 q=QueueCacheViewAuthenticPixels(polynomial_view,0,y,image->columns,1, 2374 exception); 2375 if (q == (Quantum *) NULL) 2376 { 2377 status=MagickFalse; 2378 continue; 2379 } 2380 polynomial_pixel=polynomial_pixels[id]; 2381 for (j=0; j < (ssize_t) image->columns; j++) 2382 for (i=0; i < MaxPixelChannels; i++) 2383 polynomial_pixel[j].channel[i]=0.0; 2384 next=images; 2385 for (j=0; j < (ssize_t) number_images; j++) 2386 { 2387 register const Quantum 2388 *p; 2389 2390 if (j >= (ssize_t) number_terms) 2391 continue; 2392 image_view=AcquireVirtualCacheView(next,exception); 2393 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception); 2394 if (p == (const Quantum *) NULL) 2395 { 2396 image_view=DestroyCacheView(image_view); 2397 break; 2398 } 2399 for (x=0; x < (ssize_t) image->columns; x++) 2400 { 2401 register ssize_t 2402 i; 2403 2404 if (GetPixelReadMask(next,p) == 0) 2405 { 2406 p+=GetPixelChannels(next); 2407 continue; 2408 } 2409 for (i=0; i < (ssize_t) GetPixelChannels(next); i++) 2410 { 2411 MagickRealType 2412 coefficient, 2413 degree; 2414 2415 PixelChannel channel=GetPixelChannelChannel(image,i); 2416 PixelTrait traits=GetPixelChannelTraits(next,channel); 2417 PixelTrait polynomial_traits=GetPixelChannelTraits(image,channel); 2418 if ((traits == UndefinedPixelTrait) || 2419 (polynomial_traits == UndefinedPixelTrait)) 2420 continue; 2421 if ((traits & UpdatePixelTrait) == 0) 2422 continue; 2423 coefficient=(MagickRealType) terms[2*j]; 2424 degree=(MagickRealType) terms[(j << 1)+1]; 2425 polynomial_pixel[x].channel[i]+=coefficient* 2426 pow(QuantumScale*GetPixelChannel(image,channel,p),degree); 2427 } 2428 p+=GetPixelChannels(next); 2429 } 2430 image_view=DestroyCacheView(image_view); 2431 next=GetNextImageInList(next); 2432 } 2433 for (x=0; x < (ssize_t) image->columns; x++) 2434 { 2435 register ssize_t 2436 i; 2437 2438 if (GetPixelReadMask(image,q) == 0) 2439 { 2440 q+=GetPixelChannels(image); 2441 continue; 2442 } 2443 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 2444 { 2445 PixelChannel channel=GetPixelChannelChannel(image,i); 2446 PixelTrait traits=GetPixelChannelTraits(image,channel); 2447 if (traits == UndefinedPixelTrait) 2448 continue; 2449 if ((traits & UpdatePixelTrait) == 0) 2450 continue; 2451 q[i]=ClampToQuantum(QuantumRange*polynomial_pixel[x].channel[i]); 2452 } 2453 q+=GetPixelChannels(image); 2454 } 2455 if (SyncCacheViewAuthenticPixels(polynomial_view,exception) == MagickFalse) 2456 status=MagickFalse; 2457 if (images->progress_monitor != (MagickProgressMonitor) NULL) 2458 { 2459 MagickBooleanType 2460 proceed; 2461 2462 #if defined(MAGICKCORE_OPENMP_SUPPORT) 2463 #pragma omp critical (MagickCore_PolynomialImages) 2464 #endif 2465 proceed=SetImageProgress(images,PolynomialImageTag,progress++, 2466 image->rows); 2467 if (proceed == MagickFalse) 2468 status=MagickFalse; 2469 } 2470 } 2471 polynomial_view=DestroyCacheView(polynomial_view); 2472 polynomial_pixels=DestroyPixelThreadSet(polynomial_pixels); 2473 if (status == MagickFalse) 2474 image=DestroyImage(image); 2475 return(image); 2476 } 2477 2478 /* 2480 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2481 % % 2482 % % 2483 % % 2484 % S t a t i s t i c I m a g e % 2485 % % 2486 % % 2487 % % 2488 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2489 % 2490 % StatisticImage() makes each pixel the min / max / median / mode / etc. of 2491 % the neighborhood of the specified width and height. 2492 % 2493 % The format of the StatisticImage method is: 2494 % 2495 % Image *StatisticImage(const Image *image,const StatisticType type, 2496 % const size_t width,const size_t height,ExceptionInfo *exception) 2497 % 2498 % A description of each parameter follows: 2499 % 2500 % o image: the image. 2501 % 2502 % o type: the statistic type (median, mode, etc.). 2503 % 2504 % o width: the width of the pixel neighborhood. 2505 % 2506 % o height: the height of the pixel neighborhood. 2507 % 2508 % o exception: return any errors or warnings in this structure. 2509 % 2510 */ 2511 2512 typedef struct _SkipNode 2513 { 2514 size_t 2515 next[9], 2516 count, 2517 signature; 2518 } SkipNode; 2519 2520 typedef struct _SkipList 2521 { 2522 ssize_t 2523 level; 2524 2525 SkipNode 2526 *nodes; 2527 } SkipList; 2528 2529 typedef struct _PixelList 2530 { 2531 size_t 2532 length, 2533 seed; 2534 2535 SkipList 2536 skip_list; 2537 2538 size_t 2539 signature; 2540 } PixelList; 2541 2542 static PixelList *DestroyPixelList(PixelList *pixel_list) 2543 { 2544 if (pixel_list == (PixelList *) NULL) 2545 return((PixelList *) NULL); 2546 if (pixel_list->skip_list.nodes != (SkipNode *) NULL) 2547 pixel_list->skip_list.nodes=(SkipNode *) RelinquishMagickMemory( 2548 pixel_list->skip_list.nodes); 2549 pixel_list=(PixelList *) RelinquishMagickMemory(pixel_list); 2550 return(pixel_list); 2551 } 2552 2553 static PixelList **DestroyPixelListThreadSet(PixelList **pixel_list) 2554 { 2555 register ssize_t 2556 i; 2557 2558 assert(pixel_list != (PixelList **) NULL); 2559 for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++) 2560 if (pixel_list[i] != (PixelList *) NULL) 2561 pixel_list[i]=DestroyPixelList(pixel_list[i]); 2562 pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list); 2563 return(pixel_list); 2564 } 2565 2566 static PixelList *AcquirePixelList(const size_t width,const size_t height) 2567 { 2568 PixelList 2569 *pixel_list; 2570 2571 pixel_list=(PixelList *) AcquireMagickMemory(sizeof(*pixel_list)); 2572 if (pixel_list == (PixelList *) NULL) 2573 return(pixel_list); 2574 (void) ResetMagickMemory((void *) pixel_list,0,sizeof(*pixel_list)); 2575 pixel_list->length=width*height; 2576 pixel_list->skip_list.nodes=(SkipNode *) AcquireQuantumMemory(65537UL, 2577 sizeof(*pixel_list->skip_list.nodes)); 2578 if (pixel_list->skip_list.nodes == (SkipNode *) NULL) 2579 return(DestroyPixelList(pixel_list)); 2580 (void) ResetMagickMemory(pixel_list->skip_list.nodes,0,65537UL* 2581 sizeof(*pixel_list->skip_list.nodes)); 2582 pixel_list->signature=MagickCoreSignature; 2583 return(pixel_list); 2584 } 2585 2586 static PixelList **AcquirePixelListThreadSet(const size_t width, 2587 const size_t height) 2588 { 2589 PixelList 2590 **pixel_list; 2591 2592 register ssize_t 2593 i; 2594 2595 size_t 2596 number_threads; 2597 2598 number_threads=(size_t) GetMagickResourceLimit(ThreadResource); 2599 pixel_list=(PixelList **) AcquireQuantumMemory(number_threads, 2600 sizeof(*pixel_list)); 2601 if (pixel_list == (PixelList **) NULL) 2602 return((PixelList **) NULL); 2603 (void) ResetMagickMemory(pixel_list,0,number_threads*sizeof(*pixel_list)); 2604 for (i=0; i < (ssize_t) number_threads; i++) 2605 { 2606 pixel_list[i]=AcquirePixelList(width,height); 2607 if (pixel_list[i] == (PixelList *) NULL) 2608 return(DestroyPixelListThreadSet(pixel_list)); 2609 } 2610 return(pixel_list); 2611 } 2612 2613 static void AddNodePixelList(PixelList *pixel_list,const size_t color) 2614 { 2615 register SkipList 2616 *p; 2617 2618 register ssize_t 2619 level; 2620 2621 size_t 2622 search, 2623 update[9]; 2624 2625 /* 2626 Initialize the node. 2627 */ 2628 p=(&pixel_list->skip_list); 2629 p->nodes[color].signature=pixel_list->signature; 2630 p->nodes[color].count=1; 2631 /* 2632 Determine where it belongs in the list. 2633 */ 2634 search=65536UL; 2635 for (level=p->level; level >= 0; level--) 2636 { 2637 while (p->nodes[search].next[level] < color) 2638 search=p->nodes[search].next[level]; 2639 update[level]=search; 2640 } 2641 /* 2642 Generate a pseudo-random level for this node. 2643 */ 2644 for (level=0; ; level++) 2645 { 2646 pixel_list->seed=(pixel_list->seed*42893621L)+1L; 2647 if ((pixel_list->seed & 0x300) != 0x300) 2648 break; 2649 } 2650 if (level > 8) 2651 level=8; 2652 if (level > (p->level+2)) 2653 level=p->level+2; 2654 /* 2655 If we're raising the list's level, link back to the root node. 2656 */ 2657 while (level > p->level) 2658 { 2659 p->level++; 2660 update[p->level]=65536UL; 2661 } 2662 /* 2663 Link the node into the skip-list. 2664 */ 2665 do 2666 { 2667 p->nodes[color].next[level]=p->nodes[update[level]].next[level]; 2668 p->nodes[update[level]].next[level]=color; 2669 } while (level-- > 0); 2670 } 2671 2672 static inline void GetMaximumPixelList(PixelList *pixel_list,Quantum *pixel) 2673 { 2674 register SkipList 2675 *p; 2676 2677 size_t 2678 color, 2679 maximum; 2680 2681 ssize_t 2682 count; 2683 2684 /* 2685 Find the maximum value for each of the color. 2686 */ 2687 p=(&pixel_list->skip_list); 2688 color=65536L; 2689 count=0; 2690 maximum=p->nodes[color].next[0]; 2691 do 2692 { 2693 color=p->nodes[color].next[0]; 2694 if (color > maximum) 2695 maximum=color; 2696 count+=p->nodes[color].count; 2697 } while (count < (ssize_t) pixel_list->length); 2698 *pixel=ScaleShortToQuantum((unsigned short) maximum); 2699 } 2700 2701 static inline void GetMeanPixelList(PixelList *pixel_list,Quantum *pixel) 2702 { 2703 double 2704 sum; 2705 2706 register SkipList 2707 *p; 2708 2709 size_t 2710 color; 2711 2712 ssize_t 2713 count; 2714 2715 /* 2716 Find the mean value for each of the color. 2717 */ 2718 p=(&pixel_list->skip_list); 2719 color=65536L; 2720 count=0; 2721 sum=0.0; 2722 do 2723 { 2724 color=p->nodes[color].next[0]; 2725 sum+=(double) p->nodes[color].count*color; 2726 count+=p->nodes[color].count; 2727 } while (count < (ssize_t) pixel_list->length); 2728 sum/=pixel_list->length; 2729 *pixel=ScaleShortToQuantum((unsigned short) sum); 2730 } 2731 2732 static inline void GetMedianPixelList(PixelList *pixel_list,Quantum *pixel) 2733 { 2734 register SkipList 2735 *p; 2736 2737 size_t 2738 color; 2739 2740 ssize_t 2741 count; 2742 2743 /* 2744 Find the median value for each of the color. 2745 */ 2746 p=(&pixel_list->skip_list); 2747 color=65536L; 2748 count=0; 2749 do 2750 { 2751 color=p->nodes[color].next[0]; 2752 count+=p->nodes[color].count; 2753 } while (count <= (ssize_t) (pixel_list->length >> 1)); 2754 *pixel=ScaleShortToQuantum((unsigned short) color); 2755 } 2756 2757 static inline void GetMinimumPixelList(PixelList *pixel_list,Quantum *pixel) 2758 { 2759 register SkipList 2760 *p; 2761 2762 size_t 2763 color, 2764 minimum; 2765 2766 ssize_t 2767 count; 2768 2769 /* 2770 Find the minimum value for each of the color. 2771 */ 2772 p=(&pixel_list->skip_list); 2773 count=0; 2774 color=65536UL; 2775 minimum=p->nodes[color].next[0]; 2776 do 2777 { 2778 color=p->nodes[color].next[0]; 2779 if (color < minimum) 2780 minimum=color; 2781 count+=p->nodes[color].count; 2782 } while (count < (ssize_t) pixel_list->length); 2783 *pixel=ScaleShortToQuantum((unsigned short) minimum); 2784 } 2785 2786 static inline void GetModePixelList(PixelList *pixel_list,Quantum *pixel) 2787 { 2788 register SkipList 2789 *p; 2790 2791 size_t 2792 color, 2793 max_count, 2794 mode; 2795 2796 ssize_t 2797 count; 2798 2799 /* 2800 Make each pixel the 'predominant color' of the specified neighborhood. 2801 */ 2802 p=(&pixel_list->skip_list); 2803 color=65536L; 2804 mode=color; 2805 max_count=p->nodes[mode].count; 2806 count=0; 2807 do 2808 { 2809 color=p->nodes[color].next[0]; 2810 if (p->nodes[color].count > max_count) 2811 { 2812 mode=color; 2813 max_count=p->nodes[mode].count; 2814 } 2815 count+=p->nodes[color].count; 2816 } while (count < (ssize_t) pixel_list->length); 2817 *pixel=ScaleShortToQuantum((unsigned short) mode); 2818 } 2819 2820 static inline void GetNonpeakPixelList(PixelList *pixel_list,Quantum *pixel) 2821 { 2822 register SkipList 2823 *p; 2824 2825 size_t 2826 color, 2827 next, 2828 previous; 2829 2830 ssize_t 2831 count; 2832 2833 /* 2834 Finds the non peak value for each of the colors. 2835 */ 2836 p=(&pixel_list->skip_list); 2837 color=65536L; 2838 next=p->nodes[color].next[0]; 2839 count=0; 2840 do 2841 { 2842 previous=color; 2843 color=next; 2844 next=p->nodes[color].next[0]; 2845 count+=p->nodes[color].count; 2846 } while (count <= (ssize_t) (pixel_list->length >> 1)); 2847 if ((previous == 65536UL) && (next != 65536UL)) 2848 color=next; 2849 else 2850 if ((previous != 65536UL) && (next == 65536UL)) 2851 color=previous; 2852 *pixel=ScaleShortToQuantum((unsigned short) color); 2853 } 2854 2855 static inline void GetRootMeanSquarePixelList(PixelList *pixel_list, 2856 Quantum *pixel) 2857 { 2858 double 2859 sum; 2860 2861 register SkipList 2862 *p; 2863 2864 size_t 2865 color; 2866 2867 ssize_t 2868 count; 2869 2870 /* 2871 Find the root mean square value for each of the color. 2872 */ 2873 p=(&pixel_list->skip_list); 2874 color=65536L; 2875 count=0; 2876 sum=0.0; 2877 do 2878 { 2879 color=p->nodes[color].next[0]; 2880 sum+=(double) (p->nodes[color].count*color*color); 2881 count+=p->nodes[color].count; 2882 } while (count < (ssize_t) pixel_list->length); 2883 sum/=pixel_list->length; 2884 *pixel=ScaleShortToQuantum((unsigned short) sqrt(sum)); 2885 } 2886 2887 static inline void GetStandardDeviationPixelList(PixelList *pixel_list, 2888 Quantum *pixel) 2889 { 2890 double 2891 sum, 2892 sum_squared; 2893 2894 register SkipList 2895 *p; 2896 2897 size_t 2898 color; 2899 2900 ssize_t 2901 count; 2902 2903 /* 2904 Find the standard-deviation value for each of the color. 2905 */ 2906 p=(&pixel_list->skip_list); 2907 color=65536L; 2908 count=0; 2909 sum=0.0; 2910 sum_squared=0.0; 2911 do 2912 { 2913 register ssize_t 2914 i; 2915 2916 color=p->nodes[color].next[0]; 2917 sum+=(double) p->nodes[color].count*color; 2918 for (i=0; i < (ssize_t) p->nodes[color].count; i++) 2919 sum_squared+=((double) color)*((double) color); 2920 count+=p->nodes[color].count; 2921 } while (count < (ssize_t) pixel_list->length); 2922 sum/=pixel_list->length; 2923 sum_squared/=pixel_list->length; 2924 *pixel=ScaleShortToQuantum((unsigned short) sqrt(sum_squared-(sum*sum))); 2925 } 2926 2927 static inline void InsertPixelList(const Quantum pixel,PixelList *pixel_list) 2928 { 2929 size_t 2930 signature; 2931 2932 unsigned short 2933 index; 2934 2935 index=ScaleQuantumToShort(pixel); 2936 signature=pixel_list->skip_list.nodes[index].signature; 2937 if (signature == pixel_list->signature) 2938 { 2939 pixel_list->skip_list.nodes[index].count++; 2940 return; 2941 } 2942 AddNodePixelList(pixel_list,index); 2943 } 2944 2945 static void ResetPixelList(PixelList *pixel_list) 2946 { 2947 int 2948 level; 2949 2950 register SkipNode 2951 *root; 2952 2953 register SkipList 2954 *p; 2955 2956 /* 2957 Reset the skip-list. 2958 */ 2959 p=(&pixel_list->skip_list); 2960 root=p->nodes+65536UL; 2961 p->level=0; 2962 for (level=0; level < 9; level++) 2963 root->next[level]=65536UL; 2964 pixel_list->seed=pixel_list->signature++; 2965 } 2966 2967 MagickExport Image *StatisticImage(const Image *image,const StatisticType type, 2968 const size_t width,const size_t height,ExceptionInfo *exception) 2969 { 2970 #define StatisticImageTag "Statistic/Image" 2971 2972 CacheView 2973 *image_view, 2974 *statistic_view; 2975 2976 Image 2977 *statistic_image; 2978 2979 MagickBooleanType 2980 status; 2981 2982 MagickOffsetType 2983 progress; 2984 2985 PixelList 2986 **magick_restrict pixel_list; 2987 2988 ssize_t 2989 center, 2990 y; 2991 2992 /* 2993 Initialize statistics image attributes. 2994 */ 2995 assert(image != (Image *) NULL); 2996 assert(image->signature == MagickCoreSignature); 2997 if (image->debug != MagickFalse) 2998 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 2999 assert(exception != (ExceptionInfo *) NULL); 3000 assert(exception->signature == MagickCoreSignature); 3001 statistic_image=CloneImage(image,image->columns,image->rows,MagickTrue, 3002 exception); 3003 if (statistic_image == (Image *) NULL) 3004 return((Image *) NULL); 3005 status=SetImageStorageClass(statistic_image,DirectClass,exception); 3006 if (status == MagickFalse) 3007 { 3008 statistic_image=DestroyImage(statistic_image); 3009 return((Image *) NULL); 3010 } 3011 pixel_list=AcquirePixelListThreadSet(MagickMax(width,1),MagickMax(height,1)); 3012 if (pixel_list == (PixelList **) NULL) 3013 { 3014 statistic_image=DestroyImage(statistic_image); 3015 ThrowImageException(ResourceLimitError,"MemoryAllocationFailed"); 3016 } 3017 /* 3018 Make each pixel the min / max / median / mode / etc. of the neighborhood. 3019 */ 3020 center=(ssize_t) GetPixelChannels(image)*(image->columns+MagickMax(width,1))* 3021 (MagickMax(height,1)/2L)+GetPixelChannels(image)*(MagickMax(width,1)/2L); 3022 status=MagickTrue; 3023 progress=0; 3024 image_view=AcquireVirtualCacheView(image,exception); 3025 statistic_view=AcquireAuthenticCacheView(statistic_image,exception); 3026 #if defined(MAGICKCORE_OPENMP_SUPPORT) 3027 #pragma omp parallel for schedule(static,4) shared(progress,status) \ 3028 magick_threads(image,statistic_image,statistic_image->rows,1) 3029 #endif 3030 for (y=0; y < (ssize_t) statistic_image->rows; y++) 3031 { 3032 const int 3033 id = GetOpenMPThreadId(); 3034 3035 register const Quantum 3036 *magick_restrict p; 3037 3038 register Quantum 3039 *magick_restrict q; 3040 3041 register ssize_t 3042 x; 3043 3044 if (status == MagickFalse) 3045 continue; 3046 p=GetCacheViewVirtualPixels(image_view,-((ssize_t) MagickMax(width,1)/2L),y- 3047 (ssize_t) (MagickMax(height,1)/2L),image->columns+MagickMax(width,1), 3048 MagickMax(height,1),exception); 3049 q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns, 1,exception); 3050 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL)) 3051 { 3052 status=MagickFalse; 3053 continue; 3054 } 3055 for (x=0; x < (ssize_t) statistic_image->columns; x++) 3056 { 3057 register ssize_t 3058 i; 3059 3060 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 3061 { 3062 Quantum 3063 pixel; 3064 3065 register const Quantum 3066 *magick_restrict pixels; 3067 3068 register ssize_t 3069 u; 3070 3071 ssize_t 3072 v; 3073 3074 PixelChannel channel=GetPixelChannelChannel(image,i); 3075 PixelTrait traits=GetPixelChannelTraits(image,channel); 3076 PixelTrait statistic_traits=GetPixelChannelTraits(statistic_image, 3077 channel); 3078 if ((traits == UndefinedPixelTrait) || 3079 (statistic_traits == UndefinedPixelTrait)) 3080 continue; 3081 if (((statistic_traits & CopyPixelTrait) != 0) || 3082 (GetPixelReadMask(image,p) == 0)) 3083 { 3084 SetPixelChannel(statistic_image,channel,p[center+i],q); 3085 continue; 3086 } 3087 pixels=p; 3088 ResetPixelList(pixel_list[id]); 3089 for (v=0; v < (ssize_t) MagickMax(height,1); v++) 3090 { 3091 for (u=0; u < (ssize_t) MagickMax(width,1); u++) 3092 { 3093 InsertPixelList(pixels[i],pixel_list[id]); 3094 pixels+=GetPixelChannels(image); 3095 } 3096 pixels+=GetPixelChannels(image)*image->columns; 3097 } 3098 switch (type) 3099 { 3100 case GradientStatistic: 3101 { 3102 double 3103 maximum, 3104 minimum; 3105 3106 GetMinimumPixelList(pixel_list[id],&pixel); 3107 minimum=(double) pixel; 3108 GetMaximumPixelList(pixel_list[id],&pixel); 3109 maximum=(double) pixel; 3110 pixel=ClampToQuantum(MagickAbsoluteValue(maximum-minimum)); 3111 break; 3112 } 3113 case MaximumStatistic: 3114 { 3115 GetMaximumPixelList(pixel_list[id],&pixel); 3116 break; 3117 } 3118 case MeanStatistic: 3119 { 3120 GetMeanPixelList(pixel_list[id],&pixel); 3121 break; 3122 } 3123 case MedianStatistic: 3124 default: 3125 { 3126 GetMedianPixelList(pixel_list[id],&pixel); 3127 break; 3128 } 3129 case MinimumStatistic: 3130 { 3131 GetMinimumPixelList(pixel_list[id],&pixel); 3132 break; 3133 } 3134 case ModeStatistic: 3135 { 3136 GetModePixelList(pixel_list[id],&pixel); 3137 break; 3138 } 3139 case NonpeakStatistic: 3140 { 3141 GetNonpeakPixelList(pixel_list[id],&pixel); 3142 break; 3143 } 3144 case RootMeanSquareStatistic: 3145 { 3146 GetRootMeanSquarePixelList(pixel_list[id],&pixel); 3147 break; 3148 } 3149 case StandardDeviationStatistic: 3150 { 3151 GetStandardDeviationPixelList(pixel_list[id],&pixel); 3152 break; 3153 } 3154 } 3155 SetPixelChannel(statistic_image,channel,pixel,q); 3156 } 3157 p+=GetPixelChannels(image); 3158 q+=GetPixelChannels(statistic_image); 3159 } 3160 if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse) 3161 status=MagickFalse; 3162 if (image->progress_monitor != (MagickProgressMonitor) NULL) 3163 { 3164 MagickBooleanType 3165 proceed; 3166 3167 #if defined(MAGICKCORE_OPENMP_SUPPORT) 3168 #pragma omp critical (MagickCore_StatisticImage) 3169 #endif 3170 proceed=SetImageProgress(image,StatisticImageTag,progress++, 3171 image->rows); 3172 if (proceed == MagickFalse) 3173 status=MagickFalse; 3174 } 3175 } 3176 statistic_view=DestroyCacheView(statistic_view); 3177 image_view=DestroyCacheView(image_view); 3178 pixel_list=DestroyPixelListThreadSet(pixel_list); 3179 if (status == MagickFalse) 3180 statistic_image=DestroyImage(statistic_image); 3181 return(statistic_image); 3182 } 3183
