Merge pull request #859 from radioxoma/threshold

Add ISODATA thresholding

CONTRIBUTORS.txt

@@ -151,7 +151,7 @@
   Color difference functions
 
 - Eugene Dvoretsky
-  Yen threshold implementation.
+  Yen, Ridler-Calvard (ISODATA) threshold implementations.
 
 - Riaan van den Dool
   skimage.io plugin: GDAL

skimage/filter/__init__.py

@@ -9,7 +9,8 @@ from ._denoise import denoise_tv_chambolle
 from ._denoise_cy import denoise_bilateral, denoise_tv_bregman
 from ._rank_order import rank_order
 from ._gabor import gabor_kernel, gabor_filter
-from .thresholding import threshold_otsu, threshold_adaptive
+from .thresholding import (threshold_adaptive, threshold_otsu, threshold_yen,
+                           threshold_isodata)
 from . import rank
 
 
@@ -37,6 +38,8 @@ __all__ = ['inverse',
            'rank_order',
            'gabor_kernel',
            'gabor_filter',
-           'threshold_otsu',
            'threshold_adaptive',
+           'threshold_otsu',
+           'threshold_yen',
+           'threshold_isodata',
            'rank']

skimage/filter/tests/test_thresholding.py

@@ -5,7 +5,8 @@ import skimage
 from skimage import data
 from skimage.filter.thresholding import (threshold_adaptive,
                                          threshold_otsu,
-                                         threshold_yen)
+                                         threshold_yen,
+                                         threshold_isodata)
 
 
 class TestSimpleImage():
@@ -43,10 +44,29 @@ class TestSimpleImage():
         assert threshold_yen(image) == 127
 
     def test_yen_binary(self):
-        image = np.zeros([2,256], dtype='uint8')
+        image = np.zeros([2,256], dtype=np.uint8)
         image[0] = 255
         assert threshold_yen(image) < 1
 
+    def test_yen_blank_zero(self):
+        image = np.zeros((5, 5), dtype=np.uint8)
+        assert threshold_yen(image) == 0
+
+    def test_yen_blank_max(self):
+        image = np.empty((5, 5), dtype=np.uint8)
+        image.fill(255)
+        assert threshold_yen(image) == 255
+
+    def test_isodata(self):
+        assert threshold_isodata(self.image) == 2
+
+    def test_isodata_blank_zero(self):
+        image = np.zeros((5, 5), dtype=np.uint8)
+        assert threshold_isodata(image) == 0
+
+    def test_isodata_linspace(self):
+        assert -63.8 < threshold_isodata(np.linspace(-127, 0, 256)) < -63.6
+
     def test_threshold_adaptive_generic(self):
         def func(arr):
             return arr.sum() / arr.shape[0]
@@ -114,6 +134,11 @@ def test_otsu_lena_image():
     assert 140 < threshold_otsu(lena) < 142
 
 
+def test_yen_camera_image():
+    camera = skimage.img_as_ubyte(data.camera())
+    assert 197 < threshold_yen(camera) < 199
+
+
 def test_yen_coins_image():
     coins = skimage.img_as_ubyte(data.coins())
     assert 109 < threshold_yen(coins) < 111
@@ -124,5 +149,32 @@ def test_yen_coins_image_as_float():
     assert 0.43 < threshold_yen(coins) < 0.44
 
 
+def test_isodata_camera_image():
+    camera = skimage.img_as_ubyte(data.camera())
+    assert threshold_isodata(camera) == 88
+
+
+def test_isodata_coins_image():
+    coins = skimage.img_as_ubyte(data.coins())
+    assert threshold_isodata(coins) == 107
+
+
+def test_isodata_moon_image():
+    moon = skimage.img_as_ubyte(data.moon())
+    assert threshold_isodata(moon) == 87
+
+
+def test_isodata_moon_image_negative_int():
+    moon = skimage.img_as_ubyte(data.moon()).astype(np.int32)
+    moon -= 100
+    assert threshold_isodata(moon) == -13
+
+
+def test_isodata_moon_image_negative_float():
+    moon = skimage.img_as_ubyte(data.moon()).astype(np.float64)
+    moon -= 100
+    assert -13 < threshold_isodata(moon) < -12
+
+
 if __name__ == '__main__':
     np.testing.run_module_suite()

skimage/filter/thresholding.py

@@ -1,4 +1,7 @@
-__all__ = ['threshold_adaptive', 'threshold_otsu', 'threshold_yen']
+__all__ = ['threshold_adaptive',
+           'threshold_otsu',
+           'threshold_yen',
+           'threshold_isodata']
 
 import numpy as np
 import scipy.ndimage
