Merge pull request #2104 from sciunto/threshold-minimum

Add threshold minimum algorithm

CONTRIBUTORS.txt

@@ -227,3 +227,6 @@
 
 - Alex Izvorski
   Color spaces for YUV and related spaces
+
+- Jeff Hemmelgarn
+  Minimum threshold

doc/examples/filters/plot_threshold_minimum.py

@@ -0,0 +1,36 @@
+"""
+==================================
+Minimum Algorithm For Thresholding
+==================================
+
+The minimum algorithm takes a histogram of the image and smooths it
+repeatedly until there are only two peaks in the histogram.  Then it
+finds the minimum value between the two peaks.  After smoothing the
+histogram, there can be multiple pixel values with the minimum histogram
+count, so you can pick the 'min', 'mid', or 'max' of these values.
+
+"""
+import matplotlib.pyplot as plt
+
+from skimage import data
+from skimage.filters.thresholding import threshold_minimum
+
+image = data.camera()
+
+threshold = threshold_minimum(image, bias='min')
+binarized = image > threshold
+
+fig, axes = plt.subplots(nrows=2, figsize=(7, 8))
+ax0, ax1 = axes
+plt.gray()
+
+ax0.imshow(image)
+ax0.set_title('Original image')
+
+ax1.imshow(binarized)
+ax1.set_title('Result')
+
+for ax in axes:
+    ax.axis('off')
+
+plt.show()

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
 from .thresholding import (threshold_adaptive, threshold_otsu, threshold_yen,
-                           threshold_isodata, threshold_li)
+                           threshold_isodata, threshold_li, threshold_minimum)
 from . import rank
 from .rank import median
 
@@ -57,5 +57,6 @@ __all__ = ['inverse',
            'threshold_otsu',
            'threshold_yen',
            'threshold_isodata',
-           'threshold_li',
+           'threshold_li', 
+           'threshold_minimum',
            'rank']

skimage/filters/tests/test_thresholding.py

@@ -10,7 +10,8 @@ from skimage.filters.thresholding import (threshold_adaptive,
                                           threshold_otsu,
                                           threshold_li,
                                           threshold_yen,
-                                          threshold_isodata)
+                                          threshold_isodata,
+                                          threshold_minimum)
 
 
 class TestSimpleImage():
@@ -59,7 +60,7 @@ class TestSimpleImage():
         assert threshold_yen(image) == 127
 
     def test_yen_binary(self):
-        image = np.zeros([2,256], dtype=np.uint8)
+        image = np.zeros([2, 256], dtype=np.uint8)
         image[0] = 255
         assert threshold_yen(image) < 1
 
@@ -118,7 +119,8 @@ class TestSimpleImage():
         out = threshold_adaptive(self.image, 3, method='gaussian')
         assert_equal(ref, out)
 
-        out = threshold_adaptive(self.image, 3, method='gaussian', param=1.0 / 3.0)
+        out = threshold_adaptive(self.image, 3, method='gaussian',
+                                 param=1.0 / 3.0)
         assert_equal(ref, out)
 
     def test_threshold_adaptive_mean(self):
@@ -169,6 +171,7 @@ def test_otsu_one_color_image():
     img = np.ones((10, 10), dtype=np.uint8)
     assert_raises(ValueError, threshold_otsu, img)
 
+
 def test_li_camera_image():
     camera = skimage.img_as_ubyte(data.camera())
     assert 63 < threshold_li(camera) < 65
@@ -188,6 +191,7 @@ 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
@@ -273,5 +277,39 @@ def test_isodata_moon_image_negative_float():
                          23.01757812, 24.01367188, 38.95507812, 39.95117188])
 
 
+def test_threshold_minimum():
+    camera = skimage.img_as_ubyte(data.camera())
+
+    threshold = threshold_minimum(camera)
+    assert threshold == 76
+
+    threshold = threshold_minimum(camera, bias='max')
+    assert threshold == 77
+
+    astronaut = skimage.img_as_ubyte(data.astronaut())
+    threshold = threshold_minimum(astronaut)
+    assert threshold == 117
+
+
+def test_threshold_minimum_synthetic():
+    img = np.arange(25*25, dtype=np.uint8).reshape((25, 25))
+    img[0:9, :] = 50
+    img[14:25, :] = 250
+
+    threshold = threshold_minimum(img, bias='min')
+    assert threshold == 93
+
+    threshold = threshold_minimum(img, bias='mid')
+    assert threshold == 159
+
+    threshold = threshold_minimum(img, bias='max')
+    assert threshold == 225
+
+
+def test_threshold_minimum_failure():
+    img = np.zeros((16*16), dtype=np.uint8)
+    assert_raises(RuntimeError, threshold_minimum, img)
+
+
 if __name__ == '__main__':
     np.testing.run_module_suite()

skimage/filters/thresholding.py

