From 7847287f4619e73896ebdfe039db89e8cc2c6e4d Mon Sep 17 00:00:00 2001 From: "Josh Warner (Mac)" Date: Fri, 18 Apr 2014 18:47:20 -0500 Subject: [PATCH] FIX: Shortcut output or catch trivial cases in random_walker Also: * New tests to cover these new checks * Improvements to docstrings and user warnings * Generalize handling of `sampling` in accordance with docstring * Some extra whitespace to improve readability --- .../random_walker_segmentation.py | 76 +++++++++++++------ .../segmentation/tests/test_random_walker.py | 46 ++++++++++- 2 files changed, 95 insertions(+), 27 deletions(-) diff --git a/skimage/segmentation/random_walker_segmentation.py b/skimage/segmentation/random_walker_segmentation.py index 8665dd37..8debd41d 100644 --- a/skimage/segmentation/random_walker_segmentation.py +++ b/skimage/segmentation/random_walker_segmentation.py @@ -9,7 +9,7 @@ significantly the performance. """ import warnings - +from numbers import Number import numpy as np from scipy import sparse, ndimage @@ -216,14 +216,14 @@ def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, beta : float Penalization coefficient for the random walker motion (the greater `beta`, the more difficult the diffusion). - mode : {'bf', 'cg_mg', 'cg'} (default: 'bf') - Mode for solving the linear system in the random walker - algorithm. + mode : string, available options {'cg_mg', 'cg', 'bf'} + Mode for solving the linear system in the random walker algorithm. + If no preference given, automatically attempt to use the fastest + option available ('cg_mg' from pyamg >> 'cg' with UMFPACK > 'bf'). - - 'bf' (brute force, default): an LU factorization of the Laplacian is + - 'bf' (brute force): an LU factorization of the Laplacian is computed. This is fast for small images (<1024x1024), but very slow - (due to the memory cost) and memory-consuming for big images (in 3-D - for example). + and memory-intensive for large images (e.g., 3-D volumes). - 'cg' (conjugate gradient): the linear system is solved iteratively using the Conjugate Gradient method from scipy.sparse.linalg. This is less memory-consuming than the brute force method for large images, @@ -316,11 +316,11 @@ def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, Examples -------- - >>> a = np.zeros((10, 10)) + 0.2*np.random.random((10, 10)) + >>> a = np.zeros((10, 10)) + 0.2 * np.random.random((10, 10)) >>> a[5:8, 5:8] += 1 >>> b = np.zeros_like(a) - >>> b[3,3] = 1 #Marker for first phase - >>> b[6,6] = 2 #Marker for second phase + >>> b[3, 3] = 1 # Marker for first phase + >>> b[6, 6] = 2 # Marker for second phase >>> random_walker(a, b) array([[1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], @@ -334,7 +334,7 @@ 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]], dtype=int32) """ - + # Parse input data if mode is None: if amg_loaded: mode = 'cg_mg' @@ -351,40 +351,64 @@ def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, 'You may also install pyamg and run the random_walker ' 'function in "cg_mg" mode (see docstring).') - # Spacing kwarg checks - if spacing is None: - spacing = (1., 1., 1.) - elif len(spacing) == 3: - pass - else: - raise ValueError('Input argument `spacing` incorrect, should be an ' - 'iterable of length 3.') + if (labels == 0).sum() == 0: + warnings.warn('Random walker only segments unlabeled areas, where ' + 'labels == 0. No zero valued areas in labels were ' + 'found. Returning provided labels.') + + if return_full_prob: + # Find and iterate over valid labels + unique_labels = np.unique(labels) + unique_labels = unique_labels[unique_labels > 0] + + out_labels = np.empty(labels.shape + (0, )) + for i in unique_labels: + out_labels = np.concatenate( + (out_labels, (labels == i)[..., np.newaxis]), axis=-1) + else: + out_labels = labels + return out_labels - # Parse input data if not multichannel: # We work with 4-D arrays of floats - if data.ndim < 1 or data.ndim > 4: + if data.ndim < 2 or data.ndim > 3: raise ValueError('For non-multichannel input, data must be of ' 'dimension 2 or 3.') - dims = data.shape + dims = data.shape # To reshape final labeled result data = np.atleast_3d(img_as_float(data))[..., np.newaxis] else: - if data.ndim < 2: - raise ValueError('For multichannel input, data must have >= 3 ' + if data.ndim < 3: + raise ValueError('For multichannel input, data must have 3 or 4 ' 'dimensions.') - dims = data[..., 0].shape + dims = data[..., 0].shape # To reshape final labeled result data = img_as_float(data) if data.ndim == 3: data = data[..., np.newaxis].transpose((0, 1, 3, 2)) + # Spacing kwarg checks + if spacing is None: + spacing = (1., ) * 3 + elif len(spacing) == len(dims): + for i in spacing: + if not isinstance(i, Number): + raise ValueError('Input `spacing` contained %s, which is not ' + 'a number.' % (i)) + if len(spacing) == 2: + spacing += (1., ) + else: + raise ValueError('Input argument `spacing` incorrect, should be an ' + 'iterable of length 3.') + if copy: labels = np.copy(labels) label_values = np.unique(labels) + # Reorder label values to have consecutive integers (no gaps) if np.any(np.diff(label_values) != 1): mask = labels >= 0 labels[mask] = rank_order(labels[mask])[0].astype(labels.dtype) labels = labels.astype(np.int32) + # If the array has pruned zones, be sure that no isolated pixels # exist between pruned zones (they could not be determined) if np.any(labels < 0): @@ -399,6 +423,7 @@ def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, lap_sparse = _build_laplacian(data, spacing, beta=beta, multichannel=multichannel) lap_sparse, B = _buildAB(lap_sparse, labels) + # We solve the linear system # lap_sparse X = B # where X[i, j] is the probability that a marker of label i arrives @@ -420,6 +445,7 @@ def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, if mode == 'bf': X = _solve_bf(lap_sparse, B, return_full_prob=return_full_prob) + # Clean up results if return_full_prob: labels = labels.astype(np.float) diff --git a/skimage/segmentation/tests/test_random_walker.py b/skimage/segmentation/tests/test_random_walker.py index 54e0d832..2fa0f121 100644 --- a/skimage/segmentation/tests/test_random_walker.py +++ b/skimage/segmentation/tests/test_random_walker.py @@ -255,6 +255,48 @@ def test_spacing_1(): assert (labels_aniso2[26:34, 13:17, 13:17] == 2).all() +def test_trivial_cases(): + # When all voxels are labeled + img = np.ones((10, 10)) + labels = np.ones((10, 10)) + pass_through = random_walker(img, labels) + np.testing.assert_array_equal(pass_through, labels) + + # When all voxels are labeled AND return_full_prob is True + labels[:, :5] = 3 + expected = np.concatenate(((labels == 1)[..., np.newaxis], + (labels == 3)[..., np.newaxis]), axis=2) + test = random_walker(img, labels, return_full_prob=True) + np.testing.assert_array_equal(test, expected) + + +def test_bad_inputs(): + # Too few dimensions + img = np.ones(10) + labels = np.arange(10) + np.testing.assert_raises(ValueError, random_walker, img, labels) + np.testing.assert_raises(ValueError, + random_walker, img, labels, multichannel=True) + + # Too many dimensions + img = np.random.normal(size=(3, 3, 3, 3, 3)) + labels = np.arange(3 ** 5).reshape(img.shape) + np.testing.assert_raises(ValueError, random_walker, img, labels) + np.testing.assert_raises(ValueError, + random_walker, img, labels, multichannel=True) + + # Spacing incorrect length + img = np.random.normal(size=(10, 10)) + labels = np.zeros((10, 10)) + labels[2, 4] = 2 + labels[6, 8] = 5 + np.testing.assert_raises(ValueError, + random_walker, img, labels, spacing=(1,)) + + # Spacing contains unacceptable information + np.testing.assert_raises( + ValueError, random_walker, img, labels, spacing=(1, 'chickens')) + + if __name__ == '__main__': - from numpy import testing - testing.run_module_suite() + np.testing.run_module_suite()