Add threshold minimum algorithm

doc/examples/plot_threshold_minimum.py

@@ -0,0 +1,37 @@
+"""
+==================================
+Minimum Algorithm For Thresholding
+==================================
+
+The minimum algorithm takes a histogram of the image and smooths it
+repeatedly unitl there are only two peaks in the histogram.  Then it
+finds the minimum value between the two peaks.  With the smoothing
+there can be multiple pixel values with the minimum histogram count,
+so you can pick '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('Image')
+
+ax1.imshow(binarized)
+ax1.set_title('Thresholded')
+
+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():
@@ -273,5 +274,45 @@ 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.zeros((25*25), dtype=np.uint8)
+    for i in range(25*25) :
+        img[i] = i % 256
+    img = np.reshape(img, (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)
+    for i in range(16*16) :
+        img[i] = i % 256
+    img = np.reshape(img, (16,16))
+    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
 
@@ -389,3 +391,91 @@ def threshold_li(image):
             new_thresh = temp + tolerance
 
     return threshold + immin
+
+
+def threshold_minimum(image, nbins=256, bias='min'):
+    """Return threshold value based on minimum method.
+
+    A histogram is computed and smoothed until there are only two maximums.
+    Then the minimum between these is found.
+
+    Parameters
+    ----------
+    image : array
+        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.
+
+    Returns
+    -------
+    threshold : float
+        Computed threshold value.
+        May return 0 if the algorithm fails.
+
+    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):
+        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
+
+    hist, bin_centers = histogram(image.ravel(), nbins)
+    # On blank images (e.g. filled with 0) with int dtype, `histogram()`
+    # returns `bin_centers` containing only one value. Speed up with it.
+    if bin_centers.size == 1:
+        return bin_centers[0]
+
+    threshold = 0
+    smooth_hist = np.copy(hist)
+    smooth_hist = ndif.uniform_filter1d(smooth_hist, 3)
+    maximums = find_local_maxima(smooth_hist)
+    while len(maximums) > 2:
+        smooth_hist = ndif.uniform_filter1d(smooth_hist, 3)
+        maximums = find_local_maxima(smooth_hist)
+
+    # If we don't have 2 maxima the algorithm failed
+    if len(maximums) != 2:
+        raise RuntimeError('Failure: Unable to find two maxima in histogram')
+
+    # 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]):
+        if smooth_hist[i] < minimum:
+            minimum = smooth_hist[i]
+            threshold = i
+
+    if bias == 'min':
+        return bin_centers[threshold]
+    else:
+        i = threshold
+        while smooth_hist[i] == smooth_hist[threshold]:
+            i += 1
+        i -= 1
+        if bias == 'max':
+            return bin_centers[i]
+        else:
+            return bin_centers[(threshold + i) // 2]