From a315ae718d50da67d40fa4c26e99605d74ed363c Mon Sep 17 00:00:00 2001 From: radioxoma Date: Mon, 12 Aug 2013 22:47:13 +0300 Subject: [PATCH 1/6] Add Yen threshold method --- skimage/filter/thresholding.py | 48 ++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/skimage/filter/thresholding.py b/skimage/filter/thresholding.py index 40e62d12..3258c468 100644 --- a/skimage/filter/thresholding.py +++ b/skimage/filter/thresholding.py @@ -133,3 +133,51 @@ def threshold_otsu(image, nbins=256): idx = np.argmax(variance12) threshold = bin_centers[:-1][idx] return threshold + + +def threshold_yen(image, nbins=256): + """Return threshold value based on Yen's method. + + Parameters + ---------- + image : array + Input image. + nbins : int + Number of bins used to calculate histogram. This value is ignored for + integer arrays. + + Returns + ------- + threshold : float + Threshold value. + + References + ---------- + .. [1] Yen J.C., Chang F.J., and Chang S. (1995) "A New Criterion + for Automatic Multilevel Thresholding" IEEE Trans. on Image + Processing, 4(3): 370-378 + .. [2] 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 + .. [3] ImageJ AutoThresholder code, http://fiji.sc/wiki/index.php/Auto_Threshold + + Examples + -------- + >>> from skimage.data import camera + >>> image = camera() + >>> thresh = threshold_yen(image) + >>> binary = image > thresh + """ + hist, bin_centers = histogram(img, nbins) + hist = hist.astype(float) + norm_histo = hist / hist.sum() # Probability mass function + P1 = np.cumsum(norm_histo) # Cumulative normalized histogram + P1_sq = np.cumsum(norm_histo ** 2) + 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. + # In ImageJ Yen implementation, all invalid values replaced by zero. + crit = -1*np.log(P1_sq[:-1]*P2_sq[1:]) + 2.0*np.log(P1[:-1]*(1.0-P1[:-1])) + max_crit = np.argmax(crit) + threshold = bin_centers[:-1][max_crit] + return threshold \ No newline at end of file From e41887c66d05ecd8c1dd247900d9621d4ae71ebe Mon Sep 17 00:00:00 2001 From: radioxoma Date: Tue, 13 Aug 2013 11:52:27 +0300 Subject: [PATCH 2/6] Codestyle improvement --- CONTRIBUTORS.txt | 3 +++ skimage/filter/thresholding.py | 6 +++--- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 28203676..33f763ba 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -148,3 +148,6 @@ - Matt Terry Color difference functions + +- Eugene Dvoretsky + Yen threshold implementation. diff --git a/skimage/filter/thresholding.py b/skimage/filter/thresholding.py index 3258c468..ce82500a 100644 --- a/skimage/filter/thresholding.py +++ b/skimage/filter/thresholding.py @@ -1,4 +1,4 @@ -__all__ = ['threshold_otsu', 'threshold_adaptive'] +__all__ = ['threshold_otsu', 'threshold_adaptive', 'threshold_yen'] import numpy as np import scipy.ndimage @@ -170,10 +170,10 @@ def threshold_yen(image, nbins=256): >>> binary = image > thresh """ hist, bin_centers = histogram(img, nbins) - hist = hist.astype(float) - norm_histo = hist / hist.sum() # Probability mass function + 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. # In ImageJ Yen implementation, all invalid values replaced by zero. From 8beb2ae605cbbfa7205610bee0ba6de13677cbfd Mon Sep 17 00:00:00 2001 From: radioxoma Date: Thu, 15 Aug 2013 23:03:54 +0300 Subject: [PATCH 3/6] Yen PEP8 and docstrings --- skimage/filter/thresholding.py | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/skimage/filter/thresholding.py b/skimage/filter/thresholding.py index ce82500a..10a3ab41 100644 --- a/skimage/filter/thresholding.py +++ b/skimage/filter/thresholding.py @@ -95,14 +95,15 @@ def threshold_otsu(image, nbins=256): ---------- image : array Input image. - nbins : int + nbins : int, optional Number of bins used to calculate histogram. This value is ignored for integer arrays. Returns ------- threshold : float - Threshold value. + Upper threshold value. All pixels intensities that less or equal of + this value assumed as foreground. References ---------- @@ -113,7 +114,7 @@ def threshold_otsu(image, nbins=256): >>> from skimage.data import camera >>> image = camera() >>> thresh = threshold_otsu(image) - >>> binary = image > thresh + >>> binary = image <= thresh """ hist, bin_centers = histogram(image, nbins) hist = hist.astype(float) @@ -142,14 +143,15 @@ def threshold_yen(image, nbins=256): ---------- image : array Input image. - nbins : int + nbins : int, optional Number of bins used to calculate histogram. This value is ignored for integer arrays. Returns ------- threshold : float - Threshold value. + Upper threshold value. All pixels intensities that less or equal of + this value assumed as foreground. References ---------- @@ -167,7 +169,7 @@ def threshold_yen(image, nbins=256): >>> from skimage.data import camera >>> image = camera() >>> thresh = threshold_yen(image) - >>> binary = image > thresh + >>> binary = image <= thresh """ hist, bin_centers = histogram(img, nbins) norm_histo = hist.astype(float) / hist.sum() # Probability mass function @@ -175,9 +177,10 @@ def threshold_yen(image, nbins=256): 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. - # In ImageJ Yen implementation, all invalid values replaced by zero. - crit = -1*np.log(P1_sq[:-1]*P2_sq[1:]) + 2.0*np.log(P1[:-1]*(1.0-P1[:-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 = -1 * np.log(P1_sq[:-1] * P2_sq[1:]) + \ + 2.0 * np.log(P1[:-1] * (1.0 - P1[:-1])) max_crit = np.argmax(crit) threshold = bin_centers[:-1][max_crit] - return threshold \ No newline at end of file + return threshold From 9e74f4eb70b368acf51aeaeea15482f1f73f12ab Mon Sep 17 00:00:00 2001 From: radioxoma Date: Fri, 16 Aug 2013 19:58:35 +0300 Subject: [PATCH 4/6] Added Yen tests, typo fixed --- skimage/filter/tests/test_thresholding.py | 34 ++++++++++++++++++++++- skimage/filter/thresholding.py | 4 +-- 2 files changed, 35 insertions(+), 3 deletions(-) diff --git a/skimage/filter/tests/test_thresholding.py b/skimage/filter/tests/test_thresholding.py index 97d3d9e3..0edfe4e7 100644 --- a/skimage/filter/tests/test_thresholding.py +++ b/skimage/filter/tests/test_thresholding.py @@ -3,7 +3,9 @@ from numpy.testing import assert_array_equal import skimage from skimage import data -from skimage.filter.thresholding import threshold_otsu, threshold_adaptive +from skimage.filter.thresholding import (threshold_adaptive, + threshold_otsu, + threshold_yen) class TestSimpleImage(): @@ -25,6 +27,26 @@ class TestSimpleImage(): image = np.float64(self.image) assert 2 <= threshold_otsu(image) < 3 + def test_yen(self): + assert threshold_yen(self.image) == 2 + + def test_yen_negative_int(self): + image = self.image - 2 + assert threshold_yen(image) == 0 + + def test_yen_float_image(self): + image = np.float64(self.image) + assert 2 <= threshold_yen(image) < 3 + + def test_yen_arange(self): + image = np.arange(256) + assert threshold_yen(image) == 127 + + def test_yen_binary(self): + image = np.zeros([2,256], dtype='uint8') + image[0] = 255 + assert threshold_yen(image) < 1 + def test_threshold_adaptive_generic(self): def func(arr): return arr.sum() / arr.shape[0] @@ -92,5 +114,15 @@ def test_otsu_lena_image(): assert 140 < threshold_otsu(lena) < 142 +def test_yen_coins_image(): + coins = skimage.img_as_ubyte(data.coins()) + assert 109 < threshold_yen(coins) < 111 + + +def test_yen_coins_image_as_float(): + coins = skimage.img_as_float(data.coins()) + assert 0.43 < threshold_yen(coins) < 0.44 + + if __name__ == '__main__': np.testing.run_module_suite() diff --git a/skimage/filter/thresholding.py b/skimage/filter/thresholding.py index 10a3ab41..8c78427b 100644 --- a/skimage/filter/thresholding.py +++ b/skimage/filter/thresholding.py @@ -1,4 +1,4 @@ -__all__ = ['threshold_otsu', 'threshold_adaptive', 'threshold_yen'] +__all__ = ['threshold_adaptive', 'threshold_otsu', 'threshold_yen'] import numpy as np import scipy.ndimage @@ -171,7 +171,7 @@ def threshold_yen(image, nbins=256): >>> thresh = threshold_yen(image) >>> binary = image <= thresh """ - hist, bin_centers = histogram(img, nbins) + 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) From 04804cc81ccc390e7e64e9a1fe3b214f11151d0d Mon Sep 17 00:00:00 2001 From: radioxoma Date: Sun, 18 Aug 2013 15:42:58 +0300 Subject: [PATCH 5/6] Logarithm optimisation --- skimage/filter/thresholding.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/skimage/filter/thresholding.py b/skimage/filter/thresholding.py index 8c78427b..ede3c76f 100644 --- a/skimage/filter/thresholding.py +++ b/skimage/filter/thresholding.py @@ -1,4 +1,4 @@ -__all__ = ['threshold_adaptive', 'threshold_otsu', 'threshold_yen'] +__all__ = ['threshold_otsu', 'threshold_adaptive', 'threshold_yen'] import numpy as np import scipy.ndimage @@ -179,8 +179,8 @@ def threshold_yen(image, nbins=256): 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 = -1 * np.log(P1_sq[:-1] * P2_sq[1:]) + \ - 2.0 * np.log(P1[:-1] * (1.0 - P1[:-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 From c8c95c0e633680e56c553a2c0ad9bf4ac3ee1538 Mon Sep 17 00:00:00 2001 From: radioxoma Date: Wed, 21 Aug 2013 20:18:34 +0300 Subject: [PATCH 6/6] Fix import list sequence --- skimage/filter/thresholding.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/filter/thresholding.py b/skimage/filter/thresholding.py index ede3c76f..bc052d8a 100644 --- a/skimage/filter/thresholding.py +++ b/skimage/filter/thresholding.py @@ -1,4 +1,4 @@ -__all__ = ['threshold_otsu', 'threshold_adaptive', 'threshold_yen'] +__all__ = ['threshold_adaptive', 'threshold_otsu', 'threshold_yen'] import numpy as np import scipy.ndimage