From b83ece8b0583093ea1692af4bcf84f0fb8f53066 Mon Sep 17 00:00:00 2001 From: Stefan van der Walt Date: Fri, 16 Dec 2011 23:47:32 -0800 Subject: [PATCH 1/3] ENH: Allow choice between 4 and 8 neighbor mode. Improve documentation. --- skimage/morphology/ccomp.pyx | 55 +++++++++++++++++++++----- skimage/morphology/tests/test_ccomp.py | 10 +++++ 2 files changed, 55 insertions(+), 10 deletions(-) diff --git a/skimage/morphology/ccomp.pyx b/skimage/morphology/ccomp.pyx index 54585f3a..b1fdc828 100644 --- a/skimage/morphology/ccomp.pyx +++ b/skimage/morphology/ccomp.pyx @@ -66,16 +66,27 @@ cdef join_trees(np.int_t *work, np.int_t n, np.int_t m): # 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, + int neighbors=8): """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. Returns ------- @@ -83,6 +94,24 @@ 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]] + """ cdef np.int_t rows = input.shape[0] cdef np.int_t cols = input.shape[1] @@ -97,6 +126,9 @@ def label(np.ndarray[DTYPE_t, ndim=2] input): cdef np.int_t i, j + 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] == data[0, j-1]: @@ -107,19 +139,22 @@ def label(np.ndarray[DTYPE_t, ndim=2] input): if data[i, 0] == data[i-1, 0]: join_trees(work_p, i*cols, (i-1)*cols) - if data[i, 0] == data[i-1, 1]: - join_trees(work_p, i*cols, (i-1)*cols + 1) + if neighbors == 8: + if data[i, 0] == data[i-1, 1]: + join_trees(work_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 neighbors == 8: + 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] == data[i-1, j]: join_trees(work_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(work_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) diff --git a/skimage/morphology/tests/test_ccomp.py b/skimage/morphology/tests/test_ccomp.py index 5bf4cfdc..04c21feb 100644 --- a/skimage/morphology/tests/test_ccomp.py +++ b/skimage/morphology/tests/test_ccomp.py @@ -37,5 +37,15 @@ 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]]) + if __name__ == "__main__": run_module_suite() From c78ea107aa2afd63bd4febe3ac304ef9d29f3fab Mon Sep 17 00:00:00 2001 From: Stefan van der Walt Date: Thu, 22 Dec 2011 04:47:24 -0800 Subject: [PATCH 2/3] ENH: Add background labelling. --- skimage/morphology/ccomp.pyx | 39 ++++++++++++++++++++++++-- skimage/morphology/tests/test_ccomp.py | 14 +++++++++ 2 files changed, 51 insertions(+), 2 deletions(-) diff --git a/skimage/morphology/ccomp.pyx b/skimage/morphology/ccomp.pyx index b1fdc828..84e124f8 100644 --- a/skimage/morphology/ccomp.pyx +++ b/skimage/morphology/ccomp.pyx @@ -64,10 +64,20 @@ cdef join_trees(np.int_t *work, np.int_t n, np.int_t m): set_root(work, n, root) set_root(work, m, root) +cdef link_bg(np.int_t *work, np.int_t n, np.int_t *background_node): + """ + Link a node to the background node. + + """ + if background_node[0] == -1: + background_node[0] = n + + join_trees(work, n, background_node[0]) + # Connected components search as described in Fiorio et al. def label(np.ndarray[DTYPE_t, ndim=2] input, - int neighbors=8): + np.int_t neighbors=8, np.int_t background=-1): """Label connected regions of an integer array. Two pixels are connected when they are neighbors and have the same value. @@ -87,6 +97,9 @@ def label(np.ndarray[DTYPE_t, ndim=2] input, 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 ------- @@ -112,6 +125,15 @@ def label(np.ndarray[DTYPE_t, ndim=2] input, [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] @@ -126,16 +148,24 @@ def label(np.ndarray[DTYPE_t, ndim=2] input, cdef np.int_t i, j + cdef np.int_t background_node = -1 + 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(work_p, j, &background_node) + if data[0, j] == data[0, j-1]: join_trees(work_p, j, j-1) for i in range(1, rows): # Handle the first column + if data[i, 0] == background: + link_bg(work_p, i * cols, &background_node) + if data[i, 0] == data[i-1, 0]: join_trees(work_p, i*cols, (i-1)*cols) @@ -144,6 +174,9 @@ def label(np.ndarray[DTYPE_t, ndim=2] input, join_trees(work_p, i*cols, (i-1)*cols + 1) for j in range(1, cols): + if data[i, j] == background: + link_bg(work_p, i * cols + j, &background_node) + if neighbors == 8: if data[i, j] == data[i-1, j-1]: join_trees(work_p, i*cols + j, (i-1)*cols + j - 1) @@ -164,7 +197,9 @@ def label(np.ndarray[DTYPE_t, ndim=2] input, 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) == work[i, j]: data[i, j] = ctr ctr = ctr + 1 else: diff --git a/skimage/morphology/tests/test_ccomp.py b/skimage/morphology/tests/test_ccomp.py index 04c21feb..8fb98009 100644 --- a/skimage/morphology/tests/test_ccomp.py +++ b/skimage/morphology/tests/test_ccomp.py @@ -47,5 +47,19 @@ class TestConnectedComponents: [[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() From e18fb678ec886dd3bc831bac7120c7fc81617da8 Mon Sep 17 00:00:00 2001 From: Stefan van der Walt Date: Fri, 3 Feb 2012 20:31:02 -0800 Subject: [PATCH 3/3] ENH: Cleanups as suggested by Tony. --- skimage/morphology/ccomp.pyx | 67 ++++++++++++++++++------------------ 1 file changed, 34 insertions(+), 33 deletions(-) diff --git a/skimage/morphology/ccomp.pyx b/skimage/morphology/ccomp.pyx index 84e124f8..590c5589 100644 --- a/skimage/morphology/ccomp.pyx +++ b/skimage/morphology/ccomp.pyx @@ -16,63 +16,64 @@ 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 *work, np.int_t n, np.int_t *background_node): +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] == -1: + if background_node[0] == -999: background_node[0] = n - join_trees(work, n, background_node[0]) + join_trees(forest, n, background_node[0]) # Connected components search as described in Fiorio et al. @@ -139,16 +140,16 @@ def label(np.ndarray[DTYPE_t, ndim=2] input, 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 = -1 + cdef np.int_t background_node = -999 if neighbors != 4 and neighbors != 8: raise ValueError('Neighbors must be either 4 or 8.') @@ -156,41 +157,41 @@ def label(np.ndarray[DTYPE_t, ndim=2] input, # Initialize the first row for j in range(1, cols): if data[0, j] == background: - link_bg(work_p, j, &background_node) + 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] == background: - link_bg(work_p, i * cols, &background_node) + link_bg(forest_p, i * cols, &background_node) if data[i, 0] == data[i-1, 0]: - join_trees(work_p, i*cols, (i-1)*cols) + join_trees(forest_p, i*cols, (i-1)*cols) if neighbors == 8: if data[i, 0] == data[i-1, 1]: - join_trees(work_p, i*cols, (i-1)*cols + 1) + join_trees(forest_p, i*cols, (i-1)*cols + 1) for j in range(1, cols): if data[i, j] == background: - link_bg(work_p, i * cols + j, &background_node) + link_bg(forest_p, i * cols + j, &background_node) if neighbors == 8: if data[i, j] == data[i-1, j-1]: - join_trees(work_p, i*cols + j, (i-1)*cols + 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 neighbors == 8: 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) + 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 @@ -199,10 +200,10 @@ def label(np.ndarray[DTYPE_t, ndim=2] input, for j in range(cols): if (i*cols + j) == background_node: data[i, j] = -1 - elif (i*cols + j) == work[i, j]: + 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