@@ -2,10 +2,12 @@ __all__ = ['threshold_adaptive',
            'threshold_otsu',
            'threshold_yen',
            'threshold_isodata',
-           'threshold_li', ]
+           'threshold_li',
+           'threshold_minimum', ]
 
 import numpy as np
 from scipy import ndimage as ndi
+from scipy.ndimage import filters as ndif
 from ..exposure import histogram
 from .._shared.utils import assert_nD, warn
 
@@ -99,7 +101,7 @@ def threshold_otsu(image, nbins=256):
 
     Parameters
     ----------
-    image : array
+    image : (M, N) ndarray
         Grayscale input image.
     nbins : int, optional
         Number of bins used to calculate histogram. This value is ignored for
@@ -110,7 +112,7 @@ def threshold_otsu(image, nbins=256):
     threshold : float
         Upper threshold value. All pixels intensities that less or equal of
         this value assumed as foreground.
-        
+
     Raises
     ------
     ValueError
@@ -209,8 +211,9 @@ def threshold_yen(image, nbins=256):
     P1_sq = np.cumsum(pmf ** 2)
     # Get cumsum calculated from end of squared array:
     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.
+    # P2_sq indexes is shifted +1.
+    # I assume, with P1[:-1] it helps to 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)
     return bin_centers[crit.argmax()]
@@ -343,7 +346,8 @@ def threshold_li(image):
     .. [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
+           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
@@ -389,3 +393,100 @@ def threshold_li(image):
             new_thresh = temp + tolerance
 
     return threshold + immin
+
+
+def threshold_minimum(image, nbins=256, bias='min', max_iter=10000):
+    """Return threshold value based on minimum method.
+
+    The histogram of the input `image` is computed and smoothed until there are
+    only two maxima. Then the minimum in between is the threshold value.
+
+    Parameters
+    ----------
+    image : (M, N) ndarray
+        Input image.
+    nbins : int, optional
+        Number of bins used to calculate histogram. This value is ignored for
+        integer arrays.
+    bias : {'min', 'mid', 'max'}, optional
+        'min', 'mid', 'max' return lowest, middle, or highest pixel value
+        with minimum histogram value.
+    max_iter: int, optional
+        Maximum number of iterations to smooth the histogram.
+
+    Returns
+    -------
+    threshold : float
+        Upper threshold value. All pixels with an intensity higher than
+        this value are assumed to be foreground.
+
+    Raises
+    ------
+    RuntimeError
+        If unable to find two local maxima in the histogram or if the
+        smoothing takes more than 1e4 iterations.
+
+    References
+    ----------
+    .. [1] Prewitt, JMS & Mendelsohn, ML (1966), "The analysis of cell images",
+           Annals of the New York Academy of Sciences 128: 1035-1053
+
+    Examples
+    --------
+    >>> from skimage.data import camera
+    >>> image = camera()
+    >>> thresh = threshold_minimum(image)
+    >>> binary = image > thresh
+    """
+
+    def find_local_maxima(hist):
+        # We can't use scipy.signal.argrelmax
+        # as it fails on plateaus
+        maximums = list()
+        direction = 1
+        for i in range(hist.shape[0] - 1):
+            if direction > 0:
+                if hist[i + 1] < hist[i]:
+                    direction = -1
+                    maximums.append(i)
+            else:
+                if hist[i + 1] > hist[i]:
+                    direction = 1
+        return maximums
+
+    if bias not in ('min', 'mid', 'max'):
+        raise ValueError("Unknown bias: {0}".format(bias))
+
+    hist, bin_centers = histogram(image.ravel(), nbins)
+
+    smooth_hist = np.copy(hist)
+    for counter in range(max_iter):
+        smooth_hist = ndif.uniform_filter1d(smooth_hist, 3)
+        maximums = find_local_maxima(smooth_hist)
+        if len(maximums) < 3:
+            break
+    if len(maximums) != 2:
+        raise RuntimeError('Unable to find two maxima in histogram')
+    elif counter == max_iter - 1:
+        raise RuntimeError('Maximum iteration reached for histogram'
+                           'smoothing')
+
+    # Find lowest point between the maxima, biased to the low end (min)
+    minimum = smooth_hist[maximums[0]]
+    threshold = maximums[0]
+    for i in range(maximums[0], maximums[1]+1):
+        if smooth_hist[i] < minimum:
+            minimum = smooth_hist[i]
+            threshold = i
+
+    if bias == 'min':
+        return bin_centers[threshold]
+    else:
+        upper_bound = threshold
+        while smooth_hist[upper_bound] == smooth_hist[threshold]:
+            upper_bound += 1
+        upper_bound -= 1
+        if bias == 'max':
+            return bin_centers[upper_bound]
+        elif bias == 'mid':
+            return bin_centers[(threshold + upper_bound) // 2]