Demo iaotsudemo

Illustrate the Otsu Thresholding Selection Method

More information can be seen at wikipedia:Otsu's method

Script

Image read and histograming

Gray scale image and its histogram

1. from ia636 import iahistogram
2. 
3. f = adread('woodlog.pgm');
4. adshow(f)
5. H = iahistogram(f)
6. x = arange(len(H))
7. k = x[0:-1]
8. mmplot([[x, H]])

Normalized histogram

If the histogram is divided by the number of pixels in the image, it can be seen as a probability distribution. The sum of each values gives one. The mean gray level and the variance can be computed from the normalized histogram.

01. from ia636 import iaerror
02. import numpy
03. 
04. h = 1.*H/product(f.shape)
05. print sum(h)
06. mt = sum(x * h)
07. st2 = sum((x-mt)**2 * h)
08. 
09. if abs(mt - numpy.mean(ravel(f))) > 0.01: iaerror('error in computing mean')
10. 
11. if abs(st2 - numpy.std(ravel(f))**2) > 0.0001: iaerror('error in computing var')
12. 
13. print 'mean is %2.2f, var2 is %2.2f' % (mt, st2)
14. 
15. maux = 0 * h
16. maux[int(mt)] = max(h)
17. mmplot([[x,h]])
18. mmplot([[x,maux]])
1.0
mean is 91.03, var2 is 2873.86

Two classes

Suppose the pixels are categorized in two classes and : smaller or equal than the gray scale t and larger than t. The probability of the first class occurrence is the cummulative normalized histogram. The other class is the complementary class.

01. import numpy
02. 
03. w0 = numpy.cumsum(h[0:-1])
04. aux = h[1::]
05. w1aux = numpy.cumsum(aux[::-1])[::-1]
06. w1 = 1 - w0
07. if max(abs(w1-w1aux)) > 0.0001: iaerror('error in computing w1')
08. 
09. mmplot([[k,w0]])
10. mmplot([[k,w1]])

Mean gray level of each class

The mean gray level as a function of the thresholding t is computed and displayed below.

01. m0 = numpy.cumsum(k * h[0:-1]) / (1.*w0)
02. m1 = (mt - m0*w0)/w1
03. aux = (k+1) * h[1::]
04. m1x = numpy.cumsum(aux[::-1])[::-1] / (1.*w1)
05. mm = w0 * m0 + w1 * m1
06. if max(abs(m1-m1x)) > 0.0001: iaerror('error in computing m1')
07. 
08. mmplot([[k,m0]])
09. mmplot([[k,m1]])
10. mmplot([[k,mm]])

Variance of each class

The gray level variance as a function of the thresholding t is computed and displayed below.

1. s02 = numpy.cumsum((k-m0)**2 * h[0:-1]) / (1.*w0)
2. aux = ((k+1)-m1)**2 * h[1::]
3. s12 = numpy.cumsum(aux[::-1])[::-1] / (1.*w1)
4. mmplot([[k, s02]])
5. mmplot([[k, s12]])

Class separability

The variance between class is a good measure of class separability. As higher this variance, the better the class clustering.

1. sB2 = w0 * ((m0 - mt)**2) + w1 * ((m1 - mt)**2)
2. sBaux2 = w0 * w1 * ((m0 - m1)**2)
3. if max(sB2-sBaux2) > 0.0001: iaerror('error in computing sB')
4. 
5. t = numpy.argmax(sB2)
6. eta = 1.*sBaux2[t]/st2
7. print 'Optimum threshold at %f quality factor %f' % (t, eta)
8. mmplot([[k, sB2]])
9. mmplot([[k, H[0:-1]]])
Optimum threshold at 93.000000 quality factor 0.694320

Thresholding

The thresholded image is displayed to illustrate the result of the binarization using the Otsu method.

1. x_ = adshow(f > t)

See Also

  • iaotsu - Threshold using Otsu method.