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]