Merge pull request #98 from tonysyu/thresholding

Add thresholding function using Otsu's method to filter module.

doc/examples/plot_otsu.py

@@ -0,0 +1,45 @@
+"""
+============
+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=(10, 3.5))
+plt.subplot(1, 3, 1)
+plt.imshow(image)
+plt.title('original')
+plt.axis('off')
+
+plt.subplot(1, 3, 2)
+plt.hist(image)
+plt.title('histogram')
+plt.axvline(thresh, color='r')
+
+plt.subplot(1, 3, 3)
+plt.imshow(binary)
+plt.title('thresholded')
+plt.axis('off')
+
+plt.show()
+

skimage/filter/__init__.py

@@ -4,3 +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

skimage/filter/tests/test_thresholding.py

@@ -0,0 +1,44 @@
+import numpy as np
+
+import skimage
+from skimage import data
+from skimage.filter.thresholding import threshold_otsu
+
+
+class TestSimpleImage():
+    def setup(self):
+        self.image = np.array([[0, 0, 1, 3, 5],
+                               [0, 1, 4, 3, 4],
+                               [1, 2, 5, 4, 1],
+                               [2, 4, 5, 2, 1],
+                               [4, 5, 1, 0, 0]], dtype=int)
+
+    def test_otsu(self):
+        assert threshold_otsu(self.image) == 2
+
+    def test_otsu_negative_int(self):
+        image = self.image - 2
+        assert threshold_otsu(image) == 0
+
+    def test_otsu_float_image(self):
+        image = np.float64(self.image)
+        assert 2 <= threshold_otsu(image) < 3
+
+
+def test_otsu_camera_image():
+    assert threshold_otsu(data.camera()) == 87
+
+def test_otsu_coins_image():
+    assert threshold_otsu(data.coins()) == 107
+
+def test_otsu_coins_image_as_float():
+    coins = skimage.img_as_float(data.coins())
+    assert 0.41 < threshold_otsu(coins) < 0.42
+
+def test_otsu_lena_image():
+    assert threshold_otsu(data.lena()) == 141
+
+
+if __name__ == '__main__':
+    np.testing.run_module_suite()
+

skimage/filter/thresholding.py

@@ -0,0 +1,89 @@
+import numpy as np
+
+
+__all__ = ['threshold_otsu']
+
+
+def threshold_otsu(image, nbins=256):
+    """Return threshold value based on Otsu's method.
+
+    Parameters
+    ----------
+    image : array
+        Input image.
+    nbins : int
+        Number of bins used to calculate histogram. This value is ignored for
+        integer arrays.
+
+    Returns
+    -------
+    threshold : float
+        Threshold value.
+
+    References
+    ----------
+    .. [1] Wikipedia, http://en.wikipedia.org/wiki/Otsu's_Method
+
+    Examples
+    --------
+    >>> from skimage.data import camera
+    >>> image = camera()
+    >>> thresh = threshold_otsu(image)
+    >>> binary = image > thresh
+    """
+    hist, bin_centers = histogram(image, nbins)
+    hist = hist.astype(float)
+
+    # class probabilities for all possible thresholds
+    weight1 = np.cumsum(hist)
+    weight2 = np.cumsum(hist[::-1])[::-1]
+    # class means for all possible thresholds
+    mean1 = np.cumsum(hist * bin_centers) / weight1
+    mean2 = (np.cumsum((hist * bin_centers)[::-1]) / weight2[::-1])[::-1]
+
+    # Clip ends to align class 1 and class 2 variables:
+    # The last value of `weight1`/`mean1` should pair with zero values in
+    # `weight2`/`mean2`, which do not exist.
+    variance12 = weight1[:-1] * weight2[1:] * (mean1[:-1] - mean2[1:])**2
+
+    idx = np.argmax(variance12)
+    threshold = bin_centers[:-1][idx]
+    return threshold
+
+
+def histogram(image, nbins):
+    """Return histogram of image.
+
+    Unlike `numpy.histogram`, this function returns the centers of bins and
+    does not rebin integer arrays.
+
+    Parameters
+    ----------
+    image : array
+        Input image.
+    nbins : int
+        Number of bins used to calculate histogram. This value is ignored for
+        integer arrays.
+
+    Returns
+    -------
+    hist : array
+        The values of the histogram.
+    bin_centers : array
+        The values at the center of the bins.
+    """
+    if np.issubdtype(image.dtype, np.integer):
+        offset = 0
+        if np.min(image) < 0:
+            offset = np.min(image)
+        hist = np.bincount(image.ravel() - offset)
+        bin_centers = np.arange(len(hist)) + offset
+
+        # clip histogram to return only non-zero bins
+        idx = np.nonzero(hist)[0][0]
+        return hist[idx:], bin_centers[idx:]
+    else:
+        hist, bin_edges = np.histogram(image, bins=nbins)
+        bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2.
+        return hist, bin_centers
+