From e5e1918a2bc671f9fc5b58a5bed8c6c732b531fe Mon Sep 17 00:00:00 2001 From: "Josh Warner (Mac)" Date: Sun, 13 Oct 2013 12:38:35 -0500 Subject: [PATCH] REBASE: Resolve first conflict --- .../random_walker_segmentation.py | 59 ++++++++------- .../segmentation/tests/test_random_walker.py | 72 +++++++++++++++++++ 2 files changed, 106 insertions(+), 25 deletions(-) diff --git a/skimage/segmentation/random_walker_segmentation.py b/skimage/segmentation/random_walker_segmentation.py index c475dcb7..1fbec7ae 100644 --- a/skimage/segmentation/random_walker_segmentation.py +++ b/skimage/segmentation/random_walker_segmentation.py @@ -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,11 +97,11 @@ 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() - return np.r_[gr_deep, gr_right, gr_down] +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_down, gr_right, gr_deep] def _make_laplacian_sparse(edges, weights): @@ -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() @@ -172,10 +173,11 @@ def _mask_edges_weights(edges, weights, mask): 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] * spacing[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 dmensions. + 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,11 @@ 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 version. Use `spacing` instead.') + if spacing is None: + spacing = (1., 1.) + (depth, ) # Parse input data if not multichannel: @@ -384,10 +393,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 diff --git a/skimage/segmentation/tests/test_random_walker.py b/skimage/segmentation/tests/test_random_walker.py index 1cc0a1ee..050a8590 100644 --- a/skimage/segmentation/tests/test_random_walker.py +++ b/skimage/segmentation/tests/test_random_walker.py @@ -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): @@ -181,6 +182,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 * 1.5, 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 - small_l // 2, + ly // 2 + small_l // 4, + lz // 2 - small_l // 4] = 2 + + # Anisotropic along X + labels_aniso2 = random_walker(np.rollaxis(data_aniso, 1).copy(), + np.rollaxis(labels_aniso, 1).copy(), + 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()