statistic.c revision aa8634f731b35f3924daf3aa1a8eb52d2cba0c44
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-2011 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/segment.h" 85#include "MagickCore/semaphore.h" 86#include "MagickCore/signature-private.h" 87#include "MagickCore/statistic.h" 88#include "MagickCore/string_.h" 89#include "MagickCore/thread-private.h" 90#include "MagickCore/timer.h" 91#include "MagickCore/utility.h" 92#include "MagickCore/version.h" 93 94/* 95%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 96% % 97% % 98% % 99% E v a l u a t e I m a g e % 100% % 101% % 102% % 103%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 104% 105% EvaluateImage() applies a value to the image with an arithmetic, relational, 106% or logical operator to an image. Use these operations to lighten or darken 107% an image, to increase or decrease contrast in an image, or to produce the 108% "negative" of an image. 109% 110% The format of the EvaluateImage method is: 111% 112% MagickBooleanType EvaluateImage(Image *image, 113% const MagickEvaluateOperator op,const double value, 114% ExceptionInfo *exception) 115% MagickBooleanType EvaluateImages(Image *images, 116% const MagickEvaluateOperator op,const double value, 117% ExceptionInfo *exception) 118% 119% A description of each parameter follows: 120% 121% o image: the image. 122% 123% o op: A channel op. 124% 125% o value: A value value. 126% 127% o exception: return any errors or warnings in this structure. 128% 129*/ 130 131typedef struct _WenusInfo 132{ 133 MagickRealType 134 channel[MaxPixelChannels]; 135} WenusInfo; 136 137static WenusInfo **DestroyPixelThreadSet(WenusInfo **pixels) 138{ 139 register ssize_t 140 i; 141 142 assert(pixels != (WenusInfo **) NULL); 143 for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++) 144 if (pixels[i] != (WenusInfo *) NULL) 145 pixels[i]=(WenusInfo *) RelinquishMagickMemory(pixels[i]); 146 pixels=(WenusInfo **) RelinquishMagickMemory(pixels); 147 return(pixels); 148} 149 150static WenusInfo **AcquirePixelThreadSet(const Image *image, 151 const size_t number_images) 152{ 153 register ssize_t 154 i; 155 156 WenusInfo 157 **pixels; 158 159 size_t 160 length, 161 number_threads; 162 163 number_threads=GetOpenMPMaximumThreads(); 164 pixels=(WenusInfo **) AcquireQuantumMemory(number_threads,sizeof(*pixels)); 165 if (pixels == (WenusInfo **) NULL) 166 return((WenusInfo **) NULL); 167 (void) ResetMagickMemory(pixels,0,number_threads*sizeof(*pixels)); 168 for (i=0; i < (ssize_t) number_threads; i++) 169 { 170 register ssize_t 171 j; 172 173 length=image->columns; 174 if (length < number_images) 175 length=number_images; 176 pixels[i]=(WenusInfo *) AcquireQuantumMemory(length,sizeof(**pixels)); 177 if (pixels[i] == (WenusInfo *) NULL) 178 return(DestroyPixelThreadSet(pixels)); 179 for (j=0; j < (ssize_t) length; j++) 180 { 181 register ssize_t 182 k; 183 184 for (k=0; k < MaxPixelChannels; k++) 185 pixels[i][j].channel[k]=0.0; 186 } 187 } 188 return(pixels); 189} 190 191static inline double MagickMax(const double x,const double y) 192{ 193 if (x > y) 194 return(x); 195 return(y); 196} 197 198#if defined(__cplusplus) || defined(c_plusplus) 199extern "C" { 200#endif 201 202static int IntensityCompare(const void *x,const void *y) 203{ 204 const WenusInfo 205 *color_1, 206 *color_2; 207 208 MagickRealType 209 distance; 210 211 register ssize_t 212 i; 213 214 color_1=(const WenusInfo *) x; 215 color_2=(const WenusInfo *) y; 216 distance=0.0; 217 for (i=0; i < MaxPixelChannels; i++) 218 distance+=color_1->channel[i]-(MagickRealType) color_2->channel[i]; 219 return(distance < 0 ? -1 : distance > 0 ? 1 : 0); 220} 221 222#if defined(__cplusplus) || defined(c_plusplus) 223} 224#endif 225 226static inline double MagickMin(const double x,const double y) 227{ 228 if (x < y) 229 return(x); 230 return(y); 231} 232 233static MagickRealType ApplyEvaluateOperator(RandomInfo *random_info, 234 Quantum pixel,const MagickEvaluateOperator op,const MagickRealType value) 235{ 236 MagickRealType 237 result; 238 239 result=0.0; 240 switch (op) 241 { 242 case UndefinedEvaluateOperator: 243 break; 244 case AbsEvaluateOperator: 245 { 246 result=(MagickRealType) fabs((double) (pixel+value)); 247 break; 248 } 249 case AddEvaluateOperator: 250 { 251 result=(MagickRealType) (pixel+value); 252 break; 253 } 254 case AddModulusEvaluateOperator: 255 { 256 /* 257 This returns a 'floored modulus' of the addition which is a 258 positive result. It differs from % or fmod() that returns a 259 'truncated modulus' result, where floor() is replaced by trunc() and 260 could return a negative result (which is clipped). 261 */ 262 result=pixel+value; 263 result-=(QuantumRange+1.0)*floor((double) result/(QuantumRange+1.0)); 264 break; 265 } 266 case AndEvaluateOperator: 267 { 268 result=(MagickRealType) ((size_t) pixel & (size_t) (value+0.5)); 269 break; 270 } 271 case CosineEvaluateOperator: 272 { 273 result=(MagickRealType) (QuantumRange*(0.5*cos((double) (2.0*MagickPI* 274 QuantumScale*pixel*value))+0.5)); 275 break; 276 } 277 case DivideEvaluateOperator: 278 { 279 result=pixel/(value == 0.0 ? 1.0 : value); 280 break; 281 } 282 case ExponentialEvaluateOperator: 283 { 284 result=(MagickRealType) (QuantumRange*exp((double) (value*QuantumScale* 285 pixel))); 286 break; 287 } 288 case GaussianNoiseEvaluateOperator: 289 { 290 result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel, 291 GaussianNoise,value); 292 break; 293 } 294 case ImpulseNoiseEvaluateOperator: 295 { 296 result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel, 297 ImpulseNoise,value); 298 break; 299 } 300 case LaplacianNoiseEvaluateOperator: 301 { 302 result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel, 303 LaplacianNoise,value); 304 break; 305 } 306 case LeftShiftEvaluateOperator: 307 { 308 result=(MagickRealType) ((size_t) pixel << (size_t) (value+0.5)); 309 break; 310 } 311 case LogEvaluateOperator: 312 { 313 result=(MagickRealType) (QuantumRange*log((double) (QuantumScale*value* 314 pixel+1.0))/log((double) (value+1.0))); 315 break; 316 } 317 case MaxEvaluateOperator: 318 { 319 result=(MagickRealType) MagickMax((double) pixel,value); 320 break; 321 } 322 case MeanEvaluateOperator: 323 { 324 result=(MagickRealType) (pixel+value); 325 break; 326 } 327 case MedianEvaluateOperator: 328 { 329 result=(MagickRealType) (pixel+value); 330 break; 331 } 332 case MinEvaluateOperator: 333 { 334 result=(MagickRealType) MagickMin((double) pixel,value); 335 break; 336 } 337 case MultiplicativeNoiseEvaluateOperator: 338 { 339 result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel, 340 MultiplicativeGaussianNoise,value); 341 break; 342 } 343 case MultiplyEvaluateOperator: 344 { 345 result=(MagickRealType) (value*pixel); 346 break; 347 } 348 case OrEvaluateOperator: 349 { 350 result=(MagickRealType) ((size_t) pixel | (size_t) (value+0.5)); 351 break; 352 } 353 case PoissonNoiseEvaluateOperator: 354 { 355 result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel, 356 PoissonNoise,value); 357 break; 358 } 359 case PowEvaluateOperator: 360 { 361 result=(MagickRealType) (QuantumRange*pow((double) (QuantumScale*pixel), 362 (double) value)); 363 break; 364 } 365 case RightShiftEvaluateOperator: 366 { 367 result=(MagickRealType) ((size_t) pixel >> (size_t) (value+0.5)); 368 break; 369 } 370 case SetEvaluateOperator: 371 { 372 result=value; 373 break; 374 } 375 case SineEvaluateOperator: 376 { 377 result=(MagickRealType) (QuantumRange*(0.5*sin((double) (2.0*MagickPI* 378 QuantumScale*pixel*value))+0.5)); 379 break; 380 } 381 case SubtractEvaluateOperator: 382 { 383 result=(MagickRealType) (pixel-value); 384 break; 385 } 386 case ThresholdEvaluateOperator: 387 { 388 result=(MagickRealType) (((MagickRealType) pixel <= value) ? 0 : 389 QuantumRange); 390 break; 391 } 392 case ThresholdBlackEvaluateOperator: 393 { 394 result=(MagickRealType) (((MagickRealType) pixel <= value) ? 0 : pixel); 395 break; 396 } 397 case ThresholdWhiteEvaluateOperator: 398 { 399 result=(MagickRealType) (((MagickRealType) pixel > value) ? QuantumRange : 400 pixel); 401 break; 402 } 403 case UniformNoiseEvaluateOperator: 404 { 405 result=(MagickRealType) GenerateDifferentialNoise(random_info,pixel, 406 UniformNoise,value); 407 break; 408 } 409 case XorEvaluateOperator: 410 { 411 result=(MagickRealType) ((size_t) pixel ^ (size_t) (value+0.5)); 412 break; 413 } 414 } 415 return(result); 416} 417 418MagickExport Image *EvaluateImages(const Image *images, 419 const MagickEvaluateOperator op,ExceptionInfo *exception) 420{ 421#define EvaluateImageTag "Evaluate/Image" 422 423 CacheView 424 *evaluate_view; 425 426 const Image 427 *next; 428 429 Image 430 *evaluate_image; 431 432 MagickBooleanType 433 status; 434 435 MagickOffsetType 436 progress; 437 438 WenusInfo 439 **restrict evaluate_pixels; 440 441 RandomInfo 442 **restrict random_info; 443 444 size_t 445 number_images; 446 447 ssize_t 448 y; 449 450 /* 451 Ensure the image are the same size. 452 */ 453 assert(images != (Image *) NULL); 454 assert(images->signature == MagickSignature); 455 if (images->debug != MagickFalse) 456 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename); 457 assert(exception != (ExceptionInfo *) NULL); 458 assert(exception->signature == MagickSignature); 459 for (next=images; next != (Image *) NULL; next=GetNextImageInList(next)) 460 if ((next->columns != images->columns) || (next->rows != images->rows)) 461 { 462 (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 463 "ImageWidthsOrHeightsDiffer","`%s'",images->filename); 464 return((Image *) NULL); 465 } 466 /* 467 Initialize evaluate next attributes. 468 */ 469 evaluate_image=CloneImage(images,images->columns,images->rows,MagickTrue, 470 exception); 471 if (evaluate_image == (Image *) NULL) 472 return((Image *) NULL); 473 if (SetImageStorageClass(evaluate_image,DirectClass,exception) == MagickFalse) 474 { 475 evaluate_image=DestroyImage(evaluate_image); 476 return((Image *) NULL); 477 } 478 number_images=GetImageListLength(images); 479 evaluate_pixels=AcquirePixelThreadSet(images,number_images); 480 if (evaluate_pixels == (WenusInfo **) NULL) 481 { 482 evaluate_image=DestroyImage(evaluate_image); 483 (void) ThrowMagickException(exception,GetMagickModule(), 484 ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename); 485 return((Image *) NULL); 486 } 487 /* 488 Evaluate image pixels. 489 */ 490 status=MagickTrue; 491 progress=0; 492 random_info=AcquireRandomInfoThreadSet(); 493 evaluate_view=AcquireCacheView(evaluate_image); 494 if (op == MedianEvaluateOperator) 495#if defined(MAGICKCORE_OPENMP_SUPPORT) 496 #pragma omp parallel for schedule(dynamic) shared(progress,status) 497#endif 498 for (y=0; y < (ssize_t) evaluate_image->rows; y++) 499 { 500 CacheView 501 *image_view; 502 503 const Image 504 *next; 505 506 const int 507 id = GetOpenMPThreadId(); 508 509 register WenusInfo 510 *evaluate_pixel; 511 512 register Quantum 513 *restrict q; 514 515 register ssize_t 516 x; 517 518 if (status == MagickFalse) 519 continue; 520 q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,evaluate_image->columns, 521 1,exception); 522 if (q == (Quantum *) NULL) 523 { 524 status=MagickFalse; 525 continue; 526 } 527 evaluate_pixel=evaluate_pixels[id]; 528 for (x=0; x < (ssize_t) evaluate_image->columns; x++) 529 { 530 register ssize_t 531 j, 532 k; 533 534 for (j=0; j < (ssize_t) number_images; j++) 535 for (k=0; k < MaxPixelChannels; k++) 536 evaluate_pixel[j].channel[k]=0.0; 537 next=images; 538 for (j=0; j < (ssize_t) number_images; j++) 539 { 540 register const Quantum 541 *p; 542 543 register ssize_t 544 i; 545 546 image_view=AcquireCacheView(next); 547 p=GetCacheViewVirtualPixels(image_view,x,y,1,1,exception); 548 if (p == (const Quantum *) NULL) 549 { 550 image_view=DestroyCacheView(image_view); 551 break; 552 } 553 for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++) 554 { 555 PixelChannel 556 channel; 557 558 PixelTrait 559 evaluate_traits, 560 traits; 561 562 evaluate_traits=GetPixelChannelMapTraits(evaluate_image, 563 (PixelChannel) i); 564 channel=GetPixelChannelMapChannel(evaluate_image,(PixelChannel) i); 565 traits=GetPixelChannelMapTraits(next,channel); 566 if ((traits == UndefinedPixelTrait) || 567 (evaluate_traits == UndefinedPixelTrait)) 568 continue; 569 if ((evaluate_traits & UpdatePixelTrait) == 0) 570 continue; 571 evaluate_pixel[j].channel[i]=ApplyEvaluateOperator(random_info[id], 572 GetPixelChannel(evaluate_image,channel,p),op, 573 evaluate_pixel[j].channel[i]); 574 } 575 image_view=DestroyCacheView(image_view); 576 next=GetNextImageInList(next); 577 } 578 qsort((void *) evaluate_pixel,number_images,sizeof(*evaluate_pixel), 579 IntensityCompare); 580 for (k=0; k < (ssize_t) GetPixelChannels(evaluate_image); k++) 581 q[k]=ClampToQuantum(evaluate_pixel[j/2].channel[k]); 582 q+=GetPixelChannels(evaluate_image); 583 } 584 if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse) 585 status=MagickFalse; 586 if (images->progress_monitor != (MagickProgressMonitor) NULL) 587 { 588 MagickBooleanType 589 proceed; 590 591#if defined(MAGICKCORE_OPENMP_SUPPORT) 592 #pragma omp critical (MagickCore_EvaluateImages) 593#endif 594 proceed=SetImageProgress(images,EvaluateImageTag,progress++, 595 evaluate_image->rows); 596 if (proceed == MagickFalse) 597 status=MagickFalse; 598 } 599 } 600 else 601#if defined(MAGICKCORE_OPENMP_SUPPORT) 602 #pragma omp parallel for schedule(dynamic) shared(progress,status) 603#endif 604 for (y=0; y < (ssize_t) evaluate_image->rows; y++) 605 { 606 CacheView 607 *image_view; 608 609 const Image 610 *next; 611 612 const int 613 id = GetOpenMPThreadId(); 614 615 register ssize_t 616 i, 617 x; 618 619 register WenusInfo 620 *evaluate_pixel; 621 622 register Quantum 623 *restrict q; 624 625 ssize_t 626 j; 627 628 if (status == MagickFalse) 629 continue; 630 q=QueueCacheViewAuthenticPixels(evaluate_view,0,y,evaluate_image->columns, 631 1,exception); 632 if (q == (Quantum *) NULL) 633 { 634 status=MagickFalse; 635 continue; 636 } 637 evaluate_pixel=evaluate_pixels[id]; 638 for (j=0; j < (ssize_t) evaluate_image->columns; j++) 639 for (i=0; i < MaxPixelChannels; i++) 640 evaluate_pixel[j].channel[i]=0.0; 641 next=images; 642 for (j=0; j < (ssize_t) number_images; j++) 643 { 644 register const Quantum 645 *p; 646 647 image_view=AcquireCacheView(next); 648 p=GetCacheViewVirtualPixels(image_view,0,y,next->columns,1,exception); 649 if (p == (const Quantum *) NULL) 650 { 651 image_view=DestroyCacheView(image_view); 652 break; 653 } 654 for (x=0; x < (ssize_t) next->columns; x++) 655 { 656 register ssize_t 657 i; 658 659 for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++) 660 { 661 PixelChannel 662 channel; 663 664 PixelTrait 665 evaluate_traits, 666 traits; 667 668 evaluate_traits=GetPixelChannelMapTraits(evaluate_image, 669 (PixelChannel) i); 670 channel=GetPixelChannelMapChannel(evaluate_image,(PixelChannel) i); 671 traits=GetPixelChannelMapTraits(next,channel); 672 if ((traits == UndefinedPixelTrait) || 673 (evaluate_traits == UndefinedPixelTrait)) 674 continue; 675 if ((traits & UpdatePixelTrait) == 0) 676 continue; 677 evaluate_pixel[x].channel[i]=ApplyEvaluateOperator(random_info[id], 678 GetPixelChannel(evaluate_image,channel,p),j == 0 ? 679 AddEvaluateOperator : op,evaluate_pixel[x].channel[i]); 680 } 681 p+=GetPixelChannels(next); 682 } 683 image_view=DestroyCacheView(image_view); 684 next=GetNextImageInList(next); 685 } 686 for (x=0; x < (ssize_t) evaluate_image->columns; x++) 687 { 688 register ssize_t 689 i; 690 691 switch (op) 692 { 693 case MeanEvaluateOperator: 694 { 695 for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++) 696 evaluate_pixel[x].channel[i]/=(MagickRealType) number_images; 697 break; 698 } 699 case MultiplyEvaluateOperator: 700 { 701 for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++) 702 { 703 register ssize_t 704 j; 705 706 for (j=0; j < (ssize_t) (number_images-1); j++) 707 evaluate_pixel[x].channel[i]*=QuantumScale; 708 } 709 break; 710 } 711 default: 712 break; 713 } 714 } 715 for (x=0; x < (ssize_t) evaluate_image->columns; x++) 716 { 717 register ssize_t 718 i; 719 720 for (i=0; i < (ssize_t) GetPixelChannels(evaluate_image); i++) 721 { 722 PixelTrait 723 traits; 724 725 traits=GetPixelChannelMapTraits(evaluate_image,(PixelChannel) i); 726 if (traits == UndefinedPixelTrait) 727 continue; 728 if ((traits & UpdatePixelTrait) == 0) 729 continue; 730 q[i]=ClampToQuantum(evaluate_pixel[x].channel[i]); 731 } 732 q+=GetPixelChannels(evaluate_image); 733 } 734 if (SyncCacheViewAuthenticPixels(evaluate_view,exception) == MagickFalse) 735 status=MagickFalse; 736 if (images->progress_monitor != (MagickProgressMonitor) NULL) 737 { 738 MagickBooleanType 739 proceed; 740 741#if defined(MAGICKCORE_OPENMP_SUPPORT) 742 #pragma omp critical (MagickCore_EvaluateImages) 743#endif 744 proceed=SetImageProgress(images,EvaluateImageTag,progress++, 745 evaluate_image->rows); 746 if (proceed == MagickFalse) 747 status=MagickFalse; 748 } 749 } 750 evaluate_view=DestroyCacheView(evaluate_view); 751 evaluate_pixels=DestroyPixelThreadSet(evaluate_pixels); 752 random_info=DestroyRandomInfoThreadSet(random_info); 753 if (status == MagickFalse) 754 evaluate_image=DestroyImage(evaluate_image); 755 return(evaluate_image); 756} 757 758MagickExport MagickBooleanType EvaluateImage(Image *image, 759 const MagickEvaluateOperator op,const double value,ExceptionInfo *exception) 760{ 761 CacheView 762 *image_view; 763 764 MagickBooleanType 765 status; 766 767 MagickOffsetType 768 progress; 769 770 RandomInfo 771 **restrict random_info; 772 773 ssize_t 774 y; 775 776 assert(image != (Image *) NULL); 777 assert(image->signature == MagickSignature); 778 if (image->debug != MagickFalse) 779 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 780 assert(exception != (ExceptionInfo *) NULL); 781 assert(exception->signature == MagickSignature); 782 if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse) 783 return(MagickFalse); 784 status=MagickTrue; 785 progress=0; 786 random_info=AcquireRandomInfoThreadSet(); 787 image_view=AcquireCacheView(image); 788#if defined(MAGICKCORE_OPENMP_SUPPORT) 789 #pragma omp parallel for schedule(dynamic,4) shared(progress,status) 790#endif 791 for (y=0; y < (ssize_t) image->rows; y++) 792 { 793 const int 794 id = GetOpenMPThreadId(); 795 796 register Quantum 797 *restrict q; 798 799 register ssize_t 800 x; 801 802 if (status == MagickFalse) 803 continue; 804 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception); 805 if (q == (Quantum *) NULL) 806 { 807 status=MagickFalse; 808 continue; 809 } 810 for (x=0; x < (ssize_t) image->columns; x++) 811 { 812 register ssize_t 813 i; 814 815 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 816 { 817 PixelTrait 818 traits; 819 820 traits=GetPixelChannelMapTraits(image,(PixelChannel) i); 821 if (traits == UndefinedPixelTrait) 822 continue; 823 q[i]=ClampToQuantum(ApplyEvaluateOperator(random_info[id],q[i],op, 824 value)); 825 } 826 q+=GetPixelChannels(image); 827 } 828 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse) 829 status=MagickFalse; 830 if (image->progress_monitor != (MagickProgressMonitor) NULL) 831 { 832 MagickBooleanType 833 proceed; 834 835#if defined(MAGICKCORE_OPENMP_SUPPORT) 836 #pragma omp critical (MagickCore_EvaluateImage) 837#endif 838 proceed=SetImageProgress(image,EvaluateImageTag,progress++,image->rows); 839 if (proceed == MagickFalse) 840 status=MagickFalse; 841 } 842 } 843 image_view=DestroyCacheView(image_view); 844 random_info=DestroyRandomInfoThreadSet(random_info); 845 return(status); 846} 847 848/* 849%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 850% % 851% % 852% % 853% F u n c t i o n I m a g e % 854% % 855% % 856% % 857%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 858% 859% FunctionImage() applies a value to the image with an arithmetic, relational, 860% or logical operator to an image. Use these operations to lighten or darken 861% an image, to increase or decrease contrast in an image, or to produce the 862% "negative" of an image. 863% 864% The format of the FunctionImage method is: 865% 866% MagickBooleanType FunctionImage(Image *image, 867% const MagickFunction function,const ssize_t number_parameters, 868% const double *parameters,ExceptionInfo *exception) 869% 870% A description of each parameter follows: 871% 872% o image: the image. 873% 874% o function: A channel function. 875% 876% o parameters: one or more parameters. 877% 878% o exception: return any errors or warnings in this structure. 879% 880*/ 881 882static Quantum ApplyFunction(Quantum pixel,const MagickFunction function, 883 const size_t number_parameters,const double *parameters, 884 ExceptionInfo *exception) 885{ 886 MagickRealType 887 result; 888 889 register ssize_t 890 i; 891 892 (void) exception; 893 result=0.0; 894 switch (function) 895 { 896 case PolynomialFunction: 897 { 898 /* 899 Polynomial: polynomial constants, highest to lowest order 900 (e.g. c0*x^3 + c1*x^2 + c2*x + c3). 901 */ 902 result=0.0; 903 for (i=0; i < (ssize_t) number_parameters; i++) 904 result=result*QuantumScale*pixel+parameters[i]; 905 result*=QuantumRange; 906 break; 907 } 908 case SinusoidFunction: 909 { 910 MagickRealType 911 amplitude, 912 bias, 913 frequency, 914 phase; 915 916 /* 917 Sinusoid: frequency, phase, amplitude, bias. 918 */ 919 frequency=(number_parameters >= 1) ? parameters[0] : 1.0; 920 phase=(number_parameters >= 2) ? parameters[1] : 0.0; 921 amplitude=(number_parameters >= 3) ? parameters[2] : 0.5; 922 bias=(number_parameters >= 4) ? parameters[3] : 0.5; 923 result=(MagickRealType) (QuantumRange*(amplitude*sin((double) (2.0* 924 MagickPI*(frequency*QuantumScale*pixel+phase/360.0)))+bias)); 925 break; 926 } 927 case ArcsinFunction: 928 { 929 MagickRealType 930 bias, 931 center, 932 range, 933 width; 934 935 /* 936 Arcsin (peged at range limits for invalid results): 937 width, center, range, and bias. 938 */ 939 width=(number_parameters >= 1) ? parameters[0] : 1.0; 940 center=(number_parameters >= 2) ? parameters[1] : 0.5; 941 range=(number_parameters >= 3) ? parameters[2] : 1.0; 942 bias=(number_parameters >= 4) ? parameters[3] : 0.5; 943 result=2.0/width*(QuantumScale*pixel-center); 944 if ( result <= -1.0 ) 945 result=bias-range/2.0; 946 else 947 if (result >= 1.0) 948 result=bias+range/2.0; 949 else 950 result=(MagickRealType) (range/MagickPI*asin((double) result)+bias); 951 result*=QuantumRange; 952 break; 953 } 954 case ArctanFunction: 955 { 956 MagickRealType 957 center, 958 bias, 959 range, 960 slope; 961 962 /* 963 Arctan: slope, center, range, and bias. 964 */ 965 slope=(number_parameters >= 1) ? parameters[0] : 1.0; 966 center=(number_parameters >= 2) ? parameters[1] : 0.5; 967 range=(number_parameters >= 3) ? parameters[2] : 1.0; 968 bias=(number_parameters >= 4) ? parameters[3] : 0.5; 969 result=(MagickRealType) (MagickPI*slope*(QuantumScale*pixel-center)); 970 result=(MagickRealType) (QuantumRange*(range/MagickPI*atan((double) 971 result)+bias)); 972 break; 973 } 974 case UndefinedFunction: 975 break; 976 } 977 return(ClampToQuantum(result)); 978} 979 980MagickExport MagickBooleanType FunctionImage(Image *image, 981 const MagickFunction function,const size_t number_parameters, 982 const double *parameters,ExceptionInfo *exception) 983{ 984#define FunctionImageTag "Function/Image " 985 986 CacheView 987 *image_view; 988 989 MagickBooleanType 990 status; 991 992 MagickOffsetType 993 progress; 994 995 ssize_t 996 y; 997 998 assert(image != (Image *) NULL); 999 assert(image->signature == MagickSignature); 1000 if (image->debug != MagickFalse) 1001 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 1002 assert(exception != (ExceptionInfo *) NULL); 1003 assert(exception->signature == MagickSignature); 1004 if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse) 1005 return(MagickFalse); 1006 status=MagickTrue; 1007 progress=0; 1008 image_view=AcquireCacheView(image); 1009#if defined(MAGICKCORE_OPENMP_SUPPORT) 1010 #pragma omp parallel for schedule(dynamic,4) shared(progress,status) 1011#endif 1012 for (y=0; y < (ssize_t) image->rows; y++) 1013 { 1014 register Quantum 1015 *restrict q; 1016 1017 register ssize_t 1018 x; 1019 1020 if (status == MagickFalse) 1021 continue; 1022 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception); 1023 if (q == (Quantum *) NULL) 1024 { 1025 status=MagickFalse; 1026 continue; 1027 } 1028 for (x=0; x < (ssize_t) image->columns; x++) 1029 { 1030 register ssize_t 1031 i; 1032 1033 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 1034 { 1035 PixelTrait 1036 traits; 1037 1038 traits=GetPixelChannelMapTraits(image,(PixelChannel) i); 1039 if (traits == UndefinedPixelTrait) 1040 continue; 1041 if ((traits & UpdatePixelTrait) == 0) 1042 continue; 1043 q[i]=ApplyFunction(q[i],function,number_parameters,parameters, 1044 exception); 1045 } 1046 q+=GetPixelChannels(image); 1047 } 1048 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse) 1049 status=MagickFalse; 1050 if (image->progress_monitor != (MagickProgressMonitor) NULL) 1051 { 1052 MagickBooleanType 1053 proceed; 1054 1055#if defined(MAGICKCORE_OPENMP_SUPPORT) 1056 #pragma omp critical (MagickCore_FunctionImage) 1057#endif 1058 proceed=SetImageProgress(image,FunctionImageTag,progress++,image->rows); 1059 if (proceed == MagickFalse) 1060 status=MagickFalse; 1061 } 1062 } 1063 image_view=DestroyCacheView(image_view); 1064 return(status); 1065} 1066 1067/* 1068%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1069% % 1070% % 1071% % 1072% G e t I m a g e E x t r e m a % 1073% % 1074% % 1075% % 1076%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1077% 1078% GetImageExtrema() returns the extrema of one or more image channels. 1079% 1080% The format of the GetImageExtrema method is: 1081% 1082% MagickBooleanType GetImageExtrema(const Image *image,size_t *minima, 1083% size_t *maxima,ExceptionInfo *exception) 1084% 1085% A description of each parameter follows: 1086% 1087% o image: the image. 1088% 1089% o minima: the minimum value in the channel. 1090% 1091% o maxima: the maximum value in the channel. 1092% 1093% o exception: return any errors or warnings in this structure. 1094% 1095*/ 1096MagickExport MagickBooleanType GetImageExtrema(const Image *image, 1097 size_t *minima,size_t *maxima,ExceptionInfo *exception) 1098{ 1099 double 1100 max, 1101 min; 1102 1103 MagickBooleanType 1104 status; 1105 1106 assert(image != (Image *) NULL); 1107 assert(image->signature == MagickSignature); 1108 if (image->debug != MagickFalse) 1109 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 1110 status=GetImageRange(image,&min,&max,exception); 1111 *minima=(size_t) ceil(min-0.5); 1112 *maxima=(size_t) floor(max+0.5); 1113 return(status); 1114} 1115 1116/* 1117%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1118% % 1119% % 1120% % 1121% G e t I m a g e M e a n % 1122% % 1123% % 1124% % 1125%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1126% 1127% GetImageMean() returns the mean and standard deviation of one or more 1128% image channels. 1129% 1130% The format of the GetImageMean method is: 1131% 1132% MagickBooleanType GetImageMean(const Image *image,double *mean, 1133% double *standard_deviation,ExceptionInfo *exception) 1134% 1135% A description of each parameter follows: 1136% 1137% o image: the image. 1138% 1139% o mean: the average value in the channel. 1140% 1141% o standard_deviation: the standard deviation of the channel. 1142% 1143% o exception: return any errors or warnings in this structure. 1144% 1145*/ 1146MagickExport MagickBooleanType GetImageMean(const Image *image,double *mean, 1147 double *standard_deviation,ExceptionInfo *exception) 1148{ 1149 ChannelStatistics 1150 *channel_statistics; 1151 1152 register ssize_t 1153 i; 1154 1155 size_t 1156 area; 1157 1158 assert(image != (Image *) NULL); 1159 assert(image->signature == MagickSignature); 1160 if (image->debug != MagickFalse) 1161 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 1162 channel_statistics=GetImageStatistics(image,exception); 1163 if (channel_statistics == (ChannelStatistics *) NULL) 1164 return(MagickFalse); 1165 area=0; 1166 channel_statistics[MaxPixelChannels].mean=0.0; 1167 channel_statistics[MaxPixelChannels].standard_deviation=0.0; 1168 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 1169 { 1170 PixelTrait 1171 traits; 1172 1173 traits=GetPixelChannelMapTraits(image,(PixelChannel) i); 1174 if (traits == UndefinedPixelTrait) 1175 continue; 1176 if ((traits & UpdatePixelTrait) == 0) 1177 continue; 1178 channel_statistics[MaxPixelChannels].mean+=channel_statistics[i].mean; 1179 channel_statistics[MaxPixelChannels].standard_deviation+= 1180 channel_statistics[i].variance-channel_statistics[i].mean* 1181 channel_statistics[i].mean; 1182 area++; 1183 } 1184 channel_statistics[MaxPixelChannels].mean/=area; 1185 channel_statistics[MaxPixelChannels].standard_deviation= 1186 sqrt(channel_statistics[MaxPixelChannels].standard_deviation/area); 1187 *mean=channel_statistics[MaxPixelChannels].mean; 1188 *standard_deviation=channel_statistics[MaxPixelChannels].standard_deviation; 1189 channel_statistics=(ChannelStatistics *) RelinquishMagickMemory( 1190 channel_statistics); 1191 return(MagickTrue); 1192} 1193 1194/* 1195%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1196% % 1197% % 1198% % 1199% G e t I m a g e K u r t o s i s % 1200% % 1201% % 1202% % 1203%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1204% 1205% GetImageKurtosis() returns the kurtosis and skewness of one or more 1206% image channels. 1207% 1208% The format of the GetImageKurtosis method is: 1209% 1210% MagickBooleanType GetImageKurtosis(const Image *image,double *kurtosis, 1211% double *skewness,ExceptionInfo *exception) 1212% 1213% A description of each parameter follows: 1214% 1215% o image: the image. 1216% 1217% o kurtosis: the kurtosis of the channel. 1218% 1219% o skewness: the skewness of the channel. 1220% 1221% o exception: return any errors or warnings in this structure. 1222% 1223*/ 1224MagickExport MagickBooleanType GetImageKurtosis(const Image *image, 1225 double *kurtosis,double *skewness,ExceptionInfo *exception) 1226{ 1227 CacheView 1228 *image_view; 1229 1230 double 1231 area, 1232 mean, 1233 standard_deviation, 1234 sum_squares, 1235 sum_cubes, 1236 sum_fourth_power; 1237 1238 MagickBooleanType 1239 status; 1240 1241 ssize_t 1242 y; 1243 1244 assert(image != (Image *) NULL); 1245 assert(image->signature == MagickSignature); 1246 if (image->debug != MagickFalse) 1247 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 1248 status=MagickTrue; 1249 *kurtosis=0.0; 1250 *skewness=0.0; 1251 area=0.0; 1252 mean=0.0; 1253 standard_deviation=0.0; 1254 sum_squares=0.0; 1255 sum_cubes=0.0; 1256 sum_fourth_power=0.0; 1257 image_view=AcquireCacheView(image); 1258#if defined(MAGICKCORE_OPENMP_SUPPORT) 1259 #pragma omp parallel for schedule(dynamic) shared(status) omp_throttle(1) 1260#endif 1261 for (y=0; y < (ssize_t) image->rows; y++) 1262 { 1263 register const Quantum 1264 *restrict p; 1265 1266 register ssize_t 1267 x; 1268 1269 if (status == MagickFalse) 1270 continue; 1271 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception); 1272 if (p == (const Quantum *) NULL) 1273 { 1274 status=MagickFalse; 1275 continue; 1276 } 1277 for (x=0; x < (ssize_t) image->columns; x++) 1278 { 1279 register ssize_t 1280 i; 1281 1282 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 1283 { 1284 PixelTrait 1285 traits; 1286 1287 traits=GetPixelChannelMapTraits(image,(PixelChannel) i); 1288 if (traits == UndefinedPixelTrait) 1289 continue; 1290 if ((traits & UpdatePixelTrait) == 0) 1291 continue; 1292#if defined(MAGICKCORE_OPENMP_SUPPORT) 1293 #pragma omp critical (MagickCore_GetImageKurtosis) 1294#endif 1295 { 1296 mean+=p[i]; 1297 sum_squares+=(double) p[i]*p[i]; 1298 sum_cubes+=(double) p[i]*p[i]*p[i]; 1299 sum_fourth_power+=(double) p[i]*p[i]*p[i]*p[i]; 1300 area++; 1301 } 1302 } 1303 p+=GetPixelChannels(image); 1304 } 1305 } 1306 image_view=DestroyCacheView(image_view); 1307 if (area != 0.0) 1308 { 1309 mean/=area; 1310 sum_squares/=area; 1311 sum_cubes/=area; 1312 sum_fourth_power/=area; 1313 } 1314 standard_deviation=sqrt(sum_squares-(mean*mean)); 1315 if (standard_deviation != 0.0) 1316 { 1317 *kurtosis=sum_fourth_power-4.0*mean*sum_cubes+6.0*mean*mean*sum_squares- 1318 3.0*mean*mean*mean*mean; 1319 *kurtosis/=standard_deviation*standard_deviation*standard_deviation* 1320 standard_deviation; 1321 *kurtosis-=3.0; 1322 *skewness=sum_cubes-3.0*mean*sum_squares+2.0*mean*mean*mean; 1323 *skewness/=standard_deviation*standard_deviation*standard_deviation; 1324 } 1325 return(status); 1326} 1327 1328/* 1329%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1330% % 1331% % 1332% % 1333% G e t I m a g e R a n g e % 1334% % 1335% % 1336% % 1337%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1338% 1339% GetImageRange() returns the range of one or more image channels. 1340% 1341% The format of the GetImageRange method is: 1342% 1343% MagickBooleanType GetImageRange(const Image *image,double *minima, 1344% double *maxima,ExceptionInfo *exception) 1345% 1346% A description of each parameter follows: 1347% 1348% o image: the image. 1349% 1350% o minima: the minimum value in the channel. 1351% 1352% o maxima: the maximum value in the channel. 1353% 1354% o exception: return any errors or warnings in this structure. 1355% 1356*/ 1357MagickExport MagickBooleanType GetImageRange(const Image *image,double *minima, 1358 double *maxima,ExceptionInfo *exception) 1359{ 1360 CacheView 1361 *image_view; 1362 1363 MagickBooleanType 1364 status; 1365 1366 ssize_t 1367 y; 1368 1369 assert(image != (Image *) NULL); 1370 assert(image->signature == MagickSignature); 1371 if (image->debug != MagickFalse) 1372 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 1373 status=MagickTrue; 1374 *maxima=(-MagickHuge); 1375 *minima=MagickHuge; 1376 image_view=AcquireCacheView(image); 1377#if defined(MAGICKCORE_OPENMP_SUPPORT) 1378 #pragma omp parallel for schedule(dynamic) shared(status) omp_throttle(1) 1379#endif 1380 for (y=0; y < (ssize_t) image->rows; y++) 1381 { 1382 register const Quantum 1383 *restrict p; 1384 1385 register ssize_t 1386 x; 1387 1388 if (status == MagickFalse) 1389 continue; 1390 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception); 1391 if (p == (const Quantum *) NULL) 1392 { 1393 status=MagickFalse; 1394 continue; 1395 } 1396 for (x=0; x < (ssize_t) image->columns; x++) 1397 { 1398 register ssize_t 1399 i; 1400 1401 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 1402 { 1403 PixelTrait 1404 traits; 1405 1406 traits=GetPixelChannelMapTraits(image,(PixelChannel) i); 1407 if (traits == UndefinedPixelTrait) 1408 continue; 1409 if ((traits & UpdatePixelTrait) == 0) 1410 continue; 1411#if defined(MAGICKCORE_OPENMP_SUPPORT) 1412 #pragma omp critical (MagickCore_GetImageRange) 1413#endif 1414 { 1415 if (p[i] < *minima) 1416 *minima=(double) p[i]; 1417 if (p[i] > *maxima) 1418 *maxima=(double) p[i]; 1419 } 1420 } 1421 p+=GetPixelChannels(image); 1422 } 1423 } 1424 image_view=DestroyCacheView(image_view); 1425 return(status); 1426} 1427 1428/* 1429%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1430% % 1431% % 1432% % 1433% G e t I m a g e S t a t i s t i c s % 1434% % 1435% % 1436% % 1437%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1438% 1439% GetImageStatistics() returns statistics for each channel in the 1440% image. The statistics include the channel depth, its minima, maxima, mean, 1441% standard deviation, kurtosis and skewness. You can access the red channel 1442% mean, for example, like this: 1443% 1444% channel_statistics=GetImageStatistics(image,exception); 1445% red_mean=channel_statistics[RedPixelChannel].mean; 1446% 1447% Use MagickRelinquishMemory() to free the statistics buffer. 1448% 1449% The format of the GetImageStatistics method is: 1450% 1451% ChannelStatistics *GetImageStatistics(const Image *image, 1452% ExceptionInfo *exception) 1453% 1454% A description of each parameter follows: 1455% 1456% o image: the image. 1457% 1458% o exception: return any errors or warnings in this structure. 1459% 1460*/ 1461 1462static size_t GetImageChannels(const Image *image) 1463{ 1464 register ssize_t 1465 i; 1466 1467 size_t 1468 channels; 1469 1470 channels=0; 1471 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 1472 { 1473 PixelTrait 1474 traits; 1475 1476 traits=GetPixelChannelMapTraits(image,(PixelChannel) i); 1477 if ((traits & UpdatePixelTrait) != 0) 1478 channels++; 1479 } 1480 return(channels); 1481} 1482 1483MagickExport ChannelStatistics *GetImageStatistics(const Image *image, 1484 ExceptionInfo *exception) 1485{ 1486 ChannelStatistics 1487 *channel_statistics; 1488 1489 double 1490 area; 1491 1492 MagickStatusType 1493 status; 1494 1495 QuantumAny 1496 range; 1497 1498 register ssize_t 1499 i; 1500 1501 size_t 1502 channels, 1503 depth; 1504 1505 ssize_t 1506 y; 1507 1508 assert(image != (Image *) NULL); 1509 assert(image->signature == MagickSignature); 1510 if (image->debug != MagickFalse) 1511 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 1512 channel_statistics=(ChannelStatistics *) AcquireQuantumMemory( 1513 MaxPixelChannels+1,sizeof(*channel_statistics)); 1514 if (channel_statistics == (ChannelStatistics *) NULL) 1515 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed"); 1516 (void) ResetMagickMemory(channel_statistics,0,(MaxPixelChannels+1)* 1517 sizeof(*channel_statistics)); 1518 for (i=0; i <= (ssize_t) MaxPixelChannels; i++) 1519 { 1520 channel_statistics[i].depth=1; 1521 channel_statistics[i].maxima=(-MagickHuge); 1522 channel_statistics[i].minima=MagickHuge; 1523 } 1524 for (y=0; y < (ssize_t) image->rows; y++) 1525 { 1526 register const Quantum 1527 *restrict p; 1528 1529 register ssize_t 1530 x; 1531 1532 p=GetVirtualPixels(image,0,y,image->columns,1,exception); 1533 if (p == (const Quantum *) NULL) 1534 break; 1535 for (x=0; x < (ssize_t) image->columns; x++) 1536 { 1537 register ssize_t 1538 i; 1539 1540 for (i=0; i < (ssize_t) GetPixelChannels(image); i++) 1541 { 1542 PixelTrait 1543 traits; 1544 1545 traits=GetPixelChannelMapTraits(image,(PixelChannel) i); 1546 if (traits == UndefinedPixelTrait) 1547 continue; 1548 if (channel_statistics[i].depth != MAGICKCORE_QUANTUM_DEPTH) 1549 { 1550 depth=channel_statistics[i].depth; 1551 range=GetQuantumRange(depth); 1552 status=p[i] != ScaleAnyToQuantum(ScaleQuantumToAny(p[i],range), 1553 range) ? MagickTrue : MagickFalse; 1554 if (status != MagickFalse) 1555 { 1556 channel_statistics[i].depth++; 1557 continue; 1558 } 1559 } 1560 if ((double) p[i] < channel_statistics[i].minima) 1561 channel_statistics[i].minima=(double) p[i]; 1562 if ((double) p[i] > channel_statistics[i].maxima) 1563 channel_statistics[i].maxima=(double) p[i]; 1564 channel_statistics[i].sum+=p[i]; 1565 channel_statistics[i].sum_squared+=(double) p[i]*p[i]; 1566 channel_statistics[i].sum_cubed+=(double) p[i]*p[i]*p[i]; 1567 channel_statistics[i].sum_fourth_power+=(double) p[i]*p[i]*p[i]*p[i]; 1568 } 1569 p+=GetPixelChannels(image); 1570 } 1571 } 1572 area=(double) image->columns*image->rows; 1573 for (i=0; i < (ssize_t) MaxPixelChannels; i++) 1574 { 1575 channel_statistics[i].sum/=area; 1576 channel_statistics[i].sum_squared/=area; 1577 channel_statistics[i].sum_cubed/=area; 1578 channel_statistics[i].sum_fourth_power/=area; 1579 channel_statistics[i].mean=channel_statistics[i].sum; 1580 channel_statistics[i].variance=channel_statistics[i].sum_squared; 1581 channel_statistics[i].standard_deviation=sqrt( 1582 channel_statistics[i].variance-(channel_statistics[i].mean* 1583 channel_statistics[i].mean)); 1584 } 1585 for (i=0; i < (ssize_t) MaxPixelChannels; i++) 1586 { 1587 channel_statistics[MaxPixelChannels].depth=(size_t) MagickMax((double) 1588 channel_statistics[MaxPixelChannels].depth,(double) 1589 channel_statistics[i].depth); 1590 channel_statistics[MaxPixelChannels].minima=MagickMin( 1591 channel_statistics[MaxPixelChannels].minima, 1592 channel_statistics[i].minima); 1593 channel_statistics[MaxPixelChannels].maxima=MagickMax( 1594 channel_statistics[MaxPixelChannels].maxima, 1595 channel_statistics[i].maxima); 1596 channel_statistics[MaxPixelChannels].sum+=channel_statistics[i].sum; 1597 channel_statistics[MaxPixelChannels].sum_squared+= 1598 channel_statistics[i].sum_squared; 1599 channel_statistics[MaxPixelChannels].sum_cubed+= 1600 channel_statistics[i].sum_cubed; 1601 channel_statistics[MaxPixelChannels].sum_fourth_power+= 1602 channel_statistics[i].sum_fourth_power; 1603 channel_statistics[MaxPixelChannels].mean+=channel_statistics[i].mean; 1604 channel_statistics[MaxPixelChannels].variance+= 1605 channel_statistics[i].variance-channel_statistics[i].mean* 1606 channel_statistics[i].mean; 1607 channel_statistics[MaxPixelChannels].standard_deviation+= 1608 channel_statistics[i].variance-channel_statistics[i].mean* 1609 channel_statistics[i].mean; 1610 } 1611 channels=GetImageChannels(image); 1612 channel_statistics[MaxPixelChannels].sum/=channels; 1613 channel_statistics[MaxPixelChannels].sum_squared/=channels; 1614 channel_statistics[MaxPixelChannels].sum_cubed/=channels; 1615 channel_statistics[MaxPixelChannels].sum_fourth_power/=channels; 1616 channel_statistics[MaxPixelChannels].mean/=channels; 1617 channel_statistics[MaxPixelChannels].variance/=channels; 1618 channel_statistics[MaxPixelChannels].standard_deviation= 1619 sqrt(channel_statistics[MaxPixelChannels].standard_deviation/channels); 1620 channel_statistics[MaxPixelChannels].kurtosis/=channels; 1621 channel_statistics[MaxPixelChannels].skewness/=channels; 1622 for (i=0; i <= (ssize_t) MaxPixelChannels; i++) 1623 { 1624 if (channel_statistics[i].standard_deviation == 0.0) 1625 continue; 1626 channel_statistics[i].skewness=(channel_statistics[i].sum_cubed- 1627 3.0*channel_statistics[i].mean*channel_statistics[i].sum_squared+ 1628 2.0*channel_statistics[i].mean*channel_statistics[i].mean* 1629 channel_statistics[i].mean)/(channel_statistics[i].standard_deviation* 1630 channel_statistics[i].standard_deviation* 1631 channel_statistics[i].standard_deviation); 1632 channel_statistics[i].kurtosis=(channel_statistics[i].sum_fourth_power- 1633 4.0*channel_statistics[i].mean*channel_statistics[i].sum_cubed+ 1634 6.0*channel_statistics[i].mean*channel_statistics[i].mean* 1635 channel_statistics[i].sum_squared-3.0*channel_statistics[i].mean* 1636 channel_statistics[i].mean*1.0*channel_statistics[i].mean* 1637 channel_statistics[i].mean)/(channel_statistics[i].standard_deviation* 1638 channel_statistics[i].standard_deviation* 1639 channel_statistics[i].standard_deviation* 1640 channel_statistics[i].standard_deviation)-3.0; 1641 } 1642 return(channel_statistics); 1643} 1644