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