From 8955dad32ecc2caeda2c4448b20eaf0c07d0675a Mon Sep 17 00:00:00 2001 From: JDWarner Date: Mon, 27 Aug 2012 11:42:30 -0500 Subject: [PATCH 1/6] Multispectral modifications applied to random walker. --- .../random_walker_segmentation.py | 76 +++++++++++++------ 1 file changed, 54 insertions(+), 22 deletions(-) diff --git a/skimage/segmentation/random_walker_segmentation.py b/skimage/segmentation/random_walker_segmentation.py index 46f8dfb3..49490dda 100644 --- a/skimage/segmentation/random_walker_segmentation.py +++ b/skimage/segmentation/random_walker_segmentation.py @@ -64,17 +64,30 @@ def _make_graph_edges_3d(n_x, n_y, n_z): return edges -def _compute_weights_3d(data, beta=130, eps=1.e-6): - gradients = _compute_gradients_3d(data)**2 - beta /= 10 * data.std() +def _compute_weights_3d(data, beta=130, eps=1.e-6, depth=1.): + # Weight calculation is main difference in multispectral version + # Original gradient**2 replaced with sqrt( sum of gradients**2 ) + for i, spectrum in enumerate(data): + if i == 0: + gradients = _compute_gradients_3d(spectrum, depth=depth)**2 + else: + gradients += _compute_gradients_3d(spectrum)**2 + + gradients = np.sqrt(gradients) + + # New final term in beta to give == results in trivial case where + # multiple identical spectra are passed. + # It may be faster and/or more memory efficient do an approximate + # std() combining spectrum.std() together than this 2nd term. + beta /= 10 * np.asarray(data).std() * np.sqrt( len(data) ) gradients *= beta weights = np.exp(- gradients) weights += eps return weights -def _compute_gradients_3d(data): - gr_deep = np.abs(data[:, :, :-1] - data[:, :, 1:]).ravel() +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] @@ -148,10 +161,10 @@ def _mask_edges_weights(edges, weights, mask): return edges, weights -def _build_laplacian(data, mask=None, beta=50): - l_x, l_y, l_z = data.shape +def _build_laplacian(data, mask=None, beta=50, depth=1.): + l_x, l_y, l_z = data[0].shape edges = _make_graph_edges_3d(l_x, l_y, l_z) - weights = _compute_weights_3d(data, beta=beta, eps=1.e-10) + weights = _compute_weights_3d(data, beta=beta, eps=1.e-10, depth=depth) if mask is not None: edges, weights = _mask_edges_weights(edges, weights, mask) lap = _make_laplacian_sparse(edges, weights) @@ -162,17 +175,19 @@ def _build_laplacian(data, mask=None, beta=50): #----------- Random walker algorithm -------------------------------- -def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, - return_full_prob=False): +def random_walker(data, labels, beta=130, depth=1., mode='bf', tol=1.e-3, + copy=True, return_full_prob=False): """ - Random walker algorithm for segmentation from markers. + Multispectral random walker algorithm for segmentation from markers. Parameters ---------- - data : array_like - Image to be segmented in phases. `data` can be two- or - three-dimensional. + data : array_like OR iterable of arrays + Image to be segmented in phases. Non-multispectral `data` can be + two- or three-dimensional; multispectral data is provided as an + iterable of like-sized 2D or 3D arrays. Data spacing is assumed + isotropic unless depth kwarg is used. labels : array of ints, of same shape as `data` Array of seed markers labeled with different positive integers @@ -186,6 +201,11 @@ def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, Penalization coefficient for the random walker motion (the greater `beta`, the more difficult the diffusion). + depth : float, default 1. + Correction for non-isotropic voxel depths in 3D volumes. + Default (1.) implies isotropy. This factor is derived as follows: + depth = (slice thickness) / (in-plane voxel dimension) + mode : {'bf', 'cg_mg', 'cg'} (default: 'bf') Mode for solving the linear system in the random walker algorithm. @@ -299,9 +319,19 @@ def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, [ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.]]) """ - # We work with 3-D arrays of floats - data = img_as_float(data) - data = np.atleast_3d(data) + # Parse input data + if isinstance( data, np.ndarray ): + # Wrap into single-element list + data = [ np.atleast_3d( img_as_float(data) ) ] + else: + # We work with 3-D arrays of floats + newdata = [] + for spectrum in data: + newdata.append( np.atleast_3d( img_as_float(spectrum) ) ) + del data + data = newdata + del newdata + if copy: labels = np.copy(labels) label_values = np.unique(labels) @@ -318,9 +348,10 @@ def random_walker(data, labels, beta=130, mode='bf', 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) + lap_sparse = _build_laplacian(data, mask=labels >= 0, beta=beta, + depth=depth) else: - lap_sparse = _build_laplacian(data, beta=beta) + lap_sparse = _build_laplacian(data, beta=beta, depth=depth) lap_sparse, B = _buildAB(lap_sparse, labels) # We solve the linear system # lap_sparse X = B @@ -343,18 +374,19 @@ def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, X = _solve_bf(lap_sparse, B, return_full_prob=return_full_prob) # Clean up results - data = np.squeeze(data) + for spectrum in data: + spectrum = spectrum.squeeze() if return_full_prob: labels = labels.astype(np.float) X = np.array([_clean_labels_ar(Xline, labels, - copy=True).reshape(data.shape) for Xline in X]) + copy=True).reshape(data[0].shape) for Xline in X]) for i in range(1, int(labels.max()) + 1): mask_i = np.squeeze(labels == i) X[i - 1, mask_i] = 1 X[np.setdiff1d(np.arange(0, labels.max(), dtype=np.int), [i - 1]), mask_i] = 0 else: - X = _clean_labels_ar(X + 1, labels).reshape(data.shape) + X = _clean_labels_ar(X + 1, labels).reshape(data[0].shape) return X From feca48cc49215d1813020859ddb81ad6c9096c33 Mon Sep 17 00:00:00 2001 From: JDWarner Date: Mon, 27 Aug 2012 11:48:32 -0500 Subject: [PATCH 2/6] Added return_full_prob kwarg to solve call if pyamg not present. --- skimage/segmentation/random_walker_segmentation.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/skimage/segmentation/random_walker_segmentation.py b/skimage/segmentation/random_walker_segmentation.py index 49490dda..4eb86ff7 100644 --- a/skimage/segmentation/random_walker_segmentation.py +++ b/skimage/segmentation/random_walker_segmentation.py @@ -366,7 +366,8 @@ def random_walker(data, labels, beta=130, depth=1., mode='bf', tol=1.e-3, """pyamg (http://code.google.com/p/pyamg/)) is needed to use this mode, but is not installed. The 'cg' mode will be used instead.""") - X = _solve_cg(lap_sparse, B, tol=tol) + X = _solve_cg(lap_sparse, B, tol=tol, + return_full_prob=return_full_prob) else: X = _solve_cg_mg(lap_sparse, B, tol=tol, return_full_prob=return_full_prob) From 682d0535cd31b5af65237adacf6d9fb5a12885fd Mon Sep 17 00:00:00 2001 From: JDWarner Date: Mon, 27 Aug 2012 13:41:41 -0500 Subject: [PATCH 3/6] Added multispectral random walker test. Since the multispectral path is equivalent except for gradient calcs, only one test case is needed. This test is modeled on the 3-D non-multispectral version. If deemed necessary, adding a 2-D case would be simple. --- skimage/segmentation/tests/test_random_walker.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/skimage/segmentation/tests/test_random_walker.py b/skimage/segmentation/tests/test_random_walker.py index 4e5a74c3..e0415498 100644 --- a/skimage/segmentation/tests/test_random_walker.py +++ b/skimage/segmentation/tests/test_random_walker.py @@ -138,6 +138,18 @@ def test_3d_inactive(): assert (labels.reshape(data.shape)[13:17, 13:17, 13:17] == 2).all() return data, labels, old_labels, after_labels + +def test_multispectral(): + n = 30 + lx, ly, lz = n, n, n + data, labels = make_3d_syntheticdata( lx, ly, lz ) + data = [data, data] # Result should be identical + multi_labels = random_walker(data, labels, mode='cg') + single_labels = random_walker(data[0], labels, mode='cg') + assert (multi_labels.reshape(data[0].shape)[13:17, 13:17, 13:17] == 2).all() + assert (single_labels.reshape(data[0].shape)[13:17, 13:17, 13:17] == 2).all() + return data, multi_labels, single_labels, labels + if __name__ == '__main__': from numpy import testing testing.run_module_suite() From 61320957eb1fb3c6c50094fb419bca497e47d73f Mon Sep 17 00:00:00 2001 From: JDWarner Date: Wed, 29 Aug 2012 16:33:56 -0500 Subject: [PATCH 4/6] Changes based on PR review recommendations: input format, scaling, and bugfix. In this new version, all instances of 'spectrum' have been replaced with 'channel'. The documentation also reflects this change, and the new multichannel kwarg used to indicate multichannel input is named appropriately. New boolean multichannel kwarg added, which controls if the input has multiple channels or not. Input 'data' is now array_like for both gray-level and multichannel. This kwarg is needed mainly because a 3-D array could be either 3 spatial dimensions or a set of different 2-D channels. New scaling kwarg added (may be removed in future), controlling if data scaling is applied to ALL channels or each channel individually, if multichannel=True. No effect for gray-level data. Removed np.sqrt(gradients) in _compute_weights_3d(), which was a bug. Tests now pass consistently. New method for maintaining shape from input to output, where dims = data.shape prior to np.atleast_3d(). A theoretical (70,100,1) array passed should now result in a (70,100,1) shaped output, for example. Updated and fixed multispectral test script to work with new version. TODO: Additional test(s) likely needed to cover code branches from new kwargs. --- .../random_walker_segmentation.py | 112 +++++++++++------- .../segmentation/tests/test_random_walker.py | 11 +- 2 files changed, 78 insertions(+), 45 deletions(-) diff --git a/skimage/segmentation/random_walker_segmentation.py b/skimage/segmentation/random_walker_segmentation.py index 4eb86ff7..239c4dd8 100644 --- a/skimage/segmentation/random_walker_segmentation.py +++ b/skimage/segmentation/random_walker_segmentation.py @@ -64,22 +64,29 @@ 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, beta=130, eps=1.e-6, depth=1., + multichannel=False): # Weight calculation is main difference in multispectral version # Original gradient**2 replaced with sqrt( sum of gradients**2 ) - for i, spectrum in enumerate(data): - if i == 0: - gradients = _compute_gradients_3d(spectrum, depth=depth)**2 - else: - gradients += _compute_gradients_3d(spectrum)**2 + if not multichannel: + gradients = _compute_gradients_3d( data, depth=depth )**2 + else: + for channel in range(data.shape[-1]): + if channel == 0: + gradients = _compute_gradients_3d(data[..., channel], + depth=depth)**2 + else: + gradients += _compute_gradients_3d(data[..., channel], + depth=depth)**2 - gradients = np.sqrt(gradients) + # gradients = np.sqrt(gradients) - # New final term in beta to give == results in trivial case where - # multiple identical spectra are passed. - # It may be faster and/or more memory efficient do an approximate - # std() combining spectrum.std() together than this 2nd term. - beta /= 10 * np.asarray(data).std() * np.sqrt( len(data) ) + # All channels considered together in this standard deviation + beta /= 10 * data.std() + if multichannel: + # New final term in beta to give == results in trivial case where + # multiple identical spectra are passed. + beta /= np.sqrt( data.shape[-1] ) gradients *= beta weights = np.exp(- gradients) weights += eps @@ -161,10 +168,14 @@ def _mask_edges_weights(edges, weights, mask): return edges, weights -def _build_laplacian(data, mask=None, beta=50, depth=1.): - l_x, l_y, l_z = data[0].shape +def _build_laplacian(data, mask=None, beta=50, depth=1., multichannel=False): + if not multichannel: + l_x, l_y, l_z = data.shape + else: + l_x, l_y, l_z = data.shape[0], data.shape[1], data.shape[2] 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, beta=beta, eps=1.e-10, depth=depth, + multichannel=multichannel) if mask is not None: edges, weights = _mask_edges_weights(edges, weights, mask) lap = _make_laplacian_sparse(edges, weights) @@ -175,19 +186,21 @@ def _build_laplacian(data, mask=None, beta=50, depth=1.): #----------- Random walker algorithm -------------------------------- -def random_walker(data, labels, beta=130, depth=1., mode='bf', tol=1.e-3, - copy=True, return_full_prob=False): +def random_walker(data, labels, beta=130, depth=1., mode='bf', tol=1.e-3, + copy=True, multichannel=False, scaling='all', + return_full_prob=False): """ - Multispectral random walker algorithm for segmentation from markers. + Multichannel random walker algorithm for segmentation from markers. Parameters ---------- - data : array_like OR iterable of arrays - Image to be segmented in phases. Non-multispectral `data` can be - two- or three-dimensional; multispectral data is provided as an - iterable of like-sized 2D or 3D arrays. Data spacing is assumed - isotropic unless depth kwarg is used. + data : array_like + Image to be segmented in phases. Gray-level`data` can be two- or + three-dimensional; multichannel data can be three- or four- + dimensional (requires multichannel=True) with the highest + dimension denoting channels. Data spacing is assumed isotropic + unless depth keyword argument is used. labels : array of ints, of same shape as `data` Array of seed markers labeled with different positive integers @@ -236,6 +249,21 @@ def random_walker(data, labels, beta=130, depth=1., mode='bf', tol=1.e-3, the result of the segmentation. Use copy=False if you want to save on memory. + multichannel : bool, default False + If True, input data is parsed as multichannel data (see 'data' above + for proper input format in this case) + + scaling : string, default 'all' + Controls input scaling if multichannel=True (otherwise no effect). + + - 'all' (default): Data from all channels is combined when scaling + input data to the range [0,1] as type np.float64. Recommended + option for RGB(A) inputs. + + - 'separate': Each channel is scaled individually, separate from the + others, to the range [0,1]. Select this if the channels are very + different, for example if one were x-ray CT and another MRI data. + 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. @@ -320,17 +348,22 @@ def random_walker(data, labels, beta=130, depth=1., mode='bf', tol=1.e-3, """ # Parse input data - if isinstance( data, np.ndarray ): - # Wrap into single-element list - data = [ np.atleast_3d( img_as_float(data) ) ] - else: + if not multichannel: # We work with 3-D arrays of floats - newdata = [] - for spectrum in data: - newdata.append( np.atleast_3d( img_as_float(spectrum) ) ) - del data - data = newdata - del newdata + dims = data.shape + data = np.atleast_3d( img_as_float(data) ) + else: + dims = data[..., 0].shape + data = np.atleast_3d( data ) # Should never be needed + if scaling.lower().strip() == 'all': + data = img_as_float( data ) + else: + newdata = np.zeros(data.shape, dtype=np.float64) + for channel in range( data.shape[-1] ): + newdata[..., channel] = img_as_float( data[..., channel] ) + del data + data = newdata + del newdata if copy: labels = np.copy(labels) @@ -349,9 +382,10 @@ def random_walker(data, labels, beta=130, depth=1., mode='bf', tol=1.e-3, labels = np.atleast_3d(labels) if np.any(labels < 0): lap_sparse = _build_laplacian(data, mask=labels >= 0, beta=beta, - depth=depth) + depth=depth, multichannel=multichannel) else: - lap_sparse = _build_laplacian(data, beta=beta, depth=depth) + lap_sparse = _build_laplacian(data, beta=beta, depth=depth, + multichannel=multichannel) lap_sparse, B = _buildAB(lap_sparse, labels) # We solve the linear system # lap_sparse X = B @@ -366,7 +400,7 @@ def random_walker(data, labels, beta=130, depth=1., mode='bf', tol=1.e-3, """pyamg (http://code.google.com/p/pyamg/)) is needed to use this mode, but is not installed. The 'cg' mode will be used instead.""") - X = _solve_cg(lap_sparse, B, tol=tol, + X = _solve_cg(lap_sparse, B, tol=tol, return_full_prob=return_full_prob) else: X = _solve_cg_mg(lap_sparse, B, tol=tol, @@ -375,19 +409,17 @@ def random_walker(data, labels, beta=130, depth=1., mode='bf', tol=1.e-3, X = _solve_bf(lap_sparse, B, return_full_prob=return_full_prob) # Clean up results - for spectrum in data: - spectrum = spectrum.squeeze() if return_full_prob: labels = labels.astype(np.float) X = np.array([_clean_labels_ar(Xline, labels, - copy=True).reshape(data[0].shape) for Xline in X]) + copy=True).reshape(dims) for Xline in X]) for i in range(1, int(labels.max()) + 1): mask_i = np.squeeze(labels == i) X[i - 1, mask_i] = 1 X[np.setdiff1d(np.arange(0, labels.max(), dtype=np.int), [i - 1]), mask_i] = 0 else: - X = _clean_labels_ar(X + 1, labels).reshape(data[0].shape) + X = _clean_labels_ar(X + 1, labels).reshape(dims) return X diff --git a/skimage/segmentation/tests/test_random_walker.py b/skimage/segmentation/tests/test_random_walker.py index e0415498..74397dbe 100644 --- a/skimage/segmentation/tests/test_random_walker.py +++ b/skimage/segmentation/tests/test_random_walker.py @@ -143,11 +143,12 @@ def test_multispectral(): n = 30 lx, ly, lz = n, n, n data, labels = make_3d_syntheticdata( lx, ly, lz ) - data = [data, data] # Result should be identical - multi_labels = random_walker(data, labels, mode='cg') - single_labels = random_walker(data[0], labels, mode='cg') - assert (multi_labels.reshape(data[0].shape)[13:17, 13:17, 13:17] == 2).all() - assert (single_labels.reshape(data[0].shape)[13:17, 13:17, 13:17] == 2).all() + data.shape += (1,) + data = data.repeat(2, axis=3) # Result should be identical + multi_labels = random_walker(data, labels, mode='cg', multichannel=True) + single_labels = random_walker(data[:,:,:,0], labels, mode='cg') + assert (multi_labels.reshape(labels.shape)[13:17, 13:17, 13:17] == 2).all() + assert (single_labels.reshape(labels.shape)[13:17, 13:17, 13:17] == 2).all() return data, multi_labels, single_labels, labels if __name__ == '__main__': From 99238c44a594da3eff6666aee597ea69ec464c88 Mon Sep 17 00:00:00 2001 From: JDWarner Date: Wed, 29 Aug 2012 17:07:13 -0500 Subject: [PATCH 5/6] sqrt(gradients) line removed --- skimage/segmentation/random_walker_segmentation.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/skimage/segmentation/random_walker_segmentation.py b/skimage/segmentation/random_walker_segmentation.py index 239c4dd8..6165770d 100644 --- a/skimage/segmentation/random_walker_segmentation.py +++ b/skimage/segmentation/random_walker_segmentation.py @@ -78,9 +78,6 @@ def _compute_weights_3d(data, beta=130, eps=1.e-6, depth=1., else: gradients += _compute_gradients_3d(data[..., channel], depth=depth)**2 - - # gradients = np.sqrt(gradients) - # All channels considered together in this standard deviation beta /= 10 * data.std() if multichannel: From e8ddcefae351b7096890c608df559468e7ceb2f5 Mon Sep 17 00:00:00 2001 From: JDWarner Date: Fri, 31 Aug 2012 14:14:46 -0500 Subject: [PATCH 6/6] PEP8 compliance, removed `scaling`, different data parsing. This commit represents all recommended changes since the last commit, notably: * PEP8 compliance (in new sections; a few old ones still noncompliant w/indentations) * Moved `depth` kwarg to end of list and in docstring. Clarified `depth` docstring, and added section in Notes further explaining this parameter. * Added section in Notes warning that for multichannel inputs, all channels are combined during scaling. The user must separately normalize each channel prior to calling random_walker() * New method for parsing data, allowing more elegant gradient calculation code. Probably also more extensible. The 2D multispectral case forced this change. * New test: `test_multispectral_2d()` --- .../random_walker_segmentation.py | 97 ++++++++----------- .../segmentation/tests/test_random_walker.py | 23 ++++- 2 files changed, 61 insertions(+), 59 deletions(-) diff --git a/skimage/segmentation/random_walker_segmentation.py b/skimage/segmentation/random_walker_segmentation.py index 6165770d..52bf4425 100644 --- a/skimage/segmentation/random_walker_segmentation.py +++ b/skimage/segmentation/random_walker_segmentation.py @@ -67,23 +67,17 @@ def _make_graph_edges_3d(n_x, n_y, n_z): def _compute_weights_3d(data, beta=130, eps=1.e-6, depth=1., multichannel=False): # Weight calculation is main difference in multispectral version - # Original gradient**2 replaced with sqrt( sum of gradients**2 ) - if not multichannel: - gradients = _compute_gradients_3d( data, depth=depth )**2 - else: - for channel in range(data.shape[-1]): - if channel == 0: - gradients = _compute_gradients_3d(data[..., channel], - depth=depth)**2 - else: - gradients += _compute_gradients_3d(data[..., channel], - depth=depth)**2 + # 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 # All channels considered together in this standard deviation beta /= 10 * data.std() if multichannel: # New final term in beta to give == results in trivial case where # multiple identical spectra are passed. - beta /= np.sqrt( data.shape[-1] ) + beta /= np.sqrt(data.shape[-1]) gradients *= beta weights = np.exp(- gradients) weights += eps @@ -166,10 +160,7 @@ def _mask_edges_weights(edges, weights, mask): def _build_laplacian(data, mask=None, beta=50, depth=1., multichannel=False): - if not multichannel: - l_x, l_y, l_z = data.shape - else: - l_x, l_y, l_z = data.shape[0], data.shape[1], data.shape[2] + l_x, l_y, l_z = data.shape[: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, multichannel=multichannel) @@ -183,21 +174,21 @@ def _build_laplacian(data, mask=None, beta=50, depth=1., multichannel=False): #----------- Random walker algorithm -------------------------------- -def random_walker(data, labels, beta=130, depth=1., mode='bf', tol=1.e-3, - copy=True, multichannel=False, scaling='all', - return_full_prob=False): +def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, + multichannel=False, return_full_prob=False, depth=1.): """ - Multichannel random walker algorithm for segmentation from markers. + Random walker algorithm for segmentation from markers, for gray-level or + multichannel images. Parameters ---------- data : array_like - Image to be segmented in phases. Gray-level`data` can be two- or + Image to be segmented in phases. Gray-level `data` can be two- or three-dimensional; multichannel data can be three- or four- - dimensional (requires multichannel=True) with the highest - dimension denoting channels. Data spacing is assumed isotropic - unless depth keyword argument is used. + dimensional (multichannel=True) with the highest dimension denoting + channels. Data spacing is assumed isotropic unless depth keyword + argument is used. labels : array of ints, of same shape as `data` Array of seed markers labeled with different positive integers @@ -211,11 +202,6 @@ def random_walker(data, labels, beta=130, depth=1., mode='bf', tol=1.e-3, Penalization coefficient for the random walker motion (the greater `beta`, the more difficult the diffusion). - depth : float, default 1. - Correction for non-isotropic voxel depths in 3D volumes. - Default (1.) implies isotropy. This factor is derived as follows: - depth = (slice thickness) / (in-plane voxel dimension) - mode : {'bf', 'cg_mg', 'cg'} (default: 'bf') Mode for solving the linear system in the random walker algorithm. @@ -250,21 +236,17 @@ def random_walker(data, labels, beta=130, depth=1., mode='bf', tol=1.e-3, If True, input data is parsed as multichannel data (see 'data' above for proper input format in this case) - scaling : string, default 'all' - Controls input scaling if multichannel=True (otherwise no effect). - - - 'all' (default): Data from all channels is combined when scaling - input data to the range [0,1] as type np.float64. Recommended - option for RGB(A) inputs. - - - 'separate': Each channel is scaled individually, separate from the - others, to the range [0,1]. Select this if the channels are very - different, for example if one were x-ray CT and another MRI data. - 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. + 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. + Returns ------- @@ -286,6 +268,16 @@ def random_walker(data, labels, beta=130, depth=1., mode='bf', tol=1.e-3, Notes ----- + 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 algorithm was first proposed in *Random walks for image segmentation*, Leo Grady, IEEE Trans Pattern Anal Mach Intell. 2006 Nov;28(11):1768-83. @@ -346,21 +338,18 @@ def random_walker(data, labels, beta=130, depth=1., mode='bf', tol=1.e-3, """ # Parse input data if not multichannel: - # We work with 3-D arrays of floats + # We work with 4-D arrays of floats dims = data.shape - data = np.atleast_3d( img_as_float(data) ) + data = np.atleast_3d(img_as_float(data)) + data.shape += (1,) else: dims = data[..., 0].shape - data = np.atleast_3d( data ) # Should never be needed - if scaling.lower().strip() == 'all': - data = img_as_float( data ) - else: - newdata = np.zeros(data.shape, dtype=np.float64) - for channel in range( data.shape[-1] ): - newdata[..., channel] = img_as_float( data[..., channel] ) - del data - data = newdata - del newdata + assert multichannel and data.ndim > 2, 'For multichannel input, data \ + must have >= 3 dimensions.' + data = img_as_float(data) + if data.ndim == 3: + data.shape += (1,) + data = data.transpose((0, 1, 3, 2)) if copy: labels = np.copy(labels) @@ -409,7 +398,7 @@ def random_walker(data, labels, beta=130, depth=1., mode='bf', tol=1.e-3, if return_full_prob: labels = labels.astype(np.float) X = np.array([_clean_labels_ar(Xline, labels, - copy=True).reshape(dims) for Xline in X]) + copy=True).reshape(dims) for Xline in X]) for i in range(1, int(labels.max()) + 1): mask_i = np.squeeze(labels == i) X[i - 1, mask_i] = 1 @@ -429,7 +418,7 @@ def _solve_bf(lap_sparse, B, return_full_prob=False): lap_sparse = lap_sparse.tocsc() solver = sparse.linalg.factorized(lap_sparse.astype(np.double)) X = np.array([solver(np.array((-B[i]).todense()).ravel())\ - for i in range(len(B))]) + for i in range(len(B))]) if not return_full_prob: X = np.argmax(X, axis=0) return X diff --git a/skimage/segmentation/tests/test_random_walker.py b/skimage/segmentation/tests/test_random_walker.py index 74397dbe..ecf59e99 100644 --- a/skimage/segmentation/tests/test_random_walker.py +++ b/skimage/segmentation/tests/test_random_walker.py @@ -86,6 +86,7 @@ def test_2d_cg_mg(): full_prob[0, 25:45, 40:60]).all() return data, labels_cg_mg + def test_types(): lx = 70 ly = 100 @@ -96,6 +97,7 @@ def test_types(): assert (labels_cg_mg[25:45, 40:60] == 2).all() return data, labels_cg_mg + def test_reorder_labels(): lx = 70 ly = 100 @@ -106,7 +108,6 @@ def test_reorder_labels(): return data, labels_bf - def test_2d_inactive(): lx = 70 ly = 100 @@ -139,14 +140,26 @@ def test_3d_inactive(): return data, labels, old_labels, after_labels -def test_multispectral(): +def test_multispectral_2d(): + lx, ly = 70, 100 + data, labels = make_2d_syntheticdata(lx, ly) + data2 = data.copy() + data.shape += (1,) + data = data.repeat(2, axis=2) # Result should be identical + multi_labels = random_walker(data, labels, mode='cg', multichannel=True) + single_labels = random_walker(data2, labels, mode='cg') + assert (multi_labels.reshape(labels.shape)[25:45, 40:60] == 2).all() + return data, multi_labels, single_labels, labels + + +def test_multispectral_3d(): n = 30 lx, ly, lz = n, n, n - data, labels = make_3d_syntheticdata( lx, ly, lz ) + data, labels = make_3d_syntheticdata(lx, ly, lz) data.shape += (1,) - data = data.repeat(2, axis=3) # Result should be identical + data = data.repeat(2, axis=3) # Result should be identical multi_labels = random_walker(data, labels, mode='cg', multichannel=True) - single_labels = random_walker(data[:,:,:,0], labels, mode='cg') + single_labels = random_walker(data[..., 0], labels, mode='cg') assert (multi_labels.reshape(labels.shape)[13:17, 13:17, 13:17] == 2).all() assert (single_labels.reshape(labels.shape)[13:17, 13:17, 13:17] == 2).all() return data, multi_labels, single_labels, labels