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