diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx index 9cb96ab9..b3a00a52 100644 --- a/skimage/segmentation/_slic.pyx +++ b/skimage/segmentation/_slic.pyx @@ -9,11 +9,107 @@ cimport numpy as cnp from skimage.util import regular_grid +def _enforce_label_connectivity_cython(Py_ssize_t[:, :, ::1] nearest_segments, + Py_ssize_t n_segments, + int min_size): + """ Helper function to remove small disconnected regions from the labels + + Parameters + ---------- + nearest_segments : 3D array of int, shape (Z, Y, X) + The label field/superpixels found by SLIC. + n_segments: int + number of specified segments + min_size: int + minimum size of the segment + + Returns + ------- + connected_nearest_segments : 3D array of int, shape (Z, Y, X) + A label field with connected labels starting at label=1 + """ + + #get image dimensions + cdef Py_ssize_t depth, height, width + depth = nearest_segments.shape[0] + height = nearest_segments.shape[1] + width = nearest_segments.shape[2] + + #neighborhood arrays + cdef Py_ssize_t[:] ddx = np.array((1,-1,0,0,0,0)) + cdef Py_ssize_t[:] ddy = np.array((0,0,1,-1,0,0)) + cdef Py_ssize_t[:] ddz = np.array((0,0,0,0,1,-1)) + + cdef double size = height*width*depth / n_segments + cdef double max_size = 3*size + + #new object with connected segments + cdef Py_ssize_t[:, :, ::1] new_nearest_segments \ + = np.zeros((depth, height, width), dtype=np.intp) + + cdef Py_ssize_t current_new_label = 0 + cdef Py_ssize_t label = 0 + + #variables for the breadth first search + cdef Py_ssize_t count = 1 + cdef Py_ssize_t p = 0 + cdef Py_ssize_t adjacent + + cdef Py_ssize_t zz,yy,xx + + cdef Py_ssize_t[:, :] coord_list \ + = np.zeros((max_size,3), dtype=np.intp) + + #loop through all image + for z in range(depth): + for y in range(height): + for x in range(width): + if (new_nearest_segments[z,y,x] > 0): + continue + #find the component size + adjacent = 0 + label = nearest_segments[z,y,x] + current_new_label = current_new_label+1 + new_nearest_segments[z,y,x] = current_new_label + + count = 1 + p = 0 + coord_list[p,0] = z + coord_list[p,1] = y + coord_list[p,2] = x + + #perform a breadth first search to find the size of the connected component + while (p != count): + for i in range(6): + zz = coord_list[p,0] + ddz[i] + yy = coord_list[p,1] + ddy[i] + xx = coord_list[p,2] + ddx[i] + if (xx >= 0 and xx < width and yy >= 0 and yy < height and zz >= 0 and zz < depth): + if (nearest_segments[zz,yy,xx] == label and new_nearest_segments[zz,yy,xx] == 0): + new_nearest_segments[zz,yy,xx] = current_new_label + coord_list[count,0] = zz + coord_list[count,1] = yy + coord_list[count,2] = xx + count = count + 1 + elif (new_nearest_segments[zz,yy,xx] > 0 and new_nearest_segments[zz,yy,xx] != current_new_label): + adjacent = new_nearest_segments[zz,yy,xx] + p = p + 1 + + + #change to an adjacent one, like in the original paper + if (count < min_size): + #print("Changing segment {0} label {1} ".format(current_new_label, label)) + for i in range(count): + new_nearest_segments[coord_list[i,0],coord_list[i,1],coord_list[i,2]] = adjacent + + return np.asarray(new_nearest_segments) def _slic_cython(double[:, :, :, ::1] image_zyx, double[:, ::1] segments, Py_ssize_t max_iter, - double[::1] spacing): + double[::1] spacing, + int enforce_connectivity, + Py_ssize_t min_size = True): """Helper function for SLIC segmentation. Parameters @@ -28,6 +124,8 @@ def _slic_cython(double[:, :, :, ::1] image_zyx, The voxel spacing along each image dimension. This parameter controls the weights of the distances along z, y, and x during k-means clustering. + enforce_connectivity: int indicating whether the returned label + field must have connected labels Returns ------- @@ -145,73 +243,4 @@ def _slic_cython(double[:, :, :, ::1] image_zyx, for c in range(n_features): segments[k, c] /= n_segment_elems[k] - #return np.asarray(nearest_segments) - - #enforce segment connectivity - cdef Py_ssize_t[:] ddx = np.array((1,-1,0,0,0,0)) - cdef Py_ssize_t[:] ddy = np.array((0,0,1,-1,0,0)) - cdef Py_ssize_t[:] ddz = np.array((0,0,0,0,1,-1)) - - cdef double factor = 0.25 - cdef double size = height*width*depth / n_segments - cdef double min_size = factor * size - cdef double max_size = 3*size - - #new object with connected segments - cdef Py_ssize_t[:, :, ::1] new_nearest_segments \ - = np.zeros((depth, height, width), dtype=np.intp) - - cdef Py_ssize_t current_new_label = 0 - cdef Py_ssize_t label = 0 - - #variables for the breadth first search - cdef Py_ssize_t count = 1 - cdef Py_ssize_t p = 0 - cdef Py_ssize_t adjacent - - cdef Py_ssize_t zz,yy,xx - - cdef Py_ssize_t[:, :] coord_list \ - = np.zeros((max_size,3), dtype=np.intp) - - #loop through all image - for z in range(depth): - for y in range(height): - for x in range(width): - if (new_nearest_segments[z,y,x] > 0): - continue - #find the component size - adjacent = 0 - label = nearest_segments[z,y,x] - current_new_label = current_new_label+1 - new_nearest_segments[z,y,x] = current_new_label - - count = 1 - p = 0 - coord_list[p,0] = z - coord_list[p,1] = y - coord_list[p,2] = x - - #perform a breadth first search to find the size of the connected component - while (p != count): - for i in range(6): - zz = coord_list[p,0] + ddz[i] - yy = coord_list[p,1] + ddy[i] - xx = coord_list[p,2] + ddx[i] - if (xx >= 0 and xx < width and yy >= 0 and yy < height and zz >= 0 and zz < depth): - if (nearest_segments[zz,yy,xx] == label and new_nearest_segments[zz,yy,xx] == 0): - new_nearest_segments[zz,yy,xx] = current_new_label - coord_list[count,0] = zz - coord_list[count,1] = yy - coord_list[count,2] = xx - count = count + 1 - elif (new_nearest_segments[zz,yy,xx] > 0 and new_nearest_segments[zz,yy,xx] != current_new_label): - adjacent = new_nearest_segments[zz,yy,xx] - p = p + 1 - - #change to an adjacent one, like in the original paper - if (count < min_size): - for i in range(count): - new_nearest_segments[coord_list[i,0],coord_list[i,1],coord_list[i,2]] = adjacent - - return np.asarray(new_nearest_segments) + return np.asarray(nearest_segments) diff --git a/skimage/segmentation/slic_superpixels.py b/skimage/segmentation/slic_superpixels.py index 21df45fd..920d0e48 100644 --- a/skimage/segmentation/slic_superpixels.py +++ b/skimage/segmentation/slic_superpixels.py @@ -6,12 +6,13 @@ from scipy import ndimage import warnings from skimage.util import img_as_float, regular_grid -from skimage.segmentation._slic import _slic_cython +from skimage.segmentation._slic import _slic_cython, _enforce_label_connectivity_cython from skimage.color import rgb2lab def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None, - spacing=None, multichannel=True, convert2lab=True, ratio=None): + spacing=None, multichannel=True, convert2lab=True, ratio=None, + enforce_connectivity=True, min_size_factor=0.5): """Segments image using k-means clustering in Color-(x,y,z) space. Parameters @@ -47,6 +48,11 @@ def slic(image, n_segments=100, compactness=10., max_iter=10, sigma=None, recommended. ratio : float, optional Synonym for `compactness`. This keyword is deprecated. + enforce_connectivity: bool, optional + Whether the generated segments are connected or not + min_size_factor: float + proportion of the minimum segment size to be removed with respect + to the supposed segment size (depth*width*height/n_segments) Returns ------- @@ -161,7 +167,10 @@ 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, spacing) + labels = _slic_cython(image, segments, max_iter, spacing, enforce_connectivity) + + if (enforce_connectivity): + labels = _enforce_label_connectivity_cython(labels, n_segments, min_size_factor*depth*height*width/n_segments) if is_2d: labels = labels[0]