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/filters/plot_threshold_minimum.py b/doc/examples/filters/plot_threshold_minimum.py new file mode 100644 index 00000000..950d34d2 --- /dev/null +++ b/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() 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..4891228e 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(): @@ -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() diff --git a/skimage/filters/thresholding.py b/skimage/filters/thresholding.py index 2f9023d6..241a1a85 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 @@ -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]