Merge pull request #1376 from jmetz/threshold_li

Li thresholding port from imagej

skimage/filters/LICENSE.txt

@@ -0,0 +1,28 @@
+threshold_li based on imagej code:
+
+Copyright (c) 2009 - 2015, Board of Regents of the University of
+Wisconsin-Madison, Broad Institute of MIT and Harvard, and Max Planck
+Institute of Molecular Cell Biology and Genetics.
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification,
+are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this
+   list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the documentation
+   and/or other materials provided with the distribution.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGE.

skimage/filters/__init__.py

@@ -9,7 +9,7 @@ from .edges import (sobel, hsobel, vsobel, sobel_h, sobel_v,
 from ._rank_order import rank_order
 from ._gabor import gabor_kernel, gabor_filter
 from .thresholding import (threshold_adaptive, threshold_otsu, threshold_yen,
-                           threshold_isodata)
+                           threshold_isodata, threshold_li)
 from . import rank
 from .rank import median
 
@@ -66,4 +66,5 @@ __all__ = ['inverse',
            'threshold_otsu',
            'threshold_yen',
            'threshold_isodata',
+           'threshold_li', 
            'rank']

skimage/filters/tests/test_thresholding.py

@@ -5,6 +5,7 @@ import skimage
 from skimage import data
 from skimage.filters.thresholding import (threshold_adaptive,
                                           threshold_otsu,
+                                          threshold_li,
                                           threshold_yen,
                                           threshold_isodata)
 
@@ -28,6 +29,17 @@ class TestSimpleImage():
         image = np.float64(self.image)
         assert 2 <= threshold_otsu(image) < 3
 
+    def test_li(self):
+        assert int(threshold_li(self.image)) == 2
+
+    def test_li_negative_int(self):
+        image = self.image - 2
+        assert int(threshold_li(image)) == 0
+
+    def test_li_float_image(self):
+        image = np.float64(self.image)
+        assert 2 <= threshold_li(image) < 3
+
     def test_yen(self):
         assert threshold_yen(self.image) == 2
 
@@ -152,6 +164,25 @@ def test_otsu_astro_image():
     img = skimage.img_as_ubyte(data.astronaut())
     assert 109 < threshold_otsu(img) < 111
 
+def test_li_camera_image():
+    camera = skimage.img_as_ubyte(data.camera())
+    assert 63 < threshold_li(camera) < 65
+
+
+def test_li_coins_image():
+    coins = skimage.img_as_ubyte(data.coins())
+    assert 95 < threshold_li(coins) < 97
+
+
+def test_li_coins_image_as_float():
+    coins = skimage.img_as_float(data.coins())
+    assert 0.37 < threshold_li(coins) < 0.38
+
+
+def test_li_astro_image():
+    img = skimage.img_as_ubyte(data.astronaut())
+    assert 66 < threshold_li(img) < 68
+
 def test_yen_camera_image():
     camera = skimage.img_as_ubyte(data.camera())
     assert 197 < threshold_yen(camera) < 199

skimage/filters/thresholding.py

@@ -1,7 +1,8 @@
 __all__ = ['threshold_adaptive',
            'threshold_otsu',
            'threshold_yen',
-           'threshold_isodata']
+           'threshold_isodata',
+           'threshold_li', ]
 
 import numpy as np
 import scipy.ndimage
@@ -302,3 +303,71 @@ def threshold_isodata(image, nbins=256, return_all=False):
         return thresholds
     else:
         return thresholds[0]
+
+
+def threshold_li(image):
+    """Return threshold value based on adaptation of Li's Minimum Cross Entropy method.
+
+    Parameters
+    ----------
+    image : array
+        Input image.
+
+    Returns
+    -------
+    threshold : float
+        Upper threshold value. All pixels intensities that less or equal of
+        this value assumed as foreground.
+
+    References
+    ----------
+    .. [1] Li C.H. and Lee C.K. (1993) "Minimum Cross Entropy Thresholding"
+           Pattern Recognition, 26(4): 617-625
+    .. [2] Li C.H. and Tam P.K.S. (1998) "An Iterative Algorithm for Minimum
+           Cross Entropy Thresholding" Pattern Recognition Letters, 18(8): 771-776
+    .. [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://citeseer.ist.psu.edu/sezgin04survey.html
+    .. [4] ImageJ AutoThresholder code, http://fiji.sc/wiki/index.php/Auto_Threshold
+
+    Adapted for skimage by J. Metz from ImageJ plugin by G.Landini
+
+    Examples
+    --------
+    >>> from skimage.data import camera
+    >>> image = camera()
+    >>> thresh = threshold_li(image)
+    >>> binary = image <= thresh
+    """
+    # Requires positive image (because of log(mean))
+    offset = image.min()
+    # Can not use fixed tolerance for float image
+    imrange = image.max() - offset
+    image -= offset
+
+    tolerance = 0.5 * imrange / 256.0
+    # Calculate the mean gray-level
+    mean = image.mean()
+
+    # Initial estimate
+    new_thresh = mean
+    old_thresh = new_thresh + 2 * tolerance
+
+    # Stop the iterations when the difference between the
+    # new and old threshold values is less than the tolerance
+    while abs(new_thresh - old_thresh) > tolerance:
+        old_thresh = new_thresh
+        threshold = old_thresh + tolerance   # range
+        # Calculate the means of background and object pixels
+        mean_back = image[image <= threshold].mean()
+        mean_obj = image[image > threshold].mean()
+
+        temp = (mean_back - mean_obj) / (np.log(mean_back) - np.log(mean_obj))
+
+        if temp < 0:
+            new_thresh = temp - tolerance
+        else:
+            new_thresh = temp + tolerance
+
+    return threshold + offset