diff --git a/skimage/morphology/ccomp.pyx b/skimage/morphology/ccomp.pyx index 54585f3a..590c5589 100644 --- a/skimage/morphology/ccomp.pyx +++ b/skimage/morphology/ccomp.pyx @@ -16,66 +16,91 @@ See also: Paper LBNL-56864, 2005, Lawrence Berkeley National Laboratory (University of California), - http://repositories.cdlib.org/lbnl/LBNL-56864. + http://repositories.cdlib.org/lbnl/LBNL-56864 """ # Tree operations implemented by an array as described in Wu et al. +# The term "forest" is used to indicate an array that stores one or more trees DTYPE = np.int ctypedef np.int_t DTYPE_t -cdef DTYPE_t find_root(np.int_t *work, np.int_t n): +cdef DTYPE_t find_root(np.int_t *forest, np.int_t n): """Find the root of node n. """ cdef np.int_t root = n - while (work[root] < root): - root = work[root] + while (forest[root] < root): + root = forest[root] return root -cdef set_root(np.int_t *work, np.int_t n, np.int_t root): +cdef set_root(np.int_t *forest, np.int_t n, np.int_t root): """ Set all nodes on a path to point to new_root. """ cdef np.int_t j - while (work[n] < n): - j = work[n] - work[n] = root + while (forest[n] < n): + j = forest[n] + forest[n] = root n = j - work[n] = root + forest[n] = root -cdef join_trees(np.int_t *work, np.int_t n, np.int_t m): +cdef join_trees(np.int_t *forest, np.int_t n, np.int_t m): """Join two trees containing nodes n and m. """ - cdef np.int_t root = find_root(work, n) + cdef np.int_t root = find_root(forest, n) cdef np.int_t root_m if (n != m): - root_m = find_root(work, m) + root_m = find_root(forest, m) if (root > root_m): root = root_m - set_root(work, n, root) - set_root(work, m, root) + set_root(forest, n, root) + set_root(forest, m, root) + +cdef link_bg(np.int_t *forest, np.int_t n, np.int_t *background_node): + """ + Link a node to the background node. + + """ + if background_node[0] == -999: + background_node[0] = n + + join_trees(forest, n, background_node[0]) # Connected components search as described in Fiorio et al. -def label(np.ndarray[DTYPE_t, ndim=2] input): +def label(np.ndarray[DTYPE_t, ndim=2] input, + np.int_t neighbors=8, np.int_t background=-1): """Label connected regions of an integer array. - Connectivity is defined as two (8-connected) neighboring entries - having equal value. + Two pixels are connected when they are neighbors and have the same value. + They can be neighbors either in a 4- or 8-connected sense:: + + 4-connectivity 8-connectivity + + [ ] [ ] [ ] [ ] + | \ | / + [ ]--[ ]--[ ] [ ]--[ ]--[ ] + | / | \ + [ ] [ ] [ ] [ ] Parameters ---------- input : ndarray of dtype int Image to label. + neighbors : {4, 8}, int + Whether to use 4- or 8-connectivity. + background : int + Consider all pixels with this value as background pixels, and label + them as -1. Returns ------- @@ -83,56 +108,102 @@ def label(np.ndarray[DTYPE_t, ndim=2] input): Labeled array, where all connected regions are assigned the same integer value. + Examples + -------- + >>> x = np.eye(3).astype(int) + >>> print x + [[1 0 0] + [0 1 0] + [0 0 1]] + + >>> print m.label(x, neighbors=4) + [[0 1 1] + [2 3 1] + [2 2 4]] + + >>> print m.label(x, neighbors=8) + [[0 1 1] + [1 0 1] + [1 1 0]] + + >>> x = np.array([[1, 0, 0], + ... [1, 1, 5], + ... [0, 0, 0]]) + + >>> print m.label(x, background=0) + [[ 0 -1 -1] + [ 0 0 1] + [-1 -1 -1]] + """ cdef np.int_t rows = input.shape[0] cdef np.int_t cols = input.shape[1] cdef np.ndarray[DTYPE_t, ndim=2] data = input.copy() - cdef np.ndarray[DTYPE_t, ndim=2] work + cdef np.ndarray[DTYPE_t, ndim=2] forest - work = np.arange(data.size, dtype=DTYPE).reshape((rows, cols)) + forest = np.arange(data.size, dtype=DTYPE).reshape((rows, cols)) - cdef np.int_t *work_p = work.data + cdef np.int_t *forest_p = forest.data cdef np.int_t *data_p = data.data cdef np.int_t i, j + cdef np.int_t background_node = -999 + + if neighbors != 4 and neighbors != 8: + raise ValueError('Neighbors must be either 4 or 8.') + # Initialize the first row for j in range(1, cols): + if data[0, j] == background: + link_bg(forest_p, j, &background_node) + if data[0, j] == data[0, j-1]: - join_trees(work_p, j, j-1) + join_trees(forest_p, j, j-1) for i in range(1, rows): # Handle the first column - if data[i, 0] == data[i-1, 0]: - join_trees(work_p, i*cols, (i-1)*cols) + if data[i, 0] == background: + link_bg(forest_p, i * cols, &background_node) - if data[i, 0] == data[i-1, 1]: - join_trees(work_p, i*cols, (i-1)*cols + 1) + if data[i, 0] == data[i-1, 0]: + join_trees(forest_p, i*cols, (i-1)*cols) + + if neighbors == 8: + if data[i, 0] == data[i-1, 1]: + join_trees(forest_p, i*cols, (i-1)*cols + 1) for j in range(1, cols): - if data[i, j] == data[i-1, j-1]: - join_trees(work_p, i*cols + j, (i-1)*cols + j - 1) + if data[i, j] == background: + link_bg(forest_p, i * cols + j, &background_node) + + if neighbors == 8: + if data[i, j] == data[i-1, j-1]: + join_trees(forest_p, i*cols + j, (i-1)*cols + j - 1) if data[i, j] == data[i-1, j]: - join_trees(work_p, i*cols + j, (i-1)*cols + j) + join_trees(forest_p, i*cols + j, (i-1)*cols + j) - if j < cols - 1: - if data[i, j] == data[i - 1, j + 1]: - join_trees(work_p, i*cols + j, (i-1)*cols + j + 1) + if neighbors == 8: + if j < cols - 1: + if data[i, j] == data[i - 1, j + 1]: + join_trees(forest_p, i*cols + j, (i-1)*cols + j + 1) if data[i, j] == data[i, j-1]: - join_trees(work_p, i*cols + j, i*cols + j - 1) + join_trees(forest_p, i*cols + j, i*cols + j - 1) # Label output cdef np.int_t ctr = 0 for i in range(rows): for j in range(cols): - if (i*cols + j) == work[i, j]: + if (i*cols + j) == background_node: + data[i, j] = -1 + elif (i*cols + j) == forest[i, j]: data[i, j] = ctr ctr = ctr + 1 else: - data[i, j] = data_p[work[i, j]] + data[i, j] = data_p[forest[i, j]] return data diff --git a/skimage/morphology/tests/test_ccomp.py b/skimage/morphology/tests/test_ccomp.py index 5bf4cfdc..8fb98009 100644 --- a/skimage/morphology/tests/test_ccomp.py +++ b/skimage/morphology/tests/test_ccomp.py @@ -37,5 +37,29 @@ class TestConnectedComponents: assert_array_equal(label(x), x) + def test_4_vs_8(self): + x = np.array([[0, 1], + [1, 0]], dtype=int) + assert_array_equal(label(x, 4), + [[0, 1], + [2, 3]]) + assert_array_equal(label(x, 8), + [[0, 1], + [1, 0]]) + + def test_background(self): + x = np.array([[1, 0, 0], + [1, 1, 5], + [0, 0, 0]]) + + assert_array_equal(label(x), [[0, 1, 1], + [0, 0, 2], + [3, 3, 3]]) + + assert_array_equal(label(x, background=0), + [[0, -1, -1], + [0, 0, 1], + [-1, -1, -1]]) + if __name__ == "__main__": run_module_suite()