Import libaries, read image, define common functions, use real division.


In [1]:
%run common.ipynb
from __future__ import division


Populating the interactive namespace from numpy and matplotlib

Work out an algoritm with test data.


In [2]:
a = np.array([[0, 1, 0],
              [2, 1, 1],
              [4, 2, 5]])
h = bincount(a.flatten())
pdf = h/a.size
cdf = pdf.cumsum()
plot(h, label='h')
plot(pdf, label='pdf')
plot(cdf, label='cdf')
title('Original image histogram')
legend()
b = np.zeros(a.shape, dtype=np.uint8)
for level in range(a.max()):
    mask = (a == level+1) # exclude level=0 (zero->zero)
    b[mask] = cdf[level+1] * 10 # 10 is highest new level
#print a, '\n', b
bh = bincount(b.flatten())
bpdf = bh/b.size
bcdf = bpdf.cumsum()
figure()
title('Converted image histogram')
plot(bh, label='bh')
plot(bpdf, label='bpdf')
plot(bcdf, label='bcdf')
legend()


Out[2]:
<matplotlib.legend.Legend at 0x7f37f4139dd0>

In [3]:
print pdf


[ 0.22222222  0.33333333  0.22222222  0.          0.11111111  0.11111111]

Now we do the same with the image.


In [4]:
hist = np.bincount(img.flatten())
pdf = hist / img.size
cdf = pdf.cumsum()
plot(pdf, label='pdf')
plot(cdf, label='cdf')
legend()


Out[4]:
<matplotlib.legend.Legend at 0x7f37f3f4bc50>

In [5]:
equal_img = np.zeros(img.shape, dtype=np.uint8)
for level in range(1,256): # exclude zeroes
    if hist[level] == 0: continue # short cut
    mask = (img == level)
    equal_img[mask] = cdf[level] * 255

In [6]:
equal_hist = np.bincount(equal_img.flatten())
plot(equal_hist)


Out[6]:
[<matplotlib.lines.Line2D at 0x7f37f3e84250>]

In [7]:
figure(figsize=(12,12))
imshow(img, cmap='gray')
figure(figsize=(12,12))
imshow(equal_img, cmap='gray')


Out[7]:
<matplotlib.image.AxesImage at 0x7f37f3d3ba50>

Quite a lot of noise, lets compare reduction of noise before / after equalization.


In [8]:
re_img = smooth(equal_img)
er_img = equalization(smooth(img))
gimshow(re_img, title='Noise reduction before equalization')
showhist(re_img)
gimshow(er_img, title='Noise reduction after equalization')
showhist(er_img)


Lets also try noise reduction both before and after.


In [9]:
rr_img = smooth(equalization(smooth(img)))
gimshow(rr_img)
showhist(rr_img)