diff --git a/skimage/filters/__init__.py b/skimage/filters/__init__.py index ca0c12c3..3759136a 100644 --- a/skimage/filters/__init__.py +++ b/skimage/filters/__init__.py @@ -57,6 +57,8 @@ __all__ = ['inverse', 'threshold_otsu', 'threshold_yen', 'threshold_isodata', - 'threshold_li', + 'threshold_li', 'threshold_minimum', + 'threshold_mean', + 'threshold_triangle', 'rank'] diff --git a/skimage/filters/tests/test_thresholding.py b/skimage/filters/tests/test_thresholding.py index 4891228e..dd08e36b 100644 --- a/skimage/filters/tests/test_thresholding.py +++ b/skimage/filters/tests/test_thresholding.py @@ -11,6 +11,8 @@ from skimage.filters.thresholding import (threshold_adaptive, threshold_li, threshold_yen, threshold_isodata, + threshold_mean, + threshold_triangle, threshold_minimum) @@ -120,7 +122,7 @@ class TestSimpleImage(): assert_equal(ref, out) out = threshold_adaptive(self.image, 3, method='gaussian', - param=1.0 / 3.0) + param=1./3.) assert_equal(ref, out) def test_threshold_adaptive_mean(self): @@ -311,5 +313,37 @@ def test_threshold_minimum_failure(): assert_raises(RuntimeError, threshold_minimum, img) +def test_mean(): + img = np.zeros((2, 6)) + img[:, 2:4] = 1 + img[:, 4:] = 2 + assert(threshold_mean(img) == 1.) + + +def test_triangle_images(): + assert(threshold_triangle(np.invert(data.text())) == 151) + assert(threshold_triangle(data.text()) == 104) + assert(threshold_triangle(data.coins()) == 80) + assert(threshold_triangle(np.invert(data.coins())) == 175) + + +def test_triangle_flip(): + # Depending on the skewness, the algorithm flips the histogram. + # We check that the flip doesn't affect too much the result. + img = data.camera() + inv_img = np.invert(img) + t = threshold_triangle(inv_img) + t_inv_img = inv_img > t + t_inv_inv_img = np.invert(t_inv_img) + + t = threshold_triangle(img) + t_img = img > t + + # Check that most of the pixels are identical + # See numpy #7685 for a future np.testing API + unequal_pos = np.where(t_img.ravel() != t_inv_inv_img.ravel()) + assert(len(unequal_pos[0]) / t_img.size < 1e-2) + + if __name__ == '__main__': np.testing.run_module_suite() diff --git a/skimage/filters/thresholding.py b/skimage/filters/thresholding.py index 241a1a85..2896c982 100644 --- a/skimage/filters/thresholding.py +++ b/skimage/filters/thresholding.py @@ -1,16 +1,18 @@ -__all__ = ['threshold_adaptive', - 'threshold_otsu', - 'threshold_yen', - 'threshold_isodata', - '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 +__all__ = ['threshold_adaptive', + 'threshold_otsu', + 'threshold_yen', + 'threshold_isodata', + 'threshold_li', + 'threshold_minimum', + 'threshold_mean', + 'threshold_triangle'] + def threshold_adaptive(image, block_size, method='gaussian', offset=0, mode='reflect', param=None): @@ -101,7 +103,7 @@ def threshold_otsu(image, nbins=256): Parameters ---------- - image : (M, N) ndarray + image : (N, M) ndarray Grayscale input image. nbins : int, optional Number of bins used to calculate histogram. This value is ignored for @@ -110,8 +112,8 @@ def threshold_otsu(image, nbins=256): Returns ------- threshold : float - Upper threshold value. All pixels intensities that less or equal of - this value assumed as foreground. + Upper threshold value. All pixels with an intensity higher than + this value are assumed to be foreground. Raises ------ @@ -169,7 +171,7 @@ def threshold_yen(image, nbins=256): Parameters ---------- - image : array + image : (N, M) ndarray Input image. nbins : int, optional Number of bins used to calculate histogram. This value is ignored for @@ -178,8 +180,8 @@ def threshold_yen(image, nbins=256): Returns ------- threshold : float - Upper threshold value. All pixels intensities that less or equal of - this value assumed as foreground. + Upper threshold value. All pixels with an intensity higher than + this value are assumed to be foreground. References ---------- @@ -211,9 +213,8 @@ 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 helps to avoid '-inf' in crit. - # ImageJ Yen implementation replaces those values by zero. + # 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. crit = np.log(((P1_sq[:-1] * P2_sq[1:]) ** -1) * (P1[:-1] * (1.0 - P1[:-1])) ** 2) return bin_centers[crit.argmax()] @@ -237,7 +238,7 @@ def threshold_isodata(image, nbins=256, return_all=False): Parameters ---------- - image : array + image : (N, M) ndarray Input image. nbins : int, optional Number of bins used to calculate histogram. This value is ignored for @@ -332,13 +333,13 @@ def threshold_li(image): Parameters ---------- - image : array + image : (N, M) ndarray Input image. Returns ------- threshold : float - Upper threshold value. All pixels intensities more than + Upper threshold value. All pixels with an intensity higher than this value are assumed to be foreground. References @@ -346,8 +347,7 @@ 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 @@ -490,3 +490,108 @@ def threshold_minimum(image, nbins=256, bias='min', max_iter=10000): return bin_centers[upper_bound] elif bias == 'mid': return bin_centers[(threshold + upper_bound) // 2] + + +def threshold_mean(image): + """Return threshold value based on the mean of grayscale values. + + Parameters + ---------- + image : (N, M[, ..., P]) ndarray + Grayscale input image. + + Returns + ------- + threshold : float + Upper threshold value. All pixels with an intensity higher than + this value are assumed to be foreground. + + References + ---------- + .. [1] C. A. Glasbey, "An analysis of histogram-based thresholding + algorithms," CVGIP: Graphical Models and Image Processing, + vol. 55, pp. 532-537, 1993. + + Examples + -------- + >>> from skimage.data import camera + >>> image = camera() + >>> thresh = threshold_mean(image) + >>> binary = image > thresh + """ + return np.mean(image) + + +def threshold_triangle(image, nbins=256): + """Return threshold value based on the triangle algorithm. + + Parameters + ---------- + image : (N, M[, ..., P]) ndarray + Grayscale input image. + nbins : int, optional + Number of bins used to calculate histogram. This value is ignored for + integer arrays. + + Returns + ------- + threshold : float + Upper threshold value. All pixels with an intensity higher than + this value are assumed to be foreground. + + References + ---------- + .. [1] Zack, G. W., Rogers, W. E. and Latt, S. A., 1977, + Automatic Measurement of Sister Chromatid Exchange Frequency, + Journal of Histochemistry and Cytochemistry 25 (7), pp. 741-753 + + Examples + -------- + >>> from skimage.data import camera + >>> image = camera() + >>> thresh = threshold_triangle(image) + >>> binary = image > thresh + """ + # nbins is ignored for interger arrays + # so, we recalculate the effective nbins. + hist, bin_centers = histogram(image.ravel(), nbins) + nbins = bin_centers[-1] - bin_centers[0] + 1 + + # Find peak, lowest and highest gray levels. + arg_peak_height = np.argmax(hist) + peak_height = hist[arg_peak_height] + arg_low_level, arg_high_level = np.where(hist>0)[0][[0, -1]] + + # Flip is True if left tail is shorter. + flip = arg_peak_height - arg_low_level < arg_high_level - arg_peak_height + if flip: + hist = hist[::-1] + arg_low_level = nbins - arg_high_level - 1 + arg_peak_height = nbins - arg_peak_height - 1 + + # If flip == True, arg_high_level becomes incorrect + # but we don't need it anymore. + del(arg_high_level) + + # Set up the coordinate system. + width = arg_peak_height - arg_low_level + x1 = np.arange(width) + y1 = hist[x1 + arg_low_level] + + # Normalize. + norm = np.sqrt(peak_height**2 + width**2) + peak_height /= norm + width /= norm + + # Maximize the length. + d = peak_height * arg_low_level - width * hist[arg_low_level] + length = peak_height * x1 - width * y1 - d + level = np.argmax(length) + arg_low_level + + if flip: + level = nbins - level - 1 + + # The histogram doesn't start at zero, shift it. + level += bin_centers[0] + + return level