1 # include "bitmapConfig.h" 2 3 # include <math.h> 4 5 # include "bmGrayHisto.h" 6 7 # include <appDebugon.h> 8 9 /************************************************************************/ 10 /* */ 11 /* Various gray scale histogram base Thresholding algorithms. */ 12 /* */ 13 /* Most of the approaches and some of the implementation is derived */ 14 /* from the following book: */ 15 /* */ 16 /* PARKER, J.R.: Algorithms for Image Processing and Computer Vision, */ 17 /* Wiley Computer Publishing, New York, 1996, ISBN 0-471-14056-2 */ 18 /* */ 19 /************************************************************************/ 20 bmInitThresholderHistogram(ThresholderHistogram * th)21void bmInitThresholderHistogram( ThresholderHistogram * th ) 22 { 23 int i; 24 25 for ( i= 0; i < 256; i++ ) 26 { th->thHistogram[i]= 0; } 27 28 th->thHistogramSize= 0; 29 th->thPixelCount= 0; 30 th->thThreshold= -1; 31 } 32 33 /************************************************************************/ 34 /* */ 35 /* Determine the threshold using the Mean pixel value. */ 36 /* */ 37 /************************************************************************/ 38 bmThresholdMean(ThresholderHistogram * th)39void bmThresholdMean( ThresholderHistogram * th ) 40 { 41 int i; 42 43 double sr; 44 45 sr= 0; 46 for ( i= 0; i < th->thHistogramSize; i++ ) 47 { sr += th->thHistogram[i]* i; } 48 49 th->thThreshold= sr/ th->thPixelCount; 50 51 return; 52 } 53 54 /************************************************************************/ 55 /* */ 56 /* Determine the threshold using a fraction. */ 57 /* */ 58 /************************************************************************/ 59 bmThresholdQuantile(ThresholderHistogram * th,double frac)60void bmThresholdQuantile( ThresholderHistogram * th, 61 double frac ) 62 { 63 int i; 64 65 double sr; 66 67 sr= 0; 68 for ( i= 0; i < th->thHistogramSize; i++ ) 69 { 70 if ( ( sr + th->thHistogram[i] )/ th->thPixelCount > frac ) 71 { break; } 72 73 sr += th->thHistogram[i]; 74 } 75 76 th->thThreshold= i; 77 78 return; 79 } 80 81 /************************************************************************/ 82 /* */ 83 /* Determine the threshold using the two peak method. */ 84 /* */ 85 /* 1) Determine the primary peak. */ 86 /* 2) Determine the second peak. */ 87 /* 3) Assume lowest value is background. */ 88 /* 4) Guess the threshold value is exactly between the peaks. */ 89 /* 5) But look for a valley between the two peaks. */ 90 /* */ 91 /************************************************************************/ 92 bmThreshold2Peaks(ThresholderHistogram * th)93void bmThreshold2Peaks( ThresholderHistogram * th ) 94 { 95 int i0; 96 int i1; 97 98 int p1= 0; 99 int p2= 1; 100 101 long v; 102 long v1= 0; 103 long v2= 0; 104 105 int i; 106 int p; 107 108 /* 1 */ 109 i0= i1= 0; 110 for ( i= 0; i < th->thHistogramSize; i++ ) 111 { 112 v= th->thHistogram[i]; 113 114 if ( v1 < v ) 115 { v1= v; i0= i1= i; } 116 if ( v1 == v ) 117 { i1= i; } 118 } 119 120 p1= ( i0+ i1 )/ 2; 121 122 /* 2 */ 123 for ( i= 0; i < th->thHistogramSize; i++ ) 124 { 125 v= ( i- p1 ); 126 v= v* v* th->thHistogram[i]; 127 128 if ( v2 < v ) 129 { v2= v; i0= i1= i; } 130 if ( v2 == v ) 131 { i1= i; } 132 } 133 134 p2= ( i0+ i1 )/ 2; 135 136 /* 3 */ 137 if ( p1 > p2 ) 138 { i= p1; p1= p2; p2= i; } 139 140 /* 4 */ 141 th->thThreshold= ( p1+ p2 )/ 2; 142 143 /* 5 */ 144 i0= i1= th->thThreshold; 145 v1= th->thHistogram[th->thThreshold]; 146 147 i0= i1= p1; 148 v1= th->thHistogram[p1]; 149 p= ( p1+ p2 )/ 2; 150 151 for ( i= p1; i <= p2; i++ ) 152 { 153 v= ( i- p1 )* ( p2- i ); 154 v= v* ( p- th->thHistogram[i] ); 155 156 if ( v1 < v ) 157 { v1= v; i0= i1= i; } 158 if ( v1 == v ) 159 { i1= i; } 160 } 161 162 th->thThreshold= ( i0+ i1 )/ 2; 163 164 return; 165 } 166 167 /************************************************************************/ 168 /* */ 169 /* Determine the threshold using the Ridler Iterative Selection */ 170 /* method. */ 171 /* */ 172 /* 1) First guess: the average pixel value; */ 173 /* 2) Determine the mean gray level of pixels not over the threshold */ 174 /* 3) Determine the mean gray level of pixels over the threshold */ 175 /* 4) A new estimate of the threshold is the average of the means */ 176 /* calculated above. */ 177 /* 5) Iteration stops if the new threshold is equal to the previous */ 178 /* one. */ 179 /* */ 180 /************************************************************************/ 181 bmThresholdRidler(ThresholderHistogram * th)182void bmThresholdRidler( ThresholderHistogram * th ) 183 { 184 int i; 185 186 double sr; 187 int tr; 188 189 /* 1 */ 190 sr= 0; 191 for ( i= 0; i < th->thHistogramSize; i++ ) 192 { sr += th->thHistogram[i]* i; } 193 tr= sr/ th->thPixelCount; 194 195 for (;;) 196 { 197 int no= 0; 198 int nb= 0; 199 double so= 0; 200 double sb= 0; 201 202 int tn; 203 204 /* 2 */ 205 for ( i= 0; i <= tr; i++ ) 206 { 207 nb += th->thHistogram[i]; 208 sb += th->thHistogram[i]* i; 209 } 210 211 /* 3 */ 212 for ( i= tr+ 1; i < th->thHistogramSize; i++ ) 213 { 214 no += th->thHistogram[i]; 215 so += th->thHistogram[i]* i; 216 } 217 218 if ( no == 0 ) 219 { no= 1; } 220 if ( nb == 0 ) 221 { nb= 1; } 222 223 /* 4 */ 224 tn= ( so/no+ sb/nb )/ 2.0; 225 226 /* 5 */ 227 if ( tn == tr ) 228 { 229 th->thThreshold= tr; 230 return; 231 } 232 233 tr= tn; 234 } 235 } 236 237 /************************************************************************/ 238 /* */ 239 /* Determine the threshold minimizing the ration of between- class and */ 240 /* overall variance. */ 241 /* */ 242 /* 1) Overall mean and variance. */ 243 /* 2) Below threshold. */ 244 /* 3) Above threshold. */ 245 /* */ 246 /************************************************************************/ 247 bmThresholdVariance(ThresholderHistogram * th)248void bmThresholdVariance( ThresholderHistogram * th ) 249 { 250 int i; 251 252 double sum_b; 253 double ssq_b; 254 double ub; 255 double vb; 256 int nb; 257 258 double sum_o; 259 double ssq_o; 260 double uo; 261 double vo; 262 int no; 263 264 int tr; 265 double vmin; 266 267 int tr0; 268 int tr1; 269 270 # if 0 271 272 double sum_t; 273 double ssq_t; 274 double ut; 275 double vt; 276 277 /* 1 */ 278 sum_t= 0; ssq_t= 0; 279 for ( i= 0; i < th->thHistogramSize; i++ ) 280 { 281 sum_t += (double)th->thHistogram[i]* i; 282 ssq_t += (double)th->thHistogram[i]* i* i; 283 } 284 ut= sum_t/ th->thPixelCount; 285 vt= ssq_t/ th->thPixelCount- ut* ut; 286 287 # endif 288 289 tr= 0; 290 291 /* 2 */ 292 vb= 0; ub= 0; 293 sum_b= 0; ssq_b= 0; nb= 0; 294 for ( i= 0; i <= tr; i++ ) 295 { 296 nb += th->thHistogram[i]; 297 sum_b += (double)th->thHistogram[i]* i; 298 ssq_b += (double)th->thHistogram[i]* i* i; 299 } 300 if ( nb > 0 ) 301 { 302 ub= sum_b/ nb; 303 vb= ssq_b/ nb- ub* ub; 304 } 305 306 /* 3 */ 307 vo= 0; uo= 0; 308 sum_o= 0; ssq_o= 0; no= 0; 309 for ( i= tr+ 1; i < th->thHistogramSize; i++ ) 310 { 311 no += th->thHistogram[i]; 312 sum_o += (double)th->thHistogram[i]* i; 313 ssq_o += (double)th->thHistogram[i]* i* i; 314 } 315 if ( no > 0 ) 316 { 317 uo= sum_o/ no; 318 vo= ssq_o/ no- uo* uo; 319 } 320 321 vmin= vo+ vb; 322 tr0= tr1= tr; 323 324 for ( tr= 1; tr < th->thHistogramSize; tr++ ) 325 { 326 /* 2 */ 327 vb= 0; ub= 0; 328 nb += th->thHistogram[tr]; 329 sum_b += (double)th->thHistogram[tr]* tr; 330 ssq_b += (double)th->thHistogram[tr]* tr* tr; 331 332 if ( nb > 0 ) 333 { 334 ub= sum_b/ nb; 335 vb= ssq_b/ nb- ub* ub; 336 } 337 338 /* 3 */ 339 vo= 0; uo= 0; 340 no -= th->thHistogram[tr]; 341 sum_o -= (double)th->thHistogram[tr]* tr; 342 ssq_o -= (double)th->thHistogram[tr]* tr* tr; 343 344 if ( no > 0 ) 345 { 346 uo= sum_o/ no; 347 vo= ssq_o/ no- uo* uo; 348 } 349 350 if ( vb+ vo < vmin ) 351 { tr0= tr1= tr; vmin= vb+ vo; } 352 if ( vb+ vo == vmin ) 353 { tr1= tr; } 354 } 355 356 th->thThreshold= ( tr0+ tr1 )/ 2; 357 358 return; 359 } 360 361 /************************************************************************/ 362 /* */ 363 /* Determine the threshold minimizing the error function. */ 364 /* */ 365 /* 1) Overall mean and variance. */ 366 /* 2) Below threshold. */ 367 /* 3) Above threshold. */ 368 /* */ 369 /************************************************************************/ 370 bmThresholdMinimumError(ThresholderHistogram * th)371void bmThresholdMinimumError( ThresholderHistogram * th ) 372 { 373 int i; 374 375 double sum_b; 376 double ssq_b; 377 double ub; 378 double vb; 379 double sb; 380 int nb; 381 382 double sum_o; 383 double ssq_o; 384 double uo; 385 double vo; 386 double so; 387 int no; 388 389 int tr; 390 double j; 391 double jmin= 1e+20; 392 393 int tr0; 394 int tr1; 395 396 int s; 397 int p; 398 399 s= 0; while( th->thHistogram[s] == 0 ) 400 { s++; } 401 p= th->thHistogramSize; while( th->thHistogram[p-1] == 0 ) 402 { p--; } 403 404 # if 0 405 406 double sum_t; 407 double ssq_t; 408 double ut; 409 double vt; 410 411 /* 1 */ 412 sum_t= 0; ssq_t= 0; 413 for ( i= 0; i < th->thHistogramSize; i++ ) 414 { 415 sum_t += (double)th->thHistogram[i]* i; 416 ssq_t += (double)th->thHistogram[i]* i* i; 417 } 418 ut= sum_t/ th->thPixelCount; 419 vt= ssq_t/ th->thPixelCount- ut* ut; 420 421 # endif 422 423 tr= s+ 1; 424 425 /* 2 */ 426 sb= 0.0; vb= 0; ub= 0; 427 sum_b= 0; ssq_b= 0; nb= 0; 428 for ( i= s; i <= tr; i++ ) 429 { 430 nb += th->thHistogram[i]; 431 sum_b += (double)th->thHistogram[i]* i; 432 ssq_b += (double)th->thHistogram[i]* i* i; 433 } 434 435 /* 3 */ 436 so= 0.0; vo= 0; uo= 0; 437 sum_o= 0; ssq_o= 0; no= 0; 438 for ( i= tr+ 1; i < p; i++ ) 439 { 440 no += th->thHistogram[i]; 441 sum_o += (double)th->thHistogram[i]* i; 442 ssq_o += (double)th->thHistogram[i]* i* i; 443 } 444 445 jmin= 1.0e+20; 446 tr0= tr1= tr; 447 448 for ( tr= tr+ 1; tr < p- 1; tr++ ) 449 { 450 /* 2 */ 451 sb= 0.0; vb= 0; ub= 0; 452 nb += th->thHistogram[tr]; 453 sum_b += (double)th->thHistogram[tr]* tr; 454 ssq_b += (double)th->thHistogram[tr]* tr* tr; 455 456 if ( nb > 0 ) 457 { 458 ub= sum_b/ nb; 459 vb= ssq_b/ nb- ub* ub; 460 sb= sqrt( vb ); 461 } 462 463 /* 3 */ 464 so= 0.0; vo= 0; uo= 0; 465 no -= th->thHistogram[tr]; 466 sum_o -= (double)th->thHistogram[tr]* tr; 467 ssq_o -= (double)th->thHistogram[tr]* tr* tr; 468 469 if ( no == 0 ) 470 { continue; } 471 472 if ( no > 0 ) 473 { 474 uo= sum_o/ no; 475 vo= ssq_o/ no- uo* uo; 476 so= sqrt( vo ); 477 } 478 479 if ( nb == 0 || vb == 0 || vo == 0 ) 480 { continue; } 481 482 j= 1.0+ 483 2* ( nb* log( sb )+ no* log( so ) )- 484 2* ( nb* log( 1.0*nb )+ no* log( 1.0*no ) ); 485 486 if ( j < jmin ) 487 { tr0= tr1= tr; jmin= j; } 488 if ( j == jmin ) 489 { tr1= tr; } 490 } 491 492 th->thThreshold= ( tr0+ tr1 )/ 2; 493 494 return; 495 } 496 497 498 /************************************************************************/ 499 /* */ 500 /* Determine the threshold using the Kapur maximal entropy method. */ 501 /* */ 502 /************************************************************************/ 503 bmThresholdKapur(ThresholderHistogram * const th)504void bmThresholdKapur( ThresholderHistogram * const th ) 505 { 506 int i; 507 508 double p; 509 double Pt[256]; 510 511 # if 0 512 double Ht= 0; 513 # endif 514 515 int tr; 516 int tr0; 517 int tr1; 518 double N= th->thPixelCount; 519 520 double Hmax; 521 522 double Hb; 523 double Ho; 524 525 p= 0; 526 for ( i= 0; i < th->thHistogramSize; i++ ) 527 { 528 p += (double)th->thHistogram[i]; 529 Pt[i]= p/ N; 530 531 # if 0 532 if ( th->thHistogram[i] > 0 && Pt[i] > 0 ) 533 { 534 double r= th->thHistogram[i]/ ( N* Pt[i] ); 535 536 Ht -= r* log( r ); 537 } 538 # endif 539 } 540 541 tr= 0; 542 543 Hb= 0; 544 for ( i= 0; i <= tr; i++ ) 545 { 546 if ( th->thHistogram[i] > 0 && Pt[i] > 0 ) 547 { 548 double r= th->thHistogram[i]/ ( N* Pt[i] ); 549 550 Hb -= r* log( r ); 551 } 552 } 553 554 Ho= 0; 555 for ( i= tr+ 1; i < th->thHistogramSize; i++ ) 556 { 557 if ( th->thHistogram[i] > 0 && Pt[i] < 1 ) 558 { 559 double r= th->thHistogram[i]/ ( N* ( 1- Pt[i] ) ); 560 561 Ho -= r* log( r ); 562 } 563 } 564 565 Hmax= Hb+ Ho; 566 tr0= tr1= tr; 567 568 for ( tr= 1; tr < th->thHistogramSize; tr++ ) 569 { 570 if ( th->thHistogram[tr] > 0 && Pt[tr] > 0 ) 571 { 572 double r= th->thHistogram[tr]/ ( N* Pt[tr] ); 573 574 Hb -= r* log( r ); 575 } 576 577 if ( th->thHistogram[tr] > 0 && Pt[tr] < 1 ) 578 { 579 double r= th->thHistogram[tr]/ ( N* ( 1- Pt[tr] ) ); 580 581 Ho += r* log( r ); 582 } 583 584 if ( Hb+ Ho > Hmax ) 585 { tr0= tr1= tr; Hmax= Hb+ Ho; } 586 if ( Hb+ Ho == Hmax ) 587 { tr1= tr; } 588 } 589 590 th->thThreshold= ( tr0+ tr1 )/ 2; 591 592 return; 593 } 594 595 /************************************************************************/ 596 /* */ 597 /* Determine the threshold using the Johannsen minimal entropy method. */ 598 /* */ 599 /************************************************************************/ 600 bmThresholdJohannsen(ThresholderHistogram * th)601void bmThresholdJohannsen( ThresholderHistogram * th ) 602 { 603 double Pt[256]; 604 double Pq[256]; 605 606 int tr; 607 double N= th->thPixelCount; 608 609 int e; 610 int s; 611 612 double Smin= N; 613 614 { 615 int i; 616 double p; 617 618 p= 0; 619 Pt[0]= (double)th->thHistogram[0]/N; 620 Pq[0]= 1- Pt[0]; 621 for ( i= 1; i < th->thHistogramSize; i++ ) 622 { 623 p += (double)th->thHistogram[i]; 624 Pt[i]= p/ N; 625 Pq[i]= 1.0- Pt[i-1]; 626 } 627 } 628 629 s= 0; e= 255; 630 while( th->thHistogram[s] == 0 ) 631 { s++; } 632 while( th->thHistogram[e] == 0 ) 633 { e--; } 634 635 s++; e--; 636 637 for ( tr= s; tr <= e; tr++ ) 638 { 639 double Sb; 640 double So; 641 double p; 642 double Ep; 643 644 double Epr= 0; 645 double Enx= 0; 646 647 if ( th->thHistogram[tr] == 0 ) 648 { continue; } 649 650 p= (double)th->thHistogram[tr]/ N; 651 Ep= -p* log( p ); 652 653 if ( Pt[tr-1] > 0.0 ) 654 { Epr= -Pt[tr-1]* log( Pt[tr-1] ); } 655 Sb= log( Pt[tr] )+ ( 1.0/ Pt[tr] )* ( Ep+ Epr ); 656 657 if ( Pq[tr+1] > 0.0 ) 658 { Enx= -Pq[tr+1]* log( Pq[tr+1] ); } 659 So= log( Pq[tr] )+ ( 1.0/ Pq[tr] )* ( Ep+ Enx ); 660 661 if ( Sb+ So < Smin ) 662 { 663 th->thThreshold= tr; 664 Smin= Sb+ So; 665 } 666 } 667 668 return; 669 } 670 671