added adaptive threshold

CONTRIBUTORS.txt

@@ -101,3 +101,4 @@
 
 - Johannes Schönberger
   Polygon, circle and ellipse drawing functions
+  Adaptive thresholding

doc/examples/plot_otsu.py

@@ -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()
-

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()
+
+

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

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')

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
 

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

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
-