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