From 80e85440673f79b632e198f7d71cb2255506c554 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20Boulogne?= Date: Mon, 23 May 2016 14:25:28 +0200 Subject: [PATCH] Review threshold minimum --- CONTRIBUTORS.txt | 3 + .../{ => filters}/plot_threshold_minimum.py | 13 ++-- skimage/filters/tests/test_thresholding.py | 18 ++--- skimage/filters/thresholding.py | 73 +++++++++++-------- 4 files changed, 60 insertions(+), 47 deletions(-) rename doc/examples/{ => filters}/plot_threshold_minimum.py (64%) diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 7318049a..5900487c 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -227,3 +227,6 @@ - Alex Izvorski Color spaces for YUV and related spaces + +- Jeff Hemmelgarn + Minimum threshold diff --git a/doc/examples/plot_threshold_minimum.py b/doc/examples/filters/plot_threshold_minimum.py similarity index 64% rename from doc/examples/plot_threshold_minimum.py rename to doc/examples/filters/plot_threshold_minimum.py index 67b23342..dbbb7324 100644 --- a/doc/examples/plot_threshold_minimum.py +++ b/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') diff --git a/skimage/filters/tests/test_thresholding.py b/skimage/filters/tests/test_thresholding.py index 3e772b0b..dae347fa 100644 --- a/skimage/filters/tests/test_thresholding.py +++ b/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) diff --git a/skimage/filters/thresholding.py b/skimage/filters/thresholding.py index 12c3206b..a5f21de9 100644 --- a/skimage/filters/thresholding.py +++ b/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]