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.

 1 from ia636 import iaerror
 2 import numpy
 3 
 4 h = 1.*H/product(f.shape)
 5 print sum(h)
 6 mt = sum(x * h)
 7 st2 = sum((x-mt)**2 * h)
 8 
 9 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.

 1 import numpy
 2 
 3 w0 = numpy.cumsum(h[0:-1])
 4 aux = h[1::]
 5 w1aux = numpy.cumsum(aux[::-1])[::-1]
 6 w1 = 1 - w0
 7 if max(abs(w1-w1aux)) > 0.0001: iaerror('error in computing w1')
 8 
 9 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.

 1 m0 = numpy.cumsum(k * h[0:-1]) / (1.*w0)
 2 m1 = (mt - m0*w0)/w1
 3 aux = (k+1) * h[1::]
 4 m1x = numpy.cumsum(aux[::-1])[::-1] / (1.*w1)
 5 mm = w0 * m0 + w1 * m1
 6 if max(abs(m1-m1x)) > 0.0001: iaerror('error in computing m1')
 7 
 8 mmplot([[k,m0]])
 9 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, sBaux2]])
10 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.