Merge pull request #178 from ahojnnes/adaptive-thresholding

ENH: Adaptive thresholding.

CONTRIBUTORS.txt

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

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, threshold_adaptive

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, threshold_adaptive
 
 
 class TestSimpleImage():
@@ -24,6 +25,52 @@ class TestSimpleImage():
         image = np.float64(self.image)
         assert 2 <= threshold_otsu(image) < 3
 
+    def test_threshold_adaptive_generic(self):
+        def func(arr):
+            return arr.sum() / arr.shape[0]
+        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 = threshold_adaptive(self.image, 3, method='generic', param=func)
+        assert_array_equal(ref, out)
+
+    def test_threshold_adaptive_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 = threshold_adaptive(self.image, 3, method='gaussian')
+        assert_array_equal(ref, out)
+
+    def test_threshold_adaptive_mean(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 = threshold_adaptive(self.image, 3, method='mean')
+        assert_array_equal(ref, out)
+
+    def test_threshold_adaptive_median(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 = threshold_adaptive(self.image, 3, method='median')
+        assert_array_equal(ref, out)
+
 
 def test_otsu_camera_image():
     assert threshold_otsu(data.camera()) == 87

skimage/filter/thresholding.py

@@ -1,11 +1,93 @@
 import numpy as np
-
+import scipy.ndimage
 from skimage.exposure import histogram
+from ._thresholding import _threshold_adaptive
 
 
-__all__ = ['threshold_otsu']
+__all__ = ['threshold_otsu', 'threshold_adaptive']
 
 
+def threshold_adaptive(image, block_size, method='gaussian', offset=0,
+                       mode='reflect', param=None):
+    """Applies an adaptive threshold to an array.
+
+    Also known as local or dynamic thresholding where the threshold value is the
+    weighted mean for the local neighborhood of a pixel subtracted by a
+    constant. Alternatively the threshold can be determined dynamically by a
+    a given function using the 'generic' method.
+
+    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, ...).
+    method : {'generic', 'gaussian', 'mean', 'median'}, optional
+        Method used to determine adaptive threshold for local neighbourhood in
+        weighted mean image.
+         * 'generic': use custom function (see `param` parameter)
+         * 'gaussian': apply gaussian filter (see `param` parameter for custom
+                       sigma value)
+         * 'mean': apply arithmetic mean filter
+         * 'median' apply median rank filter
+        By default the 'gaussian' method is used.
+    offset : float, optional
+        Constant subtracted from weighted mean of neighborhood to calculate
+        the local threshold value. Default offset is 0.
+    mode : {'reflect', 'constant', 'nearest', 'mirror', 'wrap'}, optional
+        The mode parameter determines how the array borders are handled, where
+        cval is the value when mode is equal to 'constant'.
+        Default is 'reflect'.
+    param : {int, function}, optional
+        Either specify sigma for 'gaussian' method or function object for
+        'generic' method. This functions takes the flat array of local
+        neighbourhood as a single argument and returns the calculated threshold
+        for the centre pixel.
+
+    Returns
+    -------
+    threshold : NxM ndarray
+        Thresholded binary image
+
+    References
+    ----------
+    http://docs.opencv.org/modules/imgproc/doc/miscellaneous_transformations
+        .html?highlight=threshold#adaptivethreshold
+
+    Examples
+    --------
+    >>> from skimage.data import camera
+    >>> image = camera()
+    >>> binary_image1 = threshold_adaptive(image, 15, 'mean')
+    >>> func = lambda arr: arr.mean()
+    >>> binary_image2 = threshold_adaptive(image, 15, 'generic', param=func)
+    """
+    thresh_image = np.zeros(image.shape, 'double')
+    if method == 'generic':
+        scipy.ndimage.generic_filter(image, param, block_size,
+            output=thresh_image, mode=mode)
+    elif method == 'gaussian':
+        if param is None:
+            # automatically determine sigma which covers > 99% of distribution
+            sigma = (block_size - 1) / 6.0
+        else:
+            sigma = param
+        scipy.ndimage.gaussian_filter(image, sigma, output=thresh_image,
+            mode=mode)
+    elif method == 'mean':
+        mask = 1. / block_size * np.ones((block_size,))
+        # separation of filters to speedup convolution
+        scipy.ndimage.convolve1d(image, mask, axis=0, output=thresh_image,
+            mode=mode)
+        scipy.ndimage.convolve1d(thresh_image, mask, axis=1,
+            output=thresh_image, mode=mode)
+    elif method == 'median':
+        scipy.ndimage.median_filter(image, block_size, output=thresh_image,
+            mode=mode)
+
+    return image > (thresh_image - offset)
+
 def threshold_otsu(image, nbins=256):
     """Return threshold value based on Otsu's method.
 
@@ -51,4 +133,3 @@ def threshold_otsu(image, nbins=256):
     idx = np.argmax(variance12)
     threshold = bin_centers[:-1][idx]
     return threshold
-