Add Yen threshold method

skimage/filter/thresholding.py

@@ -133,3 +133,51 @@ 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
+        Number of bins used to calculate histogram. This value is ignored for
+        integer arrays.
+    
+    Returns
+    -------
+    threshold : float
+        Threshold value.
+
+    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(img, nbins)
+    hist = hist.astype(float)
+    norm_histo = hist / hist.sum() # Probability mass function
+    P1 = np.cumsum(norm_histo) # Cumulative normalized histogram
+    P1_sq = np.cumsum(norm_histo ** 2)
+    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.
+    # In ImageJ Yen implementation, all invalid values replaced by zero.
+    crit = -1*np.log(P1_sq[:-1]*P2_sq[1:]) + 2.0*np.log(P1[:-1]*(1.0-P1[:-1]))
+    max_crit = np.argmax(crit)
+    threshold = bin_centers[:-1][max_crit]
+    return threshold