From 31db1c607e00fe997c0f9e2c35b716892097fecb Mon Sep 17 00:00:00 2001 From: Jeff Hemmelgarn Date: Sat, 18 Jul 2015 09:06:50 -0400 Subject: [PATCH 1/5] Add threshold minimum algorithm --- doc/examples/plot_threshold_minimum.py | 37 +++++++++ skimage/filters/__init__.py | 5 +- skimage/filters/tests/test_thresholding.py | 43 +++++++++- skimage/filters/thresholding.py | 92 +++++++++++++++++++++- 4 files changed, 173 insertions(+), 4 deletions(-) create mode 100644 doc/examples/plot_threshold_minimum.py diff --git a/doc/examples/plot_threshold_minimum.py b/doc/examples/plot_threshold_minimum.py new file mode 100644 index 00000000..67b23342 --- /dev/null +++ b/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() diff --git a/skimage/filters/__init__.py b/skimage/filters/__init__.py index a6245886..ca0c12c3 100644 --- a/skimage/filters/__init__.py +++ b/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'] diff --git a/skimage/filters/tests/test_thresholding.py b/skimage/filters/tests/test_thresholding.py index c5d583b2..3e772b0b 100644 --- a/skimage/filters/tests/test_thresholding.py +++ b/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() diff --git a/skimage/filters/thresholding.py b/skimage/filters/thresholding.py index 2f9023d6..12c3206b 100644 --- a/skimage/filters/thresholding.py +++ b/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] 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 2/5] 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] From aea36e89d19f551a0c6a0f12e6ce1e1055ba14ef Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20Boulogne?= Date: Sat, 28 May 2016 08:36:14 +0200 Subject: [PATCH 3/5] use break instead of exception --- skimage/filters/thresholding.py | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/skimage/filters/thresholding.py b/skimage/filters/thresholding.py index a5f21de9..ffab8748 100644 --- a/skimage/filters/thresholding.py +++ b/skimage/filters/thresholding.py @@ -458,16 +458,14 @@ def threshold_minimum(image, nbins=256, bias='min'): max_iter = 10000 smooth_hist = np.copy(hist) - 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: + 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 histogram' 'smoothing') From 8c9d8926c881e55a2bcb0fcbc9cb52b066aa8e72 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20Boulogne?= Date: Thu, 2 Jun 2016 10:31:40 +0200 Subject: [PATCH 4/5] Minor enhancements --- doc/examples/filters/plot_threshold_minimum.py | 2 +- skimage/filters/thresholding.py | 16 +++++++++------- 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/doc/examples/filters/plot_threshold_minimum.py b/doc/examples/filters/plot_threshold_minimum.py index dbbb7324..950d34d2 100644 --- a/doc/examples/filters/plot_threshold_minimum.py +++ b/doc/examples/filters/plot_threshold_minimum.py @@ -7,7 +7,7 @@ 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 'min', 'mid', or 'max' of these values. +count, so you can pick the 'min', 'mid', or 'max' of these values. """ import matplotlib.pyplot as plt diff --git a/skimage/filters/thresholding.py b/skimage/filters/thresholding.py index ffab8748..241a1a85 100644 --- a/skimage/filters/thresholding.py +++ b/skimage/filters/thresholding.py @@ -112,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 @@ -395,11 +395,11 @@ def threshold_li(image): return threshold + immin -def threshold_minimum(image, nbins=256, bias='min'): +def threshold_minimum(image, nbins=256, bias='min', max_iter=10000): """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. + 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 ---------- @@ -411,11 +411,14 @@ def threshold_minimum(image, nbins=256, bias='min'): 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 - Computed threshold value. + Upper threshold value. All pixels with an intensity higher than + this value are assumed to be foreground. Raises ------ @@ -456,7 +459,6 @@ def threshold_minimum(image, nbins=256, bias='min'): hist, bin_centers = histogram(image.ravel(), nbins) - max_iter = 10000 smooth_hist = np.copy(hist) for counter in range(max_iter): smooth_hist = ndif.uniform_filter1d(smooth_hist, 3) @@ -466,7 +468,7 @@ def threshold_minimum(image, nbins=256, bias='min'): 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 histogram' + raise RuntimeError('Maximum iteration reached for histogram' 'smoothing') # Find lowest point between the maxima, biased to the low end (min) From ea9180f7548766d8552db6caa8ea4bb64318984e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20Boulogne?= Date: Thu, 2 Jun 2016 13:38:32 +0200 Subject: [PATCH 5/5] more concise writting --- skimage/filters/tests/test_thresholding.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/skimage/filters/tests/test_thresholding.py b/skimage/filters/tests/test_thresholding.py index dae347fa..4891228e 100644 --- a/skimage/filters/tests/test_thresholding.py +++ b/skimage/filters/tests/test_thresholding.py @@ -292,10 +292,7 @@ def test_threshold_minimum(): 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 = np.arange(25*25, dtype=np.uint8).reshape((25, 25)) img[0:9, :] = 50 img[14:25, :] = 250