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