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

Merge pull request #775 from JDWarner/spacing_random_walker

Browse Source 
Add generalized anisotropic spacing support to random_walker
This commit is contained in:
Johannes Schönberger
2013-10-13 11:01:57 -07:00
parent 4b0cf8d27a 961c47e801
commit 0a8493805e
3 changed files with 140 additions and 48 deletions
Download Patch File Download Diff File Expand all files Collapse all files
TODO.txt
+1
View File
@@ -7,6 +7,7 @@ Version 0.10
* Change default mode of random_walker segmentation to 'cg_mg' > 'cg' > 'bf',
depending on which optional dependencies are available.
* Remove deprecated `out` parameter of `skimage.morphology.binary_*`
* Remove deprecated parameter `depth` in `skimage.segmentation.random_walker`
Version 0.9
-----------
skimage/segmentation/random_walker_segmentation.py
+42 -25
View File
@@ -77,14 +77,14 @@ def _make_graph_edges_3d(n_x, n_y, n_z):
return edges
def _compute_weights_3d(data, beta=130, eps=1.e-6, depth=1.,
def _compute_weights_3d(data, spacing, beta=130, eps=1.e-6,
multichannel=False):
# Weight calculation is main difference in multispectral version
# Original gradient**2 replaced with sum of gradients ** 2
gradients = 0
for channel in range(0, data.shape[-1]):
gradients += _compute_gradients_3d(data[..., channel],
depth=depth) ** 2
spacing) ** 2
# All channels considered together in this standard deviation
beta /= 10 * data.std()
if multichannel:
@@ -97,10 +97,10 @@ def _compute_weights_3d(data, beta=130, eps=1.e-6, depth=1.,
return weights
def _compute_gradients_3d(data, depth=1.):
gr_deep = np.abs(data[:, :, :-1] - data[:, :, 1:]).ravel() / depth
gr_right = np.abs(data[:, :-1] - data[:, 1:]).ravel()
gr_down = np.abs(data[:-1] - data[1:]).ravel()
def _compute_gradients_3d(data, spacing):
gr_deep = np.abs(data[:, :, :-1] - data[:, :, 1:]).ravel() / spacing[2]
gr_right = np.abs(data[:, :-1] - data[:, 1:]).ravel() / spacing[1]
gr_down = np.abs(data[:-1] - data[1:]).ravel() / spacing[0]
return np.r_[gr_deep, gr_right, gr_down]
@@ -116,9 +116,10 @@ def _make_laplacian_sparse(edges, weights):
lap = sparse.coo_matrix((data, (i_indices, j_indices)),
shape=(pixel_nb, pixel_nb))
connect = - np.ravel(lap.sum(axis=1))
lap = sparse.coo_matrix((np.hstack((data, connect)),
(np.hstack((i_indices, diag)), np.hstack((j_indices, diag)))),
shape=(pixel_nb, pixel_nb))
lap = sparse.coo_matrix(
(np.hstack((data, connect)), (np.hstack((i_indices, diag)),
np.hstack((j_indices, diag)))),
shape=(pixel_nb, pixel_nb))
return lap.tocsr()
@@ -168,14 +169,15 @@ def _mask_edges_weights(edges, weights, mask):
# Reassign edges labels to 0, 1, ... edges_number - 1
order = np.searchsorted(np.unique(edges.ravel()),
np.arange(max_node_index + 1))
edges = order[edges]
edges = order[edges.astype(np.int64)]
return edges, weights
def _build_laplacian(data, mask=None, beta=50, depth=1., multichannel=False):
l_x, l_y, l_z = data.shape[:3]
def _build_laplacian(data, spacing, mask=None, beta=50,
multichannel=False):
l_x, l_y, l_z = tuple(data.shape[i] for i in range(3))
edges = _make_graph_edges_3d(l_x, l_y, l_z)
weights = _compute_weights_3d(data, beta=beta, eps=1.e-10, depth=depth,
weights = _compute_weights_3d(data, spacing, beta=beta, eps=1.e-10,
multichannel=multichannel)
if mask is not None:
edges, weights = _mask_edges_weights(edges, weights, mask)
@@ -187,8 +189,9 @@ def _build_laplacian(data, mask=None, beta=50, depth=1., multichannel=False):
#----------- Random walker algorithm --------------------------------
def random_walker(data, labels, beta=130, mode=None, tol=1.e-3, copy=True,
multichannel=False, return_full_prob=False, depth=1.):
def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True,
multichannel=False, return_full_prob=False, depth=1.,
spacing=None):
"""Random walker algorithm for segmentation from markers.
Random walker algorithm is implemented for gray-level or multichannel
@@ -246,12 +249,16 @@ def random_walker(data, labels, beta=130, mode=None, tol=1.e-3, copy=True,
return_full_prob : bool, default False
If True, the probability that a pixel belongs to each of the labels
will be returned, instead of only the most likely label.
depth : float, default 1.
depth : float, default 1. [DEPRECATED]
Correction for non-isotropic voxel depths in 3D volumes.
Default (1.) implies isotropy. This factor is derived as follows:
depth = (out-of-plane voxel spacing) / (in-plane voxel spacing), where
in-plane voxel spacing represents the first two spatial dimensions and
out-of-plane voxel spacing represents the third spatial dimension.
`depth` is deprecated as of 0.9, in favor of `spacing`.
spacing : iterable of floats
Spacing between voxels in each spatial dimension. If `None`, then
the spacing between pixels/voxels in each dimension is assumed 1.
Returns
-------
@@ -274,12 +281,9 @@ def random_walker(data, labels, beta=130, mode=None, tol=1.e-3, copy=True,
Multichannel inputs are scaled with all channel data combined. Ensure all
channels are separately normalized prior to running this algorithm.
The `depth` argument is specifically for certain types of 3-dimensional
volumes which, due to how they were acquired, have different spacing
along in-plane and out-of-plane dimensions. This is commonly encountered
in medical imaging. The `depth` argument corrects gradients calculated
along the third spatial dimension for the otherwise inherent assumption
that all points are equally spaced.
The `spacing` argument is specifically for anisotropic datasets, where
data points are spaced differently in one or more spatial dimensions.
Anisotropic data is commonly encountered in medical imaging.
The algorithm was first proposed in *Random walks for image
segmentation*, Leo Grady, IEEE Trans Pattern Anal Mach Intell.
@@ -351,6 +355,19 @@ def random_walker(data, labels, beta=130, mode=None, tol=1.e-3, copy=True,
'random walker functions. You may also install pyamg '
'and run the random walker function in cg_mg mode '
'(see the docstrings)')
if depth != 1.:
warnings.warn('`depth` kwarg is deprecated, and will be removed in the'
' next major release. Use `spacing` instead.')
# Spacing kwarg checks
if spacing is None:
spacing = (1., 1.) + (depth, )
elif len(spacing) == 2:
spacing = tuple(spacing) + (depth, )
elif len(spacing) == 3:
pass
else:
raise ValueError('Input argument `spacing` incorrect, see docstring.')
# Parse input data
if not multichannel:
@@ -384,10 +401,10 @@ def random_walker(data, labels, beta=130, mode=None, tol=1.e-3, copy=True,
del filled
labels = np.atleast_3d(labels)
if np.any(labels < 0):
lap_sparse = _build_laplacian(data, mask=labels >= 0, beta=beta,
depth=depth, multichannel=multichannel)
lap_sparse = _build_laplacian(data, spacing, mask=labels >= 0,
beta=beta, multichannel=multichannel)
else:
lap_sparse = _build_laplacian(data, beta=beta, depth=depth,
lap_sparse = _build_laplacian(data, spacing, beta=beta,
multichannel=multichannel)
lap_sparse, B = _buildAB(lap_sparse, labels)
# We solve the linear system
skimage/segmentation/tests/test_random_walker.py
+97 -23
View File
@@ -1,5 +1,6 @@
import numpy as np
from skimage.segmentation import random_walker
from skimage.transform import resize
def make_2d_syntheticdata(lx, ly=None):
@@ -7,16 +8,16 @@ def make_2d_syntheticdata(lx, ly=None):
ly = lx
np.random.seed(1234)
data = np.zeros((lx, ly)) + 0.1 * np.random.randn(lx, ly)
small_l = int(lx / 5)
data[lx / 2 - small_l:lx / 2 + small_l,
ly / 2 - small_l:ly / 2 + small_l] = 1
data[lx / 2 - small_l + 1:lx / 2 + small_l - 1,
ly / 2 - small_l + 1:ly / 2 + small_l - 1] = \
0.1 * np.random.randn(2 * small_l - 2, 2 * small_l - 2)
data[lx / 2 - small_l, ly / 2 - small_l / 8:ly / 2 + small_l / 8] = 0
small_l = int(lx // 5)
data[lx // 2 - small_l:lx // 2 + small_l,
ly // 2 - small_l:ly // 2 + small_l] = 1
data[lx // 2 - small_l + 1:lx // 2 + small_l - 1,
ly // 2 - small_l + 1:ly // 2 + small_l - 1] = (
0.1 * np.random.randn(2 * small_l - 2, 2 * small_l - 2))
data[lx // 2 - small_l, ly // 2 - small_l // 8:ly // 2 + small_l // 8] = 0
seeds = np.zeros_like(data)
seeds[lx / 5, ly / 5] = 1
seeds[lx / 2 + small_l / 4, ly / 2 - small_l / 4] = 2
seeds[lx // 5, ly // 5] = 1
seeds[lx // 2 + small_l // 4, ly // 2 - small_l // 4] = 2
return data, seeds
@@ -27,21 +28,23 @@ def make_3d_syntheticdata(lx, ly=None, lz=None):
lz = lx
np.random.seed(1234)
data = np.zeros((lx, ly, lz)) + 0.1 * np.random.randn(lx, ly, lz)
small_l = int(lx / 5)
data[lx / 2 - small_l:lx / 2 + small_l,
ly / 2 - small_l:ly / 2 + small_l,
lz / 2 - small_l:lz / 2 + small_l] = 1
data[lx / 2 - small_l + 1:lx / 2 + small_l - 1,
ly / 2 - small_l + 1:ly / 2 + small_l - 1,
lz / 2 - small_l + 1:lz / 2 + small_l - 1] = 0
small_l = int(lx // 5)
data[lx // 2 - small_l:lx // 2 + small_l,
ly // 2 - small_l:ly // 2 + small_l,
lz // 2 - small_l:lz // 2 + small_l] = 1
data[lx // 2 - small_l + 1:lx // 2 + small_l - 1,
ly // 2 - small_l + 1:ly // 2 + small_l - 1,
lz // 2 - small_l + 1:lz // 2 + small_l - 1] = 0
# make a hole
hole_size = np.max([1, small_l / 8])
data[lx / 2 - small_l,
ly / 2 - hole_size:ly / 2 + hole_size,
lz / 2 - hole_size:lz / 2 + hole_size] = 0
hole_size = np.max([1, small_l // 8])
data[lx // 2 - small_l,
ly // 2 - hole_size:ly // 2 + hole_size,
lz // 2 - hole_size:lz // 2 + hole_size] = 0
seeds = np.zeros_like(data)
seeds[lx / 5, ly / 5, lz / 5] = 1
seeds[lx / 2 + small_l / 4, ly / 2 - small_l / 4, lz / 2 - small_l / 4] = 2
seeds[lx // 5, ly // 5, lz // 5] = 1
seeds[lx // 2 + small_l // 4,
ly // 2 - small_l // 4,
lz // 2 - small_l // 4] = 2
return data, seeds
@@ -101,7 +104,7 @@ def test_types():
lx = 70
ly = 100
data, labels = make_2d_syntheticdata(lx, ly)
data = 255 * (data - data.min()) / (data.max() - data.min())
data = 255 * (data - data.min()) // (data.max() - data.min())
data = data.astype(np.uint8)
labels_cg_mg = random_walker(data, labels, beta=90, mode='cg_mg')
assert (labels_cg_mg[25:45, 40:60] == 2).all()
@@ -181,6 +184,77 @@ def test_multispectral_3d():
return data, multi_labels, single_labels, labels
def test_depth():
n = 30
lx, ly, lz = n, n, n
data, _ = make_3d_syntheticdata(lx, ly, lz)
# Rescale `data` along Z axis
data_aniso = np.zeros((n, n, n // 2))
for i, yz in enumerate(data):
data_aniso[i, :, :] = resize(yz, (n, n // 2))
# Generate new labels
small_l = int(lx // 5)
labels_aniso = np.zeros_like(data_aniso)
labels_aniso[lx // 5, ly // 5, lz // 5] = 1
labels_aniso[lx // 2 + small_l // 4,
ly // 2 - small_l // 4,
lz // 4 - small_l // 8] = 2
# Test with `depth` kwarg
labels_aniso = random_walker(data_aniso, labels_aniso, mode='cg',
depth=0.5)
assert (labels_aniso[13:17, 13:17, 7:9] == 2).all()
def test_spacing():
n = 30
lx, ly, lz = n, n, n
data, _ = make_3d_syntheticdata(lx, ly, lz)
# Rescale `data` along Y axis
# `resize` is not yet 3D capable, so this must be done by looping in 2D.
data_aniso = np.zeros((n, n * 2, n))
for i, yz in enumerate(data):
data_aniso[i, :, :] = resize(yz, (n * 2, n))
# Generate new labels
small_l = int(lx // 5)
labels_aniso = np.zeros_like(data_aniso)
labels_aniso[lx // 5, ly // 5, lz // 5] = 1
labels_aniso[lx // 2 + small_l // 4,
ly - small_l // 2,
lz // 2 - small_l // 4] = 2
# Test with `spacing` kwarg
# First, anisotropic along Y
labels_aniso = random_walker(data_aniso, labels_aniso, mode='cg',
spacing=(1., 2., 1.))
assert (labels_aniso[13:17, 26:34, 13:17] == 2).all()
# Rescale `data` along X axis
# `resize` is not yet 3D capable, so this must be done by looping in 2D.
data_aniso = np.zeros((n, n * 2, n))
for i in range(data.shape[1]):
data_aniso[i, :, :] = resize(data[:, 1, :], (n * 2, n))
# Generate new labels
small_l = int(lx // 5)
labels_aniso2 = np.zeros_like(data_aniso)
labels_aniso2[lx // 5, ly // 5, lz // 5] = 1
labels_aniso2[lx - small_l // 2,
ly // 2 + small_l // 4,
lz // 2 - small_l // 4] = 2
# Anisotropic along X
labels_aniso2 = random_walker(data_aniso,
labels_aniso2,
mode='cg', spacing=(2., 1., 1.))
assert (labels_aniso2[26:34, 13:17, 13:17] == 2).all()
if __name__ == '__main__':
from numpy import testing
testing.run_module_suite()