diff --git a/TODO.txt b/TODO.txt index 7ad40875..502961c8 100644 --- a/TODO.txt +++ b/TODO.txt @@ -1,15 +1,15 @@ Version 0.10 ------------ -* Remove deprecated functions: - - ``skimage.filter.rank.*`` -* Remove deprecated parameter ``epsilon`` of ``skimage.viewer.LineProfile`` -* Remove backwards-compatability of ``skimage.measure.regionprops`` +* Remove deprecated functions in `skimage.filter.rank.*` +* Remove deprecated parameter `epsilon` of `skimage.viewer.LineProfile` +* Remove backwards-compatability of `skimage.measure.regionprops` +* Remove {`ratio`, `sigma`} deprecation warnings of `skimage.segmentation.slic` Version 0.9 ----------- * Remove deprecated functions - - ``skimage.filter.denoise_tv_chambolle`` - - ``skimage.morphology.is_local_maximum`` - - ``skimage.transform.hough`` - - ``skimage.transform.probabilistic_hough`` - - ``skimage.transform.hough_peaks`` + - `skimage.filter.denoise_tv_chambolle` + - `skimage.morphology.is_local_maximum` + - `skimage.transform.hough` + - `skimage.transform.probabilistic_hough` + - `skimage.transform.hough_peaks` diff --git a/doc/source/api_changes.txt b/doc/source/api_changes.txt index c363288d..bf9ddb60 100644 --- a/doc/source/api_changes.txt +++ b/doc/source/api_changes.txt @@ -1,6 +1,8 @@ Version 0.9 ----------- -- No longer wrap ``imread`` output in an ``Image`` class. +- No longer wrap ``imread`` output in an ``Image`` class +- Change default value of `sigma` parameter in ``skimage.segmentation.slic`` + to 0 Version 0.4 ----------- diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx index 5df747ad..c3d95ee0 100644 --- a/skimage/segmentation/_slic.pyx +++ b/skimage/segmentation/_slic.pyx @@ -2,93 +2,147 @@ #cython: boundscheck=False #cython: nonecheck=False #cython: wraparound=False -import numpy as np -from scipy import ndimage +from libc.float cimport DBL_MAX -from skimage.util import img_as_float, regular_grid -from skimage.color import rgb2lab, gray2rgb +import numpy as np +cimport numpy as cnp + +from skimage.util import regular_grid def _slic_cython(double[:, :, :, ::1] image_zyx, - Py_ssize_t[:, :, ::1] nearest_mean, - double[:, :, ::1] distance, - double[:, ::1] means, - Py_ssize_t max_iter, Py_ssize_t n_segments): + double[:, ::1] segments, + Py_ssize_t max_iter, + double[::1] spacing): """Helper function for SLIC segmentation. Parameters ---------- - image_zyx : 4D array of double, shape (Z, Y, X, 6) - The image with embedded coordinates, that is, `image_zyx[i, j, k]` is - `array([i, j, k, r, g, b])` or `array([i, j, k, L, a, b])`, depending - on the colorspace. - nearest_mean : 3D array of int, shape (Z, Y, X) - The (initially empty) label field. - distance : 3D array of double, shape (Z, Y, X) - The (initially infinity) array of distances to the nearest centroid. - means : 2D array of double, shape (n_segments, 6) - The centroids obtained by SLIC. + image_zyx : 4D array of double, shape (Z, Y, X, C) + The input image. + segments : 2D array of double, shape (N, 3 + C) + The initial centroids obtained by SLIC as [Z, Y, X, C...]. max_iter : int The maximum number of k-means iterations. - n_segments : int - The approximate/desired number of segments. + spacing : 1D array of double, shape (3,) + The voxel spacing along each image dimension. This parameter + controls the weights of the distances along z, y, and x during + k-means clustering. Returns ------- - nearest_mean : 3D array of int, shape (Z, Y, X) + nearest_segments : 3D array of int, shape (Z, Y, X) The label field/superpixels found by SLIC. + + Notes + ----- + The image is considered to be in (z, y, x) order, which can be + surprising. More commonly, the order (x, y, z) is used. However, + in 3D image analysis, 'z' is usually the "special" dimension, with, + for example, a different effective resolution than the other two + axes. Therefore, x and y are often processed together, or viewed as + a cut-plane through the volume. So, if the order was (x, y, z) and + we wanted to look at the 5th cut plane, we would write:: + + my_z_plane = img3d[:, :, 5] + + but, assuming a C-contiguous array, this would grab a discontiguous + slice of memory, which is bad for performance. In contrast, if we + see the image as (z, y, x) ordered, we would do:: + + my_z_plane = img3d[5] + + and get back a contiguous block of memory. This is better both for + performance and for readability. """ - # initialize on grid: + # initialize on grid cdef Py_ssize_t depth, height, width - depth, height, width = (image_zyx.shape[0], image_zyx.shape[1], - image_zyx.shape[2]) + depth = image_zyx.shape[0] + height = image_zyx.shape[1] + width = image_zyx.shape[2] + + cdef Py_ssize_t n_segments = segments.shape[0] + # number of features [X, Y, Z, ...] + cdef Py_ssize_t n_features = segments.shape[1] + # approximate grid size for desired n_segments cdef Py_ssize_t step_z, step_y, step_x slices = regular_grid((depth, height, width), n_segments) step_z, step_y, step_x = [int(s.step) for s in slices] - n_means = means.shape[0] - cdef Py_ssize_t i, k, x, y, z, x_min, x_max, y_min, y_max, z_min, z_max, \ - changes - cdef double dist_mean - cdef double tmp + cdef Py_ssize_t[:, :, ::1] nearest_segments \ + = np.empty((depth, height, width), dtype=np.intp) + cdef double[:, :, ::1] distance \ + = np.empty((depth, height, width), dtype=np.double) + cdef Py_ssize_t[::1] n_segment_elems = np.zeros(n_segments, dtype=np.intp) + + cdef Py_ssize_t i, c, k, x, y, z, x_min, x_max, y_min, y_max, z_min, z_max + cdef char change + cdef double dist_center, cx, cy, cz, dy, dz + + cdef double sz, sy, sx + sz = spacing[0] + sy = spacing[1] + sx = spacing[2] + for i in range(max_iter): - changes = 0 - distance[:, :, :] = np.inf - # assign pixels to means - for k in range(n_means): - # compute windows: - z_min = int(max(means[k, 0] - 2 * step_z, 0)) - z_max = int(min(means[k, 0] + 2 * step_z, depth)) - y_min = int(max(means[k, 1] - 2 * step_y, 0)) - y_max = int(min(means[k, 1] + 2 * step_y, height)) - x_min = int(max(means[k, 2] - 2 * step_x, 0)) - x_max = int(min(means[k, 2] + 2 * step_x, width)) + change = 0 + distance[:, :, :] = DBL_MAX + + # assign pixels to segments + for k in range(n_segments): + + # segment coordinate centers + cz = segments[k, 0] + cy = segments[k, 1] + cx = segments[k, 2] + + # compute windows + z_min = max(cz - 2 * step_z, 0) + z_max = min(cz + 2 * step_z + 1, depth) + y_min = max(cy - 2 * step_y, 0) + y_max = min(cy + 2 * step_y + 1, height) + x_min = max(cx - 2 * step_x, 0) + x_max = min(cx + 2 * step_x + 1, width) + for z in range(z_min, z_max): + dz = (sz * (cz - z)) ** 2 for y in range(y_min, y_max): + dy = (sy * (cy - y)) ** 2 for x in range(x_min, x_max): - dist_mean = 0 - for c in range(6): - # you would think the compiler can optimize the - # squaring itself. mine can't (with O2) - tmp = image_zyx[z, y, x, c] - means[k, c] - dist_mean += tmp * tmp - if distance[z, y, x] > dist_mean: - nearest_mean[z, y, x] = k - distance[z, y, x] = dist_mean - changes = 1 - if changes == 0: + dist_center = dz + dy + (sx * (cx - x)) ** 2 + for c in range(3, n_features): + dist_center += (image_zyx[z, y, x, c - 3] + - segments[k, c]) ** 2 + if distance[z, y, x] > dist_center: + nearest_segments[z, y, x] = k + distance[z, y, x] = dist_center + change = 1 + + # stop if no pixel changed its segment + if change == 0: break - # recompute means: - nearest_mean_ravel = np.asarray(nearest_mean).ravel() - means_list = [] - for j in range(6): - image_zyx_ravel = ( - np.ascontiguousarray(image_zyx[:, :, :, j]).ravel()) - means_list.append(np.bincount(nearest_mean_ravel, - image_zyx_ravel)) - in_mean = np.bincount(nearest_mean_ravel) - in_mean[in_mean == 0] = 1 - means = (np.vstack(means_list) / in_mean).T.copy("C") - return np.ascontiguousarray(nearest_mean) + + # recompute segment centers + + # sum features for all segments + n_segment_elems[:] = 0 + segments[:, :] = 0 + for z in range(depth): + for y in range(height): + for x in range(width): + k = nearest_segments[z, y, x] + n_segment_elems[k] += 1 + segments[k, 0] += z + segments[k, 1] += y + segments[k, 2] += x + for c in range(3, n_features): + segments[k, c] += image_zyx[z, y, x, c - 3] + + # divide by number of elements per segment to obtain mean + for k in range(n_segments): + for c in range(n_features): + segments[k, c] /= n_segment_elems[k] + + return np.asarray(nearest_segments) diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py index b905f199..8276b1a8 100644 --- a/skimage/segmentation/slic_superpixels.py +++ b/skimage/segmentation/slic_superpixels.py @@ -5,46 +5,52 @@ import numpy as np from scipy import ndimage import warnings -from ..util import img_as_float, regular_grid -from ..color import rgb2lab, gray2rgb, guess_spatial_dimensions -from ._slic import _slic_cython +from skimage.util import img_as_float, regular_grid +from skimage.segmentation._slic import _slic_cython +from skimage.color import rgb2lab -def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=1, - multichannel=None, convert2lab=True, ratio=None): - """Segments image using k-means clustering in Color-(x,y) space. +def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None, + spacing=None, multichannel=True, convert2lab=True, ratio=None): + """Segments image using k-means clustering in Color-(x,y,z) space. Parameters ---------- - image : (width, height [, depth] [, 3]) ndarray - Input image, which can be 2D or 3D, and grayscale or multi-channel + image : 2D, 3D or 4D ndarray + Input image, which can be 2D or 3D, and grayscale or multichannel (see `multichannel` parameter). - n_segments : int, optional (default: 100) + n_segments : int, optional The (approximate) number of labels in the segmented output image. - compactness: float, optional (default: 10) + compactness : float, optional Balances color-space proximity and image-space proximity. Higher values give more weight to image-space. As `compactness` tends to infinity, superpixel shapes become square/cubic. - max_iter : int, optional (default: 10) + max_iter : int, optional Maximum number of iterations of k-means. - sigma : float, optional (default: 1) - Width of Gaussian smoothing kernel for preprocessing. Zero means no - smoothing. - multichannel : bool, optional (default: None) + sigma : float or (3,) array-like of floats, optional + Width of Gaussian smoothing kernel for pre-processing for each + dimension of the image. The same sigma is applied to each dimension in + case of a scalar value. Zero means no smoothing. + Note, that `sigma` is automatically scaled if it is scalar and a + manual voxel spacing is provided (see Notes section). + spacing : (3,) array-like of floats, optional + The voxel spacing along each image dimension. By default, `slic` + assumes uniform spacing (same voxel resolution along z, y and x). + This parameter controls the weights of the distances along z, y, + and x during k-means clustering. + multichannel : bool, optional Whether the last axis of the image is to be interpreted as multiple - channels. Only 3 channels are supported. If `None`, the function will - attempt to guess this, and raise a warning if ambiguous, when the - array has shape (M, N, 3). - convert2lab : bool, optional (default: True) + channels or another spatial dimension. + convert2lab : bool, optional Whether the input should be converted to Lab colorspace prior to - segmentation. For this purpose, the input is assumed to be RGB. Highly + segmentation. For this purpose, the input is assumed to be RGB. Highly recommended. ratio : float, optional Synonym for `compactness`. This keyword is deprecated. Returns ------- - segment_mask : (width, height) ndarray + labels : 2D or 3D array Integer mask indicating segment labels. Raises @@ -52,20 +58,23 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=1, ValueError If: - the image dimension is not 2 or 3 and `multichannel == False`, OR - - the image dimension is not 3 or 4 and `multichannel == True`, OR - - `multichannel == True` and the length of the last dimension of - the image is not 3, OR + - the image dimension is not 3 or 4 and `multichannel == True` Notes ----- - If `sigma > 0` as is default, the image is smoothed using a Gaussian kernel - prior to segmentation. + * If `sigma > 0`, the image is smoothed using a Gaussian kernel prior to + segmentation. - The image is rescaled to be in [0, 1] prior to processing. + * If `sigma` is scalar and `spacing` is provided, the kernel width is + divided along each dimension by the spacing. For example, if ``sigma=1`` + and ``spacing=[5, 1, 1]``, the effective `sigma` is ``[0.2, 1, 1]``. This + ensures sensible smoothing for anisotropic images. - Images of shape (M, N, 3) are interpreted as 2D RGB images by default. To - interpret them as 3D with the last dimension having length 3, use - `multichannel=False`. + * The image is rescaled to be in [0, 1] prior to processing. + + * Images of shape (M, N, 3) are interpreted as 2D RGB images by default. To + interpret them as 3D with the last dimension having length 3, use + `multichannel=False`. References ---------- @@ -78,66 +87,80 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=1, >>> from skimage.segmentation import slic >>> from skimage.data import lena >>> img = lena() - >>> segments = slic(img, n_segments=100, ratio=10) - >>> # Increasing the ratio parameter yields more square regions - >>> segments = slic(img, n_segments=100, ratio=20) + >>> segments = slic(img, n_segments=100, compactness=10) + >>> # Increasing the compactness parameter yields more square regions + >>> segments = slic(img, n_segments=100, compactness=20) """ + + if sigma is None: + warnings.warn('Default value of keyword `sigma` changed from ``1`` ' + 'to ``0``.') + sigma = 0 if ratio is not None: - msg = 'Keyword `ratio` is deprecated. Use `compactness` instead.' - warnings.warn(msg) + warnings.warn('Keyword `ratio` is deprecated. Use `compactness` ' + 'instead.') compactness = ratio - spatial_dims = guess_spatial_dimensions(image) - if spatial_dims is None and multichannel is None: - msg = ("Images with dimensions (M, N, 3) are interpreted as 2D+RGB" + - " by default. Use `multichannel=False` to interpret as " + - " 3D image with last dimension of length 3.") - warnings.warn(RuntimeWarning(msg)) - multichannel = True - elif multichannel is None: - multichannel = (spatial_dims + 1 == image.ndim) - if ((not multichannel and image.ndim not in [2, 3]) or - (multichannel and image.ndim not in [3, 4]) or - (multichannel and image.shape[-1] != 3)): - ValueError("Only 1- or 3-channel 2- or 3-D images are supported.") + image = img_as_float(image) - if not multichannel: - image = gray2rgb(image) - if image.ndim == 3: - # See 2D RGB image as 3D RGB image with Z = 1 + is_2d = False + if image.ndim == 2: + # 2D grayscale image + image = image[np.newaxis, ..., np.newaxis] + is_2d = True + elif image.ndim == 3 and multichannel: + # Make 2D multichannel image 3D with depth = 1 image = image[np.newaxis, ...] + is_2d = True + elif image.ndim == 3 and not multichannel: + # Add channel as single last dimension + image = image[..., np.newaxis] + + if spacing is None: + spacing = np.ones(3) + elif isinstance(spacing, (list, tuple)): + spacing = np.array(spacing, dtype=np.double) + if not isinstance(sigma, coll.Iterable): - sigma = np.array([sigma, sigma, sigma, 0]) + sigma = np.array([sigma, sigma, sigma], dtype=np.double) + sigma /= spacing.astype(np.double) + elif isinstance(sigma, (list, tuple)): + sigma = np.array(sigma, dtype=np.double) if (sigma > 0).any(): + # add zero smoothing for multichannel dimension + sigma = list(sigma) + [0] image = ndimage.gaussian_filter(image, sigma) - if convert2lab: + + if convert2lab and multichannel: + if image.shape[3] != 3: + raise ValueError("Lab colorspace conversion requires a RGB image.") image = rgb2lab(image) - # initialize on grid: depth, height, width = image.shape[:3] - # approximate grid size for desired n_segments + + # initialize cluster centroids for desired number of segments grid_z, grid_y, grid_x = np.mgrid[:depth, :height, :width] slices = regular_grid(image.shape[:3], n_segments) step_z, step_y, step_x = [int(s.step) for s in slices] - means_z = grid_z[slices] - means_y = grid_y[slices] - means_x = grid_x[slices] + segments_z = grid_z[slices] + segments_y = grid_y[slices] + segments_x = grid_x[slices] + + segments_color = np.zeros(segments_z.shape + (image.shape[3],)) + segments = np.concatenate([segments_z[..., np.newaxis], + segments_y[..., np.newaxis], + segments_x[..., np.newaxis], + segments_color + ], axis=-1).reshape(-1, 3 + image.shape[3]) + segments = np.ascontiguousarray(segments) - means_color = np.zeros(means_z.shape + (3,)) - means = np.concatenate([means_z[..., np.newaxis], means_y[..., np.newaxis], - means_x[..., np.newaxis], means_color - ], axis=-1).reshape(-1, 6) - means = np.ascontiguousarray(means) # we do the scaling of ratio in the same way as in the SLIC paper # so the values have the same meaning ratio = float(max((step_z, step_y, step_x))) / compactness - image_zyx = np.concatenate([grid_z[..., np.newaxis], - grid_y[..., np.newaxis], - grid_x[..., np.newaxis], - image * ratio], axis=-1).copy("C") - nearest_mean = np.zeros((depth, height, width), dtype=np.intp) - distance = np.empty((depth, height, width), dtype=np.float) - segment_map = _slic_cython(image_zyx, nearest_mean, distance, means, - max_iter, n_segments) - if segment_map.shape[0] == 1: - segment_map = segment_map[0] - return segment_map + image = np.ascontiguousarray(image * ratio) + + labels = _slic_cython(image, segments, max_iter, spacing) + + if is_2d: + labels = labels[0] + + return labels diff --git a/skimage/segmentation/tests/test_slic.py b/skimage/segmentation/tests/test_slic.py index 9e6a39e2..a4657785 100644 --- a/skimage/segmentation/tests/test_slic.py +++ b/skimage/segmentation/tests/test_slic.py @@ -20,6 +20,7 @@ def test_color_2d(): # we expect 4 segments assert_equal(len(np.unique(seg)), 4) + assert_equal(seg.shape, img.shape[:-1]) assert_array_equal(seg[:10, :10], 0) assert_array_equal(seg[10:, :10], 2) assert_array_equal(seg[:10, 10:], 1) @@ -35,10 +36,11 @@ def test_gray_2d(): img += 0.0033 * rnd.normal(size=img.shape) img[img > 1] = 1 img[img < 0] = 0 - seg = slic(img, sigma=0, n_segments=4, compactness=20.0, - multichannel=False) + seg = slic(img, sigma=0, n_segments=4, compactness=1, + multichannel=False, convert2lab=False) assert_equal(len(np.unique(seg)), 4) + assert_equal(seg.shape, img.shape) assert_array_equal(seg[:10, :10], 0) assert_array_equal(seg[10:, :10], 2) assert_array_equal(seg[:10, 10:], 1) @@ -80,14 +82,43 @@ def test_gray_3d(): img += 0.001 * rnd.normal(size=img.shape) img[img > 1] = 1 img[img < 0] = 0 - seg = slic(img, sigma=0, n_segments=8, compactness=20.0, - multichannel=False) + seg = slic(img, sigma=0, n_segments=8, compactness=1, + multichannel=False, convert2lab=False) assert_equal(len(np.unique(seg)), 8) for s, c in zip(slices, range(8)): assert_array_equal(seg[s], c) +def test_list_sigma(): + rnd = np.random.RandomState(0) + img = np.array([[1, 1, 1, 0, 0, 0], + [0, 0, 0, 1, 1, 1]], np.float) + img += 0.1 * rnd.normal(size=img.shape) + result_sigma = np.array([[0, 0, 0, 1, 1, 1], + [0, 0, 0, 1, 1, 1]], np.int) + seg_sigma = slic(img, n_segments=2, sigma=[1, 50, 1], multichannel=False) + assert_equal(seg_sigma, result_sigma) + + +def test_spacing(): + rnd = np.random.RandomState(0) + img = np.array([[1, 1, 1, 0, 0], + [1, 1, 0, 0, 0]], np.float) + result_non_spaced = np.array([[0, 0, 0, 1, 1], + [0, 0, 1, 1, 1]], np.int) + result_spaced = np.array([[0, 0, 0, 0, 0], + [1, 1, 1, 1, 1]], np.int) + img += 0.1 * rnd.normal(size=img.shape) + seg_non_spaced = slic(img, n_segments=2, sigma=0, multichannel=False, + compactness=1.0) + seg_spaced = slic(img, n_segments=2, sigma=0, spacing=[1, 500, 1], + compactness=1.0, multichannel=False) + assert_equal(seg_non_spaced, result_non_spaced) + assert_equal(seg_spaced, result_spaced) + + + if __name__ == '__main__': from numpy import testing testing.run_module_suite()