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)21 void 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)39 void 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)60 void 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)93 void 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)182 void 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)248 void 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)371 void 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)504 void 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)601 void 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