From 31cf1acf24e5ed141bb45137dbe9a3aa75595831 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sat, 2 Aug 2014 02:12:15 -0400 Subject: [PATCH 1/7] Fix isodata-thresholding according to Zachary Pincus --- skimage/filter/thresholding.py | 87 +++++++++++++++++++++++----------- 1 file changed, 60 insertions(+), 27 deletions(-) diff --git a/skimage/filter/thresholding.py b/skimage/filter/thresholding.py index da74fe47..8a80ddb0 100644 --- a/skimage/filter/thresholding.py +++ b/skimage/filter/thresholding.py @@ -193,10 +193,19 @@ def threshold_yen(image, nbins=256): return bin_centers[crit.argmax()] -def threshold_isodata(image, nbins=256): - """Return threshold value based on ISODATA method. +def isodata(image, nbins=256, return_all=False): + """Return threshold value(s) based on ISODATA method. - Histogram-based threshold, known as Ridler-Calvard method or intermeans. + Histogram-based threshold, known as Ridler-Calvard method or inter-means. + Threshold values returned satisfy the following equality: + threshold = (image[image <= threshold].mean() + + image[image > threshold].mean()) / 2.0 + That is, returned thresholds are intensities that separate the image into + two groups of pixels, where the threshold intensity is midway between the + mean intensities of these groups. + + For integer images, the above equality holds to within one; for floating- + point images, the equality holds to within the histogram bin-width. Parameters ---------- @@ -205,12 +214,14 @@ def threshold_isodata(image, nbins=256): nbins : int, optional Number of bins used to calculate histogram. This value is ignored for integer arrays. + return_all: bool, optional + If False (default), return only the lowest threshold that satisfies + the above equality. If True, return all valid thresholds. 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. + threshold : float, int, array + Threshold value(s). References ---------- @@ -232,27 +243,49 @@ def threshold_isodata(image, nbins=256): >>> 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] + hist = hist.astype(np.float32) + # csuml and csumh contain the count of pixels in that bin or lower, and + # in all bins strictly higher than that bin, respectively + csuml = np.cumsum(hist) + csumh = np.cumsum(hist[::-1])[::-1] - hist - binnums = np.arange(pmf.size, dtype=np.min_scalar_type(nbins)) - # 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) + # intensity_sum contains the total pixel intensity from each bin + intensity_sum = hist * bin_centers - 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 + # l and h contain average value of all pixels in that bin or lower, and + # in all bins strictly higher than that bin, respectively. + # Note that since exp.histogram does not include empty bins at the low or + # high end of the range, csuml and csumh are strictly > 0, except in the + # last bin of csumh, which is zero by construction. + # So no worries about division by zero in the following lines, except + # for the last bin, but we can ignore that because no valid threshold + # can be in the top bin. So we just patch up csumh[-1] to not cause 0/0 + # errors. + csumh[-1] = 1 + l = np.cumsum(intensity_sum) / csuml + h = (np.cumsum(intensity_sum[::-1])[::-1] - intensity_sum) / csumh + + # isodata finds threshold values that meet the criterion t = (l + m)/2 + # where l is the mean of all pixels <= t and h is the mean of all pixels + # > t, as calculated above. So we are looking for places where + # (l + m) / 2 equals the intensity value for which those l and m figures + # were calculated -- which is, of course, the histogram bin centers. + # We only require this equality to be within the precision of the bin + # width, of course. + all_mean = (l + h) / 2.0 + bin_width = bin_centers[1] - bin_centers[0] + + # Look only at thresholds that are below the actual all_mean value, + # for consistency with the threshold being included in the lower pixel + # group. Otherwise can get thresholds that are not actually fixed-points + # of the isodata algorithm. For float images, this matters less, since + # there really can't be any guarantees anymore anyway. + distances = all_mean - bin_centers + thresholds = bin_centers[(distances >= 0) & (distances < bin_width)] + + if return_all: + return thresholds + else: + return thresholds[0] From 8aa5ef8697885d0a84dc0de28b44c3fce188bcef Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sat, 2 Aug 2014 09:10:50 -0400 Subject: [PATCH 2/7] Return to old function name threshold_isodata --- 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 8a80ddb0..3426e0a8 100644 --- a/skimage/filter/thresholding.py +++ b/skimage/filter/thresholding.py @@ -193,7 +193,7 @@ def threshold_yen(image, nbins=256): return bin_centers[crit.argmax()] -def isodata(image, nbins=256, return_all=False): +def threshold_isodata(image, nbins=256, return_all=False): """Return threshold value(s) based on ISODATA method. Histogram-based threshold, known as Ridler-Calvard method or inter-means. From fc085ec779707145cdba9c888c608b524ce6a952 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sat, 2 Aug 2014 09:42:28 -0400 Subject: [PATCH 3/7] Account for case when image only contains one unique value --- skimage/filter/thresholding.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/skimage/filter/thresholding.py b/skimage/filter/thresholding.py index 3426e0a8..0e908604 100644 --- a/skimage/filter/thresholding.py +++ b/skimage/filter/thresholding.py @@ -245,7 +245,16 @@ def threshold_isodata(image, nbins=256, return_all=False): """ hist, bin_centers = histogram(image, nbins) + + # image only contains one unique value + if len(bin_centers) == 1: + if return_all: + return bin_centers + else: + return bin_centers[0] + hist = hist.astype(np.float32) + # csuml and csumh contain the count of pixels in that bin or lower, and # in all bins strictly higher than that bin, respectively csuml = np.cumsum(hist) From ce0e9174c96e793474f7c44889d22155093e4b85 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sat, 2 Aug 2014 09:45:36 -0400 Subject: [PATCH 4/7] Fix test cases which are now wrong due to previous binning error --- skimage/filter/tests/test_thresholding.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/skimage/filter/tests/test_thresholding.py b/skimage/filter/tests/test_thresholding.py index 8464cb21..11e15575 100644 --- a/skimage/filter/tests/test_thresholding.py +++ b/skimage/filter/tests/test_thresholding.py @@ -157,7 +157,7 @@ def test_yen_coins_image_as_float(): def test_isodata_camera_image(): camera = skimage.img_as_ubyte(data.camera()) - assert threshold_isodata(camera) == 88 + assert threshold_isodata(camera) == 87 def test_isodata_coins_image(): @@ -167,19 +167,19 @@ def test_isodata_coins_image(): def test_isodata_moon_image(): moon = skimage.img_as_ubyte(data.moon()) - assert threshold_isodata(moon) == 87 + assert threshold_isodata(moon) == 86 def test_isodata_moon_image_negative_int(): moon = skimage.img_as_ubyte(data.moon()).astype(np.int32) moon -= 100 - assert threshold_isodata(moon) == -13 + assert threshold_isodata(moon) == -14 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 + assert -14 < threshold_isodata(moon) < -13 if __name__ == '__main__': From 23fd7b89aff9d8b2637ab39245d8617500ea64be Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sat, 2 Aug 2014 10:13:53 -0400 Subject: [PATCH 5/7] Improve test coverage --- skimage/filter/tests/test_thresholding.py | 68 +++++++++++++++++++---- 1 file changed, 56 insertions(+), 12 deletions(-) diff --git a/skimage/filter/tests/test_thresholding.py b/skimage/filter/tests/test_thresholding.py index 11e15575..5b103e6d 100644 --- a/skimage/filter/tests/test_thresholding.py +++ b/skimage/filter/tests/test_thresholding.py @@ -1,5 +1,5 @@ import numpy as np -from numpy.testing import assert_array_equal +from numpy.testing import assert_equal, assert_almost_equal import skimage from skimage import data @@ -59,19 +59,25 @@ class TestSimpleImage(): def test_isodata(self): assert threshold_isodata(self.image) == 2 + assert threshold_isodata(self.image, return_all=True) == [2] def test_isodata_blank_zero(self): image = np.zeros((5, 5), dtype=np.uint8) assert threshold_isodata(image) == 0 + assert threshold_isodata(image, return_all=True) == [0] def test_isodata_linspace(self): - assert -63.8 < threshold_isodata(np.linspace(-127, 0, 256)) < -63.6 + image = np.linspace(-127, 0, 256) + assert -63.8 < threshold_isodata(image) < -63.6 + assert_almost_equal(threshold_isodata(image, return_all=True), + [-63.74804688, -63.25195312]) def test_isodata_16bit(self): np.random.seed(0) imfloat = np.random.rand(256, 256) - t = threshold_isodata(imfloat, nbins=1024) - assert 0.49 < t < 0.51 + assert 0.49 < threshold_isodata(imfloat, nbins=1024) < 0.51 + assert all(0.49 < threshold_isodata(imfloat, nbins=1024, + return_all=True)) def test_threshold_adaptive_generic(self): def func(arr): @@ -84,7 +90,7 @@ class TestSimpleImage(): [ True, True, False, False, False]] ) out = threshold_adaptive(self.image, 3, method='generic', param=func) - assert_array_equal(ref, out) + assert_equal(ref, out) def test_threshold_adaptive_gaussian(self): ref = np.array( @@ -95,7 +101,7 @@ class TestSimpleImage(): [ True, True, False, False, False]] ) out = threshold_adaptive(self.image, 3, method='gaussian') - assert_array_equal(ref, out) + assert_equal(ref, out) def test_threshold_adaptive_mean(self): ref = np.array( @@ -106,7 +112,7 @@ class TestSimpleImage(): [ True, True, False, False, False]] ) out = threshold_adaptive(self.image, 3, method='mean') - assert_array_equal(ref, out) + assert_equal(ref, out) def test_threshold_adaptive_median(self): ref = np.array( @@ -117,7 +123,7 @@ class TestSimpleImage(): [False, True, False, False, False]] ) out = threshold_adaptive(self.image, 3, method='median') - assert_array_equal(ref, out) + assert_equal(ref, out) def test_otsu_camera_image(): @@ -157,30 +163,68 @@ def test_yen_coins_image_as_float(): def test_isodata_camera_image(): camera = skimage.img_as_ubyte(data.camera()) - assert threshold_isodata(camera) == 87 + + threshold = threshold_isodata(camera) + assert np.floor((camera[camera <= threshold].mean() + + camera[camera > threshold].mean()) / 2.0) == threshold + assert threshold == 87 + + assert threshold_isodata(camera, return_all=True) == [87] def test_isodata_coins_image(): coins = skimage.img_as_ubyte(data.coins()) - assert threshold_isodata(coins) == 107 + + threshold = threshold_isodata(coins) + assert np.floor((coins[coins <= threshold].mean() + + coins[coins > threshold].mean()) / 2.0) == threshold + assert threshold == 107 + + assert threshold_isodata(coins, return_all=True) == [107] def test_isodata_moon_image(): moon = skimage.img_as_ubyte(data.moon()) - assert threshold_isodata(moon) == 86 + + threshold = threshold_isodata(moon) + assert np.floor((moon[moon <= threshold].mean() + + moon[moon > threshold].mean()) / 2.0) == threshold + assert threshold == 86 + + thresholds = threshold_isodata(moon, return_all=True) + for threshold in thresholds: + assert np.floor((moon[moon <= threshold].mean() + + moon[moon > threshold].mean()) / 2.0) == threshold + assert_equal(thresholds, [86, 87, 88, 122, 123, 124, 139, 140]) def test_isodata_moon_image_negative_int(): moon = skimage.img_as_ubyte(data.moon()).astype(np.int32) moon -= 100 - assert threshold_isodata(moon) == -14 + + threshold = threshold_isodata(moon) + assert np.floor((moon[moon <= threshold].mean() + + moon[moon > threshold].mean()) / 2.0) == threshold + assert threshold == -14 + + thresholds = threshold_isodata(moon, return_all=True) + for threshold in thresholds: + assert np.floor((moon[moon <= threshold].mean() + + moon[moon > threshold].mean()) / 2.0) == threshold + assert_equal(thresholds, [-14, -13, -12, 22, 23, 24, 39, 40]) def test_isodata_moon_image_negative_float(): moon = skimage.img_as_ubyte(data.moon()).astype(np.float64) moon -= 100 + assert -14 < threshold_isodata(moon) < -13 + thresholds = threshold_isodata(moon, return_all=True) + assert_almost_equal(thresholds, + [-13.83789062, -12.84179688, -11.84570312, 22.02148438, + 23.01757812, 24.01367188, 38.95507812, 39.95117188]) + if __name__ == '__main__': np.testing.run_module_suite() From 83a6b6484481712025862086bc0fc6df5afa1dfd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sat, 2 Aug 2014 10:14:19 -0400 Subject: [PATCH 6/7] Fix return value description --- 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 0e908604..3985e278 100644 --- a/skimage/filter/thresholding.py +++ b/skimage/filter/thresholding.py @@ -220,7 +220,7 @@ def threshold_isodata(image, nbins=256, return_all=False): Returns ------- - threshold : float, int, array + threshold : float or int or array Threshold value(s). References From f9fd683fcbd8a5eeee1e4424fcc0b9840452fb6b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Fri, 8 Aug 2014 19:51:10 -0400 Subject: [PATCH 7/7] Test missed line for coverage --- skimage/filter/tests/test_thresholding.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/skimage/filter/tests/test_thresholding.py b/skimage/filter/tests/test_thresholding.py index 5b103e6d..69b56cde 100644 --- a/skimage/filter/tests/test_thresholding.py +++ b/skimage/filter/tests/test_thresholding.py @@ -103,6 +103,9 @@ class TestSimpleImage(): out = threshold_adaptive(self.image, 3, method='gaussian') assert_equal(ref, out) + out = threshold_adaptive(self.image, 3, method='gaussian', param=1.0 / 3.0) + assert_equal(ref, out) + def test_threshold_adaptive_mean(self): ref = np.array( [[False, False, False, False, True],