wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-08-07 11:28:14 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1090 from ahojnnes/isodata-fixes

Browse Source 
Fix isodata-thresholding according to Zachary Pincus
This commit is contained in:
Juan Nunez-Iglesias
2014-08-09 04:17:24 -05:00
parent 33a0e91120 f9fd683fcb
commit 6eb17a4786
2 changed files with 129 additions and 40 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/filter/tests/test_thresholding.py
+60 -13
View File
@@ -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,10 @@ class TestSimpleImage():
[ True, True, False, False, False]]
)
out = threshold_adaptive(self.image, 3, method='gaussian')
assert_array_equal(ref, out)
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(
@@ -106,7 +115,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 +126,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,29 +166,67 @@ def test_yen_coins_image_as_float():
def test_isodata_camera_image():
camera = skimage.img_as_ubyte(data.camera())
assert threshold_isodata(camera) == 88
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) == 87
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) == -13
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 -13 < threshold_isodata(moon) < -12
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__':
skimage/filter/thresholding.py
+69 -27
View File
@@ -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 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 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 or int or array
Threshold value(s).
References
----------
@@ -232,27 +243,58 @@ 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]
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)
# image only contains one unique value
if len(bin_centers) == 1:
if return_all:
return bin_centers
else:
return bin_centers[0]
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
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
# intensity_sum contains the total pixel intensity from each bin
intensity_sum = hist * bin_centers
# 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]