From f736047ebb4ef61409dec231e19b96f103f0f36d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mat=C4=9Bj=20T=C3=BD=C4=8D?= Date: Sat, 25 Oct 2014 18:34:32 +0200 Subject: [PATCH 01/23] Refactored measure.label, added documentation. --- skimage/measure/_ccomp.pyx | 358 +++++++++++++++++++++++++++++-------- 1 file changed, 288 insertions(+), 70 deletions(-) diff --git a/skimage/measure/_ccomp.pyx b/skimage/measure/_ccomp.pyx index 26b9645a..4ccb932a 100644 --- a/skimage/measure/_ccomp.pyx +++ b/skimage/measure/_ccomp.pyx @@ -24,15 +24,166 @@ See also: """ +DTYPE = np.intp + +# Short int - could be more graceful to the CPU cache +ctypedef cnp.int32_t INTS_t + +cdef struct s_shpinfo + +ctypedef s_shpinfo shpinfo +ctypedef s_bginfo bginfo +ctypedef int (* fun_ravel)(int, int, int, shpinfo *) + + +cdef struct s_bginfo: + DTYPE_t background_val + DTYPE_t background_node + + +# Structure for centralised access to shape data +cdef struct s_shpinfo: + INTS_t x + INTS_t y + INTS_t z + + DTYPE_t numels + INTS_t ndim + + #INTS_t Dee + INTS_t Ded + + INTS_t Dea + INTS_t Deb + INTS_t Dec + + INTS_t Def + INTS_t Deg + INTS_t Deh + INTS_t Dei + INTS_t Dej + INTS_t Dek + INTS_t Del + INTS_t Dem + INTS_t Den + + fun_ravel ravel_index + + +cdef shpinfo get_triple(inarr_shape): + cdef shpinfo res = shpinfo() + + res.y = 1 + res.z = 1 + res.ravel_index = ravel_index2D + + res.ndim = len(inarr_shape) + if res.ndim == 1: + res.x = inarr_shape[0] + res.ravel_index = ravel_index1D + elif res.ndim == 2: + res.x = inarr_shape[1] + res.y = inarr_shape[0] + res.ravel_index = ravel_index2D + elif res.ndim == 3: + res.x = inarr_shape[2] + res.y = inarr_shape[1] + res.z = inarr_shape[0] + res.ravel_index = ravel_index3D + else: + assert "Only for images of dimension 1-3 (got %s)" % res.ndim + + res.numels = res.x * res.y * res.z + + # Our point of interest is E. + # z=1 z=0 x + # ----------------------> + # | A B C F G H + # | D E . I J K + # | . . . L M N + # | + # y V + # + # Difference between E and G is (x=0, y=-1, z=-1), E and A (-1, -1, 0) etc. + # Here, it is recalculated to linear (raveled) indices of flattened arrays + # with their last (=contiguous) dimension is x. + + # So now the 1st (needed for 1D, 2D and 3D) part, y = 1, z = 1 + res.Ded = - 1 + + # Not needed, just for illustration + # + enabling it prolongs the exec time quite considerably - why? + #res.Dee = 0 + + # So now the 2nd (needed for 2D and 3D) part, y = 0, z = 1 + res.Dea = res.ravel_index(-1, -1, 0, & res) + res.Deb = res.Dea + 1 + res.Dec = res.Deb + 1 + + # And now the 3rd (needed only for 3D) part, z = 0 + res.Def = res.ravel_index(-1, -1, -1, & res) + res.Deg = res.Def + 1 + res.Deh = res.Def + 2 + res.Dei = res.Def - res.Deb # Deb = one row up, remember? + res.Dej = res.Dei + 1 + res.Dek = res.Dei + 2 + res.Del = res.Dei - 2 * res.Deb + res.Dem = res.Del + 1 + res.Den = res.Del + 2 + + return res + + +cdef int ravel_index1D(int x, int y, int z, shpinfo * shapeinfo): + """ + Ravel index of a 1D array - trivial. y and z are ignored. + """ + return x + + +cdef int ravel_index2D(int x, int y, int z, shpinfo * shapeinfo): + """ + Ravel index of a 2D array. z is ignored + """ + cdef int ret = x + y * shapeinfo.x + return ret + + +cdef int ravel_index3D(int x, int y, int z, shpinfo * shapeinfo): + """ + Ravel index of a 3D array + """ + cdef int ret = x + y * shapeinfo.x + z * shapeinfo.y * shapeinfo.x + return ret + + # 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.intp - +# Consider a following tree: +# +# 5 ----> 3 ----> 2 ----> 1 <---- 6 <---- 7 +# | | +# 4 >----/ \----< 8 <---- 9 +# +# The vertices are a unique number, so the tree can be represented by an +# array where a the tuple (index, array[index]) represents an edge, +# so for our example, array[2] == 1, array[7] == 6 and array[1] == 1, because +# 1 is the root. +# Last but not least, one array can hold more than one tree as long as their +# indices are different. It is the case in this algorithm, so for that reason +# the array is referred to as the "forrest" = multiple trees next to each +# other. +# +# In this algorithm, there are as many indices as there are elements in the +# array to label and array[x] == x for all x. As the labelling progresses, +# equivalence between so-called provisional (i.e. not final) labels is +# discovered and trees begin to surface. +# When we found out that label 5 and 3 are the same, we assign array[5] = 3. cdef DTYPE_t find_root(DTYPE_t *forest, DTYPE_t n): """Find the root of node n. - + Given the example above, for any integer from 1 to 9, 1 is always returned """ cdef DTYPE_t root = n while (forest[root] < root): @@ -43,7 +194,10 @@ cdef DTYPE_t find_root(DTYPE_t *forest, DTYPE_t n): cdef inline void set_root(DTYPE_t *forest, DTYPE_t n, DTYPE_t root): """ Set all nodes on a path to point to new_root. - + Given the example above, given n=9, root=6, it would "reconnect" the tree. + so forest[9] = 6 and forest[8] = 6 + The ultimate goal is that all tree nodes point to the real root, + which is element 1 in this case. """ cdef DTYPE_t j while (forest[n] < n): @@ -56,12 +210,17 @@ cdef inline void set_root(DTYPE_t *forest, DTYPE_t n, DTYPE_t root): cdef inline void join_trees(DTYPE_t *forest, DTYPE_t n, DTYPE_t m): """Join two trees containing nodes n and m. - + If we imagine that in the example tree, the root 1 is not known, we + rather have two disjoint trees with roots 2 and 6. + Joining them would mean that all elements of both trees become connected + to the element 2, so forest[9] == 2, forest[6] == 2 etc. + However, when the relationship between 1 and 2 can still be discovered later. """ - cdef DTYPE_t root = find_root(forest, n) + cdef DTYPE_t root cdef DTYPE_t root_m if (n != m): + root = find_root(forest, n) root_m = find_root(forest, m) if (root > root_m): @@ -147,94 +306,153 @@ def label(input, DTYPE_t neighbors=8, background=None, return_num=False): [-1 -1 -1]] """ - cdef DTYPE_t rows = input.shape[0] - cdef DTYPE_t cols = input.shape[1] - cdef cnp.ndarray[DTYPE_t, ndim=2] data = np.array(input, copy=True, - dtype=DTYPE) - cdef cnp.ndarray[DTYPE_t, ndim=2] forest + cdef cnp.ndarray[DTYPE_t, ndim=1] data + cdef cnp.ndarray[DTYPE_t, ndim=1] forest - forest = np.arange(data.size, dtype=DTYPE).reshape((rows, cols)) + # Having data a 2D array slows down access considerably using linear + # indices even when using the data_p pointer :-( + data = input.flatten().astype(DTYPE, copy=True) + forest = np.arange(data.size, dtype=DTYPE) cdef DTYPE_t *forest_p = forest.data cdef DTYPE_t *data_p = data.data - cdef DTYPE_t i, j + cdef shpinfo shapeinfo + cdef bginfo bg - cdef DTYPE_t background_val + shapeinfo = get_triple(input.shape) + + bg.background_val = 0 + bg.background_node = -999 if background is None: - background_val = -1 + bg.background_val = -1 warnings.warn(DeprecationWarning( 'The default value for `background` will change to 0 in v0.12' )) else: - background_val = background - - cdef DTYPE_t background_node = -999 + bg.background_val = background if neighbors != 4 and neighbors != 8: raise ValueError('Neighbors must be either 4 or 8.') - # Initialize the first row - if data[0, 0] == background_val: - link_bg(forest_p, 0, &background_node) - - for j in range(1, cols): - if data[0, j] == background_val: - link_bg(forest_p, j, &background_node) - - if data[0, j] == data[0, j-1]: - join_trees(forest_p, j, j-1) - - for i in range(1, rows): - # Handle the first column - if data[i, 0] == background_val: - link_bg(forest_p, i * cols, &background_node) - - 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] == background_val: - 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(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(forest_p, i*cols + j, (i-1)*cols + j + 1) - - if data[i, j] == data[i, j-1]: - join_trees(forest_p, i*cols + j, i*cols + j - 1) + if shapeinfo.ndim == 1: + scan1D(data_p, forest_p, & shapeinfo, & bg, neighbors, 0, 0) + elif shapeinfo.ndim == 2: + scan2D(data_p, forest_p, & shapeinfo, & bg, neighbors, 0) # Label output cdef DTYPE_t ctr = 0 - for i in range(rows): - for j in range(cols): - 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[forest[i, j]] + ctr = resolve_labels(data_p, forest_p, & shapeinfo, & bg) # Work around a bug in ndimage's type checking on 32-bit platforms if data.dtype == np.int32: data = data.view(np.int32) + res = data.reshape(input.shape) + if return_num: - return data, ctr + return res, ctr else: - return data + return res + + +cdef DTYPE_t resolve_labels(DTYPE_t * data_p, DTYPE_t * forest_p, + shpinfo * shapeinfo, bginfo * bg): + """ + We iterate through the provisional labels and assign final labels based on + our knowledge of prov. labels relationship. + We also track how many distinct final labels we have. + """ + cdef DTYPE_t counter = 0 + for i in range(shapeinfo.numels): + if i == bg.background_node: + data_p[i] = -1 + elif i == forest_p[i]: + data_p[i] = counter + counter += 1 + else: + data_p[i] = data_p[forest_p[i]] + return counter + + +# Here, we work with flat arrays regardless whether the data is 1, 2 or 3D. +# The lookup to the neighbor in a 2D array is achieved by precalculating an +# offset and ading it to the index. +# The forward scan mask looks like this (the center point is actually E): +# (take a look at shpinfo docs for more exhaustive info) +# A B C +# D E +# +# So if I am in the point E and want to take a look to A, I take the index of +# E and add shapeinfo.Dea to it and teg the index of A. +# The 1D indices are "raveled" or "linear", that's where "rindex" comes from. + + +cdef scan1D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, + bginfo * bg, DTYPE_t neighbors, DTYPE_t y, DTYPE_t z): + """ + Perform forward scan on a 1D object, usually the first row of an image + """ + # Initialize the first row + cdef DTYPE_t x, rindex + rindex = shapeinfo.ravel_index(0, y, z, shapeinfo) + + if data_p[rindex] == bg.background_val: + link_bg(forest_p, rindex, & bg.background_node) + + for x in range(1, shapeinfo.x): + rindex += 1 + # Handle the first row + # First row => rindex == j + if data_p[rindex] == bg.background_val: + link_bg(forest_p, rindex, & bg.background_node) + + if data_p[rindex] == data_p[rindex + shapeinfo.Ded]: + join_trees(forest_p, rindex, rindex + shapeinfo.Ded) + + +cdef scan2D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, + bginfo * bg, DTYPE_t neighbors, DTYPE_t z): + """ + Perform forward scan on a 2D array. + """ + cdef DTYPE_t x, y, rindex + scan1D(data_p, forest_p, shapeinfo, bg, neighbors, 0, z) + for y in range(1, shapeinfo.y): + rindex = shapeinfo.ravel_index(0, y, 0, shapeinfo) + # Handle the first column + if data_p[rindex] == bg.background_val: + link_bg(forest_p, rindex, & bg.background_node) + + if data_p[rindex] == data_p[rindex + shapeinfo.Deb]: + join_trees(forest_p, rindex, rindex + shapeinfo.Deb) + + if neighbors == 8: + if data_p[rindex] == data_p[rindex + shapeinfo.Dec]: + join_trees(forest_p, rindex, rindex + shapeinfo.Dec) + + # Handle the rest of columns + for x in range(1, shapeinfo.x): + # We have just moved to another column (of the same row) + # so we increment the raveled index. It will be reset when we get + # to another row, so we don't have to worry about altering it here. + rindex += 1 + if data_p[rindex] == bg.background_val: + link_bg(forest_p, rindex, & bg.background_node) + + if neighbors == 8: + if data_p[rindex] == data_p[rindex + shapeinfo.Dea]: + join_trees(forest_p, rindex, rindex + shapeinfo.Dea) + + if data_p[rindex] == data_p[rindex + shapeinfo.Deb]: + join_trees(forest_p, rindex, rindex + shapeinfo.Deb) + + if neighbors == 8: + if x < shapeinfo.x - 1: + if data_p[rindex] == data_p[rindex + shapeinfo.Dec]: + join_trees(forest_p, rindex, rindex + shapeinfo.Dec) + + if data_p[rindex] == data_p[rindex + shapeinfo.Ded]: + join_trees(forest_p, rindex, rindex + shapeinfo.Ded) From b65ebdfb13a305c8902dfe29f009537cdd887b05 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mat=C4=9Bj=20T=C3=BD=C4=8D?= Date: Sun, 26 Oct 2014 21:30:05 +0100 Subject: [PATCH 02/23] Added support for 3D arrays --- skimage/measure/_ccomp.pyx | 150 +++++++++++++++++++++++++++++++++++++ 1 file changed, 150 insertions(+) diff --git a/skimage/measure/_ccomp.pyx b/skimage/measure/_ccomp.pyx index 4ccb932a..6ba36bf8 100644 --- a/skimage/measure/_ccomp.pyx +++ b/skimage/measure/_ccomp.pyx @@ -341,6 +341,8 @@ def label(input, DTYPE_t neighbors=8, background=None, return_num=False): scan1D(data_p, forest_p, & shapeinfo, & bg, neighbors, 0, 0) elif shapeinfo.ndim == 2: scan2D(data_p, forest_p, & shapeinfo, & bg, neighbors, 0) + elif shapeinfo.ndim == 3: + scan3D(data_p, forest_p, & shapeinfo, & bg, neighbors) # Label output cdef DTYPE_t ctr = 0 @@ -456,3 +458,151 @@ cdef scan2D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, if data_p[rindex] == data_p[rindex + shapeinfo.Ded]: join_trees(forest_p, rindex, rindex + shapeinfo.Ded) + + +cdef scan3D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, + bginfo * bg, DTYPE_t neighbors): + """ + Perform forward scan on a 2D array. + """ + cdef DTYPE_t x, y, z, rindex + # Handle first plane + scan2D(data_p, forest_p, shapeinfo, bg, neighbors, 0) + for z in range(1, shapeinfo.z): + # Handle first row in 3D manner + rindex = shapeinfo.ravel_index(0, 0, z, shapeinfo) + if data_p[rindex] == bg.background_val: + link_bg(forest_p, rindex, & bg.background_node) + + # Now we have pixels below + if data_p[rindex] == data_p[rindex + shapeinfo.Dej]: + join_trees(forest_p, rindex, rindex + shapeinfo.Dej) + + if neighbors == 8: + if data_p[rindex] == data_p[rindex + shapeinfo.Dek]: + join_trees(forest_p, rindex, rindex + shapeinfo.Dek) + + if data_p[rindex] == data_p[rindex + shapeinfo.Dem]: + join_trees(forest_p, rindex, rindex + shapeinfo.Dem) + + if data_p[rindex] == data_p[rindex + shapeinfo.Den]: + join_trees(forest_p, rindex, rindex + shapeinfo.Den) + + for x in range(1, shapeinfo.x): + rindex += 1 + # Handle the first row + if data_p[rindex] == bg.background_val: + link_bg(forest_p, rindex, & bg.background_node) + + if data_p[rindex] == data_p[rindex + shapeinfo.Ded]: + join_trees(forest_p, rindex, rindex + shapeinfo.Ded) + + if neighbors == 8: + if data_p[rindex] == data_p[rindex + shapeinfo.Dei]: + join_trees(forest_p, rindex, rindex + shapeinfo.Dei) + + if data_p[rindex] == data_p[rindex + shapeinfo.Dej]: + join_trees(forest_p, rindex, rindex + shapeinfo.Dej) + + if neighbors == 8: + if x + 1 < shapeinfo.x: + if data_p[rindex] == data_p[rindex + shapeinfo.Dek]: + join_trees(forest_p, rindex, rindex + shapeinfo.Dek) + + if data_p[rindex] == data_p[rindex + shapeinfo.Del]: + join_trees(forest_p, rindex, rindex + shapeinfo.Del) + + if data_p[rindex] == data_p[rindex + shapeinfo.Dem]: + join_trees(forest_p, rindex, rindex + shapeinfo.Dem) + + if x + 1 < shapeinfo.x: + if data_p[rindex] == data_p[rindex + shapeinfo.Den]: + join_trees(forest_p, rindex, rindex + shapeinfo.Den) + + + for y in range(1, shapeinfo.y): + rindex = shapeinfo.ravel_index(0, y, z, shapeinfo) + # Handle the first column in 3D manner + if data_p[rindex] == bg.background_val: + link_bg(forest_p, rindex, & bg.background_node) + + if data_p[rindex] == data_p[rindex + shapeinfo.Deb]: + join_trees(forest_p, rindex, rindex + shapeinfo.Deb) + + if neighbors == 8: + if data_p[rindex] == data_p[rindex + shapeinfo.Dec]: + join_trees(forest_p, rindex, rindex + shapeinfo.Dec) + + if data_p[rindex] == data_p[rindex + shapeinfo.Deg]: + join_trees(forest_p, rindex, rindex + shapeinfo.Deg) + + if data_p[rindex] == data_p[rindex + shapeinfo.Deh]: + join_trees(forest_p, rindex, rindex + shapeinfo.Deh) + + if data_p[rindex] == data_p[rindex + shapeinfo.Dej]: + join_trees(forest_p, rindex, rindex + shapeinfo.Dej) + + if neighbors == 8: + if data_p[rindex] == data_p[rindex + shapeinfo.Dek]: + join_trees(forest_p, rindex, rindex + shapeinfo.Dek) + + if y + 1 < shapeinfo.y: + if data_p[rindex] == data_p[rindex + shapeinfo.Dem]: + join_trees(forest_p, rindex, rindex + shapeinfo.Dem) + + if data_p[rindex] == data_p[rindex + shapeinfo.Den]: + join_trees(forest_p, rindex, rindex + shapeinfo.Den) + + # Handle the rest of columns + for x in range(1, shapeinfo.x): + rindex += 1 + if data_p[rindex] == bg.background_val: + link_bg(forest_p, rindex, & bg.background_node) + + if neighbors == 8: + if data_p[rindex] == data_p[rindex + shapeinfo.Dea]: + join_trees(forest_p, rindex, rindex + shapeinfo.Dea) + + if data_p[rindex] == data_p[rindex + shapeinfo.Deb]: + join_trees(forest_p, rindex, rindex + shapeinfo.Deb) + + if neighbors == 8: + if x < shapeinfo.x - 1: + if data_p[rindex] == data_p[rindex + shapeinfo.Dec]: + join_trees(forest_p, rindex, rindex + shapeinfo.Dec) + + if data_p[rindex] == data_p[rindex + shapeinfo.Ded]: + join_trees(forest_p, rindex, rindex + shapeinfo.Ded) + + # Now pixels below: + if neighbors == 8: + if data_p[rindex] == data_p[rindex + shapeinfo.Def]: + join_trees(forest_p, rindex, rindex + shapeinfo.Def) + + if data_p[rindex] == data_p[rindex + shapeinfo.Deg]: + join_trees(forest_p, rindex, rindex + shapeinfo.Deg) + + if x + 1 < shapeinfo.x: + if data_p[rindex] == data_p[rindex + shapeinfo.Deh]: + join_trees(forest_p, rindex, rindex + shapeinfo.Deh) + + if data_p[rindex] == data_p[rindex + shapeinfo.Dei]: + join_trees(forest_p, rindex, rindex + shapeinfo.Dei) + + if data_p[rindex] == data_p[rindex + shapeinfo.Dej]: + join_trees(forest_p, rindex, rindex + shapeinfo.Dej) + + if neighbors == 8: + if x + 1 < shapeinfo.x: + if data_p[rindex] == data_p[rindex + shapeinfo.Dek]: + join_trees(forest_p, rindex, rindex + shapeinfo.Dek) + + if data_p[rindex] == data_p[rindex + shapeinfo.Del]: + join_trees(forest_p, rindex, rindex + shapeinfo.Del) + + if data_p[rindex] == data_p[rindex + shapeinfo.Dem]: + join_trees(forest_p, rindex, rindex + shapeinfo.Dem) + + if x + 1 < shapeinfo.x: + if data_p[rindex] == data_p[rindex + shapeinfo.Den]: + join_trees(forest_p, rindex, rindex + shapeinfo.Den) From d3049a087cf4b729af6d3460f75bc97a144d0eb1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mat=C4=9Bj=20T=C3=BD=C4=8D?= Date: Sun, 26 Oct 2014 23:22:54 +0100 Subject: [PATCH 03/23] Refactored Cython stuff with easier code readability in mind --- skimage/measure/_ccomp.pyx | 189 ++++++++++++++----------------------- 1 file changed, 72 insertions(+), 117 deletions(-) diff --git a/skimage/measure/_ccomp.pyx b/skimage/measure/_ccomp.pyx index 6ba36bf8..e537b598 100644 --- a/skimage/measure/_ccomp.pyx +++ b/skimage/measure/_ccomp.pyx @@ -32,15 +32,22 @@ ctypedef cnp.int32_t INTS_t cdef struct s_shpinfo ctypedef s_shpinfo shpinfo -ctypedef s_bginfo bginfo ctypedef int (* fun_ravel)(int, int, int, shpinfo *) -cdef struct s_bginfo: +ctypedef struct bginfo: DTYPE_t background_val DTYPE_t background_node +cdef enum: + # D_ee, # We don't need D_ee + D_ed, + D_ea, D_eb, D_ec, + D_ef, D_eg, D_eh, D_ei, D_ej, D_ek, D_el, D_em, D_en, + D_COUNT + + # Structure for centralised access to shape data cdef struct s_shpinfo: INTS_t x @@ -50,22 +57,7 @@ cdef struct s_shpinfo: DTYPE_t numels INTS_t ndim - #INTS_t Dee - INTS_t Ded - - INTS_t Dea - INTS_t Deb - INTS_t Dec - - INTS_t Def - INTS_t Deg - INTS_t Deh - INTS_t Dei - INTS_t Dej - INTS_t Dek - INTS_t Del - INTS_t Dem - INTS_t Den + INTS_t DEX[D_COUNT] fun_ravel ravel_index @@ -109,31 +101,36 @@ cdef shpinfo get_triple(inarr_shape): # with their last (=contiguous) dimension is x. # So now the 1st (needed for 1D, 2D and 3D) part, y = 1, z = 1 - res.Ded = - 1 + res.DEX[D_ed] = -1 # Not needed, just for illustration - # + enabling it prolongs the exec time quite considerably - why? - #res.Dee = 0 + # res.DEX[D_ee] = 0 # So now the 2nd (needed for 2D and 3D) part, y = 0, z = 1 - res.Dea = res.ravel_index(-1, -1, 0, & res) - res.Deb = res.Dea + 1 - res.Dec = res.Deb + 1 + res.DEX[D_ea] = res.ravel_index(-1, -1, 0, & res) + res.DEX[D_eb] = res.DEX[D_ea] + 1 + res.DEX[D_ec] = res.DEX[D_eb] + 1 # And now the 3rd (needed only for 3D) part, z = 0 - res.Def = res.ravel_index(-1, -1, -1, & res) - res.Deg = res.Def + 1 - res.Deh = res.Def + 2 - res.Dei = res.Def - res.Deb # Deb = one row up, remember? - res.Dej = res.Dei + 1 - res.Dek = res.Dei + 2 - res.Del = res.Dei - 2 * res.Deb - res.Dem = res.Del + 1 - res.Den = res.Del + 2 + res.DEX[D_ef] = res.ravel_index(-1, -1, -1, & res) + res.DEX[D_eg] = res.DEX[D_ef] + 1 + res.DEX[D_eh] = res.DEX[D_ef] + 2 + res.DEX[D_ei] = res.DEX[D_ef] - res.DEX[D_eb] # DEX[D_eb] = one row up, remember? + res.DEX[D_ej] = res.DEX[D_ei] + 1 + res.DEX[D_ek] = res.DEX[D_ei] + 2 + res.DEX[D_el] = res.DEX[D_ei] - 2 * res.DEX[D_eb] + res.DEX[D_em] = res.DEX[D_el] + 1 + res.DEX[D_en] = res.DEX[D_el] + 2 return res +cdef inline void join_trees_wrapper(DTYPE_t * data_p, DTYPE_t * forest_p, + DTYPE_t rindex, INTS_t idxdiff): + if data_p[rindex] == data_p[rindex + idxdiff]: + join_trees(forest_p, rindex, rindex + idxdiff) + + cdef int ravel_index1D(int x, int y, int z, shpinfo * shapeinfo): """ Ravel index of a 1D array - trivial. y and z are ignored. @@ -399,6 +396,7 @@ cdef scan1D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, """ # Initialize the first row cdef DTYPE_t x, rindex + cdef INTS_t * DEX = shapeinfo.DEX rindex = shapeinfo.ravel_index(0, y, z, shapeinfo) if data_p[rindex] == bg.background_val: @@ -411,8 +409,7 @@ cdef scan1D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, if data_p[rindex] == bg.background_val: link_bg(forest_p, rindex, & bg.background_node) - if data_p[rindex] == data_p[rindex + shapeinfo.Ded]: - join_trees(forest_p, rindex, rindex + shapeinfo.Ded) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) cdef scan2D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, @@ -421,6 +418,7 @@ cdef scan2D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, Perform forward scan on a 2D array. """ cdef DTYPE_t x, y, rindex + cdef INTS_t * DEX = shapeinfo.DEX scan1D(data_p, forest_p, shapeinfo, bg, neighbors, 0, z) for y in range(1, shapeinfo.y): rindex = shapeinfo.ravel_index(0, y, 0, shapeinfo) @@ -428,12 +426,10 @@ cdef scan2D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, if data_p[rindex] == bg.background_val: link_bg(forest_p, rindex, & bg.background_node) - if data_p[rindex] == data_p[rindex + shapeinfo.Deb]: - join_trees(forest_p, rindex, rindex + shapeinfo.Deb) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) if neighbors == 8: - if data_p[rindex] == data_p[rindex + shapeinfo.Dec]: - join_trees(forest_p, rindex, rindex + shapeinfo.Dec) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) # Handle the rest of columns for x in range(1, shapeinfo.x): @@ -445,19 +441,15 @@ cdef scan2D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, link_bg(forest_p, rindex, & bg.background_node) if neighbors == 8: - if data_p[rindex] == data_p[rindex + shapeinfo.Dea]: - join_trees(forest_p, rindex, rindex + shapeinfo.Dea) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ea]) - if data_p[rindex] == data_p[rindex + shapeinfo.Deb]: - join_trees(forest_p, rindex, rindex + shapeinfo.Deb) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) if neighbors == 8: if x < shapeinfo.x - 1: - if data_p[rindex] == data_p[rindex + shapeinfo.Dec]: - join_trees(forest_p, rindex, rindex + shapeinfo.Dec) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) - if data_p[rindex] == data_p[rindex + shapeinfo.Ded]: - join_trees(forest_p, rindex, rindex + shapeinfo.Ded) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) cdef scan3D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, @@ -466,6 +458,7 @@ cdef scan3D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, Perform forward scan on a 2D array. """ cdef DTYPE_t x, y, z, rindex + cdef INTS_t * DEX = shapeinfo.DEX # Handle first plane scan2D(data_p, forest_p, shapeinfo, bg, neighbors, 0) for z in range(1, shapeinfo.z): @@ -475,18 +468,12 @@ cdef scan3D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, link_bg(forest_p, rindex, & bg.background_node) # Now we have pixels below - if data_p[rindex] == data_p[rindex + shapeinfo.Dej]: - join_trees(forest_p, rindex, rindex + shapeinfo.Dej) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) if neighbors == 8: - if data_p[rindex] == data_p[rindex + shapeinfo.Dek]: - join_trees(forest_p, rindex, rindex + shapeinfo.Dek) - - if data_p[rindex] == data_p[rindex + shapeinfo.Dem]: - join_trees(forest_p, rindex, rindex + shapeinfo.Dem) - - if data_p[rindex] == data_p[rindex + shapeinfo.Den]: - join_trees(forest_p, rindex, rindex + shapeinfo.Den) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_en]) for x in range(1, shapeinfo.x): rindex += 1 @@ -494,30 +481,22 @@ cdef scan3D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, if data_p[rindex] == bg.background_val: link_bg(forest_p, rindex, & bg.background_node) - if data_p[rindex] == data_p[rindex + shapeinfo.Ded]: - join_trees(forest_p, rindex, rindex + shapeinfo.Ded) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) if neighbors == 8: - if data_p[rindex] == data_p[rindex + shapeinfo.Dei]: - join_trees(forest_p, rindex, rindex + shapeinfo.Dei) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) - if data_p[rindex] == data_p[rindex + shapeinfo.Dej]: - join_trees(forest_p, rindex, rindex + shapeinfo.Dej) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) if neighbors == 8: if x + 1 < shapeinfo.x: - if data_p[rindex] == data_p[rindex + shapeinfo.Dek]: - join_trees(forest_p, rindex, rindex + shapeinfo.Dek) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) - if data_p[rindex] == data_p[rindex + shapeinfo.Del]: - join_trees(forest_p, rindex, rindex + shapeinfo.Del) - - if data_p[rindex] == data_p[rindex + shapeinfo.Dem]: - join_trees(forest_p, rindex, rindex + shapeinfo.Dem) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_el]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) if x + 1 < shapeinfo.x: - if data_p[rindex] == data_p[rindex + shapeinfo.Den]: - join_trees(forest_p, rindex, rindex + shapeinfo.Den) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_en]) for y in range(1, shapeinfo.y): @@ -526,32 +505,21 @@ cdef scan3D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, if data_p[rindex] == bg.background_val: link_bg(forest_p, rindex, & bg.background_node) - if data_p[rindex] == data_p[rindex + shapeinfo.Deb]: - join_trees(forest_p, rindex, rindex + shapeinfo.Deb) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) if neighbors == 8: - if data_p[rindex] == data_p[rindex + shapeinfo.Dec]: - join_trees(forest_p, rindex, rindex + shapeinfo.Dec) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eh]) - if data_p[rindex] == data_p[rindex + shapeinfo.Deg]: - join_trees(forest_p, rindex, rindex + shapeinfo.Deg) - - if data_p[rindex] == data_p[rindex + shapeinfo.Deh]: - join_trees(forest_p, rindex, rindex + shapeinfo.Deh) - - if data_p[rindex] == data_p[rindex + shapeinfo.Dej]: - join_trees(forest_p, rindex, rindex + shapeinfo.Dej) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) if neighbors == 8: - if data_p[rindex] == data_p[rindex + shapeinfo.Dek]: - join_trees(forest_p, rindex, rindex + shapeinfo.Dek) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) if y + 1 < shapeinfo.y: - if data_p[rindex] == data_p[rindex + shapeinfo.Dem]: - join_trees(forest_p, rindex, rindex + shapeinfo.Dem) - - if data_p[rindex] == data_p[rindex + shapeinfo.Den]: - join_trees(forest_p, rindex, rindex + shapeinfo.Den) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_en]) # Handle the rest of columns for x in range(1, shapeinfo.x): @@ -560,49 +528,36 @@ cdef scan3D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, link_bg(forest_p, rindex, & bg.background_node) if neighbors == 8: - if data_p[rindex] == data_p[rindex + shapeinfo.Dea]: - join_trees(forest_p, rindex, rindex + shapeinfo.Dea) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ea]) - if data_p[rindex] == data_p[rindex + shapeinfo.Deb]: - join_trees(forest_p, rindex, rindex + shapeinfo.Deb) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) if neighbors == 8: if x < shapeinfo.x - 1: - if data_p[rindex] == data_p[rindex + shapeinfo.Dec]: - join_trees(forest_p, rindex, rindex + shapeinfo.Dec) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) - if data_p[rindex] == data_p[rindex + shapeinfo.Ded]: - join_trees(forest_p, rindex, rindex + shapeinfo.Ded) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) # Now pixels below: if neighbors == 8: - if data_p[rindex] == data_p[rindex + shapeinfo.Def]: - join_trees(forest_p, rindex, rindex + shapeinfo.Def) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ef]) - if data_p[rindex] == data_p[rindex + shapeinfo.Deg]: - join_trees(forest_p, rindex, rindex + shapeinfo.Deg) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) if x + 1 < shapeinfo.x: - if data_p[rindex] == data_p[rindex + shapeinfo.Deh]: - join_trees(forest_p, rindex, rindex + shapeinfo.Deh) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eh]) - if data_p[rindex] == data_p[rindex + shapeinfo.Dei]: - join_trees(forest_p, rindex, rindex + shapeinfo.Dei) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) - if data_p[rindex] == data_p[rindex + shapeinfo.Dej]: - join_trees(forest_p, rindex, rindex + shapeinfo.Dej) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) if neighbors == 8: if x + 1 < shapeinfo.x: - if data_p[rindex] == data_p[rindex + shapeinfo.Dek]: - join_trees(forest_p, rindex, rindex + shapeinfo.Dek) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) - if data_p[rindex] == data_p[rindex + shapeinfo.Del]: - join_trees(forest_p, rindex, rindex + shapeinfo.Del) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_el]) - if data_p[rindex] == data_p[rindex + shapeinfo.Dem]: - join_trees(forest_p, rindex, rindex + shapeinfo.Dem) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) if x + 1 < shapeinfo.x: - if data_p[rindex] == data_p[rindex + shapeinfo.Den]: - join_trees(forest_p, rindex, rindex + shapeinfo.Den) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_en]) From 98047f5ce87939daf39c6c0e692829e1a04167a4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mat=C4=9Bj=20T=C3=BD=C4=8D?= Date: Sun, 26 Oct 2014 23:38:12 +0100 Subject: [PATCH 04/23] Cython minor fixes --- skimage/measure/_ccomp.pyx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/skimage/measure/_ccomp.pyx b/skimage/measure/_ccomp.pyx index e537b598..03be39d7 100644 --- a/skimage/measure/_ccomp.pyx +++ b/skimage/measure/_ccomp.pyx @@ -63,7 +63,7 @@ cdef struct s_shpinfo: cdef shpinfo get_triple(inarr_shape): - cdef shpinfo res = shpinfo() + cdef shpinfo res res.y = 1 res.z = 1 @@ -389,7 +389,7 @@ cdef DTYPE_t resolve_labels(DTYPE_t * data_p, DTYPE_t * forest_p, # The 1D indices are "raveled" or "linear", that's where "rindex" comes from. -cdef scan1D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, +cdef void scan1D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, bginfo * bg, DTYPE_t neighbors, DTYPE_t y, DTYPE_t z): """ Perform forward scan on a 1D object, usually the first row of an image @@ -412,7 +412,7 @@ cdef scan1D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) -cdef scan2D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, +cdef void scan2D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, bginfo * bg, DTYPE_t neighbors, DTYPE_t z): """ Perform forward scan on a 2D array. @@ -452,7 +452,7 @@ cdef scan2D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) -cdef scan3D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, +cdef void scan3D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, bginfo * bg, DTYPE_t neighbors): """ Perform forward scan on a 2D array. From 152a6d0de84d0ff279d755a4141065f05dafdb06 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mat=C4=9Bj=20T=C3=BD=C4=8D?= Date: Mon, 27 Oct 2014 22:37:04 +0100 Subject: [PATCH 05/23] Fixed typos + added docstrings and comments --- skimage/measure/_ccomp.pyx | 61 ++++++++++++++++++++++++-------------- 1 file changed, 39 insertions(+), 22 deletions(-) diff --git a/skimage/measure/_ccomp.pyx b/skimage/measure/_ccomp.pyx index 03be39d7..be0148cf 100644 --- a/skimage/measure/_ccomp.pyx +++ b/skimage/measure/_ccomp.pyx @@ -31,8 +31,8 @@ ctypedef cnp.int32_t INTS_t cdef struct s_shpinfo -ctypedef s_shpinfo shpinfo -ctypedef int (* fun_ravel)(int, int, int, shpinfo *) +ctypedef s_shpinfo shape_info +ctypedef int (* fun_ravel)(int, int, int, shape_info *) ctypedef struct bginfo: @@ -49,21 +49,33 @@ cdef enum: # Structure for centralised access to shape data +# Contains information related to the shape of the input array cdef struct s_shpinfo: INTS_t x INTS_t y INTS_t z + # Number of elements DTYPE_t numels + # Number of of the input array INTS_t ndim + # Offsets between elements recalculated to linear index increments + # DEX[D_ea] is offset between E and A (i.e. to the point to upper left) + # The name DEX is supposed to evoke DE., where . = A, B, C, D, F etc. INTS_t DEX[D_COUNT] + # Function pointer to a function that recalculates multiindex to linear + # index. Heavily depends on dimensions of the input array. fun_ravel ravel_index -cdef shpinfo get_triple(inarr_shape): - cdef shpinfo res +cdef shape_info get_shape_info(inarr_shape): + """ + Precalculates all the needed data from the input array shape + and stores them in the shape_info struct. + """ + cdef shape_info res res.y = 1 res.z = 1 @@ -131,14 +143,14 @@ cdef inline void join_trees_wrapper(DTYPE_t * data_p, DTYPE_t * forest_p, join_trees(forest_p, rindex, rindex + idxdiff) -cdef int ravel_index1D(int x, int y, int z, shpinfo * shapeinfo): +cdef int ravel_index1D(int x, int y, int z, shape_info * shapeinfo): """ Ravel index of a 1D array - trivial. y and z are ignored. """ return x -cdef int ravel_index2D(int x, int y, int z, shpinfo * shapeinfo): +cdef int ravel_index2D(int x, int y, int z, shape_info * shapeinfo): """ Ravel index of a 2D array. z is ignored """ @@ -146,7 +158,7 @@ cdef int ravel_index2D(int x, int y, int z, shpinfo * shapeinfo): return ret -cdef int ravel_index3D(int x, int y, int z, shpinfo * shapeinfo): +cdef int ravel_index3D(int x, int y, int z, shape_info * shapeinfo): """ Ravel index of a 3D array """ @@ -169,7 +181,7 @@ cdef int ravel_index3D(int x, int y, int z, shpinfo * shapeinfo): # 1 is the root. # Last but not least, one array can hold more than one tree as long as their # indices are different. It is the case in this algorithm, so for that reason -# the array is referred to as the "forrest" = multiple trees next to each +# the array is referred to as the "forest" = multiple trees next to each # other. # # In this algorithm, there are as many indices as there are elements in the @@ -259,6 +271,8 @@ def label(input, DTYPE_t neighbors=8, background=None, return_num=False): Image to label. neighbors : {4, 8}, int, optional Whether to use 4- or 8-connectivity. + In 3D, 4-connectivity means connected pixels share have to share face, + whereas with 8-connectivity, they have to share only edge or vertex. background : int, optional Consider all pixels with this value as background pixels, and label them as -1. (Note: background pixels will be labeled as 0 starting with @@ -309,16 +323,16 @@ def label(input, DTYPE_t neighbors=8, background=None, return_num=False): # Having data a 2D array slows down access considerably using linear # indices even when using the data_p pointer :-( - data = input.flatten().astype(DTYPE, copy=True) + data = np.copy(input.flatten().astype(DTYPE), order="C") forest = np.arange(data.size, dtype=DTYPE) cdef DTYPE_t *forest_p = forest.data cdef DTYPE_t *data_p = data.data - cdef shpinfo shapeinfo + cdef shape_info shapeinfo cdef bginfo bg - shapeinfo = get_triple(input.shape) + shapeinfo = get_shape_info(input.shape) bg.background_val = 0 bg.background_node = -999 @@ -342,7 +356,7 @@ def label(input, DTYPE_t neighbors=8, background=None, return_num=False): scan3D(data_p, forest_p, & shapeinfo, & bg, neighbors) # Label output - cdef DTYPE_t ctr = 0 + cdef DTYPE_t ctr ctr = resolve_labels(data_p, forest_p, & shapeinfo, & bg) # Work around a bug in ndimage's type checking on 32-bit platforms @@ -358,17 +372,20 @@ def label(input, DTYPE_t neighbors=8, background=None, return_num=False): cdef DTYPE_t resolve_labels(DTYPE_t * data_p, DTYPE_t * forest_p, - shpinfo * shapeinfo, bginfo * bg): + shape_info * shapeinfo, bginfo * bg): """ We iterate through the provisional labels and assign final labels based on our knowledge of prov. labels relationship. We also track how many distinct final labels we have. """ - cdef DTYPE_t counter = 0 + cdef DTYPE_t counter = 0, i + for i in range(shapeinfo.numels): if i == bg.background_node: data_p[i] = -1 elif i == forest_p[i]: + # We have stumbled across a root which is something new to us (root + # is the LOWEST of all prov. labels that are equivalent to it) data_p[i] = counter counter += 1 else: @@ -378,9 +395,9 @@ cdef DTYPE_t resolve_labels(DTYPE_t * data_p, DTYPE_t * forest_p, # Here, we work with flat arrays regardless whether the data is 1, 2 or 3D. # The lookup to the neighbor in a 2D array is achieved by precalculating an -# offset and ading it to the index. +# offset and adding it to the index. # The forward scan mask looks like this (the center point is actually E): -# (take a look at shpinfo docs for more exhaustive info) +# (take a look at shape_info docs for more exhaustive info) # A B C # D E # @@ -389,8 +406,8 @@ cdef DTYPE_t resolve_labels(DTYPE_t * data_p, DTYPE_t * forest_p, # The 1D indices are "raveled" or "linear", that's where "rindex" comes from. -cdef void scan1D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, - bginfo * bg, DTYPE_t neighbors, DTYPE_t y, DTYPE_t z): +cdef void scan1D(DTYPE_t * data_p, DTYPE_t * forest_p, shape_info * shapeinfo, + bginfo * bg, DTYPE_t neighbors, DTYPE_t y, DTYPE_t z): """ Perform forward scan on a 1D object, usually the first row of an image """ @@ -412,8 +429,8 @@ cdef void scan1D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) -cdef void scan2D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, - bginfo * bg, DTYPE_t neighbors, DTYPE_t z): +cdef void scan2D(DTYPE_t * data_p, DTYPE_t * forest_p, shape_info * shapeinfo, + bginfo * bg, DTYPE_t neighbors, DTYPE_t z): """ Perform forward scan on a 2D array. """ @@ -452,8 +469,8 @@ cdef void scan2D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) -cdef void scan3D(DTYPE_t * data_p, DTYPE_t * forest_p, shpinfo * shapeinfo, - bginfo * bg, DTYPE_t neighbors): +cdef void scan3D(DTYPE_t * data_p, DTYPE_t * forest_p, shape_info * shapeinfo, + bginfo * bg, DTYPE_t neighbors): """ Perform forward scan on a 2D array. """ From 416f3a2eb55c8815717cb8ba3cd1d78d55fb7250 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mat=C4=9Bj=20T=C3=BD=C4=8D?= Date: Mon, 27 Oct 2014 23:33:59 +0100 Subject: [PATCH 06/23] numpy 1.6 compat fix --- skimage/measure/_ccomp.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/measure/_ccomp.pyx b/skimage/measure/_ccomp.pyx index be0148cf..bed57069 100644 --- a/skimage/measure/_ccomp.pyx +++ b/skimage/measure/_ccomp.pyx @@ -323,7 +323,7 @@ def label(input, DTYPE_t neighbors=8, background=None, return_num=False): # Having data a 2D array slows down access considerably using linear # indices even when using the data_p pointer :-( - data = np.copy(input.flatten().astype(DTYPE), order="C") + data = np.copy(input.flatten().astype(DTYPE)) forest = np.arange(data.size, dtype=DTYPE) cdef DTYPE_t *forest_p = forest.data From 815ea29d26dad62f2e88620f2127ca22f42580e0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mat=C4=9Bj=20T=C3=BD=C4=8D?= Date: Wed, 29 Oct 2014 00:58:01 +0100 Subject: [PATCH 07/23] Added some comments + requested reformatting --- skimage/measure/_ccomp.pyx | 37 ++++++++++++++++++++++--------------- 1 file changed, 22 insertions(+), 15 deletions(-) diff --git a/skimage/measure/_ccomp.pyx b/skimage/measure/_ccomp.pyx index bed57069..cc4a5320 100644 --- a/skimage/measure/_ccomp.pyx +++ b/skimage/measure/_ccomp.pyx @@ -35,15 +35,22 @@ ctypedef s_shpinfo shape_info ctypedef int (* fun_ravel)(int, int, int, shape_info *) +# For having stuff concerning background in one place ctypedef struct bginfo: DTYPE_t background_val DTYPE_t background_node +# A pixel has neighbors that have already been scanned. +# In the paper, the pixel is denoted by E and its neighbors +# by A, B, C and D (in 2D) - see doc for function get_shape_info() cdef enum: # D_ee, # We don't need D_ee + # the 1D neighbor D_ed, + # 2D neighbors D_ea, D_eb, D_ec, + # 3D neighbors D_ef, D_eg, D_eh, D_ei, D_ej, D_ek, D_el, D_em, D_en, D_COUNT @@ -137,20 +144,20 @@ cdef shape_info get_shape_info(inarr_shape): return res -cdef inline void join_trees_wrapper(DTYPE_t * data_p, DTYPE_t * forest_p, +cdef inline void join_trees_wrapper(DTYPE_t *data_p, DTYPE_t *forest_p, DTYPE_t rindex, INTS_t idxdiff): if data_p[rindex] == data_p[rindex + idxdiff]: join_trees(forest_p, rindex, rindex + idxdiff) -cdef int ravel_index1D(int x, int y, int z, shape_info * shapeinfo): +cdef int ravel_index1D(int x, int y, int z, shape_info *shapeinfo): """ Ravel index of a 1D array - trivial. y and z are ignored. """ return x -cdef int ravel_index2D(int x, int y, int z, shape_info * shapeinfo): +cdef int ravel_index2D(int x, int y, int z, shape_info *shapeinfo): """ Ravel index of a 2D array. z is ignored """ @@ -158,7 +165,7 @@ cdef int ravel_index2D(int x, int y, int z, shape_info * shapeinfo): return ret -cdef int ravel_index3D(int x, int y, int z, shape_info * shapeinfo): +cdef int ravel_index3D(int x, int y, int z, shape_info *shapeinfo): """ Ravel index of a 3D array """ @@ -371,8 +378,8 @@ def label(input, DTYPE_t neighbors=8, background=None, return_num=False): return res -cdef DTYPE_t resolve_labels(DTYPE_t * data_p, DTYPE_t * forest_p, - shape_info * shapeinfo, bginfo * bg): +cdef DTYPE_t resolve_labels(DTYPE_t *data_p, DTYPE_t *forest_p, + shape_info *shapeinfo, bginfo *bg): """ We iterate through the provisional labels and assign final labels based on our knowledge of prov. labels relationship. @@ -406,14 +413,14 @@ cdef DTYPE_t resolve_labels(DTYPE_t * data_p, DTYPE_t * forest_p, # The 1D indices are "raveled" or "linear", that's where "rindex" comes from. -cdef void scan1D(DTYPE_t * data_p, DTYPE_t * forest_p, shape_info * shapeinfo, - bginfo * bg, DTYPE_t neighbors, DTYPE_t y, DTYPE_t z): +cdef void scan1D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, + bginfo *bg, DTYPE_t neighbors, DTYPE_t y, DTYPE_t z): """ Perform forward scan on a 1D object, usually the first row of an image """ # Initialize the first row cdef DTYPE_t x, rindex - cdef INTS_t * DEX = shapeinfo.DEX + cdef INTS_t *DEX = shapeinfo.DEX rindex = shapeinfo.ravel_index(0, y, z, shapeinfo) if data_p[rindex] == bg.background_val: @@ -429,13 +436,13 @@ cdef void scan1D(DTYPE_t * data_p, DTYPE_t * forest_p, shape_info * shapeinfo, join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) -cdef void scan2D(DTYPE_t * data_p, DTYPE_t * forest_p, shape_info * shapeinfo, - bginfo * bg, DTYPE_t neighbors, DTYPE_t z): +cdef void scan2D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, + bginfo *bg, DTYPE_t neighbors, DTYPE_t z): """ Perform forward scan on a 2D array. """ cdef DTYPE_t x, y, rindex - cdef INTS_t * DEX = shapeinfo.DEX + cdef INTS_t *DEX = shapeinfo.DEX scan1D(data_p, forest_p, shapeinfo, bg, neighbors, 0, z) for y in range(1, shapeinfo.y): rindex = shapeinfo.ravel_index(0, y, 0, shapeinfo) @@ -469,13 +476,13 @@ cdef void scan2D(DTYPE_t * data_p, DTYPE_t * forest_p, shape_info * shapeinfo, join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) -cdef void scan3D(DTYPE_t * data_p, DTYPE_t * forest_p, shape_info * shapeinfo, - bginfo * bg, DTYPE_t neighbors): +cdef void scan3D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, + bginfo *bg, DTYPE_t neighbors): """ Perform forward scan on a 2D array. """ cdef DTYPE_t x, y, z, rindex - cdef INTS_t * DEX = shapeinfo.DEX + cdef INTS_t *DEX = shapeinfo.DEX # Handle first plane scan2D(data_p, forest_p, shapeinfo, bg, neighbors, 0) for z in range(1, shapeinfo.z): From 53f93a7f8e5f3687bef2b62bcb6a04071ba0d2e3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mat=C4=9Bj=20T=C3=BD=C4=8D?= Date: Sat, 8 Nov 2014 12:09:19 +0100 Subject: [PATCH 08/23] Added tests for 3D labelling (that fail) --- skimage/morphology/tests/test_ccomp.py | 153 ++++++++++++++++++++++++- 1 file changed, 151 insertions(+), 2 deletions(-) diff --git a/skimage/morphology/tests/test_ccomp.py b/skimage/morphology/tests/test_ccomp.py index fb596481..c7c52d96 100644 --- a/skimage/morphology/tests/test_ccomp.py +++ b/skimage/morphology/tests/test_ccomp.py @@ -6,6 +6,7 @@ from warnings import catch_warnings from skimage._shared.utils import skimage_deprecation np.random.seed(0) +BGL = -1 class TestConnectedComponents: @@ -78,7 +79,7 @@ class TestConnectedComponents: assert_array_equal(label(x, background=0), [[-1, -1, 0], [-1, -1, 0], - [ 1, 1, 1]]) + [+1, 1, 1]]) def test_background_one_region_center(self): x = np.array([[0, 0, 0], @@ -101,5 +102,153 @@ class TestConnectedComponents: assert_array_equal(label(x, background=0, return_num=True)[1], 3) +class TestConnectedComponents3d: + def setup(self): + self.x = np.zeros((3, 4, 5), int) + self.x[0] = np.array([[0, 3, 2, 1, 9], + [0, 1, 9, 2, 9], + [0, 1, 9, 9, 9], + [3, 1, 5, 3, 0]]) + + self.x[1] = np.array([[3, 3, 2, 1, 9], + [0, 3, 9, 2, 1], + [0, 3, 3, 1, 1], + [3, 1, 3, 3, 0]]) + + self.x[2] = np.array([[3, 3, 8, 8, 0], + [2, 3, 9, 8, 8], + [2, 3, 0, 8, 0], + [2, 1, 0, 0, 0]]) + + self.labels = np.zeros((3, 4, 5), int) + + self.labels[0] = np.array([[0, 1, 2, 3, 4], + [0, 5, 4, 2, 4], + [0, 5, 4, 4, 4], + [1, 5, 6, 1, 7]]) + + self.labels[1] = np.array([[1, 1, 2, 3, 4], + [0, 1, 4, 2, 3], + [0, 1, 1, 3, 3], + [1, 5, 1, 1, 7]]) + + self.labels[2] = np.array([[1, 1, 8, 8, 10], + [9, 1, 4, 8, 8], + [9, 1, 7, 8, 7], + [9, 5, 7, 7, 7]]) + + def test_basic(self): + labels = label(self.x) + assert_array_equal(labels, self.labels) + + assert self.x[0, 0, 2] == 3, \ + "Data was modified!" + + def test_random(self): + x = (np.random.rand(20, 30) * 5).astype(np.int) + + with catch_warnings(): + labels = label(x) + + n = labels.max() + for i in range(n): + values = x[labels == i] + assert np.all(values == values[0]) + + def test_diag(self): + x = np.zeros((3, 3, 3), int) + x[0, 2, 2] = 1 + x[1, 1, 1] = 1 + x[2, 0, 0] = 1 + with catch_warnings(): + assert_array_equal(label(x), x) + + def test_4_vs_8(self): + x = np.zeros((2, 2, 2), int) + x[0, 1, 1] = 1 + x[1, 0, 0] = 1 + label4 = x.copy() + label4[1, 0, 0] = 2 + with catch_warnings(): + assert_array_equal(label(x, 4), label4) + assert_array_equal(label(x, 8), x) + + def test_background(self): + x = np.zeros((2, 3, 3), int) + x[0] = np.array([[1, 0, 0], + [1, 0, 0], + [0, 0, 0]]) + x[1] = np.array([[0, 0, 0], + [0, 1, 5], + [0, 0, 0]]) + + lnb = x.copy() + lnb[0] = np.array([[0, 1, 1], + [0, 1, 1], + [1, 1, 1]]) + lnb[1] = np.array([[1, 1, 1], + [1, 0, 2], + [1, 1, 1]]) + lb = x.copy() + lb[0] = np.array([[0, BGL, BGL], + [0, BGL, BGL], + [BGL, BGL, BGL]]) + lb[1] = np.array([[BGL, BGL, BGL], + [BGL, 0, 1], + [BGL, BGL, BGL]]) + + with catch_warnings(): + assert_array_equal(label(x), lnb) + + assert_array_equal(label(x, background=0), lb) + + def test_background_two_regions(self): + x = np.zeros((2, 3, 3), int) + x[0] = np.array([[0, 0, 6], + [0, 0, 6], + [5, 5, 5]]) + x[1] = np.array([[6, 6, 0], + [5, 0, 0], + [0, 0, 0]]) + lb = x.copy() + lb[0] = np.array([[BGL, BGL, 0], + [BGL, BGL, 0], + [1, 1, 1]]) + lb[1] = np.array([[0, 0, BGL], + [1, BGL, BGL], + [BGL, BGL, BGL]]) + + assert_array_equal(label(x, background=0), lb) + + def test_background_one_region_center(self): + x = np.zeros((3, 3, 3), int) + x[1, 1, 1] = 1 + + lb = np.ones_like(x) * BGL + lb[1, 1, 1] = 0 + + assert_array_equal(label(x, neighbors=4, background=0), lb) + + def test_return_num(self): + x = np.array([[1, 0, 6], + [0, 0, 6], + [5, 5, 5]]) + + with catch_warnings(): + assert_array_equal(label(x, return_num=True)[1], 4) + + assert_array_equal(label(x, background=0, return_num=True)[1], 3) + + if __name__ == "__main__": - run_module_suite() + #run_module_suite() + lol = TestConnectedComponents3d() + lol.setup() + lol.test_basic() # 1 failiure + lol.test_random() + #lol.test_diag() # epic failiure + #lol.test_4_vs_8() # label8 epic + #lol.test_background() + #lol.test_background_two_regions() + lol.test_background_one_region_center() + lol.test_return_num() From 438a9355d6aa0085a9ed7d12afee7c8b2d92351f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mat=C4=9Bj=20T=C3=BD=C4=8D?= Date: Tue, 11 Nov 2014 00:39:49 +0100 Subject: [PATCH 09/23] Fixed bugs so test pass now --- skimage/measure/_ccomp.pyx | 16 ++++++++------- skimage/morphology/tests/test_ccomp.py | 28 +++++++++----------------- 2 files changed, 19 insertions(+), 25 deletions(-) diff --git a/skimage/measure/_ccomp.pyx b/skimage/measure/_ccomp.pyx index cc4a5320..5344aa9a 100644 --- a/skimage/measure/_ccomp.pyx +++ b/skimage/measure/_ccomp.pyx @@ -137,7 +137,7 @@ cdef shape_info get_shape_info(inarr_shape): res.DEX[D_ei] = res.DEX[D_ef] - res.DEX[D_eb] # DEX[D_eb] = one row up, remember? res.DEX[D_ej] = res.DEX[D_ei] + 1 res.DEX[D_ek] = res.DEX[D_ei] + 2 - res.DEX[D_el] = res.DEX[D_ei] - 2 * res.DEX[D_eb] + res.DEX[D_el] = res.DEX[D_ei] - res.DEX[D_eb] res.DEX[D_em] = res.DEX[D_el] + 1 res.DEX[D_en] = res.DEX[D_el] + 2 @@ -470,7 +470,7 @@ cdef void scan2D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) if neighbors == 8: - if x < shapeinfo.x - 1: + if x + 1 < shapeinfo.x: join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) @@ -557,7 +557,7 @@ cdef void scan3D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) if neighbors == 8: - if x < shapeinfo.x - 1: + if x + 1 < shapeinfo.x: join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) @@ -579,9 +579,11 @@ cdef void scan3D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, if x + 1 < shapeinfo.x: join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_el]) + if y + 1 < shapeinfo.y: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_el]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) - if x + 1 < shapeinfo.x: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_en]) + if x + 1 < shapeinfo.x: + join_trees_wrapper(data_p, forest_p, + rindex, DEX[D_en]) diff --git a/skimage/morphology/tests/test_ccomp.py b/skimage/morphology/tests/test_ccomp.py index c7c52d96..490c32d3 100644 --- a/skimage/morphology/tests/test_ccomp.py +++ b/skimage/morphology/tests/test_ccomp.py @@ -76,7 +76,8 @@ class TestConnectedComponents: [0, 0, 6], [5, 5, 5]]) - assert_array_equal(label(x, background=0), + res = label(x, background=0) + assert_array_equal(res, [[-1, -1, 0], [-1, -1, 0], [+1, 1, 1]]) @@ -132,16 +133,16 @@ class TestConnectedComponents3d: [0, 1, 1, 3, 3], [1, 5, 1, 1, 7]]) - self.labels[2] = np.array([[1, 1, 8, 8, 10], - [9, 1, 4, 8, 8], - [9, 1, 7, 8, 7], - [9, 5, 7, 7, 7]]) + self.labels[2] = np.array([[1, 1, 8, 8, 9], + [10, 1, 4, 8, 8], + [10, 1, 7, 8, 7], + [10, 5, 7, 7, 7]]) def test_basic(self): labels = label(self.x) assert_array_equal(labels, self.labels) - assert self.x[0, 0, 2] == 3, \ + assert self.x[0, 0, 2] == 2, \ "Data was modified!" def test_random(self): @@ -218,7 +219,8 @@ class TestConnectedComponents3d: [1, BGL, BGL], [BGL, BGL, BGL]]) - assert_array_equal(label(x, background=0), lb) + res = label(x, background=0) + assert_array_equal(res, lb) def test_background_one_region_center(self): x = np.zeros((3, 3, 3), int) @@ -241,14 +243,4 @@ class TestConnectedComponents3d: if __name__ == "__main__": - #run_module_suite() - lol = TestConnectedComponents3d() - lol.setup() - lol.test_basic() # 1 failiure - lol.test_random() - #lol.test_diag() # epic failiure - #lol.test_4_vs_8() # label8 epic - #lol.test_background() - #lol.test_background_two_regions() - lol.test_background_one_region_center() - lol.test_return_num() + run_module_suite() From 7c7ecc4bdc3a87a4c600cb663300d84f44f9359c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mat=C4=9Bj=20T=C3=BD=C4=8D?= Date: Tue, 11 Nov 2014 21:49:57 +0100 Subject: [PATCH 10/23] Mainly typo fixes --- skimage/measure/_ccomp.pyx | 30 ++++++++++++++----------- skimage/morphology/tests/test_ccomp.py | 31 ++++++++++++++------------ 2 files changed, 34 insertions(+), 27 deletions(-) diff --git a/skimage/measure/_ccomp.pyx b/skimage/measure/_ccomp.pyx index 5344aa9a..98502579 100644 --- a/skimage/measure/_ccomp.pyx +++ b/skimage/measure/_ccomp.pyx @@ -42,9 +42,20 @@ ctypedef struct bginfo: # A pixel has neighbors that have already been scanned. -# In the paper, the pixel is denoted by E and its neighbors -# by A, B, C and D (in 2D) - see doc for function get_shape_info() +# In the paper, the pixel is denoted by E and its neighbors: +# +# z=1 z=0 x +# ----------------------> +# | A B C F G H +# | D E . I J K +# | . . . L M N +# | +# y V +# +# D_ea represents offset of A from E etc. - see the definition of +# get_shape_info cdef enum: + # the 0D neighbor # D_ee, # We don't need D_ee # the 1D neighbor D_ed, @@ -106,15 +117,8 @@ cdef shape_info get_shape_info(inarr_shape): res.numels = res.x * res.y * res.z - # Our point of interest is E. - # z=1 z=0 x - # ----------------------> - # | A B C F G H - # | D E . I J K - # | . . . L M N - # | - # y V - # + # When reading this for the first time, look at the diagram by the enum + # definition above (keyword D_ee) # Difference between E and G is (x=0, y=-1, z=-1), E and A (-1, -1, 0) etc. # Here, it is recalculated to linear (raveled) indices of flattened arrays # with their last (=contiguous) dimension is x. @@ -278,7 +282,7 @@ def label(input, DTYPE_t neighbors=8, background=None, return_num=False): Image to label. neighbors : {4, 8}, int, optional Whether to use 4- or 8-connectivity. - In 3D, 4-connectivity means connected pixels share have to share face, + In 3D, 4-connectivity means connected pixels have to share face, whereas with 8-connectivity, they have to share only edge or vertex. background : int, optional Consider all pixels with this value as background pixels, and label @@ -389,7 +393,7 @@ cdef DTYPE_t resolve_labels(DTYPE_t *data_p, DTYPE_t *forest_p, for i in range(shapeinfo.numels): if i == bg.background_node: - data_p[i] = -1 + data_p[i] = bg.background_val elif i == forest_p[i]: # We have stumbled across a root which is something new to us (root # is the LOWEST of all prov. labels that are equivalent to it) diff --git a/skimage/morphology/tests/test_ccomp.py b/skimage/morphology/tests/test_ccomp.py index 490c32d3..b33bec4e 100644 --- a/skimage/morphology/tests/test_ccomp.py +++ b/skimage/morphology/tests/test_ccomp.py @@ -6,7 +6,10 @@ from warnings import catch_warnings from skimage._shared.utils import skimage_deprecation np.random.seed(0) -BGL = -1 + +# The background label value +# is supposed to be changed to 0 soon +BG = -1 class TestConnectedComponents: @@ -191,12 +194,12 @@ class TestConnectedComponents3d: [1, 0, 2], [1, 1, 1]]) lb = x.copy() - lb[0] = np.array([[0, BGL, BGL], - [0, BGL, BGL], - [BGL, BGL, BGL]]) - lb[1] = np.array([[BGL, BGL, BGL], - [BGL, 0, 1], - [BGL, BGL, BGL]]) + lb[0] = np.array([[0, BG, BG], + [0, BG, BG], + [BG, BG, BG]]) + lb[1] = np.array([[BG, BG, BG], + [BG, 0, 1], + [BG, BG, BG]]) with catch_warnings(): assert_array_equal(label(x), lnb) @@ -212,12 +215,12 @@ class TestConnectedComponents3d: [5, 0, 0], [0, 0, 0]]) lb = x.copy() - lb[0] = np.array([[BGL, BGL, 0], - [BGL, BGL, 0], - [1, 1, 1]]) - lb[1] = np.array([[0, 0, BGL], - [1, BGL, BGL], - [BGL, BGL, BGL]]) + lb[0] = np.array([[BG, BG, 0], + [BG, BG, 0], + [1, 1, 1]]) + lb[1] = np.array([[0, 0, BG], + [1, BG, BG], + [BG, BG, BG]]) res = label(x, background=0) assert_array_equal(res, lb) @@ -226,7 +229,7 @@ class TestConnectedComponents3d: x = np.zeros((3, 3, 3), int) x[1, 1, 1] = 1 - lb = np.ones_like(x) * BGL + lb = np.ones_like(x) * BG lb[1, 1, 1] = 0 assert_array_equal(label(x, neighbors=4, background=0), lb) From 67464e20816e942c2fee9dc581344c5b7f2de8f3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mat=C4=9Bj=20T=C3=BD=C4=8D?= Date: Sat, 15 Nov 2014 16:11:31 +0100 Subject: [PATCH 11/23] Fixed a typo + cleaned the code semantics --- skimage/measure/_ccomp.pyx | 38 +++++++++++++++++++++++++------------- 1 file changed, 25 insertions(+), 13 deletions(-) diff --git a/skimage/measure/_ccomp.pyx b/skimage/measure/_ccomp.pyx index 98502579..91475036 100644 --- a/skimage/measure/_ccomp.pyx +++ b/skimage/measure/_ccomp.pyx @@ -37,8 +37,30 @@ ctypedef int (* fun_ravel)(int, int, int, shape_info *) # For having stuff concerning background in one place ctypedef struct bginfo: + ## The value in the image (i.e. not the label!) that identifies + ## the background. DTYPE_t background_val DTYPE_t background_node + ## Identification of the background in the labelled image + DTYPE_t background_label + + +cdef bginfo get_bginfo(background_val): + cdef bginfo ret + + if background_val is None: + warnings.warn(DeprecationWarning( + 'The default value for `background` will change to 0 in v0.12' + )) + ret.background_val = -1 + else: + ret.background_val = background_val + + # The node -999 doesn't exist, it will get subsituted by a meaningful value + # upon the first background pixel occurence + ret.background_node = -999 + ret.background_label = -1 + return ret # A pixel has neighbors that have already been scanned. @@ -344,17 +366,7 @@ def label(input, DTYPE_t neighbors=8, background=None, return_num=False): cdef bginfo bg shapeinfo = get_shape_info(input.shape) - - bg.background_val = 0 - bg.background_node = -999 - - if background is None: - bg.background_val = -1 - warnings.warn(DeprecationWarning( - 'The default value for `background` will change to 0 in v0.12' - )) - else: - bg.background_val = background + bg = get_bginfo(background) if neighbors != 4 and neighbors != 8: raise ValueError('Neighbors must be either 4 or 8.') @@ -389,11 +401,11 @@ cdef DTYPE_t resolve_labels(DTYPE_t *data_p, DTYPE_t *forest_p, our knowledge of prov. labels relationship. We also track how many distinct final labels we have. """ - cdef DTYPE_t counter = 0, i + cdef DTYPE_t counter = bg.background_label + 1, i for i in range(shapeinfo.numels): if i == bg.background_node: - data_p[i] = bg.background_val + data_p[i] = bg.background_label elif i == forest_p[i]: # We have stumbled across a root which is something new to us (root # is the LOWEST of all prov. labels that are equivalent to it) From d29831882daef9b7d7691becc3669349d642b547 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mat=C4=9Bj=20T=C3=BD=C4=8D?= Date: Sat, 15 Nov 2014 16:42:22 +0100 Subject: [PATCH 12/23] Redone the background handling = considerable speedup --- skimage/measure/_ccomp.pyx | 138 ++++++++++++++++++++----------------- 1 file changed, 76 insertions(+), 62 deletions(-) diff --git a/skimage/measure/_ccomp.pyx b/skimage/measure/_ccomp.pyx index 91475036..dd7143a2 100644 --- a/skimage/measure/_ccomp.pyx +++ b/skimage/measure/_ccomp.pyx @@ -272,17 +272,6 @@ cdef inline void join_trees(DTYPE_t *forest, DTYPE_t n, DTYPE_t m): set_root(forest, m, root) -cdef inline void link_bg(DTYPE_t *forest, DTYPE_t n, DTYPE_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(input, DTYPE_t neighbors=8, background=None, return_num=False): """Label connected regions of an integer array. @@ -371,6 +360,7 @@ def label(input, DTYPE_t neighbors=8, background=None, return_num=False): if neighbors != 4 and neighbors != 8: raise ValueError('Neighbors must be either 4 or 8.') + scanBG(data_p, forest_p, & shapeinfo, & bg) if shapeinfo.ndim == 1: scan1D(data_p, forest_p, & shapeinfo, & bg, neighbors, 0, 0) elif shapeinfo.ndim == 2: @@ -416,6 +406,30 @@ cdef DTYPE_t resolve_labels(DTYPE_t *data_p, DTYPE_t *forest_p, return counter +cdef void scanBG(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, + bginfo *bg): + """ + Settle all background pixels now and don't bother with them later. + Since this only requires one linar sweep through the array, it is fast + and it makes sense to do it separately. + + The result of this function is update of forest_p and bg parameter. + """ + cdef DTYPE_t i, bgval = bg.background_val, firstbg + # We find the provisional label of the background, which is the index of + # the first background pixel + for i in range(shapeinfo.numels): + if data_p[i] == bgval: + firstbg = i + bg.background_node = firstbg + break + + # And then we apply this provisional label to the whole background + for i in range(firstbg, shapeinfo.numels): + if data_p[i] == bgval: + forest_p[i] = firstbg + + # Here, we work with flat arrays regardless whether the data is 1, 2 or 3D. # The lookup to the neighbor in a 2D array is achieved by precalculating an # offset and adding it to the index. @@ -435,19 +449,17 @@ cdef void scan1D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, Perform forward scan on a 1D object, usually the first row of an image """ # Initialize the first row - cdef DTYPE_t x, rindex + cdef DTYPE_t x, rindex, bgval = bg.background_val cdef INTS_t *DEX = shapeinfo.DEX rindex = shapeinfo.ravel_index(0, y, z, shapeinfo) - if data_p[rindex] == bg.background_val: - link_bg(forest_p, rindex, & bg.background_node) - for x in range(1, shapeinfo.x): rindex += 1 # Handle the first row # First row => rindex == j - if data_p[rindex] == bg.background_val: - link_bg(forest_p, rindex, & bg.background_node) + if data_p[rindex] == bgval: + # Nothing to do if we are background + continue join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) @@ -457,19 +469,19 @@ cdef void scan2D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, """ Perform forward scan on a 2D array. """ - cdef DTYPE_t x, y, rindex + cdef DTYPE_t x, y, rindex, bgval = bg.background_val cdef INTS_t *DEX = shapeinfo.DEX scan1D(data_p, forest_p, shapeinfo, bg, neighbors, 0, z) for y in range(1, shapeinfo.y): rindex = shapeinfo.ravel_index(0, y, 0, shapeinfo) # Handle the first column - if data_p[rindex] == bg.background_val: - link_bg(forest_p, rindex, & bg.background_node) + if data_p[rindex] != bgval: + # Nothing to do if we are background - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) - if neighbors == 8: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) + if neighbors == 8: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) # Handle the rest of columns for x in range(1, shapeinfo.x): @@ -477,8 +489,9 @@ cdef void scan2D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, # so we increment the raveled index. It will be reset when we get # to another row, so we don't have to worry about altering it here. rindex += 1 - if data_p[rindex] == bg.background_val: - link_bg(forest_p, rindex, & bg.background_node) + if data_p[rindex] == bgval: + # Nothing to do if we are background + continue if neighbors == 8: join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ea]) @@ -497,75 +510,76 @@ cdef void scan3D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, """ Perform forward scan on a 2D array. """ - cdef DTYPE_t x, y, z, rindex + cdef DTYPE_t x, y, z, rindex, bgval = bg.background_val cdef INTS_t *DEX = shapeinfo.DEX # Handle first plane scan2D(data_p, forest_p, shapeinfo, bg, neighbors, 0) for z in range(1, shapeinfo.z): # Handle first row in 3D manner rindex = shapeinfo.ravel_index(0, 0, z, shapeinfo) - if data_p[rindex] == bg.background_val: - link_bg(forest_p, rindex, & bg.background_node) + if data_p[rindex] != bgval: + # Nothing to do if we are background - # Now we have pixels below - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) + # Now we have pixels below + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) - if neighbors == 8: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_en]) + if neighbors == 8: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_en]) for x in range(1, shapeinfo.x): rindex += 1 # Handle the first row - if data_p[rindex] == bg.background_val: - link_bg(forest_p, rindex, & bg.background_node) + if data_p[rindex] != bgval: + # Nothing to do if we are background - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) - if neighbors == 8: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) + if neighbors == 8: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) - if neighbors == 8: - if x + 1 < shapeinfo.x: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) + if neighbors == 8: + if x + 1 < shapeinfo.x: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_el]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_el]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) - if x + 1 < shapeinfo.x: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_en]) + if x + 1 < shapeinfo.x: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_en]) for y in range(1, shapeinfo.y): rindex = shapeinfo.ravel_index(0, y, z, shapeinfo) # Handle the first column in 3D manner - if data_p[rindex] == bg.background_val: - link_bg(forest_p, rindex, & bg.background_node) + if data_p[rindex] != bgval: + # Nothing to do if we are background - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) - if neighbors == 8: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eh]) + if neighbors == 8: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eh]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) - if neighbors == 8: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) + if neighbors == 8: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) - if y + 1 < shapeinfo.y: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_en]) + if y + 1 < shapeinfo.y: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_en]) # Handle the rest of columns for x in range(1, shapeinfo.x): rindex += 1 - if data_p[rindex] == bg.background_val: - link_bg(forest_p, rindex, & bg.background_node) + if data_p[rindex] == bgval: + # Nothing to do if we are background + continue if neighbors == 8: join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ea]) From eb55d177c3ad136836fc91b921cba2588c1ee18b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mat=C4=9Bj=20T=C3=BD=C4=8D?= Date: Sat, 15 Nov 2014 17:11:45 +0100 Subject: [PATCH 13/23] Removed declaration of deceased function --- skimage/measure/_ccomp.pxd | 1 - 1 file changed, 1 deletion(-) diff --git a/skimage/measure/_ccomp.pxd b/skimage/measure/_ccomp.pxd index 62d21fec..98ecf15a 100644 --- a/skimage/measure/_ccomp.pxd +++ b/skimage/measure/_ccomp.pxd @@ -7,4 +7,3 @@ ctypedef cnp.intp_t DTYPE_t cdef DTYPE_t find_root(DTYPE_t *forest, DTYPE_t n) cdef void set_root(DTYPE_t *forest, DTYPE_t n, DTYPE_t root) cdef void join_trees(DTYPE_t *forest, DTYPE_t n, DTYPE_t m) -cdef void link_bg(DTYPE_t *forest, DTYPE_t n, DTYPE_t *background_node) From abf3dbc593782db05f2bfd55fe408518f72de7d0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mat=C4=9Bj=20T=C3=BD=C4=8D?= Date: Sun, 16 Nov 2014 15:33:13 +0100 Subject: [PATCH 14/23] Added and tested proper exception handling concerning array shape --- skimage/measure/_ccomp.pyx | 30 +++++++++----------------- skimage/morphology/tests/test_ccomp.py | 4 ++++ 2 files changed, 14 insertions(+), 20 deletions(-) diff --git a/skimage/measure/_ccomp.pyx b/skimage/measure/_ccomp.pyx index dd7143a2..2a7ff72d 100644 --- a/skimage/measure/_ccomp.pyx +++ b/skimage/measure/_ccomp.pyx @@ -45,9 +45,7 @@ ctypedef struct bginfo: DTYPE_t background_label -cdef bginfo get_bginfo(background_val): - cdef bginfo ret - +cdef void get_bginfo(background_val, bginfo *ret) except *: if background_val is None: warnings.warn(DeprecationWarning( 'The default value for `background` will change to 0 in v0.12' @@ -60,7 +58,6 @@ cdef bginfo get_bginfo(background_val): # upon the first background pixel occurence ret.background_node = -999 ret.background_label = -1 - return ret # A pixel has neighbors that have already been scanned. @@ -110,13 +107,11 @@ cdef struct s_shpinfo: fun_ravel ravel_index -cdef shape_info get_shape_info(inarr_shape): +cdef void get_shape_info(inarr_shape, shape_info *res) except *: """ Precalculates all the needed data from the input array shape and stores them in the shape_info struct. """ - cdef shape_info res - res.y = 1 res.z = 1 res.ravel_index = ravel_index2D @@ -135,7 +130,9 @@ cdef shape_info get_shape_info(inarr_shape): res.z = inarr_shape[0] res.ravel_index = ravel_index3D else: - assert "Only for images of dimension 1-3 (got %s)" % res.ndim + raise NotImplementedError( + "Only for images of dimension 1-3 are supported, got a %sD one" + % res.ndim) res.numels = res.x * res.y * res.z @@ -152,12 +149,12 @@ cdef shape_info get_shape_info(inarr_shape): # res.DEX[D_ee] = 0 # So now the 2nd (needed for 2D and 3D) part, y = 0, z = 1 - res.DEX[D_ea] = res.ravel_index(-1, -1, 0, & res) + res.DEX[D_ea] = res.ravel_index(-1, -1, 0, res) res.DEX[D_eb] = res.DEX[D_ea] + 1 res.DEX[D_ec] = res.DEX[D_eb] + 1 # And now the 3rd (needed only for 3D) part, z = 0 - res.DEX[D_ef] = res.ravel_index(-1, -1, -1, & res) + res.DEX[D_ef] = res.ravel_index(-1, -1, -1, res) res.DEX[D_eg] = res.DEX[D_ef] + 1 res.DEX[D_eh] = res.DEX[D_ef] + 2 res.DEX[D_ei] = res.DEX[D_ef] - res.DEX[D_eb] # DEX[D_eb] = one row up, remember? @@ -167,8 +164,6 @@ cdef shape_info get_shape_info(inarr_shape): res.DEX[D_em] = res.DEX[D_el] + 1 res.DEX[D_en] = res.DEX[D_el] + 2 - return res - cdef inline void join_trees_wrapper(DTYPE_t *data_p, DTYPE_t *forest_p, DTYPE_t rindex, INTS_t idxdiff): @@ -354,19 +349,14 @@ def label(input, DTYPE_t neighbors=8, background=None, return_num=False): cdef shape_info shapeinfo cdef bginfo bg - shapeinfo = get_shape_info(input.shape) - bg = get_bginfo(background) + get_shape_info(input.shape, &shapeinfo) + get_bginfo(background, &bg) if neighbors != 4 and neighbors != 8: raise ValueError('Neighbors must be either 4 or 8.') scanBG(data_p, forest_p, & shapeinfo, & bg) - if shapeinfo.ndim == 1: - scan1D(data_p, forest_p, & shapeinfo, & bg, neighbors, 0, 0) - elif shapeinfo.ndim == 2: - scan2D(data_p, forest_p, & shapeinfo, & bg, neighbors, 0) - elif shapeinfo.ndim == 3: - scan3D(data_p, forest_p, & shapeinfo, & bg, neighbors) + scan3D(data_p, forest_p, & shapeinfo, & bg, neighbors) # Label output cdef DTYPE_t ctr diff --git a/skimage/morphology/tests/test_ccomp.py b/skimage/morphology/tests/test_ccomp.py index b33bec4e..ea3124a0 100644 --- a/skimage/morphology/tests/test_ccomp.py +++ b/skimage/morphology/tests/test_ccomp.py @@ -244,6 +244,10 @@ class TestConnectedComponents3d: assert_array_equal(label(x, background=0, return_num=True)[1], 3) + def test_nd(self): + x = np.ones((1, 2, 3, 4)) + np.testing.assert_raises(NotImplementedError, label, x) + if __name__ == "__main__": run_module_suite() From 06c46ec51e0c193c9b9c7400a369a420c3bfcb56 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mat=C4=9Bj=20T=C3=BD=C4=8D?= Date: Tue, 18 Nov 2014 23:45:54 +0100 Subject: [PATCH 15/23] Introduced the connectivity parameter, regrouped some code --- skimage/measure/_ccomp.pyx | 186 ++++++++++++++++++++++--------------- 1 file changed, 111 insertions(+), 75 deletions(-) diff --git a/skimage/measure/_ccomp.pyx b/skimage/measure/_ccomp.pyx index 2a7ff72d..d471ef16 100644 --- a/skimage/measure/_ccomp.pyx +++ b/skimage/measure/_ccomp.pyx @@ -267,19 +267,49 @@ cdef inline void join_trees(DTYPE_t *forest, DTYPE_t n, DTYPE_t m): set_root(forest, m, root) +def _norm_connectivity(connectivity, ndim): + """ + Takes the value of the connectivity parameter, validates it and converts + it to a value that the subsequent algorithm may use as-is safely. + + Parameters + ---------- + connectivity : int + The following should be true: -ndim < connectivity < 0, or + 0 < connectivity <= ndim. + + Returns + ------- + connectivity : int + Connectivity, 0 < connectivity < ndim + """ + if connectivity == 0 or connectivity > ndim: + raise ValueError( + "Connectivity of 0 or above %d doesn't make sense" % ndim) + res = connectivity + if res < 0: + res = res % ndim + # we want -1 to be normed to ndim, -2 to (ndim - 1) etc. + res += 1 + return res + + # Connected components search as described in Fiorio et al. -def label(input, DTYPE_t neighbors=8, background=None, return_num=False): +def label(input, DTYPE_t neighbors=8, background=None, return_num=False, + connectivity=None): """Label connected regions of an integer array. 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:: + In 2D, they can be neighbors either in a 1- or 2-connected sense. + The value refers to the greatest number of orthogonal hops between the + starting point and the neighbor. - 4-connectivity 8-connectivity + 1-connectivity 2-connectivity - [ ] [ ] [ ] [ ] - | \ | / - [ ]--[ ]--[ ] [ ]--[ ]--[ ] - | / | \\ + [ ] [ ] [ ] [ ] [ ] + | \\ | / | <- hop 2 + [ ]--[x]--[ ] [ ]--[x]--[ ] [x]--[ ] + | / | \\ hop 1 [ ] [ ] [ ] [ ] Parameters @@ -290,12 +320,21 @@ def label(input, DTYPE_t neighbors=8, background=None, return_num=False): Whether to use 4- or 8-connectivity. In 3D, 4-connectivity means connected pixels have to share face, whereas with 8-connectivity, they have to share only edge or vertex. + **Depreceated, use ``connectivity`` instead.** background : int, optional Consider all pixels with this value as background pixels, and label them as -1. (Note: background pixels will be labeled as 0 starting with version 0.12). return_num : bool, optional Whether to return the number of assigned labels. + connectivity : int + Number of orthogonal hops + For the 2D case, 1 considers horizontal and vertical neighbors, whereas + 2 adds the diagonals (you hop once vertically and once horizontally). + 1 is the lowest value of connection (4 neighbors in 2D, 6 in 3D). + Moreover, the value of -1 specifies the highest connectivity available. + So for example in 2D, -1 is equivalent of 2, resulting in considering + all 8 neighbors. Returns ------- @@ -352,15 +391,32 @@ def label(input, DTYPE_t neighbors=8, background=None, return_num=False): get_shape_info(input.shape, &shapeinfo) get_bginfo(background, &bg) - if neighbors != 4 and neighbors != 8: - raise ValueError('Neighbors must be either 4 or 8.') + if neighbors is None and connectivity is None: + # Pure default + connectivity = -1 + elif neighbors is not None: + # Pure fail + if neighbors != 4 and neighbors != 8: + msg = "Neighbors must be either 4 or 8, got '%d'.\n" % neighbors + msg += "Moreover, this arg is depreceated, use 'connectivity' instead" + raise ValueError(msg) + else: + # backwards-compatible neighbors recalc to connectivity, + # depreciation warning + nei2conn = {4: 1, 8: -1} + connectivity = nei2conn[neighbors] + msg = "Argument 'neighbors' is depreceated, use 'connectivity' " + msg += "instead. Its coresponing value is likely '%d'" % connectivity + DeprecationWarning(msg) - scanBG(data_p, forest_p, & shapeinfo, & bg) - scan3D(data_p, forest_p, & shapeinfo, & bg, neighbors) + connectivity = _norm_connectivity(connectivity, shapeinfo.ndim) + + scanBG(data_p, forest_p, &shapeinfo, &bg) + scan3D(data_p, forest_p, &shapeinfo, &bg, connectivity) # Label output cdef DTYPE_t ctr - ctr = resolve_labels(data_p, forest_p, & shapeinfo, & bg) + ctr = resolve_labels(data_p, forest_p, &shapeinfo, &bg) # Work around a bug in ndimage's type checking on 32-bit platforms if data.dtype == np.int32: @@ -429,12 +485,12 @@ cdef void scanBG(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, # D E # # So if I am in the point E and want to take a look to A, I take the index of -# E and add shapeinfo.Dea to it and teg the index of A. +# E and add shapeinfo.DEX[D_ea] to it and get the index of A. # The 1D indices are "raveled" or "linear", that's where "rindex" comes from. cdef void scan1D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, - bginfo *bg, DTYPE_t neighbors, DTYPE_t y, DTYPE_t z): + bginfo *bg, DTYPE_t connectivity, DTYPE_t y, DTYPE_t z): """ Perform forward scan on a 1D object, usually the first row of an image """ @@ -455,13 +511,13 @@ cdef void scan1D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, cdef void scan2D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, - bginfo *bg, DTYPE_t neighbors, DTYPE_t z): + bginfo *bg, DTYPE_t connectivity, DTYPE_t z): """ Perform forward scan on a 2D array. """ cdef DTYPE_t x, y, rindex, bgval = bg.background_val cdef INTS_t *DEX = shapeinfo.DEX - scan1D(data_p, forest_p, shapeinfo, bg, neighbors, 0, z) + scan1D(data_p, forest_p, shapeinfo, bg, connectivity, 0, z) for y in range(1, shapeinfo.y): rindex = shapeinfo.ravel_index(0, y, 0, shapeinfo) # Handle the first column @@ -470,7 +526,7 @@ cdef void scan2D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) - if neighbors == 8: + if connectivity >= 2: join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) # Handle the rest of columns @@ -483,27 +539,24 @@ cdef void scan2D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, # Nothing to do if we are background continue - if neighbors == 8: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ea]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) - if neighbors == 8: + if connectivity >= 2: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ea]) if x + 1 < shapeinfo.x: join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) - cdef void scan3D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, - bginfo *bg, DTYPE_t neighbors): + bginfo *bg, DTYPE_t connectivity): """ Perform forward scan on a 2D array. """ cdef DTYPE_t x, y, z, rindex, bgval = bg.background_val cdef INTS_t *DEX = shapeinfo.DEX # Handle first plane - scan2D(data_p, forest_p, shapeinfo, bg, neighbors, 0) + scan2D(data_p, forest_p, shapeinfo, bg, connectivity, 0) for z in range(1, shapeinfo.z): # Handle first row in 3D manner rindex = shapeinfo.ravel_index(0, 0, z, shapeinfo) @@ -513,10 +566,11 @@ cdef void scan3D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, # Now we have pixels below join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) - if neighbors == 8: + if connectivity >= 2: join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_en]) + if connectivity >= 3: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_en]) for x in range(1, shapeinfo.x): rindex += 1 @@ -525,22 +579,18 @@ cdef void scan3D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, # Nothing to do if we are background join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) - - if neighbors == 8: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) - if neighbors == 8: + if connectivity >= 2: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) if x + 1 < shapeinfo.x: join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) - - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_el]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) - - if x + 1 < shapeinfo.x: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_en]) - + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) + if connectivity >= 3: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_el]) + if x + 1 < shapeinfo.x: + join_trees_wrapper(data_p, forest_p, rindex, + DEX[D_en]) for y in range(1, shapeinfo.y): rindex = shapeinfo.ravel_index(0, y, z, shapeinfo) @@ -549,20 +599,20 @@ cdef void scan3D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, # Nothing to do if we are background join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) - - if neighbors == 8: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eh]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) - if neighbors == 8: + if connectivity >= 2: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) - if y + 1 < shapeinfo.y: join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_en]) + if connectivity >= 3: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eh]) + if y + 1 < shapeinfo.y: + join_trees_wrapper(data_p, forest_p, rindex, + DEX[D_en]) + # Handle the rest of columns for x in range(1, shapeinfo.x): @@ -571,39 +621,25 @@ cdef void scan3D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, # Nothing to do if we are background continue - if neighbors == 8: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ea]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) - - if neighbors == 8: - if x + 1 < shapeinfo.x: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) - - # Now pixels below: - if neighbors == 8: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ef]) - - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) - - if x + 1 < shapeinfo.x: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eh]) - - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) - if neighbors == 8: + if connectivity >= 2: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ea]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) if x + 1 < shapeinfo.x: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) - if y + 1 < shapeinfo.y: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_el]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) - + if connectivity >= 3: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ef]) if x + 1 < shapeinfo.x: - join_trees_wrapper(data_p, forest_p, - rindex, DEX[D_en]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eh]) + if y + 1 < shapeinfo.y: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_el]) + if x + 1 < shapeinfo.x: + join_trees_wrapper(data_p, forest_p, + rindex, DEX[D_en]) From 20755da3fea5e198878e6a71f88f03dc16122794 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mat=C4=9Bj=20T=C3=BD=C4=8D?= Date: Wed, 19 Nov 2014 21:32:27 +0100 Subject: [PATCH 16/23] Fixed typos and made some style fixes --- skimage/measure/_ccomp.pyx | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/skimage/measure/_ccomp.pyx b/skimage/measure/_ccomp.pyx index d471ef16..612616c2 100644 --- a/skimage/measure/_ccomp.pyx +++ b/skimage/measure/_ccomp.pyx @@ -295,16 +295,16 @@ def _norm_connectivity(connectivity, ndim): # Connected components search as described in Fiorio et al. -def label(input, DTYPE_t neighbors=8, background=None, return_num=False, +def label(input, neighbors=None, background=None, return_num=False, connectivity=None): """Label connected regions of an integer array. Two pixels are connected when they are neighbors and have the same value. In 2D, they can be neighbors either in a 1- or 2-connected sense. - The value refers to the greatest number of orthogonal hops between the - starting point and the neighbor. + The value refers to the maximum number of orthogonal hops to consider a + pixel/voxel a neighbor. - 1-connectivity 2-connectivity + 1-connectivity 2-connectivity diagonal connection close-up [ ] [ ] [ ] [ ] [ ] | \\ | / | <- hop 2 @@ -320,14 +320,14 @@ def label(input, DTYPE_t neighbors=8, background=None, return_num=False, Whether to use 4- or 8-connectivity. In 3D, 4-connectivity means connected pixels have to share face, whereas with 8-connectivity, they have to share only edge or vertex. - **Depreceated, use ``connectivity`` instead.** + **Deprecated, use ``connectivity`` instead.** background : int, optional Consider all pixels with this value as background pixels, and label them as -1. (Note: background pixels will be labeled as 0 starting with version 0.12). return_num : bool, optional Whether to return the number of assigned labels. - connectivity : int + connectivity : int, optional Number of orthogonal hops For the 2D case, 1 considers horizontal and vertical neighbors, whereas 2 adds the diagonals (you hop once vertically and once horizontally). @@ -392,22 +392,23 @@ def label(input, DTYPE_t neighbors=8, background=None, return_num=False, get_bginfo(background, &bg) if neighbors is None and connectivity is None: - # Pure default + # default connectivity = -1 elif neighbors is not None: - # Pure fail + depr_msg = ("The argument 'neighbors' is deprecated, use 'connectivity'" + " instead") + # fail if neighbors != 4 and neighbors != 8: + DeprecationWarning(depr_msg) msg = "Neighbors must be either 4 or 8, got '%d'.\n" % neighbors - msg += "Moreover, this arg is depreceated, use 'connectivity' instead" raise ValueError(msg) else: # backwards-compatible neighbors recalc to connectivity, - # depreciation warning + # deprecation warning nei2conn = {4: 1, 8: -1} connectivity = nei2conn[neighbors] - msg = "Argument 'neighbors' is depreceated, use 'connectivity' " - msg += "instead. Its coresponing value is likely '%d'" % connectivity - DeprecationWarning(msg) + msg = " Its corresponing value is '%d'" % connectivity + DeprecationWarning(depr_msg + msg) connectivity = _norm_connectivity(connectivity, shapeinfo.ndim) From fe4a9b58471af2a123266eecfdb063a9a5f240e7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mat=C4=9Bj=20T=C3=BD=C4=8D?= Date: Wed, 19 Nov 2014 23:02:34 +0100 Subject: [PATCH 17/23] Made the scan code cleaner and faster (but longer) --- skimage/measure/_ccomp.pyx | 178 +++++++++++++++++++++++++++++-------- 1 file changed, 139 insertions(+), 39 deletions(-) diff --git a/skimage/measure/_ccomp.pyx b/skimage/measure/_ccomp.pyx index d471ef16..685cb338 100644 --- a/skimage/measure/_ccomp.pyx +++ b/skimage/measure/_ccomp.pyx @@ -502,7 +502,6 @@ cdef void scan1D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, for x in range(1, shapeinfo.x): rindex += 1 # Handle the first row - # First row => rindex == j if data_p[rindex] == bgval: # Nothing to do if we are background continue @@ -519,6 +518,7 @@ cdef void scan2D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, cdef INTS_t *DEX = shapeinfo.DEX scan1D(data_p, forest_p, shapeinfo, bg, connectivity, 0, z) for y in range(1, shapeinfo.y): + # BEGINNING of x = 0 rindex = shapeinfo.ravel_index(0, y, 0, shapeinfo) # Handle the first column if data_p[rindex] != bgval: @@ -528,9 +528,9 @@ cdef void scan2D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, if connectivity >= 2: join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) + # END of x = 0 - # Handle the rest of columns - for x in range(1, shapeinfo.x): + for x in range(1, shapeinfo.x - 1): # We have just moved to another column (of the same row) # so we increment the raveled index. It will be reset when we get # to another row, so we don't have to worry about altering it here. @@ -544,14 +544,27 @@ cdef void scan2D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, if connectivity >= 2: join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ea]) - if x + 1 < shapeinfo.x: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) + + # Finally, the last column + # BEGINNING of x = max + rindex += 1 + if data_p[rindex] != bgval: + # Nothing to do if we are background + + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) + + if connectivity >= 2: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ea]) + # END of x = max cdef void scan3D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, bginfo *bg, DTYPE_t connectivity): """ Perform forward scan on a 2D array. + """ cdef DTYPE_t x, y, z, rindex, bgval = bg.background_val cdef INTS_t *DEX = shapeinfo.DEX @@ -559,6 +572,8 @@ cdef void scan3D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, scan2D(data_p, forest_p, shapeinfo, bg, connectivity, 0) for z in range(1, shapeinfo.z): # Handle first row in 3D manner + # BEGINNING of y = 0 + # BEGINNING of x = 0 rindex = shapeinfo.ravel_index(0, 0, z, shapeinfo) if data_p[rindex] != bgval: # Nothing to do if we are background @@ -571,28 +586,45 @@ cdef void scan3D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) if connectivity >= 3: join_trees_wrapper(data_p, forest_p, rindex, DEX[D_en]) + # END of x = 0 - for x in range(1, shapeinfo.x): + for x in range(1, shapeinfo.x - 1): rindex += 1 # Handle the first row - if data_p[rindex] != bgval: + if data_p[rindex] == bgval: # Nothing to do if we are background + continue - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) - if connectivity >= 2: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) - if x + 1 < shapeinfo.x: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) - if connectivity >= 3: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_el]) - if x + 1 < shapeinfo.x: - join_trees_wrapper(data_p, forest_p, rindex, - DEX[D_en]) + if connectivity >= 2: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) + if connectivity >= 3: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_el]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_en]) - for y in range(1, shapeinfo.y): + # BEGINNING of x = max + rindex += 1 + # Handle the last element of the first row + if data_p[rindex] != bgval: + # Nothing to do if we are background + + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) + + if connectivity >= 2: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) + if connectivity >= 3: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_el]) + # END of x = max + # END of y = 0 + + for y in range(1, shapeinfo.y - 1): + # BEGINNING of y = ... + # BEGINNING of x = 0 rindex = shapeinfo.ravel_index(0, y, z, shapeinfo) # Handle the first column in 3D manner if data_p[rindex] != bgval: @@ -605,17 +637,14 @@ cdef void scan3D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) - if y + 1 < shapeinfo.y: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) if connectivity >= 3: join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eh]) - if y + 1 < shapeinfo.y: - join_trees_wrapper(data_p, forest_p, rindex, - DEX[D_en]) - + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_en]) + # END of x = 0 # Handle the rest of columns - for x in range(1, shapeinfo.x): + for x in range(1, shapeinfo.x - 1): rindex += 1 if data_p[rindex] == bgval: # Nothing to do if we are background @@ -629,17 +658,88 @@ cdef void scan3D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ea]) join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) - if x + 1 < shapeinfo.x: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) - if y + 1 < shapeinfo.y: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) if connectivity >= 3: join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ef]) - if x + 1 < shapeinfo.x: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eh]) - if y + 1 < shapeinfo.y: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_el]) - if x + 1 < shapeinfo.x: - join_trees_wrapper(data_p, forest_p, - rindex, DEX[D_en]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eh]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_el]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_en]) + + # BEGINNING of x = max + rindex += 1 + if data_p[rindex] != bgval: + # Nothing to do if we are background + + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) + + if connectivity >= 2: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ea]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) + if connectivity >= 3: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ef]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_el]) + # END of x = max + # END of y = ... + + # BEGINNING of y = max + # BEGINNING of x = 0 + rindex = shapeinfo.ravel_index(0, shapeinfo.y, z, shapeinfo) + # Handle the first column in 3D manner + if data_p[rindex] != bgval: + # Nothing to do if we are background + + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) + + if connectivity >= 2: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) + if connectivity >= 3: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eh]) + # END of x = 0 + + # Handle the rest of columns + for x in range(1, shapeinfo.x - 1): + rindex += 1 + if data_p[rindex] == bgval: + # Nothing to do if we are background + continue + + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) + + if connectivity >= 2: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ea]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) + if connectivity >= 3: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ef]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eh]) + + # BEGINNING of x = max + rindex += 1 + if data_p[rindex] != bgval: + # Nothing to do if we are background + + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) + + if connectivity >= 2: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ea]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) + if connectivity >= 3: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ef]) + # END of x = max + # END of y = max From a4b841daf097edbf19a774374bfede1a60733d65 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mat=C4=9Bj=20T=C3=BD=C4=8D?= Date: Wed, 19 Nov 2014 23:16:33 +0100 Subject: [PATCH 18/23] Fixed bugs --- skimage/measure/_ccomp.pyx | 105 +++++++++++++++++++------------------ 1 file changed, 53 insertions(+), 52 deletions(-) diff --git a/skimage/measure/_ccomp.pyx b/skimage/measure/_ccomp.pyx index f3aacb2e..460ff45a 100644 --- a/skimage/measure/_ccomp.pyx +++ b/skimage/measure/_ccomp.pyx @@ -618,13 +618,14 @@ cdef void scan3D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, if connectivity >= 2: join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) if connectivity >= 3: join_trees_wrapper(data_p, forest_p, rindex, DEX[D_el]) # END of x = max # END of y = 0 + # BEGINNING of y = ... for y in range(1, shapeinfo.y - 1): - # BEGINNING of y = ... # BEGINNING of x = 0 rindex = shapeinfo.ravel_index(0, y, z, shapeinfo) # Handle the first column in 3D manner @@ -657,9 +658,9 @@ cdef void scan3D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, if connectivity >= 2: join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ea]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) join_trees_wrapper(data_p, forest_p, rindex, DEX[D_em]) if connectivity >= 3: @@ -686,61 +687,61 @@ cdef void scan3D(DTYPE_t *data_p, DTYPE_t *forest_p, shape_info *shapeinfo, join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ef]) join_trees_wrapper(data_p, forest_p, rindex, DEX[D_el]) # END of x = max - # END of y = ... + # END of y = ... - # BEGINNING of y = max - # BEGINNING of x = 0 - rindex = shapeinfo.ravel_index(0, shapeinfo.y, z, shapeinfo) - # Handle the first column in 3D manner - if data_p[rindex] != bgval: - # Nothing to do if we are background + # BEGINNING of y = max + # BEGINNING of x = 0 + rindex = shapeinfo.ravel_index(0, shapeinfo.y - 1, z, shapeinfo) + # Handle the first column in 3D manner + if data_p[rindex] != bgval: + # Nothing to do if we are background - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) - if connectivity >= 2: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) - if connectivity >= 3: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eh]) - # END of x = 0 + if connectivity >= 2: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) + if connectivity >= 3: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eh]) + # END of x = 0 - # Handle the rest of columns - for x in range(1, shapeinfo.x - 1): - rindex += 1 - if data_p[rindex] == bgval: - # Nothing to do if we are background - continue - - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) - - if connectivity >= 2: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ea]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) - if connectivity >= 3: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ef]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eh]) - - # BEGINNING of x = max + # Handle the rest of columns + for x in range(1, shapeinfo.x - 1): rindex += 1 - if data_p[rindex] != bgval: + if data_p[rindex] == bgval: # Nothing to do if we are background + continue - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) - if connectivity >= 2: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ea]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) - if connectivity >= 3: - join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ef]) - # END of x = max - # END of y = max + if connectivity >= 2: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ea]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ek]) + if connectivity >= 3: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ef]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eh]) + + # BEGINNING of x = max + rindex += 1 + if data_p[rindex] != bgval: + # Nothing to do if we are background + + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ed]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej]) + + if connectivity >= 2: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ea]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eg]) + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ei]) + if connectivity >= 3: + join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ef]) + # END of x = max + # END of y = max From d07caea5a843f69a35a06d886fd4aaefe0d55f71 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mat=C4=9Bj=20T=C3=BD=C4=8D?= Date: Sun, 23 Nov 2014 23:17:01 +0100 Subject: [PATCH 19/23] Mended code style --- skimage/measure/_ccomp.pyx | 50 +++++++++++++++++++------------------- 1 file changed, 25 insertions(+), 25 deletions(-) diff --git a/skimage/measure/_ccomp.pyx b/skimage/measure/_ccomp.pyx index 460ff45a..f1d38fb7 100644 --- a/skimage/measure/_ccomp.pyx +++ b/skimage/measure/_ccomp.pyx @@ -283,21 +283,23 @@ def _norm_connectivity(connectivity, ndim): connectivity : int Connectivity, 0 < connectivity < ndim """ - if connectivity == 0 or connectivity > ndim: + if connectivity == 0: raise ValueError( "Connectivity of 0 or above %d doesn't make sense" % ndim) - res = connectivity - if res < 0: - res = res % ndim - # we want -1 to be normed to ndim, -2 to (ndim - 1) etc. - res += 1 - return res + if not -ndim <= connectivity <= ndim: + raise ValueError( + "Connectivity below %d or above %d is illegal." + % (-ndim, ndim)) + if connectivity < 0: + # Just as we say in the docs + connectivity += ndim + 1 + return connectivity # Connected components search as described in Fiorio et al. def label(input, neighbors=None, background=None, return_num=False, connectivity=None): - """Label connected regions of an integer array. + r"""Label connected regions of an integer array. Two pixels are connected when they are neighbors and have the same value. In 2D, they can be neighbors either in a 1- or 2-connected sense. @@ -307,9 +309,9 @@ def label(input, neighbors=None, background=None, return_num=False, 1-connectivity 2-connectivity diagonal connection close-up [ ] [ ] [ ] [ ] [ ] - | \\ | / | <- hop 2 + | \ | / | <- hop 2 [ ]--[x]--[ ] [ ]--[x]--[ ] [x]--[ ] - | / | \\ hop 1 + | / | \ hop 1 [ ] [ ] [ ] [ ] Parameters @@ -328,13 +330,14 @@ def label(input, neighbors=None, background=None, return_num=False, return_num : bool, optional Whether to return the number of assigned labels. connectivity : int, optional - Number of orthogonal hops - For the 2D case, 1 considers horizontal and vertical neighbors, whereas - 2 adds the diagonals (you hop once vertically and once horizontally). - 1 is the lowest value of connection (4 neighbors in 2D, 6 in 3D). - Moreover, the value of -1 specifies the highest connectivity available. - So for example in 2D, -1 is equivalent of 2, resulting in considering - all 8 neighbors. + Maximum number of orthogonal hops to consider a pixel/voxel + as a neighbor. + Accepted values are ranging from 1 to input.ndim. Negative values from + -input.ndim to -1 are also accepted. Negative connectivity value + :math:`x` is equivalent to the positive connectivity value + :math:`x + input.ndim + 1`. For example in 2D, -1 and 2 = -1 + 2 + 1 + are equivalent, values. + Returns ------- @@ -353,12 +356,12 @@ def label(input, neighbors=None, background=None, return_num=False, [0 1 0] [0 0 1]] - >>> print(m.label(x, neighbors=4)) + >>> print(m.label(x, connectivity=1)) [[0 1 1] [2 3 1] [2 2 4]] - >>> print(m.label(x, neighbors=8)) + >>> print(m.label(x, connectivity=2)) [[0 1 1] [1 0 1] [1 1 0]] @@ -395,20 +398,17 @@ def label(input, neighbors=None, background=None, return_num=False, # default connectivity = -1 elif neighbors is not None: - depr_msg = ("The argument 'neighbors' is deprecated, use 'connectivity'" - " instead") + depr_msg = ("The argument 'neighbors' is deprecated, use " + "'connectivity' instead") + DeprecationWarning(depr_msg) # fail if neighbors != 4 and neighbors != 8: - DeprecationWarning(depr_msg) msg = "Neighbors must be either 4 or 8, got '%d'.\n" % neighbors raise ValueError(msg) else: # backwards-compatible neighbors recalc to connectivity, - # deprecation warning nei2conn = {4: 1, 8: -1} connectivity = nei2conn[neighbors] - msg = " Its corresponing value is '%d'" % connectivity - DeprecationWarning(depr_msg + msg) connectivity = _norm_connectivity(connectivity, shapeinfo.ndim) From ac11b61502f78c3c0bb7c34940beaf70ac40cd01 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mat=C4=9Bj=20T=C3=BD=C4=8D?= Date: Sun, 23 Nov 2014 23:45:32 +0100 Subject: [PATCH 20/23] Fixed typos, added array shape checking --- skimage/measure/_ccomp.pyx | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/skimage/measure/_ccomp.pyx b/skimage/measure/_ccomp.pyx index f1d38fb7..00c0645e 100644 --- a/skimage/measure/_ccomp.pyx +++ b/skimage/measure/_ccomp.pyx @@ -115,6 +115,11 @@ cdef void get_shape_info(inarr_shape, shape_info *res) except *: res.y = 1 res.z = 1 res.ravel_index = ravel_index2D + # A shape (3, 1, 4) would make the algorithm crash, but the corresponding + # good_shape (i.e. the array with axis swapped) (1, 3, 4) is OK. + # Having an axis length of 1 when an axis on the left is longer than 1 + # (in this case, it has length of 3) is NOT ALLOWED. + good_shape = tuple(sorted(inarr_shape)) res.ndim = len(inarr_shape) if res.ndim == 1: @@ -124,11 +129,20 @@ cdef void get_shape_info(inarr_shape, shape_info *res) except *: res.x = inarr_shape[1] res.y = inarr_shape[0] res.ravel_index = ravel_index2D + if res.x == 1 and res.y > 1: + raise ValueError( + "Swap axis of your %s array so it has a %s shape" + % (inarr_shape, good_shape)) elif res.ndim == 3: res.x = inarr_shape[2] res.y = inarr_shape[1] res.z = inarr_shape[0] res.ravel_index = ravel_index3D + if ((res.x == 1 and res.y > 1) + or res.y == 1 and res.z > 1): + raise ValueError( + "Swap axes of your %s array so it has a %s shape" + % (inarr_shape, good_shape)) else: raise NotImplementedError( "Only for images of dimension 1-3 are supported, got a %sD one" @@ -285,7 +299,7 @@ def _norm_connectivity(connectivity, ndim): """ if connectivity == 0: raise ValueError( - "Connectivity of 0 or above %d doesn't make sense" % ndim) + "Connectivity of 0 doesn't make sense") if not -ndim <= connectivity <= ndim: raise ValueError( "Connectivity below %d or above %d is illegal." From 3b4425813b03b692a8ef95a3f205f1a6af6c4091 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mat=C4=9Bj=20T=C3=BD=C4=8D?= Date: Tue, 25 Nov 2014 00:37:19 +0100 Subject: [PATCH 21/23] Streamlined connectivity vs neighbors handling --- skimage/measure/_ccomp.pyx | 23 +++++++++++------------ 1 file changed, 11 insertions(+), 12 deletions(-) diff --git a/skimage/measure/_ccomp.pyx b/skimage/measure/_ccomp.pyx index 00c0645e..e9ee2514 100644 --- a/skimage/measure/_ccomp.pyx +++ b/skimage/measure/_ccomp.pyx @@ -409,20 +409,19 @@ def label(input, neighbors=None, background=None, return_num=False, get_bginfo(background, &bg) if neighbors is None and connectivity is None: - # default - connectivity = -1 + # use the full connectivity by default + connectivity = input.ndim elif neighbors is not None: - depr_msg = ("The argument 'neighbors' is deprecated, use " - "'connectivity' instead") - DeprecationWarning(depr_msg) - # fail - if neighbors != 4 and neighbors != 8: - msg = "Neighbors must be either 4 or 8, got '%d'.\n" % neighbors - raise ValueError(msg) + DeprecationWarning("The argument 'neighbors' is deprecated, use " + "'connectivity' instead") + # backwards-compatible neighbors recalc to connectivity, + if neighbors == 4: + connectivity = 1 + elif neighbors == 8: + connectivity = input.ndim else: - # backwards-compatible neighbors recalc to connectivity, - nei2conn = {4: 1, 8: -1} - connectivity = nei2conn[neighbors] + raise ValueError("Neighbors must be either 4 or 8, got '%d'.\n" + % neighbors) connectivity = _norm_connectivity(connectivity, shapeinfo.ndim) From e41bff898a45d688b93a8b0ec870585aa379a7bc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mat=C4=9Bj=20T=C3=BD=C4=8D?= Date: Fri, 28 Nov 2014 16:45:52 +0100 Subject: [PATCH 22/23] Added fix for the strange input array shapes + tests --- skimage/measure/_ccomp.pyx | 100 ++++++++++++++++++++++++- skimage/morphology/tests/test_ccomp.py | 32 ++++++++ 2 files changed, 128 insertions(+), 4 deletions(-) diff --git a/skimage/measure/_ccomp.pyx b/skimage/measure/_ccomp.pyx index e9ee2514..466fc910 100644 --- a/skimage/measure/_ccomp.pyx +++ b/skimage/measure/_ccomp.pyx @@ -130,6 +130,7 @@ cdef void get_shape_info(inarr_shape, shape_info *res) except *: res.y = inarr_shape[0] res.ravel_index = ravel_index2D if res.x == 1 and res.y > 1: + # Should not occur, but better be safe than sorry raise ValueError( "Swap axis of your %s array so it has a %s shape" % (inarr_shape, good_shape)) @@ -140,6 +141,7 @@ cdef void get_shape_info(inarr_shape, shape_info *res) except *: res.ravel_index = ravel_index3D if ((res.x == 1 and res.y > 1) or res.y == 1 and res.z > 1): + # Should not occur, but better be safe than sorry raise ValueError( "Swap axes of your %s array so it has a %s shape" % (inarr_shape, good_shape)) @@ -310,6 +312,83 @@ def _norm_connectivity(connectivity, ndim): return connectivity +def _get_swaps(shp): + """ + What axes to swap if we want to convert an illegal array shape + to a legal one. + + Args: + shp (tuple-like): The array shape + + Returns: + list: List of tuples to be passed to np.swapaxes + """ + shp = np.array(shp) + swaps = [] + + # Dimensions where the array is "flat" + ones = np.where(shp == 1)[0][::-1] + # The other dimensions + bigs = np.where(shp > 1)[0] + # We now go from opposite directions and swap axes if an index of a flat + # axis is higher than the one of a thick axis + for one, big in zip(ones, bigs): + if one < big: + # We are ordered already + break + else: + swaps.append((one, big)) + # TODO: Add more swaps so the array is longer along x and shorter along y + return swaps + + +def _apply_swaps(arr, swaps): + arr2 = arr + for one, two in swaps: + arr2 = arr.swapaxes(one, two) + return arr2 + + +def reshape_array(arr): + """ + "Rotates" the array so it gains a shape that the algorithm can work with. + An illegal shape is (3, 1, 4), and correct one is (1, 3, 4) or (1, 4, 3). + The point is to have all 1s of the shape at the beginning, not in the + middle of the shape tuple. + + Note: The greater-than-one shape component should go from greatest to + lowest numbers since it is more friendly to the CPU cache (so (1, 3, 4) is + less optimal than (1, 4, 3). Keyword to this is "memory spatial locality" + + Args: + arr (np.ndarray): The array we want to fix + + Returns: + tuple (result, swaps): The result is the "fixed" array and a record + of what has been done with it. + """ + swaps = _get_swaps(arr.shape) + reshaped = _apply_swaps(arr, swaps) + return reshaped, swaps + + +def undo_reshape_array(arr, swaps): + """ + Does the opposite of what :func:`reshape_array` does + + Args: + arr (np.ndarray): The array to "revert" + swaps (list): The second result of :func:`reshape_array` + + Returns: + np.ndarray: The array of the original shape + """ + # It is safer to undo axes swaps in the opposite order + # than the application order + reshaped = _apply_swaps(arr, swaps[::-1]) + return reshaped + + # Connected components search as described in Fiorio et al. def label(input, neighbors=None, background=None, return_num=False, connectivity=None): @@ -390,7 +469,23 @@ def label(input, neighbors=None, background=None, return_num=False, [-1 -1 -1]] """ + # We have to ensure that the shape of the input can be handled by the + # algorithm the input if it is the case + input_corrected, swaps = reshape_array(input) + # Do the labelling + res, ctr = _label(input_corrected, neighbors, background, connectivity) + + res_orig = undo_reshape_array(res, swaps) + + if return_num: + return res_orig, ctr + else: + return res_orig + + +# Connected components search as described in Fiorio et al. +def _label(input, neighbors=None, background=None, connectivity=None): cdef cnp.ndarray[DTYPE_t, ndim=1] data cdef cnp.ndarray[DTYPE_t, ndim=1] forest @@ -438,10 +533,7 @@ def label(input, neighbors=None, background=None, return_num=False, res = data.reshape(input.shape) - if return_num: - return res, ctr - else: - return res + return res, ctr cdef DTYPE_t resolve_labels(DTYPE_t *data_p, DTYPE_t *forest_p, diff --git a/skimage/morphology/tests/test_ccomp.py b/skimage/morphology/tests/test_ccomp.py index ea3124a0..aebb443a 100644 --- a/skimage/morphology/tests/test_ccomp.py +++ b/skimage/morphology/tests/test_ccomp.py @@ -2,6 +2,7 @@ import numpy as np from numpy.testing import assert_array_equal, run_module_suite from skimage.morphology import label +import skimage.measure._ccomp as ccomp from warnings import catch_warnings from skimage._shared.utils import skimage_deprecation @@ -244,10 +245,41 @@ class TestConnectedComponents3d: assert_array_equal(label(x, background=0, return_num=True)[1], 3) + def test_1D(self): + x = np.array((0, 1, 2, 2, 1, 1, 0, 0)) + xlen = len(x) + y = np.array((0, 1, 2, 2, 3, 3, 4, 4)) + reshapes = ((xlen,), + (1, xlen), (xlen, 1), + (1, xlen, 1), (xlen, 1, 1), (1, 1, xlen)) + for reshape in reshapes: + x2 = x.reshape(reshape) + labelled = label(x2) + assert_array_equal(y, labelled.flatten()) + def test_nd(self): x = np.ones((1, 2, 3, 4)) np.testing.assert_raises(NotImplementedError, label, x) +class TestSupport: + def test_reshape(self): + shapes_in = ((3, 1, 2), (1, 4, 5), (3, 1, 1), (2, 1), (1,)) + for shape in shapes_in: + shape = np.array(shape) + numones = sum(shape == 1) + inp = np.random.random(shape) + + fixed, swaps = ccomp.reshape_array(inp) + shape2 = fixed.shape + # now check that all ones are at the beginning + for i in range(numones): + assert shape2[i] == 1 + + back = ccomp.undo_reshape_array(fixed, swaps) + # check that the undo works as expected + assert_array_equal(inp, back) + + if __name__ == "__main__": run_module_suite() From c6ebe6c661b7368ae26842a9ecbfd1ad4b17747a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mat=C4=9Bj=20T=C3=BD=C4=8D?= Date: Tue, 2 Dec 2014 22:34:04 +0100 Subject: [PATCH 23/23] Removed negative connectivity support as requested --- skimage/measure/_ccomp.pyx | 49 ++++++++------------------------------ 1 file changed, 10 insertions(+), 39 deletions(-) diff --git a/skimage/measure/_ccomp.pyx b/skimage/measure/_ccomp.pyx index 466fc910..a9c4ab55 100644 --- a/skimage/measure/_ccomp.pyx +++ b/skimage/measure/_ccomp.pyx @@ -283,35 +283,6 @@ cdef inline void join_trees(DTYPE_t *forest, DTYPE_t n, DTYPE_t m): set_root(forest, m, root) -def _norm_connectivity(connectivity, ndim): - """ - Takes the value of the connectivity parameter, validates it and converts - it to a value that the subsequent algorithm may use as-is safely. - - Parameters - ---------- - connectivity : int - The following should be true: -ndim < connectivity < 0, or - 0 < connectivity <= ndim. - - Returns - ------- - connectivity : int - Connectivity, 0 < connectivity < ndim - """ - if connectivity == 0: - raise ValueError( - "Connectivity of 0 doesn't make sense") - if not -ndim <= connectivity <= ndim: - raise ValueError( - "Connectivity below %d or above %d is illegal." - % (-ndim, ndim)) - if connectivity < 0: - # Just as we say in the docs - connectivity += ndim + 1 - return connectivity - - def _get_swaps(shp): """ What axes to swap if we want to convert an illegal array shape @@ -412,9 +383,9 @@ def label(input, neighbors=None, background=None, return_num=False, input : ndarray of dtype int Image to label. neighbors : {4, 8}, int, optional - Whether to use 4- or 8-connectivity. - In 3D, 4-connectivity means connected pixels have to share face, - whereas with 8-connectivity, they have to share only edge or vertex. + Whether to use 4- or 8-"connectivity". + In 3D, 4-"connectivity" means connected pixels have to share face, + whereas with 8-"connectivity", they have to share only edge or vertex. **Deprecated, use ``connectivity`` instead.** background : int, optional Consider all pixels with this value as background pixels, and label @@ -425,12 +396,7 @@ def label(input, neighbors=None, background=None, return_num=False, connectivity : int, optional Maximum number of orthogonal hops to consider a pixel/voxel as a neighbor. - Accepted values are ranging from 1 to input.ndim. Negative values from - -input.ndim to -1 are also accepted. Negative connectivity value - :math:`x` is equivalent to the positive connectivity value - :math:`x + input.ndim + 1`. For example in 2D, -1 and 2 = -1 + 2 + 1 - are equivalent, values. - + Accepted values are ranging from 1 to input.ndim. Returns ------- @@ -518,9 +484,14 @@ def _label(input, neighbors=None, background=None, connectivity=None): raise ValueError("Neighbors must be either 4 or 8, got '%d'.\n" % neighbors) - connectivity = _norm_connectivity(connectivity, shapeinfo.ndim) + if not 1 <= connectivity <= input.ndim: + raise ValueError( + "Connectivity below 1 or above %d is illegal." + % input.ndim) scanBG(data_p, forest_p, &shapeinfo, &bg) + # the data are treated as degenerated 3D arrays if needed + # witout any performance sacrifice scan3D(data_p, forest_p, &shapeinfo, &bg, connectivity) # Label output