@@ -129,7 +132,7 @@ def threshold_otsu(image, nbins=256):
     # 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
+    variance12 = weight1[:-1] * weight2[1:] * (mean1[:-1] - mean2[1:]) ** 2
 
     idx = np.argmax(variance12)
     threshold = bin_centers[:-1][idx]
@@ -172,15 +175,84 @@ def threshold_yen(image, nbins=256):
     >>> binary = image <= thresh
     """
     hist, bin_centers = histogram(image, nbins)
-    norm_histo = hist.astype(float) / hist.sum() # Probability mass function
-    P1 = np.cumsum(norm_histo) # Cumulative normalized histogram
-    P1_sq = np.cumsum(norm_histo ** 2)
+    # On blank images (e.g. filled with 0) with int dtype, `histogram()`
+    # returns `bin_centers` containing only one value. Speed up with it.
+    if bin_centers.size == 1:
+        return bin_centers[0]
+
+    # Calculate probability mass function
+    pmf = hist.astype(np.float32) / hist.sum()
+    P1 = np.cumsum(pmf)  # Cumulative normalized histogram
+    P1_sq = np.cumsum(pmf ** 2)
     # Get cumsum calculated from end of squared array:
-    P2_sq = np.cumsum(norm_histo[::-1] ** 2)[::-1]
+    P2_sq = np.cumsum(pmf[::-1] ** 2)[::-1]
     # P2_sq indexes is shifted +1. I assume, with P1[:-1] it's help avoid '-inf'
     # in crit. ImageJ Yen implementation replaces those values by zero.
-    crit = np.log(((P1_sq[:-1] * P2_sq[1:]) ** -1) * \
+    crit = np.log(((P1_sq[:-1] * P2_sq[1:]) ** -1) *
                   (P1[:-1] * (1.0 - P1[:-1])) ** 2)
-    max_crit = np.argmax(crit)
-    threshold = bin_centers[:-1][max_crit]
+    return bin_centers[crit.argmax()]
+
+
+def threshold_isodata(image, nbins=256):
+    """Return threshold value based on ISODATA method.
+
+    Histogram-based threshold, known as Ridler-Calvard method or intermeans.
+
+    Parameters
+    ----------
+    image : array
+        Input image.
+    nbins : int, optional
+        Number of bins used to calculate histogram. This value is ignored for
+        integer arrays.
+
+    Returns
+    -------
+    threshold : float or int, corresponding input array dtype.
+        Upper threshold value. All pixels intensities that less or equal of
+        this value assumed as background.
+
+    References
+    ----------
+    .. [1] Ridler, TW & Calvard, S (1978), "Picture thresholding using an
+           iterative selection method"
+    .. [2] IEEE Transactions on Systems, Man and Cybernetics 8: 630-632,
+           http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=4310039
+    .. [3] Sezgin M. and Sankur B. (2004) "Survey over Image Thresholding
+           Techniques and Quantitative Performance Evaluation" Journal of
+           Electronic Imaging, 13(1): 146-165,
+           http://www.busim.ee.boun.edu.tr/~sankur/SankurFolder/Threshold_survey.pdf
+    .. [4] ImageJ AutoThresholder code,
+           http://fiji.sc/wiki/index.php/Auto_Threshold
+
+    Examples
+    --------
+    >>> from skimage.data import coins
+    >>> image = coins()
+    >>> thresh = threshold_isodata(image)
+    >>> binary = image > thresh
+    """
+    hist, bin_centers = histogram(image, nbins)
+    # On blank images (e.g. filled with 0) with int dtype, `histogram()`
+    # returns `bin_centers` containing only one value. Speed up with it.
+    if bin_centers.size == 1:
+        return bin_centers[0]
+    # It is not necessary to calculate the probability mass function here,
+    # because the l and h fractions already include the normalization.
+    pmf = hist.astype(np.float32)  # / hist.sum()
+    cpmfl = np.cumsum(pmf, dtype=np.float32)
+    cpmfh = np.cumsum(pmf[::-1], dtype=np.float32)[::-1]
+
+    binnums = np.arange(pmf.size, dtype=np.uint8)
+    # l and h contain average value of pixels in sum of bins, calculated
+    # from lower to higher and from higher to lower respectively.
+    l = np.ma.divide(np.cumsum(pmf * binnums, dtype=np.float32), cpmfl)
+    h = np.ma.divide(
+        np.cumsum((pmf[::-1] * binnums[::-1]), dtype=np.float32)[::-1],
+        cpmfh)
+
+    allmean = (l + h) / 2.0
+    threshold = bin_centers[np.nonzero(allmean.round() == binnums)[0][0]]
+    # This implementation returns threshold where
+    # `background <= threshold < foreground`.
     return threshold