# Demo iaotsudemo

Illustrate the Otsu Thresholding Selection Method

# Script

Gray scale image and its histogram

```1 from ia636 import iahistogram
2
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)
```