Local Otsu ThresholdΒΆ

This example shows how Otsu’s threshold [1] method can be applied locally. For each pixel, an “optimal” threshold is determined by maximizing the variance between two classes of pixels of the local neighborhood defined by a structuring element.

The example compares the local threshold with the global threshold.

[1]http://en.wikipedia.org/wiki/Otsu’s_method
../_images/plot_local_otsu_1.png

import matplotlib
import matplotlib.pyplot as plt

from skimage import data
from skimage.morphology import disk
from skimage.filter import threshold_otsu, rank
from skimage.util import img_as_ubyte


matplotlib.rcParams['font.size'] = 9


img = img_as_ubyte(data.page())

radius = 15
selem = disk(radius)

local_otsu = rank.otsu(img, selem)
threshold_global_otsu = threshold_otsu(img)
global_otsu = img >= threshold_global_otsu


fig, ax = plt.subplots(2, 2, figsize=(8, 5))
ax1, ax2, ax3, ax4 = ax.ravel()

fig.colorbar(ax1.imshow(img, cmap=plt.cm.gray),
           ax=ax1, orientation='horizontal')
ax1.set_title('Original')
ax1.axis('off')

fig.colorbar(ax2.imshow(local_otsu, cmap=plt.cm.gray),
           ax=ax2, orientation='horizontal')
ax2.set_title('Local Otsu (radius=%d)' % radius)
ax2.axis('off')

ax3.imshow(img >= local_otsu, cmap=plt.cm.gray)
ax3.set_title('Original >= Local Otsu' % threshold_global_otsu)
ax3.axis('off')

ax4.imshow(global_otsu, cmap=plt.cm.gray)
ax4.set_title('Global Otsu (threshold = %d)' % threshold_global_otsu)
ax4.axis('off')

plt.show()

STDOUT


        

STDERR


        

Python source code: download (generated using skimage 0.11dev)

IPython Notebook: download (generated using skimage 0.11dev)

aW1wb3J0IG1hdHBsb3RsaWIKaW1wb3J0IG1hdHBsb3RsaWIucHlwbG90IGFzIHBsdAoKZnJvbSBza2ltYWdlIGltcG9ydCBkYXRhCmZyb20gc2tpbWFnZS5tb3JwaG9sb2d5IGltcG9ydCBkaXNrCmZyb20gc2tpbWFnZS5maWx0ZXIgaW1wb3J0IHRocmVzaG9sZF9vdHN1LCByYW5rCmZyb20gc2tpbWFnZS51dGlsIGltcG9ydCBpbWdfYXNfdWJ5dGUKCgptYXRwbG90bGliLnJjUGFyYW1zWydmb250LnNpemUnXSA9IDkKCgppbWcgPSBpbWdfYXNfdWJ5dGUoZGF0YS5wYWdlKCkpCgpyYWRpdXMgPSAxNQpzZWxlbSA9IGRpc2socmFkaXVzKQoKbG9jYWxfb3RzdSA9IHJhbmsub3RzdShpbWcsIHNlbGVtKQp0aHJlc2hvbGRfZ2xvYmFsX290c3UgPSB0aHJlc2hvbGRfb3RzdShpbWcpCmdsb2JhbF9vdHN1ID0gaW1nID49IHRocmVzaG9sZF9nbG9iYWxfb3RzdQoKCmZpZywgYXggPSBwbHQuc3VicGxvdHMoMiwgMiwgZmlnc2l6ZT0oOCwgNSkpCmF4MSwgYXgyLCBheDMsIGF4NCA9IGF4LnJhdmVsKCkKCmZpZy5jb2xvcmJhcihheDEuaW1zaG93KGltZywgY21hcD1wbHQuY20uZ3JheSksCiAgICAgICAgICAgYXg9YXgxLCBvcmllbnRhdGlvbj0naG9yaXpvbnRhbCcpCmF4MS5zZXRfdGl0bGUoJ09yaWdpbmFsJykKYXgxLmF4aXMoJ29mZicpCgpmaWcuY29sb3JiYXIoYXgyLmltc2hvdyhsb2NhbF9vdHN1LCBjbWFwPXBsdC5jbS5ncmF5KSwKICAgICAgICAgICBheD1heDIsIG9yaWVudGF0aW9uPSdob3Jpem9udGFsJykKYXgyLnNldF90aXRsZSgnTG9jYWwgT3RzdSAocmFkaXVzPSVkKScgJSByYWRpdXMpCmF4Mi5heGlzKCdvZmYnKQoKYXgzLmltc2hvdyhpbWcgPj0gbG9jYWxfb3RzdSwgY21hcD1wbHQuY20uZ3JheSkKYXgzLnNldF90aXRsZSgnT3JpZ2luYWwgPj0gTG9jYWwgT3RzdScgJSB0aHJlc2hvbGRfZ2xvYmFsX290c3UpCmF4My5heGlzKCdvZmYnKQoKYXg0Lmltc2hvdyhnbG9iYWxfb3RzdSwgY21hcD1wbHQuY20uZ3JheSkKYXg0LnNldF90aXRsZSgnR2xvYmFsIE90c3UgKHRocmVzaG9sZCA9ICVkKScgJSB0aHJlc2hvbGRfZ2xvYmFsX290c3UpCmF4NC5heGlzKCdvZmYnKQoKcGx0LnNob3coKQ==