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