From 7fec88fc82982de1a8e4df959e905f77b7fa7977 Mon Sep 17 00:00:00 2001 From: radioxoma Date: Sun, 29 Dec 2013 18:29:51 +0300 Subject: [PATCH 1/8] fix threshold_yen import --- skimage/filter/__init__.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/skimage/filter/__init__.py b/skimage/filter/__init__.py index cef9b2e2..fed8a83b 100644 --- a/skimage/filter/__init__.py +++ b/skimage/filter/__init__.py @@ -9,7 +9,7 @@ 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 from . import rank @@ -37,6 +37,7 @@ __all__ = ['inverse', 'rank_order', 'gabor_kernel', 'gabor_filter', - 'threshold_otsu', 'threshold_adaptive', + 'threshold_otsu', + 'threshold_yen', 'rank'] From dfa2ba7bcd2946d2377336b4f71f95fa4926b44e Mon Sep 17 00:00:00 2001 From: radioxoma Date: Sun, 29 Dec 2013 21:21:06 +0300 Subject: [PATCH 2/8] fix yen blank image bug --- skimage/filter/tests/test_thresholding.py | 16 +++++++++++++++- skimage/filter/thresholding.py | 13 +++++++------ 2 files changed, 22 insertions(+), 7 deletions(-) diff --git a/skimage/filter/tests/test_thresholding.py b/skimage/filter/tests/test_thresholding.py index 0edfe4e7..d4fc36ed 100644 --- a/skimage/filter/tests/test_thresholding.py +++ b/skimage/filter/tests/test_thresholding.py @@ -43,10 +43,19 @@ 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_threshold_adaptive_generic(self): def func(arr): return arr.sum() / arr.shape[0] @@ -124,5 +133,10 @@ def test_yen_coins_image_as_float(): assert 0.43 < threshold_yen(coins) < 0.44 +def test_yen_camera_image(): + camera = skimage.img_as_ubyte(data.camera()) + assert 197 < threshold_yen(camera) < 199 + + if __name__ == '__main__': np.testing.run_module_suite() diff --git a/skimage/filter/thresholding.py b/skimage/filter/thresholding.py index 2ce8a3e1..c5c060e1 100644 --- a/skimage/filter/thresholding.py +++ b/skimage/filter/thresholding.py @@ -172,15 +172,16 @@ 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 + if bin_centers.size == 1: + return bin_centers[0] + + 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) # Get cumsum calculated from end of squared array: P2_sq = np.cumsum(norm_histo[::-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 threshold + return bin_centers[crit.argmax()] From 2c197846f9a7de782df909b132588837aaf82f90 Mon Sep 17 00:00:00 2001 From: radioxoma Date: Sun, 29 Dec 2013 21:58:13 +0300 Subject: [PATCH 3/8] Add ISODATA threshold with tests --- skimage/filter/__init__.py | 4 +- skimage/filter/tests/test_thresholding.py | 32 +++++++++-- skimage/filter/thresholding.py | 65 ++++++++++++++++++++++- 3 files changed, 96 insertions(+), 5 deletions(-) diff --git a/skimage/filter/__init__.py b/skimage/filter/__init__.py index fed8a83b..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_adaptive, threshold_otsu, threshold_yen +from .thresholding import (threshold_adaptive, threshold_otsu, threshold_yen, + threshold_isodata) from . import rank @@ -40,4 +41,5 @@ __all__ = ['inverse', '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 d4fc36ed..a0d39e74 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(): @@ -56,6 +57,16 @@ class TestSimpleImage(): 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] @@ -123,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 @@ -133,9 +149,19 @@ def test_yen_coins_image_as_float(): assert 0.43 < threshold_yen(coins) < 0.44 -def test_yen_camera_image(): +def test_isodata_camera_image(): camera = skimage.img_as_ubyte(data.camera()) - assert 197 < threshold_yen(camera) < 199 + 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 if __name__ == '__main__': diff --git a/skimage/filter/thresholding.py b/skimage/filter/thresholding.py index c5c060e1..f2de468f 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 @@ -185,3 +188,63 @@ def threshold_yen(image, nbins=256): crit = np.log(((P1_sq[:-1] * P2_sq[1:]) ** -1) * (P1[:-1] * (1.0 - P1[:-1])) ** 2) 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 float or int of any range. + nbins : int, optional + Number of bins used to calculate histogram. This value is ignored for + integer arrays. + + Returns + ------- + threshold : float64 or int64 + 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) + if bin_centers.size == 1: + return bin_centers[0] + # It is not necessary to calculate probability mass function in this case + # since in the l and h fractions it's reduced. + pmf = hist.astype(float)# / hist.sum() + cpmfl = np.cumsum(pmf, dtype=float) # Cumulative probability mass function + cpmfh = np.cumsum(pmf[::-1], dtype=float)[::-1] + + binnums = np.arange(pmf.size, dtype=np.uint8) + l = np.ma.divide(np.cumsum(pmf * binnums, dtype=float), cpmfl) + h = np.ma.divide(np.cumsum((pmf[::-1] * binnums[::-1]), dtype=float)[::-1], + cpmfh) + + allmean = (l + h) / 2.0 + threshold = bin_centers[np.nonzero(allmean.round() == binnums)[0][0]] + # ImageJ shows *inclusive* threshold. This implementation returns + # threshold, where `background <= threshold_value < foreground`. + return threshold From 08dae8e53c483d62651f3525b5b358be587fe621 Mon Sep 17 00:00:00 2001 From: radioxoma Date: Sun, 29 Dec 2013 22:09:14 +0300 Subject: [PATCH 4/8] contrib --- CONTRIBUTORS.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From 4d5889964bbda6f4df17449d7141e7937a7f1c08 Mon Sep 17 00:00:00 2001 From: radioxoma Date: Fri, 10 Jan 2014 22:25:48 +0300 Subject: [PATCH 5/8] Unit-tests and small fixes --- skimage/filter/tests/test_thresholding.py | 12 ++++++++++++ skimage/filter/thresholding.py | 23 ++++++++++++----------- 2 files changed, 24 insertions(+), 11 deletions(-) diff --git a/skimage/filter/tests/test_thresholding.py b/skimage/filter/tests/test_thresholding.py index a0d39e74..7d42fbcf 100644 --- a/skimage/filter/tests/test_thresholding.py +++ b/skimage/filter/tests/test_thresholding.py @@ -164,5 +164,17 @@ def test_isodata_moon_image(): assert threshold_isodata(moon) == 87 +def test_isodata_moon_image_negative_int(): + moon = data.moon().astype(np.int32) + moon -= 100 + assert threshold_isodata(moon) == -13 + + +def test_isodata_moon_image_negative_float(): + moon = 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 f2de468f..932089ba 100644 --- a/skimage/filter/thresholding.py +++ b/skimage/filter/thresholding.py @@ -198,14 +198,14 @@ def threshold_isodata(image, nbins=256): Parameters ---------- image : array - Input image float or int of any range. + Input image. nbins : int, optional Number of bins used to calculate histogram. This value is ignored for integer arrays. Returns ------- - threshold : float64 or int64 + threshold : float or int Upper threshold value. All pixels intensities that less or equal of this value assumed as background. @@ -233,18 +233,19 @@ def threshold_isodata(image, nbins=256): if bin_centers.size == 1: return bin_centers[0] # It is not necessary to calculate probability mass function in this case - # since in the l and h fractions it's reduced. - pmf = hist.astype(float)# / hist.sum() - cpmfl = np.cumsum(pmf, dtype=float) # Cumulative probability mass function - cpmfh = np.cumsum(pmf[::-1], dtype=float)[::-1] + # since the l and h fractions are reduced. + 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 = np.ma.divide(np.cumsum(pmf * binnums, dtype=float), cpmfl) - h = np.ma.divide(np.cumsum((pmf[::-1] * binnums[::-1]), dtype=float)[::-1], - cpmfh) + 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]] - # ImageJ shows *inclusive* threshold. This implementation returns - # threshold, where `background <= threshold_value < foreground`. + # This implementation returns threshold where + # `background <= threshold < foreground`. return threshold From cd01bad72134721d417591cc8cc5eb548dfabd4f Mon Sep 17 00:00:00 2001 From: radioxoma Date: Sat, 11 Jan 2014 02:00:03 +0300 Subject: [PATCH 6/8] fix tests --- skimage/filter/tests/test_thresholding.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skimage/filter/tests/test_thresholding.py b/skimage/filter/tests/test_thresholding.py index 7d42fbcf..f0a68f7d 100644 --- a/skimage/filter/tests/test_thresholding.py +++ b/skimage/filter/tests/test_thresholding.py @@ -165,13 +165,13 @@ def test_isodata_moon_image(): def test_isodata_moon_image_negative_int(): - moon = data.moon().astype(np.int32) + 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 = data.moon().astype(np.float64) + moon = skimage.img_as_ubyte(data.moon()).astype(np.float64) moon -= 100 assert -13 < threshold_isodata(moon) < -12 From 0079e8de5cf8eef57aa4cf102c7397e3861a4a1b Mon Sep 17 00:00:00 2001 From: radioxoma Date: Sun, 12 Jan 2014 18:18:50 +0300 Subject: [PATCH 7/8] Adding comments --- skimage/filter/thresholding.py | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/skimage/filter/thresholding.py b/skimage/filter/thresholding.py index 932089ba..86f380dd 100644 --- a/skimage/filter/thresholding.py +++ b/skimage/filter/thresholding.py @@ -132,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] @@ -175,14 +175,17 @@ def threshold_yen(image, nbins=256): >>> 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] - 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) + # 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) * @@ -205,7 +208,7 @@ def threshold_isodata(image, nbins=256): Returns ------- - threshold : float or int + threshold : float or int, corresponding input array dtype. Upper threshold value. All pixels intensities that less or equal of this value assumed as background. @@ -230,15 +233,19 @@ def threshold_isodata(image, nbins=256): >>> 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 probability mass function in this case # since the l and h fractions are reduced. - pmf = hist.astype(np.float32)# / hist.sum() + 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], From 38525fb4237a27fda87cacc2f523dccd202414f8 Mon Sep 17 00:00:00 2001 From: radioxoma Date: Mon, 13 Jan 2014 19:46:10 +0300 Subject: [PATCH 8/8] Explanation by ahojnnes --- skimage/filter/thresholding.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skimage/filter/thresholding.py b/skimage/filter/thresholding.py index 86f380dd..6ccf8f6b 100644 --- a/skimage/filter/thresholding.py +++ b/skimage/filter/thresholding.py @@ -237,8 +237,8 @@ def threshold_isodata(image, nbins=256): # 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 probability mass function in this case - # since the l and h fractions are reduced. + # 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]