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