From 238d4eb4adc2881f043791394a1e1285085fed45 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Tue, 24 Apr 2012 14:13:08 +0200 Subject: [PATCH] added adaptive threshold --- CONTRIBUTORS.txt | 1 + doc/examples/plot_otsu.py | 45 ------------------ doc/examples/plot_thresholding.py | 56 +++++++++++++++++++++++ skimage/filter/__init__.py | 2 +- skimage/filter/_thresholding.pyx | 28 ++++++++++++ skimage/filter/setup.py | 3 ++ skimage/filter/tests/test_thresholding.py | 25 +++++++++- skimage/filter/thresholding.py | 39 +++++++++++++++- 8 files changed, 150 insertions(+), 49 deletions(-) delete mode 100644 doc/examples/plot_otsu.py create mode 100644 doc/examples/plot_thresholding.py create mode 100644 skimage/filter/_thresholding.pyx diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 8258eada..a7a7e83e 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -101,3 +101,4 @@ - Johannes Schönberger Polygon, circle and ellipse drawing functions + Adaptive thresholding diff --git a/doc/examples/plot_otsu.py b/doc/examples/plot_otsu.py deleted file mode 100644 index f2335fc0..00000000 --- a/doc/examples/plot_otsu.py +++ /dev/null @@ -1,45 +0,0 @@ -""" -============ -Thresholding -============ - -Thresholding is used to create a binary image. This example uses Otsu's method -to calculate the threshold value. - -Otsu's method calculates an "optimal" threshold (marked by a red line in the -histogram below) by maximizing the variance between two classes of pixels, -which are separated by the threshold. Equivalently, this threshold minimizes -the intra-class variance. - -.. [1] http://en.wikipedia.org/wiki/Otsu's_method - -""" - -import matplotlib.pyplot as plt - -from skimage.data import camera -from skimage.filter import threshold_otsu - - -image = camera() -thresh = threshold_otsu(image) -binary = image > thresh - -plt.figure(figsize=(8, 2.5)) -plt.subplot(1, 3, 1) -plt.imshow(image, cmap=plt.cm.gray) -plt.title('Original') -plt.axis('off') - -plt.subplot(1, 3, 2, aspect='equal') -plt.hist(image) -plt.title('Histogram') -plt.axvline(thresh, color='r') - -plt.subplot(1, 3, 3) -plt.imshow(binary, cmap=plt.cm.gray) -plt.title('Thresholded') -plt.axis('off') - -plt.show() - diff --git a/doc/examples/plot_thresholding.py b/doc/examples/plot_thresholding.py new file mode 100644 index 00000000..e860f243 --- /dev/null +++ b/doc/examples/plot_thresholding.py @@ -0,0 +1,56 @@ +""" +============ +Thresholding +============ + +Thresholding is used to create a binary image. + +This example uses Otsu's method to calculate the threshold value. Otsu's method +calculates an "optimal" threshold (marked by a red line in the histogram below) +by maximizing the variance between two classes of pixels, which are separated by +the threshold. Equivalently, this threshold minimizes the intra-class variance. + +Additionnally an adaptive thresholding is applied. Also known as local or +dynamic thresholding where the the threshold value is the weighted mean for the +local neighborhood of a pixel subtracted by a constant. + +.. [1] http://en.wikipedia.org/wiki/Otsu's_method + +""" + +import matplotlib.pyplot as plt +import numpy as np + +from skimage.data import camera +from skimage.filter import threshold_otsu, adaptive_threshold + + +image = camera() +thresh = threshold_otsu(image) +otsu_binary = image > thresh +adaptive_binary = np.invert(adaptive_threshold(image, 9, 5)) + +plt.figure(figsize=(8, 2.5)) +plt.subplot(2, 2, 1) +plt.imshow(image, cmap=plt.cm.gray) +plt.title('Original') +plt.axis('off') + +plt.subplot(2, 2, 2, aspect='equal') +plt.hist(image) +plt.title('Histogram') +plt.axvline(thresh, color='r') + +plt.subplot(2, 2, 3) +plt.imshow(otsu_binary, cmap=plt.cm.gray) +plt.title('Thresholded with Otsu') +plt.axis('off') + +plt.subplot(2, 2, 4) +plt.imshow(adaptive_binary, cmap=plt.cm.gray) +plt.title('Adaptively thresholded') +plt.axis('off') + +plt.show() + + diff --git a/skimage/filter/__init__.py b/skimage/filter/__init__.py index 748487ac..8a9cffaa 100644 --- a/skimage/filter/__init__.py +++ b/skimage/filter/__init__.py @@ -4,4 +4,4 @@ from .canny import canny from .edges import sobel, hsobel, vsobel, hprewitt, vprewitt, prewitt from .tv_denoise import tv_denoise from .rank_order import rank_order -from .thresholding import threshold_otsu +from .thresholding import threshold_otsu, adaptive_threshold diff --git a/skimage/filter/_thresholding.pyx b/skimage/filter/_thresholding.pyx new file mode 100644 index 00000000..303c1a5f --- /dev/null +++ b/skimage/filter/_thresholding.pyx @@ -0,0 +1,28 @@ +import numpy as np +import scipy.ndimage +cimport numpy as np +cimport cython + + +@cython.boundscheck(False) +@cython.wraparound(False) +def _adaptive_threshold( + np.ndarray[np.double_t, ndim=2] image, + int block_size, + double offset, + method +): + cdef int r, c + cdef np.ndarray[np.float64_t, ndim=2] mean_image + if method == 'gaussian': + # covers > 99% of distribution + sigma = (block_size - 1) / 6.0 + mean_image = scipy.ndimage.gaussian_filter(image, sigma) + elif method == 'mean': + mean_image = scipy.ndimage.median_filter(image, block_size) + + for r in range(image.shape[0]): + for c in range(image.shape[1]): + mean_image[r,c] = image[r,c] > (mean_image[r,c] - offset) + + return mean_image.astype('bool') diff --git a/skimage/filter/setup.py b/skimage/filter/setup.py index 12cb84a7..f93f7143 100644 --- a/skimage/filter/setup.py +++ b/skimage/filter/setup.py @@ -12,9 +12,12 @@ def configuration(parent_package='', top_path=None): config.add_data_dir('tests') cython(['_ctmf.pyx'], working_path=base_path) + cython(['_thresholding.pyx'], working_path=base_path) config.add_extension('_ctmf', sources=['_ctmf.c'], include_dirs=[get_numpy_include_dirs()]) + config.add_extension('_thresholding', sources=['_thresholding.c'], + include_dirs=[get_numpy_include_dirs()]) return config diff --git a/skimage/filter/tests/test_thresholding.py b/skimage/filter/tests/test_thresholding.py index b0044a47..64eb48cc 100644 --- a/skimage/filter/tests/test_thresholding.py +++ b/skimage/filter/tests/test_thresholding.py @@ -1,8 +1,9 @@ import numpy as np +from numpy.testing import assert_array_equal import skimage from skimage import data -from skimage.filter.thresholding import threshold_otsu +from skimage.filter.thresholding import threshold_otsu, adaptive_threshold class TestSimpleImage(): @@ -24,6 +25,28 @@ class TestSimpleImage(): image = np.float64(self.image) assert 2 <= threshold_otsu(image) < 3 + def test_adaptive_threshold_gaussian(self): + ref = np.array( + [[False, False, False, False, True], + [False, False, True, False, True], + [False, False, True, True, False], + [False, True, True, False, False], + [ True, True, False, False, False]] + ) + out = adaptive_threshold(self.image, 3, 0, 'gaussian') + assert_array_equal(ref, out) + + def test_adaptive_threshold_mean(self): + ref = np.array( + [[False, False, False, False, True], + [False, False, True, False, False], + [False, False, True, False, False], + [False, False, True, True, False], + [False, True, False, False, False]] + ) + out = adaptive_threshold(self.image, 3, 0, 'mean') + assert_array_equal(ref, out) + def test_otsu_camera_image(): assert threshold_otsu(data.camera()) == 87 diff --git a/skimage/filter/thresholding.py b/skimage/filter/thresholding.py index 856bd78d..931aeaa2 100644 --- a/skimage/filter/thresholding.py +++ b/skimage/filter/thresholding.py @@ -1,11 +1,47 @@ import numpy as np from skimage.exposure import histogram +from ._thresholding import _adaptive_threshold -__all__ = ['threshold_otsu'] +__all__ = ['threshold_otsu', 'adaptive_threshold'] +def adaptive_threshold(image, block_size, offset, method='gaussian'): + """Applies an adaptive threshold to an array. + + Also known as local or dynamic thresholding where the the threshold value is + the weighted mean for the local neighborhood of a pixel subtracted by a + constant. + + Parameters + ---------- + image : NxM ndarray + Input image. + block_size : int + uneven size of pixel neighborhood which is used to calculate the + threshold value (e.g. 3, 5, 7, ..., 21, ...) + offset : float + constant subtracted from weighted mean of neighborhood to calculate + the local threshold value + method : string, optional + thresholding type which must be one of `gaussian` or `mean`. + By default the `gaussian` method is used. + + Returns + ------- + threshold : NxM ndarray + thresholded binary image + + References + ---------- + http://docs.opencv.org/modules/imgproc/doc/miscellaneous_transformations + .html?highlight=threshold#adaptivethreshold + """ + # not using img_as_float because threshold parameter wouldn't work + image = image.astype('double') + return _adaptive_threshold(image, block_size, offset, method) + def threshold_otsu(image, nbins=256): """Return threshold value based on Otsu's method. @@ -51,4 +87,3 @@ def threshold_otsu(image, nbins=256): idx = np.argmax(variance12) threshold = bin_centers[:-1][idx] return threshold -