Merge pull request #6 from jni/slic-speed-spacing

Slic speed spacing

skimage/segmentation/_slic.pyx

@@ -12,7 +12,8 @@ from skimage.util import regular_grid
 
 def _slic_cython(double[:, :, :, ::1] image_zyx,
                  double[:, ::1] segments,
-                 Py_ssize_t max_iter):
+                 Py_ssize_t max_iter,
+                 double[:] spacing):
     """Helper function for SLIC segmentation.
 
     Parameters
@@ -23,6 +24,10 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
         The initial centroids obtained by SLIC as [Z, Y, X, C...].
     max_iter : int
         The maximum number of k-means iterations.
+    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
     -------
@@ -55,6 +60,11 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
     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):
         change = 0
         distance[:, :, :] = DBL_MAX
@@ -76,11 +86,11 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
             x_max = <Py_ssize_t>min(cx + 2 * step_x + 1, width)
 
             for z in range(z_min, z_max):
-                dz = (cz - z) ** 2
+                dz = (sz * (cz - z)) ** 2
                 for y in range(y_min, y_max):
-                    dy = (cy - y) ** 2
+                    dy = (sy * (cy - y)) ** 2
                     for x in range(x_min, x_max):
-                        dist_center = dz + dy + (cx - x) ** 2
+                        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

skimage/segmentation/slic_superpixels.py

@@ -11,7 +11,7 @@ from skimage.color import rgb2lab
 
 
 def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
-         multichannel=True, convert2lab=True, ratio=None):
+         spacing=None, multichannel=True, convert2lab=True, ratio=None):
     """Segments image using k-means clustering in Color-(x,y,z) space.
 
     Parameters
@@ -27,10 +27,15 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
         infinity, superpixel shapes become square/cubic.
     max_iter : int, optional
         Maximum number of iterations of k-means.
-    sigma : float or (3,) array of floats, optional
+    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.
+    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 or another spatial dimension.
@@ -58,6 +63,11 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
     If `sigma > 0`, the image is smoothed using a Gaussian kernel prior to
     segmentation.
 
+    If `sigma > 0` 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.
+
     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
@@ -75,9 +85,9 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
     >>> 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:
@@ -103,15 +113,21 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
         # 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, np.double)
     if not isinstance(sigma, coll.Iterable):
-        sigma = np.array([sigma, sigma, sigma])
+        sigma = np.array([sigma, sigma, sigma], np.double)
+    elif isinstance(sigma, (list, tuple)):
+        sigma = np.array(sigma, np.double)
     if (sigma > 0).any():
+        sigma /= spacing.astype(np.double)
         sigma = list(sigma) + [0]
         image = ndimage.gaussian_filter(image, sigma)
 
-    if convert2lab:
-
-        if not multichannel or image.shape[3] != 3:
+    if convert2lab and multichannel:
+        if image.shape[3] != 3:
             raise ValueError("Lab colorspace conversion requires a RGB image.")
         image = rgb2lab(image)
 
@@ -138,7 +154,7 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None,
     ratio = float(max((step_z, step_y, step_x))) / compactness
     image = np.ascontiguousarray(image * ratio)
 
-    labels = _slic_cython(image, segments, max_iter)
+    labels = _slic_cython(image, segments, max_iter, spacing)
 
     if is2d:
         labels = labels[0]

skimage/segmentation/tests/test_slic.py

@@ -90,6 +90,35 @@ def test_gray_3d():
         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()