Merge pull request #686 from radioxoma/yen-threshold

Add Yen threshold method

CONTRIBUTORS.txt

@@ -148,3 +148,6 @@
 
 - Matt Terry
   Color difference functions
+
+- Eugene Dvoretsky
+  Yen threshold implementation.

skimage/filter/tests/test_thresholding.py

@@ -3,7 +3,9 @@ from numpy.testing import assert_array_equal
 
 import skimage
 from skimage import data
-from skimage.filter.thresholding import threshold_otsu, threshold_adaptive
+from skimage.filter.thresholding import (threshold_adaptive,
+                                         threshold_otsu,
+                                         threshold_yen)
 
 
 class TestSimpleImage():
@@ -25,6 +27,26 @@ class TestSimpleImage():
         image = np.float64(self.image)
         assert 2 <= threshold_otsu(image) < 3
 
+    def test_yen(self):
+        assert threshold_yen(self.image) == 2
+
+    def test_yen_negative_int(self):
+        image = self.image - 2
+        assert threshold_yen(image) == 0
+
+    def test_yen_float_image(self):
+        image = np.float64(self.image)
+        assert 2 <= threshold_yen(image) < 3
+
+    def test_yen_arange(self):
+        image = np.arange(256)
+        assert threshold_yen(image) == 127
+
+    def test_yen_binary(self):
+        image = np.zeros([2,256], dtype='uint8')
+        image[0] = 255
+        assert threshold_yen(image) < 1
+
     def test_threshold_adaptive_generic(self):
         def func(arr):
             return arr.sum() / arr.shape[0]
@@ -92,5 +114,15 @@ def test_otsu_lena_image():
     assert 140 < threshold_otsu(lena) < 142
 
 
+def test_yen_coins_image():
+    coins = skimage.img_as_ubyte(data.coins())
+    assert 109 < threshold_yen(coins) < 111
+
+
+def test_yen_coins_image_as_float():
+    coins = skimage.img_as_float(data.coins())
+    assert 0.43 < threshold_yen(coins) < 0.44
+
+
 if __name__ == '__main__':
     np.testing.run_module_suite()

skimage/filter/thresholding.py

@@ -1,4 +1,4 @@
-__all__ = ['threshold_otsu', 'threshold_adaptive']
+__all__ = ['threshold_adaptive', 'threshold_otsu', 'threshold_yen'] 
 
 import numpy as np
 import scipy.ndimage
@@ -95,14 +95,15 @@ def threshold_otsu(image, nbins=256):
     ----------
     image : array
         Input image.
-    nbins : int
+    nbins : int, optional
         Number of bins used to calculate histogram. This value is ignored for
         integer arrays.
 
     Returns
     -------
     threshold : float
-        Threshold value.
+        Upper threshold value. All pixels intensities that less or equal of
+        this value assumed as foreground.
 
     References
     ----------
@@ -113,7 +114,7 @@ def threshold_otsu(image, nbins=256):
     >>> from skimage.data import camera
     >>> image = camera()
     >>> thresh = threshold_otsu(image)
-    >>> binary = image > thresh
+    >>> binary = image <= thresh
     """
     hist, bin_centers = histogram(image, nbins)
     hist = hist.astype(float)
@@ -133,3 +134,53 @@ def threshold_otsu(image, nbins=256):
     idx = np.argmax(variance12)
     threshold = bin_centers[:-1][idx]
     return threshold
+
+
+def threshold_yen(image, nbins=256):
+    """Return threshold value based on Yen's method.
+
+    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
+        Upper threshold value. All pixels intensities that less or equal of
+        this value assumed as foreground.
+
+    References
+    ----------
+    .. [1] Yen J.C., Chang F.J., and Chang S. (1995) "A New Criterion 
+           for Automatic Multilevel Thresholding" IEEE Trans. on Image 
+           Processing, 4(3): 370-378
+    .. [2] 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
+    .. [3] ImageJ AutoThresholder code, http://fiji.sc/wiki/index.php/Auto_Threshold
+
+    Examples
+    --------
+    >>> from skimage.data import camera
+    >>> image = camera()
+    >>> thresh = threshold_yen(image)
+    >>> 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)
+    # Get cumsum calculated from end of squared array:
+    P2_sq = np.cumsum(norm_histo[::-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) * \
+                  (P1[:-1] * (1.0 - P1[:-1])) ** 2)
+    max_crit = np.argmax(crit)
+    threshold = bin_centers[:-1][max_crit]
+    return threshold