Add ISODATA threshold with tests

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_adaptive, threshold_otsu, threshold_yen
+from .thresholding import (threshold_adaptive, threshold_otsu, threshold_yen,
+                           threshold_isodata)
 from . import rank
 
 
@@ -40,4 +41,5 @@ __all__ = ['inverse',
            '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():
@@ -56,6 +57,16 @@ class TestSimpleImage():
         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]
@@ -123,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
@@ -133,9 +149,19 @@ def test_yen_coins_image_as_float():
     assert 0.43 < threshold_yen(coins) < 0.44
 
 
-def test_yen_camera_image():
+def test_isodata_camera_image():
     camera = skimage.img_as_ubyte(data.camera())
-    assert 197 < threshold_yen(camera) < 199
+    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
 
 
 if __name__ == '__main__':

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
@@ -185,3 +188,63 @@ def threshold_yen(image, nbins=256):
     crit = np.log(((P1_sq[:-1] * P2_sq[1:]) ** -1) *
                   (P1[:-1] * (1.0 - P1[:-1])) ** 2)
     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 float or int of any range.
+    nbins : int, optional
+        Number of bins used to calculate histogram. This value is ignored for
+        integer arrays.
+
+    Returns
+    -------
+    threshold : float64 or int64
+        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)
+    if bin_centers.size == 1:
+        return bin_centers[0]
+    # It is not necessary to calculate probability mass function in this case
+    # since in the l and h fractions it's reduced.
+    pmf = hist.astype(float)# / hist.sum()
+    cpmfl = np.cumsum(pmf, dtype=float)  # Cumulative probability mass function
+    cpmfh = np.cumsum(pmf[::-1], dtype=float)[::-1]
+
+    binnums = np.arange(pmf.size, dtype=np.uint8)
+    l = np.ma.divide(np.cumsum(pmf * binnums, dtype=float), cpmfl)
+    h = np.ma.divide(np.cumsum((pmf[::-1] * binnums[::-1]), dtype=float)[::-1], 
+                     cpmfh)
+
+    allmean = (l + h) / 2.0
+    threshold = bin_centers[np.nonzero(allmean.round() == binnums)[0][0]]
+    # ImageJ shows *inclusive* threshold. This implementation returns
+    # threshold, where `background <= threshold_value < foreground`.
+    return threshold