From e5e1918a2bc671f9fc5b58a5bed8c6c732b531fe Mon Sep 17 00:00:00 2001 From: "Josh Warner (Mac)" Date: Sun, 13 Oct 2013 12:38:35 -0500 Subject: [PATCH 1/5] 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() From f25ca3a835a7f90251fc28ac2e82bf55ef5fe4d5 Mon Sep 17 00:00:00 2001 From: "Josh Warner (Mac)" Date: Sat, 12 Oct 2013 00:41:33 -0500 Subject: [PATCH 2/5] ENH: `spacing` kwarg for random_walker and improved tests --- .../random_walker_segmentation.py | 18 ++++-- .../segmentation/tests/test_random_walker.py | 64 ++++++++++--------- 2 files changed, 46 insertions(+), 36 deletions(-) diff --git a/skimage/segmentation/random_walker_segmentation.py b/skimage/segmentation/random_walker_segmentation.py index 1fbec7ae..7c3916eb 100644 --- a/skimage/segmentation/random_walker_segmentation.py +++ b/skimage/segmentation/random_walker_segmentation.py @@ -169,13 +169,13 @@ 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, spacing, mask=None, beta=50, multichannel=False): - l_x, l_y, l_z = tuple(data.shape[i] * spacing[i] for i in range(3)) + 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, spacing, beta=beta, eps=1.e-10, multichannel=multichannel) @@ -257,7 +257,7 @@ def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, 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 + Spacing between voxels in each spatial dimension. If `None`, then the spacing between pixels/voxels in each dimension is assumed 1. Returns @@ -357,9 +357,17 @@ def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, '(see the docstrings)') if depth != 1.: warnings.warn('`depth` kwarg is deprecated, and will be removed in the' - ' next major version. Use `spacing` instead.') + ' next major release. Use `spacing` instead.') + + # Spacing kwarg checks if spacing is None: - spacing = (1., 1.) + (depth, ) + 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: diff --git a/skimage/segmentation/tests/test_random_walker.py b/skimage/segmentation/tests/test_random_walker.py index 050a8590..46a82d8e 100644 --- a/skimage/segmentation/tests/test_random_walker.py +++ b/skimage/segmentation/tests/test_random_walker.py @@ -8,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 @@ -28,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 @@ -102,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() @@ -236,19 +238,19 @@ def test_spacing(): # `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)) + data_aniso[i, :, :] = resize(data[:, 1, :], (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 - small_l // 2, - ly // 2 + small_l // 4, - lz // 2 - small_l // 4] = 2 + 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(np.rollaxis(data_aniso, 1).copy(), - np.rollaxis(labels_aniso, 1).copy(), + 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() From 927ba2cd8ec9c337989eda60352b6841d2495ac1 Mon Sep 17 00:00:00 2001 From: "Josh Warner (Mac)" Date: Sat, 12 Oct 2013 01:05:58 -0500 Subject: [PATCH 3/5] FIX: roll back incorrect testing change of gradient order --- skimage/segmentation/random_walker_segmentation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/segmentation/random_walker_segmentation.py b/skimage/segmentation/random_walker_segmentation.py index 7c3916eb..213cf94f 100644 --- a/skimage/segmentation/random_walker_segmentation.py +++ b/skimage/segmentation/random_walker_segmentation.py @@ -101,7 +101,7 @@ 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] + return np.r_[gr_deep, gr_right, gr_down] def _make_laplacian_sparse(edges, weights): From 03349250bcb90ad980c778d2a78786e1c4c4321e Mon Sep 17 00:00:00 2001 From: "Josh Warner (Mac)" Date: Sat, 12 Oct 2013 09:26:44 -0500 Subject: [PATCH 4/5] DOC: Add removal of `depth` from random_walker in 0.10 --- TODO.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/TODO.txt b/TODO.txt index 29d1ebac..5f5171ec 100644 --- a/TODO.txt +++ b/TODO.txt @@ -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 ----------- From 961c47e8013df3ccc7b9cb8766a32b194faa7dde Mon Sep 17 00:00:00 2001 From: "Josh Warner (Mac)" Date: Sun, 13 Oct 2013 12:34:30 -0500 Subject: [PATCH 5/5] DOC: Fix minor typo in docstring --- skimage/segmentation/random_walker_segmentation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/segmentation/random_walker_segmentation.py b/skimage/segmentation/random_walker_segmentation.py index 213cf94f..17c4c98b 100644 --- a/skimage/segmentation/random_walker_segmentation.py +++ b/skimage/segmentation/random_walker_segmentation.py @@ -282,7 +282,7 @@ def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, channels are separately normalized prior to running this algorithm. The `spacing` argument is specifically for anisotropic datasets, where - data points are spaced differently in one or more spatial dmensions. + 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