diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 92d0d10e..00a76e59 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -151,7 +151,7 @@ Color difference functions - Eugene Dvoretsky - Yen threshold implementation. + Yen, Ridler-Calvard (ISODATA) threshold implementations. - Riaan van den Dool skimage.io plugin: GDAL diff --git a/skimage/filter/__init__.py b/skimage/filter/__init__.py index cef9b2e2..f4dea16f 100644 --- a/skimage/filter/__init__.py +++ b/skimage/filter/__init__.py @@ -9,7 +9,8 @@ from ._denoise import denoise_tv_chambolle from ._denoise_cy import denoise_bilateral, denoise_tv_bregman from ._rank_order import rank_order from ._gabor import gabor_kernel, gabor_filter -from .thresholding import threshold_otsu, threshold_adaptive +from .thresholding import (threshold_adaptive, threshold_otsu, threshold_yen, + threshold_isodata) from . import rank @@ -37,6 +38,8 @@ __all__ = ['inverse', 'rank_order', 'gabor_kernel', 'gabor_filter', - 'threshold_otsu', 'threshold_adaptive', + 'threshold_otsu', + 'threshold_yen', + 'threshold_isodata', 'rank'] diff --git a/skimage/filter/tests/test_thresholding.py b/skimage/filter/tests/test_thresholding.py index 0edfe4e7..f0a68f7d 100644 --- a/skimage/filter/tests/test_thresholding.py +++ b/skimage/filter/tests/test_thresholding.py @@ -5,7 +5,8 @@ import skimage from skimage import data from skimage.filter.thresholding import (threshold_adaptive, threshold_otsu, - threshold_yen) + threshold_yen, + threshold_isodata) class TestSimpleImage(): @@ -43,10 +44,29 @@ class TestSimpleImage(): assert threshold_yen(image) == 127 def test_yen_binary(self): - image = np.zeros([2,256], dtype='uint8') + image = np.zeros([2,256], dtype=np.uint8) image[0] = 255 assert threshold_yen(image) < 1 + def test_yen_blank_zero(self): + image = np.zeros((5, 5), dtype=np.uint8) + assert threshold_yen(image) == 0 + + def test_yen_blank_max(self): + image = np.empty((5, 5), dtype=np.uint8) + image.fill(255) + assert threshold_yen(image) == 255 + + def test_isodata(self): + assert threshold_isodata(self.image) == 2 + + def test_isodata_blank_zero(self): + image = np.zeros((5, 5), dtype=np.uint8) + assert threshold_isodata(image) == 0 + + def test_isodata_linspace(self): + assert -63.8 < threshold_isodata(np.linspace(-127, 0, 256)) < -63.6 + def test_threshold_adaptive_generic(self): def func(arr): return arr.sum() / arr.shape[0] @@ -114,6 +134,11 @@ def test_otsu_lena_image(): assert 140 < threshold_otsu(lena) < 142 +def test_yen_camera_image(): + camera = skimage.img_as_ubyte(data.camera()) + assert 197 < threshold_yen(camera) < 199 + + def test_yen_coins_image(): coins = skimage.img_as_ubyte(data.coins()) assert 109 < threshold_yen(coins) < 111 @@ -124,5 +149,32 @@ def test_yen_coins_image_as_float(): assert 0.43 < threshold_yen(coins) < 0.44 +def test_isodata_camera_image(): + camera = skimage.img_as_ubyte(data.camera()) + assert threshold_isodata(camera) == 88 + + +def test_isodata_coins_image(): + coins = skimage.img_as_ubyte(data.coins()) + assert threshold_isodata(coins) == 107 + + +def test_isodata_moon_image(): + moon = skimage.img_as_ubyte(data.moon()) + assert threshold_isodata(moon) == 87 + + +def test_isodata_moon_image_negative_int(): + moon = skimage.img_as_ubyte(data.moon()).astype(np.int32) + moon -= 100 + assert threshold_isodata(moon) == -13 + + +def test_isodata_moon_image_negative_float(): + moon = skimage.img_as_ubyte(data.moon()).astype(np.float64) + moon -= 100 + assert -13 < threshold_isodata(moon) < -12 + + if __name__ == '__main__': np.testing.run_module_suite() diff --git a/skimage/filter/thresholding.py b/skimage/filter/thresholding.py index 2ce8a3e1..6ccf8f6b 100644 --- a/skimage/filter/thresholding.py +++ b/skimage/filter/thresholding.py @@ -1,4 +1,7 @@ -__all__ = ['threshold_adaptive', 'threshold_otsu', 'threshold_yen'] +__all__ = ['threshold_adaptive', + 'threshold_otsu', + 'threshold_yen', + 'threshold_isodata'] import numpy as np import scipy.ndimage @@ -129,7 +132,7 @@ def threshold_otsu(image, nbins=256): # Clip ends to align class 1 and class 2 variables: # The last value of `weight1`/`mean1` should pair with zero values in # `weight2`/`mean2`, which do not exist. - variance12 = weight1[:-1] * weight2[1:] * (mean1[:-1] - mean2[1:])**2 + variance12 = weight1[:-1] * weight2[1:] * (mean1[:-1] - mean2[1:]) ** 2 idx = np.argmax(variance12) threshold = bin_centers[:-1][idx] @@ -172,15 +175,84 @@ def threshold_yen(image, nbins=256): >>> binary = image <= thresh """ hist, bin_centers = histogram(image, nbins) - norm_histo = hist.astype(float) / hist.sum() # Probability mass function - P1 = np.cumsum(norm_histo) # Cumulative normalized histogram - P1_sq = np.cumsum(norm_histo ** 2) + # 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] + + # Calculate probability mass function + pmf = hist.astype(np.float32) / hist.sum() + P1 = np.cumsum(pmf) # Cumulative normalized histogram + P1_sq = np.cumsum(pmf ** 2) # Get cumsum calculated from end of squared array: - P2_sq = np.cumsum(norm_histo[::-1] ** 2)[::-1] + 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. - crit = np.log(((P1_sq[:-1] * P2_sq[1:]) ** -1) * \ + crit = np.log(((P1_sq[:-1] * P2_sq[1:]) ** -1) * (P1[:-1] * (1.0 - P1[:-1])) ** 2) - max_crit = np.argmax(crit) - threshold = bin_centers[:-1][max_crit] + return bin_centers[crit.argmax()] + + +def threshold_isodata(image, nbins=256): + """Return threshold value based on ISODATA method. + + Histogram-based threshold, known as Ridler-Calvard method or intermeans. + + Parameters + ---------- + image : array + Input image. + nbins : int, optional + Number of bins used to calculate histogram. This value is ignored for + integer arrays. + + Returns + ------- + threshold : float or int, corresponding input array dtype. + Upper threshold value. All pixels intensities that less or equal of + this value assumed as background. + + References + ---------- + .. [1] Ridler, TW & Calvard, S (1978), "Picture thresholding using an + iterative selection method" + .. [2] IEEE Transactions on Systems, Man and Cybernetics 8: 630-632, + http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=4310039 + .. [3] Sezgin M. and Sankur B. (2004) "Survey over Image Thresholding + Techniques and Quantitative Performance Evaluation" Journal of + Electronic Imaging, 13(1): 146-165, + http://www.busim.ee.boun.edu.tr/~sankur/SankurFolder/Threshold_survey.pdf + .. [4] ImageJ AutoThresholder code, + http://fiji.sc/wiki/index.php/Auto_Threshold + + Examples + -------- + >>> from skimage.data import coins + >>> image = coins() + >>> thresh = threshold_isodata(image) + >>> binary = image > thresh + """ + hist, bin_centers = histogram(image, 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] + # It is not necessary to calculate the probability mass function here, + # because the l and h fractions already include the normalization. + pmf = hist.astype(np.float32) # / hist.sum() + cpmfl = np.cumsum(pmf, dtype=np.float32) + cpmfh = np.cumsum(pmf[::-1], dtype=np.float32)[::-1] + + binnums = np.arange(pmf.size, dtype=np.uint8) + # l and h contain average value of pixels in sum of bins, calculated + # from lower to higher and from higher to lower respectively. + l = np.ma.divide(np.cumsum(pmf * binnums, dtype=np.float32), cpmfl) + h = np.ma.divide( + np.cumsum((pmf[::-1] * binnums[::-1]), dtype=np.float32)[::-1], + cpmfh) + + allmean = (l + h) / 2.0 + threshold = bin_centers[np.nonzero(allmean.round() == binnums)[0][0]] + # This implementation returns threshold where + # `background <= threshold < foreground`. return threshold