From 7161175e43ed16dbfcf5f9b6e547fcb58f4079b8 Mon Sep 17 00:00:00 2001 From: "Josh Warner (Mac)" Date: Fri, 18 Apr 2014 16:52:19 -0500 Subject: [PATCH 1/5] DOC: Clarify UMFPACK warning; use if statements instead of assert --- .../random_walker_segmentation.py | 24 ++++++++++--------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/skimage/segmentation/random_walker_segmentation.py b/skimage/segmentation/random_walker_segmentation.py index ffd7526d..8665dd37 100644 --- a/skimage/segmentation/random_walker_segmentation.py +++ b/skimage/segmentation/random_walker_segmentation.py @@ -202,7 +202,7 @@ def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, Image to be segmented in phases. Gray-level `data` can be two- or three-dimensional; multichannel data can be three- or four- dimensional (multichannel=True) with the highest dimension denoting - channels. Data spacing is assumed isotropic unless the `spacing` + channels. Data spacing is assumed isotropic unless the `spacing` keyword argument is used. labels : array of ints, of same shape as `data` without channels dimension Array of seed markers labeled with different positive integers @@ -344,11 +344,12 @@ def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, mode = 'bf' if UmfpackContext is None and mode == 'cg': - warnings.warn('SciPy was built without UMFPACK. Consider rebuilding ' - 'SciPy with UMFPACK, this will greatly speed up the ' - 'random walker functions. You may also install pyamg ' - 'and run the random walker function in cg_mg mode ' - '(see the docstrings)') + warnings.warn('"cg" mode will be used, but it may be slower than ' + '"bf" because SciPy was built without UMFPACK. Consider' + ' rebuilding SciPy with UMFPACK; this will greatly ' + 'accelerate the conjugate gradient ("cg") solver. ' + 'You may also install pyamg and run the random_walker ' + 'function in "cg_mg" mode (see docstring).') # Spacing kwarg checks if spacing is None: @@ -362,15 +363,16 @@ def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, # Parse input data if not multichannel: # We work with 4-D arrays of floats - assert data.ndim > 1 and data.ndim < 4, 'For non-multichannel input, \ - data must be of dimension 2 \ - or 3.' + if data.ndim < 1 or data.ndim > 4: + raise ValueError('For non-multichannel input, data must be of ' + 'dimension 2 or 3.') dims = data.shape data = np.atleast_3d(img_as_float(data))[..., np.newaxis] else: + if data.ndim < 2: + raise ValueError('For multichannel input, data must have >= 3 ' + 'dimensions.') dims = data[..., 0].shape - assert multichannel and data.ndim > 2, 'For multichannel input, data \ - must have >= 3 dimensions.' data = img_as_float(data) if data.ndim == 3: data = data[..., np.newaxis].transpose((0, 1, 3, 2)) From 7847287f4619e73896ebdfe039db89e8cc2c6e4d Mon Sep 17 00:00:00 2001 From: "Josh Warner (Mac)" Date: Fri, 18 Apr 2014 18:47:20 -0500 Subject: [PATCH 2/5] 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() From 88f08ebea2799bcad179074bcd01cb823c662a91 Mon Sep 17 00:00:00 2001 From: "Josh Warner (Mac)" Date: Fri, 18 Apr 2014 21:48:18 -0500 Subject: [PATCH 3/5] TST: Add unit test for length-2 spacing with rank-2 input --- 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 2fa0f121..a7318daa 100644 --- a/skimage/segmentation/tests/test_random_walker.py +++ b/skimage/segmentation/tests/test_random_walker.py @@ -270,6 +270,17 @@ def test_trivial_cases(): np.testing.assert_array_equal(test, expected) +def test_length2_spacing(): + # If this passes without raising an exception (warnings OK), the new + # spacing code is working properly. + np.random.seed(42) + img = np.ones((10, 10)) + 0.2 * np.random.normal(size=(10, 10)) + labels = np.zeros((10, 10), dtype=np.uint8) + labels[2, 4] = 1 + labels[6, 8] = 4 + random_walker(img, labels, spacing=(1., 2.)) + + def test_bad_inputs(): # Too few dimensions img = np.ones(10) @@ -279,6 +290,7 @@ def test_bad_inputs(): random_walker, img, labels, multichannel=True) # Too many dimensions + np.random.seed(42) 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) From 0ac86b3e96b126da9dab5c2fa33b5806937ca1a9 Mon Sep 17 00:00:00 2001 From: "Josh Warner (Mac)" Date: Tue, 6 May 2014 02:16:45 -0500 Subject: [PATCH 4/5] STY: More elegant and maintainable style per code review --- .../random_walker_segmentation.py | 33 ++++++++++++------- 1 file changed, 21 insertions(+), 12 deletions(-) diff --git a/skimage/segmentation/random_walker_segmentation.py b/skimage/segmentation/random_walker_segmentation.py index 8debd41d..c66e9c6f 100644 --- a/skimage/segmentation/random_walker_segmentation.py +++ b/skimage/segmentation/random_walker_segmentation.py @@ -22,6 +22,7 @@ from scipy import sparse, ndimage try: from scipy.sparse.linalg.dsolve import umfpack old_del = umfpack.UmfpackContext.__del__ + def new_del(self): try: old_del(self) @@ -351,7 +352,7 @@ 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).') - if (labels == 0).sum() == 0: + if (labels != 0).all(): warnings.warn('Random walker only segments unlabeled areas, where ' 'labels == 0. No zero valued areas in labels were ' 'found. Returning provided labels.') @@ -361,16 +362,22 @@ def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, 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) + out_labels = np.empty(labels.shape + (len(unique_labels),), + dtype=np.bool) + for n, i in enumerate(unique_labels): + out_labels[..., n] = (labels == i) + else: out_labels = labels return out_labels + # This algorithm expects 4-D arrays of floats, where the first three + # dimensions are spatial and the final denotes channels. 2-D images have + # a singleton placeholder dimension added for the third spatial dimension, + # and single channel images likewise have a singleton added for channels. + # The following block ensures valid input and coerces it to the correct + # form. if not multichannel: - # We work with 4-D arrays of floats if data.ndim < 2 or data.ndim > 3: raise ValueError('For non-multichannel input, data must be of ' 'dimension 2 or 3.') @@ -382,22 +389,24 @@ def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, 'dimensions.') 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)) + if data.ndim == 3: # 2D multispectral, needs singleton in 3rd axis + data = data[:, :, np.newaxis, :].transpose((0, 1, 3, 2)) # Spacing kwarg checks if spacing is None: - spacing = (1., ) * 3 + spacing = np.asarray((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., ) + if len(spacing) == 2: # Need a dummy spacing for singleton 3rd dim + spacing = np.r_[spacing, 1.] + else: # Convert to array + spacing = np.asarray(spacing) else: raise ValueError('Input argument `spacing` incorrect, should be an ' - 'iterable of length 3.') + 'iterable with one number per spatial dimension.') if copy: labels = np.copy(labels) From 0d775c3aab3bd20210b546cc5988832f6f203a65 Mon Sep 17 00:00:00 2001 From: "Josh Warner (Mac)" Date: Wed, 7 May 2014 02:42:38 -0500 Subject: [PATCH 5/5] Remove check for non-numeric values in spacing. --- skimage/segmentation/random_walker_segmentation.py | 5 ----- skimage/segmentation/tests/test_random_walker.py | 4 ---- 2 files changed, 9 deletions(-) diff --git a/skimage/segmentation/random_walker_segmentation.py b/skimage/segmentation/random_walker_segmentation.py index c66e9c6f..d722569b 100644 --- a/skimage/segmentation/random_walker_segmentation.py +++ b/skimage/segmentation/random_walker_segmentation.py @@ -9,7 +9,6 @@ significantly the performance. """ import warnings -from numbers import Number import numpy as np from scipy import sparse, ndimage @@ -396,10 +395,6 @@ def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True, if spacing is None: spacing = np.asarray((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: # Need a dummy spacing for singleton 3rd dim spacing = np.r_[spacing, 1.] else: # Convert to array diff --git a/skimage/segmentation/tests/test_random_walker.py b/skimage/segmentation/tests/test_random_walker.py index a7318daa..9d30473b 100644 --- a/skimage/segmentation/tests/test_random_walker.py +++ b/skimage/segmentation/tests/test_random_walker.py @@ -305,10 +305,6 @@ def test_bad_inputs(): 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__': np.testing.run_module_suite()