Review threshold minimum

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

@@ -4,10 +4,10 @@ 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.
+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 'min', 'mid', or 'max' of these values.
 
 """
 import matplotlib.pyplot as plt
@@ -18,7 +18,6 @@ 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))
@@ -26,10 +25,10 @@ ax0, ax1 = axes
 plt.gray()
 
 ax0.imshow(image)
-ax0.set_title('Image')
+ax0.set_title('Original image')
 
 ax1.imshow(binarized)
-ax1.set_title('Thresholded')
+ax1.set_title('Result')
 
 for ax in axes:
     ax.axis('off')

skimage/filters/tests/test_thresholding.py

@@ -60,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
 
@@ -119,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):
@@ -170,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
@@ -189,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
@@ -290,11 +293,11 @@ def test_threshold_minimum():
 
 def test_threshold_minimum_synthetic():
     img = np.zeros((25*25), dtype=np.uint8)
-    for i in range(25*25) :
+    for i in range(25*25):
         img[i] = i % 256
-    img = np.reshape(img, (25,25))
-    img[0:9,:] = 50
-    img[14:25,:] = 250
+    img = np.reshape(img, (25, 25))
+    img[0:9, :] = 50
+    img[14:25, :] = 250
 
     threshold = threshold_minimum(img, bias='min')
     assert threshold == 93
@@ -308,9 +311,6 @@ def test_threshold_minimum_synthetic():
 
 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)
 
 

skimage/filters/thresholding.py

@@ -101,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
@@ -211,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()]
@@ -345,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
@@ -401,12 +403,12 @@ def threshold_minimum(image, nbins=256, bias='min'):
 
     Parameters
     ----------
-    image : array
+    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
+    bias : {'min', 'mid', 'max'}, optional
         'min', 'mid', 'max' return lowest, middle, or highest pixel value
         with minimum histogram value.
 
@@ -414,7 +416,12 @@ def threshold_minimum(image, nbins=256, bias='min'):
     -------
     threshold : float
         Computed threshold value.
-        May return 0 if the algorithm fails.
+
+    Raises
+    ------
+    RuntimeError
+        If unable to find two local maxima in the histogram or if the
+        smoothing takes more than 1e4 iterations.
 
     References
     ----------
@@ -430,8 +437,10 @@ def threshold_minimum(image, nbins=256, bias='min'):
     """
 
     def find_local_maxima(hist):
+        # We can't use scipy.signal.argrelmax
+        # as it fails on plateaus
         maximums = list()
-        direction = +1
+        direction = 1
         for i in range(hist.shape[0] - 1):
             if direction > 0:
                 if hist[i + 1] < hist[i]:
@@ -439,31 +448,33 @@ def threshold_minimum(image, nbins=256, bias='min'):
                     maximums.append(i)
             else:
                 if hist[i + 1] > hist[i]:
-                    direction = +1
+                    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)
-    # 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
+    max_iter = 10000
     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')
+    try:
+        for counter in range(max_iter):
+            smooth_hist = ndif.uniform_filter1d(smooth_hist, 3)
+            maximums = find_local_maxima(smooth_hist)
+            if len(maximums) < 3:
+                raise StopIteration
+    except StopIteration:
+        if len(maximums) != 2:
+            raise RuntimeError('Unable to find two maxima in histogram')
+    else:
+        raise RuntimeError('Maximum iteration reached for histogram 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]):
+    for i in range(maximums[0], maximums[1]+1):
         if smooth_hist[i] < minimum:
             minimum = smooth_hist[i]
             threshold = i
@@ -471,11 +482,11 @@ def threshold_minimum(image, nbins=256, bias='min'):
     if bias == 'min':
         return bin_centers[threshold]
     else:
-        i = threshold
-        while smooth_hist[i] == smooth_hist[threshold]:
-            i += 1
-        i -= 1
+        upper_bound = threshold
+        while smooth_hist[upper_bound] == smooth_hist[threshold]:
+            upper_bound += 1
+        upper_bound -= 1
         if bias == 'max':
-            return bin_centers[i]
-        else:
-            return bin_centers[(threshold + i) // 2]
+            return bin_centers[upper_bound]
+        elif bias == 'mid':
+            return bin_centers[(threshold + upper_bound) // 2]