1 /* 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % % 4 % % 5 % % 6 % SSSSS EEEEE GGGG M M EEEEE N N TTTTT % 7 % SS E G MM MM E NN N T % 8 % SSS EEE G GGG M M M EEE N N N T % 9 % SS E G G M M E N NN T % 10 % SSSSS EEEEE GGGG M M EEEEE N N T % 11 % % 12 % % 13 % MagickCore Methods to Segment an Image with Thresholding Fuzzy c-Means % 14 % % 15 % Software Design % 16 % Cristy % 17 % April 1993 % 18 % % 19 % % 20 % Copyright 1999-2016 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 % Segment segments an image by analyzing the histograms of the color 37 % components and identifying units that are homogeneous with the fuzzy 38 % c-means technique. The scale-space filter analyzes the histograms of 39 % the three color components of the image and identifies a set of 40 % classes. The extents of each class is used to coarsely segment the 41 % image with thresholding. The color associated with each class is 42 % determined by the mean color of all pixels within the extents of a 43 % particular class. Finally, any unclassified pixels are assigned to 44 % the closest class with the fuzzy c-means technique. 45 % 46 % The fuzzy c-Means algorithm can be summarized as follows: 47 % 48 % o Build a histogram, one for each color component of the image. 49 % 50 % o For each histogram, successively apply the scale-space filter and 51 % build an interval tree of zero crossings in the second derivative 52 % at each scale. Analyze this scale-space ''fingerprint'' to 53 % determine which peaks and valleys in the histogram are most 54 % predominant. 55 % 56 % o The fingerprint defines intervals on the axis of the histogram. 57 % Each interval contains either a minima or a maxima in the original 58 % signal. If each color component lies within the maxima interval, 59 % that pixel is considered ''classified'' and is assigned an unique 60 % class number. 61 % 62 % o Any pixel that fails to be classified in the above thresholding 63 % pass is classified using the fuzzy c-Means technique. It is 64 % assigned to one of the classes discovered in the histogram analysis 65 % phase. 66 % 67 % The fuzzy c-Means technique attempts to cluster a pixel by finding 68 % the local minima of the generalized within group sum of squared error 69 % objective function. A pixel is assigned to the closest class of 70 % which the fuzzy membership has a maximum value. 71 % 72 % Segment is strongly based on software written by Andy Gallo, 73 % University of Delaware. 74 % 75 % The following reference was used in creating this program: 76 % 77 % Young Won Lim, Sang Uk Lee, "On The Color Image Segmentation 78 % Algorithm Based on the Thresholding and the Fuzzy c-Means 79 % Techniques", Pattern Recognition, Volume 23, Number 9, pages 80 % 935-952, 1990. 81 % 82 % 83 */ 84 85 #include "MagickCore/studio.h" 87 #include "MagickCore/cache.h" 88 #include "MagickCore/color.h" 89 #include "MagickCore/colormap.h" 90 #include "MagickCore/colorspace.h" 91 #include "MagickCore/colorspace-private.h" 92 #include "MagickCore/exception.h" 93 #include "MagickCore/exception-private.h" 94 #include "MagickCore/image.h" 95 #include "MagickCore/image-private.h" 96 #include "MagickCore/memory_.h" 97 #include "MagickCore/monitor.h" 98 #include "MagickCore/monitor-private.h" 99 #include "MagickCore/pixel-accessor.h" 100 #include "MagickCore/pixel-private.h" 101 #include "MagickCore/quantize.h" 102 #include "MagickCore/quantum.h" 103 #include "MagickCore/quantum-private.h" 104 #include "MagickCore/resource_.h" 105 #include "MagickCore/segment.h" 106 #include "MagickCore/string_.h" 107 #include "MagickCore/thread-private.h" 108 109 /* 111 Define declarations. 112 */ 113 #define MaxDimension 3 114 #define DeltaTau 0.5f 115 #if defined(FastClassify) 116 #define WeightingExponent 2.0 117 #define SegmentPower(ratio) (ratio) 118 #else 119 #define WeightingExponent 2.5 120 #define SegmentPower(ratio) pow(ratio,(double) (1.0/(weighting_exponent-1.0))); 121 #endif 122 #define Tau 5.2f 123 124 /* 126 Typedef declarations. 127 */ 128 typedef struct _ExtentPacket 129 { 130 double 131 center; 132 133 ssize_t 134 index, 135 left, 136 right; 137 } ExtentPacket; 138 139 typedef struct _Cluster 140 { 141 struct _Cluster 142 *next; 143 144 ExtentPacket 145 red, 146 green, 147 blue; 148 149 ssize_t 150 count, 151 id; 152 } Cluster; 153 154 typedef struct _IntervalTree 155 { 156 double 157 tau; 158 159 ssize_t 160 left, 161 right; 162 163 double 164 mean_stability, 165 stability; 166 167 struct _IntervalTree 168 *sibling, 169 *child; 170 } IntervalTree; 171 172 typedef struct _ZeroCrossing 173 { 174 double 175 tau, 176 histogram[256]; 177 178 short 179 crossings[256]; 180 } ZeroCrossing; 181 182 /* 184 Constant declarations. 185 */ 186 static const int 187 Blue = 2, 188 Green = 1, 189 Red = 0, 190 SafeMargin = 3, 191 TreeLength = 600; 192 193 /* 195 Method prototypes. 196 */ 197 static double 198 OptimalTau(const ssize_t *,const double,const double,const double, 199 const double,short *); 200 201 static ssize_t 202 DefineRegion(const short *,ExtentPacket *); 203 204 static void 205 InitializeHistogram(const Image *,ssize_t **,ExceptionInfo *), 206 ScaleSpace(const ssize_t *,const double,double *), 207 ZeroCrossHistogram(double *,const double,short *); 208 209 /* 211 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 212 % % 213 % % 214 % % 215 + C l a s s i f y % 216 % % 217 % % 218 % % 219 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 220 % 221 % Classify() defines one or more classes. Each pixel is thresholded to 222 % determine which class it belongs to. If the class is not identified it is 223 % assigned to the closest class based on the fuzzy c-Means technique. 224 % 225 % The format of the Classify method is: 226 % 227 % MagickBooleanType Classify(Image *image,short **extrema, 228 % const double cluster_threshold, 229 % const double weighting_exponent, 230 % const MagickBooleanType verbose,ExceptionInfo *exception) 231 % 232 % A description of each parameter follows. 233 % 234 % o image: the image. 235 % 236 % o extrema: Specifies a pointer to an array of integers. They 237 % represent the peaks and valleys of the histogram for each color 238 % component. 239 % 240 % o cluster_threshold: This double represents the minimum number of 241 % pixels contained in a hexahedra before it can be considered valid 242 % (expressed as a percentage). 243 % 244 % o weighting_exponent: Specifies the membership weighting exponent. 245 % 246 % o verbose: A value greater than zero prints detailed information about 247 % the identified classes. 248 % 249 % o exception: return any errors or warnings in this structure. 250 % 251 */ 252 static MagickBooleanType Classify(Image *image,short **extrema, 253 const double cluster_threshold, 254 const double weighting_exponent,const MagickBooleanType verbose, 255 ExceptionInfo *exception) 256 { 257 #define SegmentImageTag "Segment/Image" 258 259 CacheView 260 *image_view; 261 262 Cluster 263 *cluster, 264 *head, 265 *last_cluster, 266 *next_cluster; 267 268 ExtentPacket 269 blue, 270 green, 271 red; 272 273 MagickOffsetType 274 progress; 275 276 double 277 *free_squares; 278 279 MagickStatusType 280 status; 281 282 register ssize_t 283 i; 284 285 register double 286 *squares; 287 288 size_t 289 number_clusters; 290 291 ssize_t 292 count, 293 y; 294 295 /* 296 Form clusters. 297 */ 298 cluster=(Cluster *) NULL; 299 head=(Cluster *) NULL; 300 (void) ResetMagickMemory(&red,0,sizeof(red)); 301 (void) ResetMagickMemory(&green,0,sizeof(green)); 302 (void) ResetMagickMemory(&blue,0,sizeof(blue)); 303 while (DefineRegion(extrema[Red],&red) != 0) 304 { 305 green.index=0; 306 while (DefineRegion(extrema[Green],&green) != 0) 307 { 308 blue.index=0; 309 while (DefineRegion(extrema[Blue],&blue) != 0) 310 { 311 /* 312 Allocate a new class. 313 */ 314 if (head != (Cluster *) NULL) 315 { 316 cluster->next=(Cluster *) AcquireMagickMemory( 317 sizeof(*cluster->next)); 318 cluster=cluster->next; 319 } 320 else 321 { 322 cluster=(Cluster *) AcquireMagickMemory(sizeof(*cluster)); 323 head=cluster; 324 } 325 if (cluster == (Cluster *) NULL) 326 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed", 327 image->filename); 328 /* 329 Initialize a new class. 330 */ 331 cluster->count=0; 332 cluster->red=red; 333 cluster->green=green; 334 cluster->blue=blue; 335 cluster->next=(Cluster *) NULL; 336 } 337 } 338 } 339 if (head == (Cluster *) NULL) 340 { 341 /* 342 No classes were identified-- create one. 343 */ 344 cluster=(Cluster *) AcquireMagickMemory(sizeof(*cluster)); 345 if (cluster == (Cluster *) NULL) 346 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed", 347 image->filename); 348 /* 349 Initialize a new class. 350 */ 351 cluster->count=0; 352 cluster->red=red; 353 cluster->green=green; 354 cluster->blue=blue; 355 cluster->next=(Cluster *) NULL; 356 head=cluster; 357 } 358 /* 359 Count the pixels for each cluster. 360 */ 361 status=MagickTrue; 362 count=0; 363 progress=0; 364 image_view=AcquireVirtualCacheView(image,exception); 365 for (y=0; y < (ssize_t) image->rows; y++) 366 { 367 register const Quantum 368 *p; 369 370 register ssize_t 371 x; 372 373 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception); 374 if (p == (const Quantum *) NULL) 375 break; 376 for (x=0; x < (ssize_t) image->columns; x++) 377 { 378 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next) 379 if (((ssize_t) ScaleQuantumToChar(GetPixelRed(image,p)) >= 380 (cluster->red.left-SafeMargin)) && 381 ((ssize_t) ScaleQuantumToChar(GetPixelRed(image,p)) <= 382 (cluster->red.right+SafeMargin)) && 383 ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,p)) >= 384 (cluster->green.left-SafeMargin)) && 385 ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,p)) <= 386 (cluster->green.right+SafeMargin)) && 387 ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,p)) >= 388 (cluster->blue.left-SafeMargin)) && 389 ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,p)) <= 390 (cluster->blue.right+SafeMargin))) 391 { 392 /* 393 Count this pixel. 394 */ 395 count++; 396 cluster->red.center+=(double) ScaleQuantumToChar( 397 GetPixelRed(image,p)); 398 cluster->green.center+=(double) ScaleQuantumToChar( 399 GetPixelGreen(image,p)); 400 cluster->blue.center+=(double) ScaleQuantumToChar( 401 GetPixelBlue(image,p)); 402 cluster->count++; 403 break; 404 } 405 p+=GetPixelChannels(image); 406 } 407 if (image->progress_monitor != (MagickProgressMonitor) NULL) 408 { 409 MagickBooleanType 410 proceed; 411 412 #if defined(MAGICKCORE_OPENMP_SUPPORT) 413 #pragma omp critical (MagickCore_Classify) 414 #endif 415 proceed=SetImageProgress(image,SegmentImageTag,progress++,2* 416 image->rows); 417 if (proceed == MagickFalse) 418 status=MagickFalse; 419 } 420 } 421 image_view=DestroyCacheView(image_view); 422 /* 423 Remove clusters that do not meet minimum cluster threshold. 424 */ 425 count=0; 426 last_cluster=head; 427 next_cluster=head; 428 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster) 429 { 430 next_cluster=cluster->next; 431 if ((cluster->count > 0) && 432 (cluster->count >= (count*cluster_threshold/100.0))) 433 { 434 /* 435 Initialize cluster. 436 */ 437 cluster->id=count; 438 cluster->red.center/=cluster->count; 439 cluster->green.center/=cluster->count; 440 cluster->blue.center/=cluster->count; 441 count++; 442 last_cluster=cluster; 443 continue; 444 } 445 /* 446 Delete cluster. 447 */ 448 if (cluster == head) 449 head=next_cluster; 450 else 451 last_cluster->next=next_cluster; 452 cluster=(Cluster *) RelinquishMagickMemory(cluster); 453 } 454 number_clusters=(size_t) count; 455 if (verbose != MagickFalse) 456 { 457 /* 458 Print cluster statistics. 459 */ 460 (void) FormatLocaleFile(stdout,"Fuzzy C-means Statistics\n"); 461 (void) FormatLocaleFile(stdout,"===================\n\n"); 462 (void) FormatLocaleFile(stdout,"\tCluster Threshold = %g\n",(double) 463 cluster_threshold); 464 (void) FormatLocaleFile(stdout,"\tWeighting Exponent = %g\n",(double) 465 weighting_exponent); 466 (void) FormatLocaleFile(stdout,"\tTotal Number of Clusters = %.20g\n\n", 467 (double) number_clusters); 468 /* 469 Print the total number of points per cluster. 470 */ 471 (void) FormatLocaleFile(stdout,"\n\nNumber of Vectors Per Cluster\n"); 472 (void) FormatLocaleFile(stdout,"=============================\n\n"); 473 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next) 474 (void) FormatLocaleFile(stdout,"Cluster #%.20g = %.20g\n",(double) 475 cluster->id,(double) cluster->count); 476 /* 477 Print the cluster extents. 478 */ 479 (void) FormatLocaleFile(stdout, 480 "\n\n\nCluster Extents: (Vector Size: %d)\n",MaxDimension); 481 (void) FormatLocaleFile(stdout,"================"); 482 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next) 483 { 484 (void) FormatLocaleFile(stdout,"\n\nCluster #%.20g\n\n",(double) 485 cluster->id); 486 (void) FormatLocaleFile(stdout, 487 "%.20g-%.20g %.20g-%.20g %.20g-%.20g\n",(double) 488 cluster->red.left,(double) cluster->red.right,(double) 489 cluster->green.left,(double) cluster->green.right,(double) 490 cluster->blue.left,(double) cluster->blue.right); 491 } 492 /* 493 Print the cluster center values. 494 */ 495 (void) FormatLocaleFile(stdout, 496 "\n\n\nCluster Center Values: (Vector Size: %d)\n",MaxDimension); 497 (void) FormatLocaleFile(stdout,"====================="); 498 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next) 499 { 500 (void) FormatLocaleFile(stdout,"\n\nCluster #%.20g\n\n",(double) 501 cluster->id); 502 (void) FormatLocaleFile(stdout,"%g %g %g\n",(double) 503 cluster->red.center,(double) cluster->green.center,(double) 504 cluster->blue.center); 505 } 506 (void) FormatLocaleFile(stdout,"\n"); 507 } 508 if (number_clusters > 256) 509 ThrowBinaryException(ImageError,"TooManyClusters",image->filename); 510 /* 511 Speed up distance calculations. 512 */ 513 squares=(double *) AcquireQuantumMemory(513UL,sizeof(*squares)); 514 if (squares == (double *) NULL) 515 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed", 516 image->filename); 517 squares+=255; 518 for (i=(-255); i <= 255; i++) 519 squares[i]=(double) i*(double) i; 520 /* 521 Allocate image colormap. 522 */ 523 if (AcquireImageColormap(image,number_clusters,exception) == MagickFalse) 524 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed", 525 image->filename); 526 i=0; 527 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next) 528 { 529 image->colormap[i].red=(double) ScaleCharToQuantum((unsigned char) 530 (cluster->red.center+0.5)); 531 image->colormap[i].green=(double) ScaleCharToQuantum((unsigned char) 532 (cluster->green.center+0.5)); 533 image->colormap[i].blue=(double) ScaleCharToQuantum((unsigned char) 534 (cluster->blue.center+0.5)); 535 i++; 536 } 537 /* 538 Do course grain classes. 539 */ 540 image_view=AcquireAuthenticCacheView(image,exception); 541 #if defined(MAGICKCORE_OPENMP_SUPPORT) 542 #pragma omp parallel for schedule(static,4) shared(progress,status) \ 543 magick_threads(image,image,image->rows,1) 544 #endif 545 for (y=0; y < (ssize_t) image->rows; y++) 546 { 547 Cluster 548 *cluster; 549 550 register const PixelInfo 551 *magick_restrict p; 552 553 register ssize_t 554 x; 555 556 register Quantum 557 *magick_restrict q; 558 559 if (status == MagickFalse) 560 continue; 561 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception); 562 if (q == (Quantum *) NULL) 563 { 564 status=MagickFalse; 565 continue; 566 } 567 for (x=0; x < (ssize_t) image->columns; x++) 568 { 569 SetPixelIndex(image,0,q); 570 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next) 571 { 572 if (((ssize_t) ScaleQuantumToChar(GetPixelRed(image,q)) >= 573 (cluster->red.left-SafeMargin)) && 574 ((ssize_t) ScaleQuantumToChar(GetPixelRed(image,q)) <= 575 (cluster->red.right+SafeMargin)) && 576 ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,q)) >= 577 (cluster->green.left-SafeMargin)) && 578 ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,q)) <= 579 (cluster->green.right+SafeMargin)) && 580 ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,q)) >= 581 (cluster->blue.left-SafeMargin)) && 582 ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,q)) <= 583 (cluster->blue.right+SafeMargin))) 584 { 585 /* 586 Classify this pixel. 587 */ 588 SetPixelIndex(image,(Quantum) cluster->id,q); 589 break; 590 } 591 } 592 if (cluster == (Cluster *) NULL) 593 { 594 double 595 distance_squared, 596 local_minima, 597 numerator, 598 ratio, 599 sum; 600 601 register ssize_t 602 j, 603 k; 604 605 /* 606 Compute fuzzy membership. 607 */ 608 local_minima=0.0; 609 for (j=0; j < (ssize_t) image->colors; j++) 610 { 611 sum=0.0; 612 p=image->colormap+j; 613 distance_squared=squares[(ssize_t) ScaleQuantumToChar( 614 GetPixelRed(image,q))-(ssize_t) 615 ScaleQuantumToChar(ClampToQuantum(p->red))]+squares[(ssize_t) 616 ScaleQuantumToChar(GetPixelGreen(image,q))-(ssize_t) 617 ScaleQuantumToChar(ClampToQuantum(p->green))]+squares[(ssize_t) 618 ScaleQuantumToChar(GetPixelBlue(image,q))-(ssize_t) 619 ScaleQuantumToChar(ClampToQuantum(p->blue))]; 620 numerator=distance_squared; 621 for (k=0; k < (ssize_t) image->colors; k++) 622 { 623 p=image->colormap+k; 624 distance_squared=squares[(ssize_t) ScaleQuantumToChar( 625 GetPixelRed(image,q))-(ssize_t) 626 ScaleQuantumToChar(ClampToQuantum(p->red))]+squares[ 627 (ssize_t) ScaleQuantumToChar(GetPixelGreen(image,q))-(ssize_t) 628 ScaleQuantumToChar(ClampToQuantum(p->green))]+squares[ 629 (ssize_t) ScaleQuantumToChar(GetPixelBlue(image,q))-(ssize_t) 630 ScaleQuantumToChar(ClampToQuantum(p->blue))]; 631 ratio=numerator/distance_squared; 632 sum+=SegmentPower(ratio); 633 } 634 if ((sum != 0.0) && ((1.0/sum) > local_minima)) 635 { 636 /* 637 Classify this pixel. 638 */ 639 local_minima=1.0/sum; 640 SetPixelIndex(image,(Quantum) j,q); 641 } 642 } 643 } 644 q+=GetPixelChannels(image); 645 } 646 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse) 647 status=MagickFalse; 648 if (image->progress_monitor != (MagickProgressMonitor) NULL) 649 { 650 MagickBooleanType 651 proceed; 652 653 #if defined(MAGICKCORE_OPENMP_SUPPORT) 654 #pragma omp critical (MagickCore_Classify) 655 #endif 656 proceed=SetImageProgress(image,SegmentImageTag,progress++, 657 2*image->rows); 658 if (proceed == MagickFalse) 659 status=MagickFalse; 660 } 661 } 662 image_view=DestroyCacheView(image_view); 663 status&=SyncImage(image,exception); 664 /* 665 Relinquish resources. 666 */ 667 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster) 668 { 669 next_cluster=cluster->next; 670 cluster=(Cluster *) RelinquishMagickMemory(cluster); 671 } 672 squares-=255; 673 free_squares=squares; 674 free_squares=(double *) RelinquishMagickMemory(free_squares); 675 return(MagickTrue); 676 } 677 678 /* 680 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 681 % % 682 % % 683 % % 684 + C o n s o l i d a t e C r o s s i n g s % 685 % % 686 % % 687 % % 688 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 689 % 690 % ConsolidateCrossings() guarantees that an even number of zero crossings 691 % always lie between two crossings. 692 % 693 % The format of the ConsolidateCrossings method is: 694 % 695 % ConsolidateCrossings(ZeroCrossing *zero_crossing, 696 % const size_t number_crossings) 697 % 698 % A description of each parameter follows. 699 % 700 % o zero_crossing: Specifies an array of structures of type ZeroCrossing. 701 % 702 % o number_crossings: This size_t specifies the number of elements 703 % in the zero_crossing array. 704 % 705 */ 706 static void ConsolidateCrossings(ZeroCrossing *zero_crossing, 707 const size_t number_crossings) 708 { 709 register ssize_t 710 i, 711 j, 712 k, 713 l; 714 715 ssize_t 716 center, 717 correct, 718 count, 719 left, 720 right; 721 722 /* 723 Consolidate zero crossings. 724 */ 725 for (i=(ssize_t) number_crossings-1; i >= 0; i--) 726 for (j=0; j <= 255; j++) 727 { 728 if (zero_crossing[i].crossings[j] == 0) 729 continue; 730 /* 731 Find the entry that is closest to j and still preserves the 732 property that there are an even number of crossings between 733 intervals. 734 */ 735 for (k=j-1; k > 0; k--) 736 if (zero_crossing[i+1].crossings[k] != 0) 737 break; 738 left=MagickMax(k,0); 739 center=j; 740 for (k=j+1; k < 255; k++) 741 if (zero_crossing[i+1].crossings[k] != 0) 742 break; 743 right=MagickMin(k,255); 744 /* 745 K is the zero crossing just left of j. 746 */ 747 for (k=j-1; k > 0; k--) 748 if (zero_crossing[i].crossings[k] != 0) 749 break; 750 if (k < 0) 751 k=0; 752 /* 753 Check center for an even number of crossings between k and j. 754 */ 755 correct=(-1); 756 if (zero_crossing[i+1].crossings[j] != 0) 757 { 758 count=0; 759 for (l=k+1; l < center; l++) 760 if (zero_crossing[i+1].crossings[l] != 0) 761 count++; 762 if (((count % 2) == 0) && (center != k)) 763 correct=center; 764 } 765 /* 766 Check left for an even number of crossings between k and j. 767 */ 768 if (correct == -1) 769 { 770 count=0; 771 for (l=k+1; l < left; l++) 772 if (zero_crossing[i+1].crossings[l] != 0) 773 count++; 774 if (((count % 2) == 0) && (left != k)) 775 correct=left; 776 } 777 /* 778 Check right for an even number of crossings between k and j. 779 */ 780 if (correct == -1) 781 { 782 count=0; 783 for (l=k+1; l < right; l++) 784 if (zero_crossing[i+1].crossings[l] != 0) 785 count++; 786 if (((count % 2) == 0) && (right != k)) 787 correct=right; 788 } 789 l=(ssize_t) zero_crossing[i].crossings[j]; 790 zero_crossing[i].crossings[j]=0; 791 if (correct != -1) 792 zero_crossing[i].crossings[correct]=(short) l; 793 } 794 } 795 796 /* 798 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 799 % % 800 % % 801 % % 802 + D e f i n e R e g i o n % 803 % % 804 % % 805 % % 806 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 807 % 808 % DefineRegion() defines the left and right boundaries of a peak region. 809 % 810 % The format of the DefineRegion method is: 811 % 812 % ssize_t DefineRegion(const short *extrema,ExtentPacket *extents) 813 % 814 % A description of each parameter follows. 815 % 816 % o extrema: Specifies a pointer to an array of integers. They 817 % represent the peaks and valleys of the histogram for each color 818 % component. 819 % 820 % o extents: This pointer to an ExtentPacket represent the extends 821 % of a particular peak or valley of a color component. 822 % 823 */ 824 static ssize_t DefineRegion(const short *extrema,ExtentPacket *extents) 825 { 826 /* 827 Initialize to default values. 828 */ 829 extents->left=0; 830 extents->center=0.0; 831 extents->right=255; 832 /* 833 Find the left side (maxima). 834 */ 835 for ( ; extents->index <= 255; extents->index++) 836 if (extrema[extents->index] > 0) 837 break; 838 if (extents->index > 255) 839 return(MagickFalse); /* no left side - no region exists */ 840 extents->left=extents->index; 841 /* 842 Find the right side (minima). 843 */ 844 for ( ; extents->index <= 255; extents->index++) 845 if (extrema[extents->index] < 0) 846 break; 847 extents->right=extents->index-1; 848 return(MagickTrue); 849 } 850 851 /* 853 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 854 % % 855 % % 856 % % 857 + D e r i v a t i v e H i s t o g r a m % 858 % % 859 % % 860 % % 861 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 862 % 863 % DerivativeHistogram() determines the derivative of the histogram using 864 % central differencing. 865 % 866 % The format of the DerivativeHistogram method is: 867 % 868 % DerivativeHistogram(const double *histogram, 869 % double *derivative) 870 % 871 % A description of each parameter follows. 872 % 873 % o histogram: Specifies an array of doubles representing the number 874 % of pixels for each intensity of a particular color component. 875 % 876 % o derivative: This array of doubles is initialized by 877 % DerivativeHistogram to the derivative of the histogram using central 878 % differencing. 879 % 880 */ 881 static void DerivativeHistogram(const double *histogram, 882 double *derivative) 883 { 884 register ssize_t 885 i, 886 n; 887 888 /* 889 Compute endpoints using second order polynomial interpolation. 890 */ 891 n=255; 892 derivative[0]=(-1.5*histogram[0]+2.0*histogram[1]-0.5*histogram[2]); 893 derivative[n]=(0.5*histogram[n-2]-2.0*histogram[n-1]+1.5*histogram[n]); 894 /* 895 Compute derivative using central differencing. 896 */ 897 for (i=1; i < n; i++) 898 derivative[i]=(histogram[i+1]-histogram[i-1])/2.0; 899 return; 900 } 901 902 /* 904 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 905 % % 906 % % 907 % % 908 + G e t I m a g e D y n a m i c T h r e s h o l d % 909 % % 910 % % 911 % % 912 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 913 % 914 % GetImageDynamicThreshold() returns the dynamic threshold for an image. 915 % 916 % The format of the GetImageDynamicThreshold method is: 917 % 918 % MagickBooleanType GetImageDynamicThreshold(const Image *image, 919 % const double cluster_threshold,const double smooth_threshold, 920 % PixelInfo *pixel,ExceptionInfo *exception) 921 % 922 % A description of each parameter follows. 923 % 924 % o image: the image. 925 % 926 % o cluster_threshold: This double represents the minimum number of 927 % pixels contained in a hexahedra before it can be considered valid 928 % (expressed as a percentage). 929 % 930 % o smooth_threshold: the smoothing threshold eliminates noise in the second 931 % derivative of the histogram. As the value is increased, you can expect a 932 % smoother second derivative. 933 % 934 % o pixel: return the dynamic threshold here. 935 % 936 % o exception: return any errors or warnings in this structure. 937 % 938 */ 939 MagickExport MagickBooleanType GetImageDynamicThreshold(const Image *image, 940 const double cluster_threshold,const double smooth_threshold, 941 PixelInfo *pixel,ExceptionInfo *exception) 942 { 943 Cluster 944 *background, 945 *cluster, 946 *object, 947 *head, 948 *last_cluster, 949 *next_cluster; 950 951 ExtentPacket 952 blue, 953 green, 954 red; 955 956 MagickBooleanType 957 proceed; 958 959 double 960 threshold; 961 962 register const Quantum 963 *p; 964 965 register ssize_t 966 i, 967 x; 968 969 short 970 *extrema[MaxDimension]; 971 972 ssize_t 973 count, 974 *histogram[MaxDimension], 975 y; 976 977 /* 978 Allocate histogram and extrema. 979 */ 980 assert(image != (Image *) NULL); 981 assert(image->signature == MagickCoreSignature); 982 if (image->debug != MagickFalse) 983 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 984 GetPixelInfo(image,pixel); 985 for (i=0; i < MaxDimension; i++) 986 { 987 histogram[i]=(ssize_t *) AcquireQuantumMemory(256UL,sizeof(**histogram)); 988 extrema[i]=(short *) AcquireQuantumMemory(256UL,sizeof(**histogram)); 989 if ((histogram[i] == (ssize_t *) NULL) || (extrema[i] == (short *) NULL)) 990 { 991 for (i-- ; i >= 0; i--) 992 { 993 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]); 994 histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]); 995 } 996 (void) ThrowMagickException(exception,GetMagickModule(), 997 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename); 998 return(MagickFalse); 999 } 1000 } 1001 /* 1002 Initialize histogram. 1003 */ 1004 InitializeHistogram(image,histogram,exception); 1005 (void) OptimalTau(histogram[Red],Tau,0.2f,DeltaTau, 1006 (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Red]); 1007 (void) OptimalTau(histogram[Green],Tau,0.2f,DeltaTau, 1008 (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Green]); 1009 (void) OptimalTau(histogram[Blue],Tau,0.2f,DeltaTau, 1010 (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Blue]); 1011 /* 1012 Form clusters. 1013 */ 1014 cluster=(Cluster *) NULL; 1015 head=(Cluster *) NULL; 1016 (void) ResetMagickMemory(&red,0,sizeof(red)); 1017 (void) ResetMagickMemory(&green,0,sizeof(green)); 1018 (void) ResetMagickMemory(&blue,0,sizeof(blue)); 1019 while (DefineRegion(extrema[Red],&red) != 0) 1020 { 1021 green.index=0; 1022 while (DefineRegion(extrema[Green],&green) != 0) 1023 { 1024 blue.index=0; 1025 while (DefineRegion(extrema[Blue],&blue) != 0) 1026 { 1027 /* 1028 Allocate a new class. 1029 */ 1030 if (head != (Cluster *) NULL) 1031 { 1032 cluster->next=(Cluster *) AcquireMagickMemory( 1033 sizeof(*cluster->next)); 1034 cluster=cluster->next; 1035 } 1036 else 1037 { 1038 cluster=(Cluster *) AcquireMagickMemory(sizeof(*cluster)); 1039 head=cluster; 1040 } 1041 if (cluster == (Cluster *) NULL) 1042 { 1043 (void) ThrowMagickException(exception,GetMagickModule(), 1044 ResourceLimitError,"MemoryAllocationFailed","`%s'", 1045 image->filename); 1046 return(MagickFalse); 1047 } 1048 /* 1049 Initialize a new class. 1050 */ 1051 cluster->count=0; 1052 cluster->red=red; 1053 cluster->green=green; 1054 cluster->blue=blue; 1055 cluster->next=(Cluster *) NULL; 1056 } 1057 } 1058 } 1059 if (head == (Cluster *) NULL) 1060 { 1061 /* 1062 No classes were identified-- create one. 1063 */ 1064 cluster=(Cluster *) AcquireMagickMemory(sizeof(*cluster)); 1065 if (cluster == (Cluster *) NULL) 1066 { 1067 (void) ThrowMagickException(exception,GetMagickModule(), 1068 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename); 1069 return(MagickFalse); 1070 } 1071 /* 1072 Initialize a new class. 1073 */ 1074 cluster->count=0; 1075 cluster->red=red; 1076 cluster->green=green; 1077 cluster->blue=blue; 1078 cluster->next=(Cluster *) NULL; 1079 head=cluster; 1080 } 1081 /* 1082 Count the pixels for each cluster. 1083 */ 1084 count=0; 1085 for (y=0; y < (ssize_t) image->rows; y++) 1086 { 1087 p=GetVirtualPixels(image,0,y,image->columns,1,exception); 1088 if (p == (const Quantum *) NULL) 1089 break; 1090 for (x=0; x < (ssize_t) image->columns; x++) 1091 { 1092 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next) 1093 if (((ssize_t) ScaleQuantumToChar(GetPixelRed(image,p)) >= 1094 (cluster->red.left-SafeMargin)) && 1095 ((ssize_t) ScaleQuantumToChar(GetPixelRed(image,p)) <= 1096 (cluster->red.right+SafeMargin)) && 1097 ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,p)) >= 1098 (cluster->green.left-SafeMargin)) && 1099 ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,p)) <= 1100 (cluster->green.right+SafeMargin)) && 1101 ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,p)) >= 1102 (cluster->blue.left-SafeMargin)) && 1103 ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,p)) <= 1104 (cluster->blue.right+SafeMargin))) 1105 { 1106 /* 1107 Count this pixel. 1108 */ 1109 count++; 1110 cluster->red.center+=(double) ScaleQuantumToChar( 1111 GetPixelRed(image,p)); 1112 cluster->green.center+=(double) ScaleQuantumToChar( 1113 GetPixelGreen(image,p)); 1114 cluster->blue.center+=(double) ScaleQuantumToChar( 1115 GetPixelBlue(image,p)); 1116 cluster->count++; 1117 break; 1118 } 1119 p+=GetPixelChannels(image); 1120 } 1121 proceed=SetImageProgress(image,SegmentImageTag,(MagickOffsetType) y, 1122 2*image->rows); 1123 if (proceed == MagickFalse) 1124 break; 1125 } 1126 /* 1127 Remove clusters that do not meet minimum cluster threshold. 1128 */ 1129 count=0; 1130 last_cluster=head; 1131 next_cluster=head; 1132 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster) 1133 { 1134 next_cluster=cluster->next; 1135 if ((cluster->count > 0) && 1136 (cluster->count >= (count*cluster_threshold/100.0))) 1137 { 1138 /* 1139 Initialize cluster. 1140 */ 1141 cluster->id=count; 1142 cluster->red.center/=cluster->count; 1143 cluster->green.center/=cluster->count; 1144 cluster->blue.center/=cluster->count; 1145 count++; 1146 last_cluster=cluster; 1147 continue; 1148 } 1149 /* 1150 Delete cluster. 1151 */ 1152 if (cluster == head) 1153 head=next_cluster; 1154 else 1155 last_cluster->next=next_cluster; 1156 cluster=(Cluster *) RelinquishMagickMemory(cluster); 1157 } 1158 object=head; 1159 background=head; 1160 if (count > 1) 1161 { 1162 object=head->next; 1163 for (cluster=object; cluster->next != (Cluster *) NULL; ) 1164 { 1165 if (cluster->count < object->count) 1166 object=cluster; 1167 cluster=cluster->next; 1168 } 1169 background=head->next; 1170 for (cluster=background; cluster->next != (Cluster *) NULL; ) 1171 { 1172 if (cluster->count > background->count) 1173 background=cluster; 1174 cluster=cluster->next; 1175 } 1176 } 1177 if (background != (Cluster *) NULL) 1178 { 1179 threshold=(background->red.center+object->red.center)/2.0; 1180 pixel->red=(double) ScaleCharToQuantum((unsigned char) 1181 (threshold+0.5)); 1182 threshold=(background->green.center+object->green.center)/2.0; 1183 pixel->green=(double) ScaleCharToQuantum((unsigned char) 1184 (threshold+0.5)); 1185 threshold=(background->blue.center+object->blue.center)/2.0; 1186 pixel->blue=(double) ScaleCharToQuantum((unsigned char) 1187 (threshold+0.5)); 1188 } 1189 /* 1190 Relinquish resources. 1191 */ 1192 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster) 1193 { 1194 next_cluster=cluster->next; 1195 cluster=(Cluster *) RelinquishMagickMemory(cluster); 1196 } 1197 for (i=0; i < MaxDimension; i++) 1198 { 1199 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]); 1200 histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]); 1201 } 1202 return(MagickTrue); 1203 } 1204 1205 /* 1207 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1208 % % 1209 % % 1210 % % 1211 + I n i t i a l i z e H i s t o g r a m % 1212 % % 1213 % % 1214 % % 1215 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1216 % 1217 % InitializeHistogram() computes the histogram for an image. 1218 % 1219 % The format of the InitializeHistogram method is: 1220 % 1221 % InitializeHistogram(const Image *image,ssize_t **histogram) 1222 % 1223 % A description of each parameter follows. 1224 % 1225 % o image: Specifies a pointer to an Image structure; returned from 1226 % ReadImage. 1227 % 1228 % o histogram: Specifies an array of integers representing the number 1229 % of pixels for each intensity of a particular color component. 1230 % 1231 */ 1232 static void InitializeHistogram(const Image *image,ssize_t **histogram, 1233 ExceptionInfo *exception) 1234 { 1235 register const Quantum 1236 *p; 1237 1238 register ssize_t 1239 i, 1240 x; 1241 1242 ssize_t 1243 y; 1244 1245 /* 1246 Initialize histogram. 1247 */ 1248 for (i=0; i <= 255; i++) 1249 { 1250 histogram[Red][i]=0; 1251 histogram[Green][i]=0; 1252 histogram[Blue][i]=0; 1253 } 1254 for (y=0; y < (ssize_t) image->rows; y++) 1255 { 1256 p=GetVirtualPixels(image,0,y,image->columns,1,exception); 1257 if (p == (const Quantum *) NULL) 1258 break; 1259 for (x=0; x < (ssize_t) image->columns; x++) 1260 { 1261 histogram[Red][(ssize_t) ScaleQuantumToChar(GetPixelRed(image,p))]++; 1262 histogram[Green][(ssize_t) ScaleQuantumToChar(GetPixelGreen(image,p))]++; 1263 histogram[Blue][(ssize_t) ScaleQuantumToChar(GetPixelBlue(image,p))]++; 1264 p+=GetPixelChannels(image); 1265 } 1266 } 1267 } 1268 1269 /* 1271 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1272 % % 1273 % % 1274 % % 1275 + I n i t i a l i z e I n t e r v a l T r e e % 1276 % % 1277 % % 1278 % % 1279 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1280 % 1281 % InitializeIntervalTree() initializes an interval tree from the lists of 1282 % zero crossings. 1283 % 1284 % The format of the InitializeIntervalTree method is: 1285 % 1286 % InitializeIntervalTree(IntervalTree **list,ssize_t *number_nodes, 1287 % IntervalTree *node) 1288 % 1289 % A description of each parameter follows. 1290 % 1291 % o zero_crossing: Specifies an array of structures of type ZeroCrossing. 1292 % 1293 % o number_crossings: This size_t specifies the number of elements 1294 % in the zero_crossing array. 1295 % 1296 */ 1297 1298 static void InitializeList(IntervalTree **list,ssize_t *number_nodes, 1299 IntervalTree *node) 1300 { 1301 if (node == (IntervalTree *) NULL) 1302 return; 1303 if (node->child == (IntervalTree *) NULL) 1304 list[(*number_nodes)++]=node; 1305 InitializeList(list,number_nodes,node->sibling); 1306 InitializeList(list,number_nodes,node->child); 1307 } 1308 1309 static void MeanStability(IntervalTree *node) 1310 { 1311 register IntervalTree 1312 *child; 1313 1314 if (node == (IntervalTree *) NULL) 1315 return; 1316 node->mean_stability=0.0; 1317 child=node->child; 1318 if (child != (IntervalTree *) NULL) 1319 { 1320 register ssize_t 1321 count; 1322 1323 register double 1324 sum; 1325 1326 sum=0.0; 1327 count=0; 1328 for ( ; child != (IntervalTree *) NULL; child=child->sibling) 1329 { 1330 sum+=child->stability; 1331 count++; 1332 } 1333 node->mean_stability=sum/(double) count; 1334 } 1335 MeanStability(node->sibling); 1336 MeanStability(node->child); 1337 } 1338 1339 static void Stability(IntervalTree *node) 1340 { 1341 if (node == (IntervalTree *) NULL) 1342 return; 1343 if (node->child == (IntervalTree *) NULL) 1344 node->stability=0.0; 1345 else 1346 node->stability=node->tau-(node->child)->tau; 1347 Stability(node->sibling); 1348 Stability(node->child); 1349 } 1350 1351 static IntervalTree *InitializeIntervalTree(const ZeroCrossing *zero_crossing, 1352 const size_t number_crossings) 1353 { 1354 IntervalTree 1355 *head, 1356 **list, 1357 *node, 1358 *root; 1359 1360 register ssize_t 1361 i; 1362 1363 ssize_t 1364 j, 1365 k, 1366 left, 1367 number_nodes; 1368 1369 /* 1370 Allocate interval tree. 1371 */ 1372 list=(IntervalTree **) AcquireQuantumMemory((size_t) TreeLength, 1373 sizeof(*list)); 1374 if (list == (IntervalTree **) NULL) 1375 return((IntervalTree *) NULL); 1376 /* 1377 The root is the entire histogram. 1378 */ 1379 root=(IntervalTree *) AcquireMagickMemory(sizeof(*root)); 1380 root->child=(IntervalTree *) NULL; 1381 root->sibling=(IntervalTree *) NULL; 1382 root->tau=0.0; 1383 root->left=0; 1384 root->right=255; 1385 for (i=(-1); i < (ssize_t) number_crossings; i++) 1386 { 1387 /* 1388 Initialize list with all nodes with no children. 1389 */ 1390 number_nodes=0; 1391 InitializeList(list,&number_nodes,root); 1392 /* 1393 Split list. 1394 */ 1395 for (j=0; j < number_nodes; j++) 1396 { 1397 head=list[j]; 1398 left=head->left; 1399 node=head; 1400 for (k=head->left+1; k < head->right; k++) 1401 { 1402 if (zero_crossing[i+1].crossings[k] != 0) 1403 { 1404 if (node == head) 1405 { 1406 node->child=(IntervalTree *) AcquireMagickMemory( 1407 sizeof(*node->child)); 1408 node=node->child; 1409 } 1410 else 1411 { 1412 node->sibling=(IntervalTree *) AcquireMagickMemory( 1413 sizeof(*node->sibling)); 1414 node=node->sibling; 1415 } 1416 node->tau=zero_crossing[i+1].tau; 1417 node->child=(IntervalTree *) NULL; 1418 node->sibling=(IntervalTree *) NULL; 1419 node->left=left; 1420 node->right=k; 1421 left=k; 1422 } 1423 } 1424 if (left != head->left) 1425 { 1426 node->sibling=(IntervalTree *) AcquireMagickMemory( 1427 sizeof(*node->sibling)); 1428 node=node->sibling; 1429 node->tau=zero_crossing[i+1].tau; 1430 node->child=(IntervalTree *) NULL; 1431 node->sibling=(IntervalTree *) NULL; 1432 node->left=left; 1433 node->right=head->right; 1434 } 1435 } 1436 } 1437 /* 1438 Determine the stability: difference between a nodes tau and its child. 1439 */ 1440 Stability(root->child); 1441 MeanStability(root->child); 1442 list=(IntervalTree **) RelinquishMagickMemory(list); 1443 return(root); 1444 } 1445 1446 /* 1448 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1449 % % 1450 % % 1451 % % 1452 + O p t i m a l T a u % 1453 % % 1454 % % 1455 % % 1456 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1457 % 1458 % OptimalTau() finds the optimal tau for each band of the histogram. 1459 % 1460 % The format of the OptimalTau method is: 1461 % 1462 % double OptimalTau(const ssize_t *histogram,const double max_tau, 1463 % const double min_tau,const double delta_tau, 1464 % const double smooth_threshold,short *extrema) 1465 % 1466 % A description of each parameter follows. 1467 % 1468 % o histogram: Specifies an array of integers representing the number 1469 % of pixels for each intensity of a particular color component. 1470 % 1471 % o extrema: Specifies a pointer to an array of integers. They 1472 % represent the peaks and valleys of the histogram for each color 1473 % component. 1474 % 1475 */ 1476 1477 static void ActiveNodes(IntervalTree **list,ssize_t *number_nodes, 1478 IntervalTree *node) 1479 { 1480 if (node == (IntervalTree *) NULL) 1481 return; 1482 if (node->stability >= node->mean_stability) 1483 { 1484 list[(*number_nodes)++]=node; 1485 ActiveNodes(list,number_nodes,node->sibling); 1486 } 1487 else 1488 { 1489 ActiveNodes(list,number_nodes,node->sibling); 1490 ActiveNodes(list,number_nodes,node->child); 1491 } 1492 } 1493 1494 static void FreeNodes(IntervalTree *node) 1495 { 1496 if (node == (IntervalTree *) NULL) 1497 return; 1498 FreeNodes(node->sibling); 1499 FreeNodes(node->child); 1500 node=(IntervalTree *) RelinquishMagickMemory(node); 1501 } 1502 1503 static double OptimalTau(const ssize_t *histogram,const double max_tau, 1504 const double min_tau,const double delta_tau,const double smooth_threshold, 1505 short *extrema) 1506 { 1507 IntervalTree 1508 **list, 1509 *node, 1510 *root; 1511 1512 MagickBooleanType 1513 peak; 1514 1515 double 1516 average_tau, 1517 *derivative, 1518 *second_derivative, 1519 tau, 1520 value; 1521 1522 register ssize_t 1523 i, 1524 x; 1525 1526 size_t 1527 count, 1528 number_crossings; 1529 1530 ssize_t 1531 index, 1532 j, 1533 k, 1534 number_nodes; 1535 1536 ZeroCrossing 1537 *zero_crossing; 1538 1539 /* 1540 Allocate interval tree. 1541 */ 1542 list=(IntervalTree **) AcquireQuantumMemory((size_t) TreeLength, 1543 sizeof(*list)); 1544 if (list == (IntervalTree **) NULL) 1545 return(0.0); 1546 /* 1547 Allocate zero crossing list. 1548 */ 1549 count=(size_t) ((max_tau-min_tau)/delta_tau)+2; 1550 zero_crossing=(ZeroCrossing *) AcquireQuantumMemory((size_t) count, 1551 sizeof(*zero_crossing)); 1552 if (zero_crossing == (ZeroCrossing *) NULL) 1553 return(0.0); 1554 for (i=0; i < (ssize_t) count; i++) 1555 zero_crossing[i].tau=(-1.0); 1556 /* 1557 Initialize zero crossing list. 1558 */ 1559 derivative=(double *) AcquireQuantumMemory(256,sizeof(*derivative)); 1560 second_derivative=(double *) AcquireQuantumMemory(256, 1561 sizeof(*second_derivative)); 1562 if ((derivative == (double *) NULL) || 1563 (second_derivative == (double *) NULL)) 1564 ThrowFatalException(ResourceLimitFatalError, 1565 "UnableToAllocateDerivatives"); 1566 i=0; 1567 for (tau=max_tau; tau >= min_tau; tau-=delta_tau) 1568 { 1569 zero_crossing[i].tau=tau; 1570 ScaleSpace(histogram,tau,zero_crossing[i].histogram); 1571 DerivativeHistogram(zero_crossing[i].histogram,derivative); 1572 DerivativeHistogram(derivative,second_derivative); 1573 ZeroCrossHistogram(second_derivative,smooth_threshold, 1574 zero_crossing[i].crossings); 1575 i++; 1576 } 1577 /* 1578 Add an entry for the original histogram. 1579 */ 1580 zero_crossing[i].tau=0.0; 1581 for (j=0; j <= 255; j++) 1582 zero_crossing[i].histogram[j]=(double) histogram[j]; 1583 DerivativeHistogram(zero_crossing[i].histogram,derivative); 1584 DerivativeHistogram(derivative,second_derivative); 1585 ZeroCrossHistogram(second_derivative,smooth_threshold, 1586 zero_crossing[i].crossings); 1587 number_crossings=(size_t) i; 1588 derivative=(double *) RelinquishMagickMemory(derivative); 1589 second_derivative=(double *) 1590 RelinquishMagickMemory(second_derivative); 1591 /* 1592 Ensure the scale-space fingerprints form lines in scale-space, not loops. 1593 */ 1594 ConsolidateCrossings(zero_crossing,number_crossings); 1595 /* 1596 Force endpoints to be included in the interval. 1597 */ 1598 for (i=0; i <= (ssize_t) number_crossings; i++) 1599 { 1600 for (j=0; j < 255; j++) 1601 if (zero_crossing[i].crossings[j] != 0) 1602 break; 1603 zero_crossing[i].crossings[0]=(-zero_crossing[i].crossings[j]); 1604 for (j=255; j > 0; j--) 1605 if (zero_crossing[i].crossings[j] != 0) 1606 break; 1607 zero_crossing[i].crossings[255]=(-zero_crossing[i].crossings[j]); 1608 } 1609 /* 1610 Initialize interval tree. 1611 */ 1612 root=InitializeIntervalTree(zero_crossing,number_crossings); 1613 if (root == (IntervalTree *) NULL) 1614 return(0.0); 1615 /* 1616 Find active nodes: stability is greater (or equal) to the mean stability of 1617 its children. 1618 */ 1619 number_nodes=0; 1620 ActiveNodes(list,&number_nodes,root->child); 1621 /* 1622 Initialize extrema. 1623 */ 1624 for (i=0; i <= 255; i++) 1625 extrema[i]=0; 1626 for (i=0; i < number_nodes; i++) 1627 { 1628 /* 1629 Find this tau in zero crossings list. 1630 */ 1631 k=0; 1632 node=list[i]; 1633 for (j=0; j <= (ssize_t) number_crossings; j++) 1634 if (zero_crossing[j].tau == node->tau) 1635 k=j; 1636 /* 1637 Find the value of the peak. 1638 */ 1639 peak=zero_crossing[k].crossings[node->right] == -1 ? MagickTrue : 1640 MagickFalse; 1641 index=node->left; 1642 value=zero_crossing[k].histogram[index]; 1643 for (x=node->left; x <= node->right; x++) 1644 { 1645 if (peak != MagickFalse) 1646 { 1647 if (zero_crossing[k].histogram[x] > value) 1648 { 1649 value=zero_crossing[k].histogram[x]; 1650 index=x; 1651 } 1652 } 1653 else 1654 if (zero_crossing[k].histogram[x] < value) 1655 { 1656 value=zero_crossing[k].histogram[x]; 1657 index=x; 1658 } 1659 } 1660 for (x=node->left; x <= node->right; x++) 1661 { 1662 if (index == 0) 1663 index=256; 1664 if (peak != MagickFalse) 1665 extrema[x]=(short) index; 1666 else 1667 extrema[x]=(short) (-index); 1668 } 1669 } 1670 /* 1671 Determine the average tau. 1672 */ 1673 average_tau=0.0; 1674 for (i=0; i < number_nodes; i++) 1675 average_tau+=list[i]->tau; 1676 average_tau/=(double) number_nodes; 1677 /* 1678 Relinquish resources. 1679 */ 1680 FreeNodes(root); 1681 zero_crossing=(ZeroCrossing *) RelinquishMagickMemory(zero_crossing); 1682 list=(IntervalTree **) RelinquishMagickMemory(list); 1683 return(average_tau); 1684 } 1685 1686 /* 1688 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1689 % % 1690 % % 1691 % % 1692 + S c a l e S p a c e % 1693 % % 1694 % % 1695 % % 1696 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1697 % 1698 % ScaleSpace() performs a scale-space filter on the 1D histogram. 1699 % 1700 % The format of the ScaleSpace method is: 1701 % 1702 % ScaleSpace(const ssize_t *histogram,const double tau, 1703 % double *scale_histogram) 1704 % 1705 % A description of each parameter follows. 1706 % 1707 % o histogram: Specifies an array of doubles representing the number 1708 % of pixels for each intensity of a particular color component. 1709 % 1710 */ 1711 1712 static void ScaleSpace(const ssize_t *histogram,const double tau, 1713 double *scale_histogram) 1714 { 1715 double 1716 alpha, 1717 beta, 1718 *gamma, 1719 sum; 1720 1721 register ssize_t 1722 u, 1723 x; 1724 1725 gamma=(double *) AcquireQuantumMemory(256,sizeof(*gamma)); 1726 if (gamma == (double *) NULL) 1727 ThrowFatalException(ResourceLimitFatalError, 1728 "UnableToAllocateGammaMap"); 1729 alpha=PerceptibleReciprocal(tau*sqrt(2.0*MagickPI)); 1730 beta=(-1.0*PerceptibleReciprocal(2.0*tau*tau)); 1731 for (x=0; x <= 255; x++) 1732 gamma[x]=0.0; 1733 for (x=0; x <= 255; x++) 1734 { 1735 gamma[x]=exp((double) beta*x*x); 1736 if (gamma[x] < MagickEpsilon) 1737 break; 1738 } 1739 for (x=0; x <= 255; x++) 1740 { 1741 sum=0.0; 1742 for (u=0; u <= 255; u++) 1743 sum+=(double) histogram[u]*gamma[MagickAbsoluteValue(x-u)]; 1744 scale_histogram[x]=alpha*sum; 1745 } 1746 gamma=(double *) RelinquishMagickMemory(gamma); 1747 } 1748 1749 /* 1751 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1752 % % 1753 % % 1754 % % 1755 % S e g m e n t I m a g e % 1756 % % 1757 % % 1758 % % 1759 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1760 % 1761 % SegmentImage() segment an image by analyzing the histograms of the color 1762 % components and identifying units that are homogeneous with the fuzzy 1763 % C-means technique. 1764 % 1765 % The format of the SegmentImage method is: 1766 % 1767 % MagickBooleanType SegmentImage(Image *image, 1768 % const ColorspaceType colorspace,const MagickBooleanType verbose, 1769 % const double cluster_threshold,const double smooth_threshold, 1770 % ExceptionInfo *exception) 1771 % 1772 % A description of each parameter follows. 1773 % 1774 % o image: the image. 1775 % 1776 % o colorspace: Indicate the colorspace. 1777 % 1778 % o verbose: Set to MagickTrue to print detailed information about the 1779 % identified classes. 1780 % 1781 % o cluster_threshold: This represents the minimum number of pixels 1782 % contained in a hexahedra before it can be considered valid (expressed 1783 % as a percentage). 1784 % 1785 % o smooth_threshold: the smoothing threshold eliminates noise in the second 1786 % derivative of the histogram. As the value is increased, you can expect a 1787 % smoother second derivative. 1788 % 1789 % o exception: return any errors or warnings in this structure. 1790 % 1791 */ 1792 MagickExport MagickBooleanType SegmentImage(Image *image, 1793 const ColorspaceType colorspace,const MagickBooleanType verbose, 1794 const double cluster_threshold,const double smooth_threshold, 1795 ExceptionInfo *exception) 1796 { 1797 ColorspaceType 1798 previous_colorspace; 1799 1800 MagickBooleanType 1801 status; 1802 1803 register ssize_t 1804 i; 1805 1806 short 1807 *extrema[MaxDimension]; 1808 1809 ssize_t 1810 *histogram[MaxDimension]; 1811 1812 /* 1813 Allocate histogram and extrema. 1814 */ 1815 assert(image != (Image *) NULL); 1816 assert(image->signature == MagickCoreSignature); 1817 if (image->debug != MagickFalse) 1818 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); 1819 for (i=0; i < MaxDimension; i++) 1820 { 1821 histogram[i]=(ssize_t *) AcquireQuantumMemory(256,sizeof(**histogram)); 1822 extrema[i]=(short *) AcquireQuantumMemory(256,sizeof(**extrema)); 1823 if ((histogram[i] == (ssize_t *) NULL) || (extrema[i] == (short *) NULL)) 1824 { 1825 for (i-- ; i >= 0; i--) 1826 { 1827 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]); 1828 histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]); 1829 } 1830 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed", 1831 image->filename) 1832 } 1833 } 1834 /* 1835 Initialize histogram. 1836 */ 1837 previous_colorspace=image->colorspace; 1838 (void) TransformImageColorspace(image,colorspace,exception); 1839 InitializeHistogram(image,histogram,exception); 1840 (void) OptimalTau(histogram[Red],Tau,0.2,DeltaTau, 1841 smooth_threshold == 0.0 ? 1.0 : smooth_threshold,extrema[Red]); 1842 (void) OptimalTau(histogram[Green],Tau,0.2,DeltaTau, 1843 smooth_threshold == 0.0 ? 1.0 : smooth_threshold,extrema[Green]); 1844 (void) OptimalTau(histogram[Blue],Tau,0.2,DeltaTau, 1845 smooth_threshold == 0.0 ? 1.0 : smooth_threshold,extrema[Blue]); 1846 /* 1847 Classify using the fuzzy c-Means technique. 1848 */ 1849 status=Classify(image,extrema,cluster_threshold,WeightingExponent,verbose, 1850 exception); 1851 (void) TransformImageColorspace(image,previous_colorspace,exception); 1852 /* 1853 Relinquish resources. 1854 */ 1855 for (i=0; i < MaxDimension; i++) 1856 { 1857 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]); 1858 histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]); 1859 } 1860 return(status); 1861 } 1862 1863 /* 1865 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1866 % % 1867 % % 1868 % % 1869 + Z e r o C r o s s H i s t o g r a m % 1870 % % 1871 % % 1872 % % 1873 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1874 % 1875 % ZeroCrossHistogram() find the zero crossings in a histogram and marks 1876 % directions as: 1 is negative to positive; 0 is zero crossing; and -1 1877 % is positive to negative. 1878 % 1879 % The format of the ZeroCrossHistogram method is: 1880 % 1881 % ZeroCrossHistogram(double *second_derivative, 1882 % const double smooth_threshold,short *crossings) 1883 % 1884 % A description of each parameter follows. 1885 % 1886 % o second_derivative: Specifies an array of doubles representing the 1887 % second derivative of the histogram of a particular color component. 1888 % 1889 % o crossings: This array of integers is initialized with 1890 % -1, 0, or 1 representing the slope of the first derivative of the 1891 % of a particular color component. 1892 % 1893 */ 1894 static void ZeroCrossHistogram(double *second_derivative, 1895 const double smooth_threshold,short *crossings) 1896 { 1897 register ssize_t 1898 i; 1899 1900 ssize_t 1901 parity; 1902 1903 /* 1904 Merge low numbers to zero to help prevent noise. 1905 */ 1906 for (i=0; i <= 255; i++) 1907 if ((second_derivative[i] < smooth_threshold) && 1908 (second_derivative[i] >= -smooth_threshold)) 1909 second_derivative[i]=0.0; 1910 /* 1911 Mark zero crossings. 1912 */ 1913 parity=0; 1914 for (i=0; i <= 255; i++) 1915 { 1916 crossings[i]=0; 1917 if (second_derivative[i] < 0.0) 1918 { 1919 if (parity > 0) 1920 crossings[i]=(-1); 1921 parity=1; 1922 } 1923 else 1924 if (second_derivative[i] > 0.0) 1925 { 1926 if (parity < 0) 1927 crossings[i]=1; 1928 parity=(-1); 1929 } 1930 } 1931 } 1932
