mirror of
https://github.com/wassname/scikit-image.git
synced 2026-06-30 10:15:18 +08:00
Merge pull request #1254 from matejak/label3d
Refactor labeling and extend to 3D
This commit is contained in:
Stefan van der Walt
@@ -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)
|
||||
|
||||
+679
-96
@@ -24,15 +24,219 @@ 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 shape_info
|
||||
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 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'
|
||||
))
|
||||
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
|
||||
|
||||
|
||||
# A pixel has neighbors that have already been scanned.
|
||||
# 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,
|
||||
# 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
|
||||
|
||||
|
||||
# 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 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.
|
||||
"""
|
||||
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:
|
||||
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
|
||||
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))
|
||||
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):
|
||||
# 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))
|
||||
else:
|
||||
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
|
||||
|
||||
# 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.
|
||||
|
||||
# So now the 1st (needed for 1D, 2D and 3D) part, y = 1, z = 1
|
||||
res.DEX[D_ed] = -1
|
||||
|
||||
# Not needed, just for illustration
|
||||
# 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_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_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] - res.DEX[D_eb]
|
||||
res.DEX[D_em] = res.DEX[D_el] + 1
|
||||
res.DEX[D_en] = res.DEX[D_el] + 2
|
||||
|
||||
|
||||
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):
|
||||
"""
|
||||
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):
|
||||
"""
|
||||
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, shape_info *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 "forest" = 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 +247,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 +263,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):
|
||||
@@ -71,30 +283,99 @@ 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):
|
||||
def _get_swaps(shp):
|
||||
"""
|
||||
Link a node to the background node.
|
||||
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
|
||||
"""
|
||||
if background_node[0] == -999:
|
||||
background_node[0] = n
|
||||
shp = np.array(shp)
|
||||
swaps = []
|
||||
|
||||
join_trees(forest, n, background_node[0])
|
||||
# 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, DTYPE_t neighbors=8, background=None, return_num=False):
|
||||
"""Label connected regions of an integer array.
|
||||
def label(input, neighbors=None, background=None, return_num=False,
|
||||
connectivity=None):
|
||||
r"""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 maximum number of orthogonal hops to consider a
|
||||
pixel/voxel a neighbor.
|
||||
|
||||
4-connectivity 8-connectivity
|
||||
1-connectivity 2-connectivity diagonal connection close-up
|
||||
|
||||
[ ] [ ] [ ] [ ]
|
||||
| \ | /
|
||||
[ ]--[ ]--[ ] [ ]--[ ]--[ ]
|
||||
| / | \\
|
||||
[ ] [ ] [ ] [ ] [ ]
|
||||
| \ | / | <- hop 2
|
||||
[ ]--[x]--[ ] [ ]--[x]--[ ] [x]--[ ]
|
||||
| / | \ hop 1
|
||||
[ ] [ ] [ ] [ ]
|
||||
|
||||
Parameters
|
||||
@@ -102,13 +383,20 @@ def label(input, DTYPE_t neighbors=8, 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.
|
||||
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
|
||||
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, optional
|
||||
Maximum number of orthogonal hops to consider a pixel/voxel
|
||||
as a neighbor.
|
||||
Accepted values are ranging from 1 to input.ndim.
|
||||
|
||||
Returns
|
||||
-------
|
||||
@@ -127,12 +415,12 @@ def label(input, DTYPE_t neighbors=8, 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]]
|
||||
@@ -147,94 +435,389 @@ 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]
|
||||
# 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)
|
||||
|
||||
cdef cnp.ndarray[DTYPE_t, ndim=2] data = np.array(input, copy=True,
|
||||
dtype=DTYPE)
|
||||
cdef cnp.ndarray[DTYPE_t, ndim=2] forest
|
||||
# Do the labelling
|
||||
res, ctr = _label(input_corrected, neighbors, background, connectivity)
|
||||
|
||||
forest = np.arange(data.size, dtype=DTYPE).reshape((rows, cols))
|
||||
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
|
||||
|
||||
# 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))
|
||||
forest = np.arange(data.size, dtype=DTYPE)
|
||||
|
||||
cdef DTYPE_t *forest_p = <DTYPE_t*>forest.data
|
||||
cdef DTYPE_t *data_p = <DTYPE_t*>data.data
|
||||
|
||||
cdef DTYPE_t i, j
|
||||
cdef shape_info shapeinfo
|
||||
cdef bginfo bg
|
||||
|
||||
cdef DTYPE_t background_val
|
||||
get_shape_info(input.shape, &shapeinfo)
|
||||
get_bginfo(background, &bg)
|
||||
|
||||
if background is None:
|
||||
background_val = -1
|
||||
warnings.warn(DeprecationWarning(
|
||||
'The default value for `background` will change to 0 in v0.12'
|
||||
))
|
||||
else:
|
||||
background_val = background
|
||||
if neighbors is None and connectivity is None:
|
||||
# use the full connectivity by default
|
||||
connectivity = input.ndim
|
||||
elif neighbors is not None:
|
||||
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:
|
||||
raise ValueError("Neighbors must be either 4 or 8, got '%d'.\n"
|
||||
% neighbors)
|
||||
|
||||
cdef DTYPE_t background_node = -999
|
||||
if not 1 <= connectivity <= input.ndim:
|
||||
raise ValueError(
|
||||
"Connectivity below 1 or above %d is illegal."
|
||||
% input.ndim)
|
||||
|
||||
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)
|
||||
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
|
||||
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]]
|
||||
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
|
||||
if data.dtype == np.int32:
|
||||
data = data.view(np.int32)
|
||||
|
||||
if return_num:
|
||||
return data, ctr
|
||||
else:
|
||||
return data
|
||||
res = data.reshape(input.shape)
|
||||
|
||||
return res, ctr
|
||||
|
||||
|
||||
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.
|
||||
We also track how many distinct final labels we have.
|
||||
"""
|
||||
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_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)
|
||||
data_p[i] = counter
|
||||
counter += 1
|
||||
else:
|
||||
data_p[i] = data_p[forest_p[i]]
|
||||
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.
|
||||
# The forward scan mask looks like this (the center point is actually E):
|
||||
# (take a look at shape_info 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.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 connectivity, 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, bgval = bg.background_val
|
||||
cdef INTS_t *DEX = shapeinfo.DEX
|
||||
rindex = shapeinfo.ravel_index(0, y, z, shapeinfo)
|
||||
|
||||
for x in range(1, shapeinfo.x):
|
||||
rindex += 1
|
||||
# Handle the first row
|
||||
if data_p[rindex] == bgval:
|
||||
# Nothing to do if we are background
|
||||
continue
|
||||
|
||||
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 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, 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:
|
||||
# Nothing to do if we are background
|
||||
|
||||
join_trees_wrapper(data_p, forest_p, rindex, DEX[D_eb])
|
||||
|
||||
if connectivity >= 2:
|
||||
join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ec])
|
||||
# END of x = 0
|
||||
|
||||
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.
|
||||
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])
|
||||
|
||||
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])
|
||||
|
||||
# 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
|
||||
# Handle first plane
|
||||
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
|
||||
|
||||
# Now we have pixels below
|
||||
join_trees_wrapper(data_p, forest_p, rindex, DEX[D_ej])
|
||||
|
||||
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])
|
||||
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 - 1):
|
||||
rindex += 1
|
||||
# Handle the first row
|
||||
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])
|
||||
|
||||
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])
|
||||
|
||||
# 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])
|
||||
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 x = 0
|
||||
rindex = shapeinfo.ravel_index(0, 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])
|
||||
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])
|
||||
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 - 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_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])
|
||||
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_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 - 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])
|
||||
|
||||
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_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
|
||||
|
||||
@@ -2,11 +2,16 @@ 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
|
||||
|
||||
np.random.seed(0)
|
||||
|
||||
# The background label value
|
||||
# is supposed to be changed to 0 soon
|
||||
BG = -1
|
||||
|
||||
|
||||
class TestConnectedComponents:
|
||||
def setup(self):
|
||||
@@ -75,10 +80,11 @@ 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]])
|
||||
[+1, 1, 1]])
|
||||
|
||||
def test_background_one_region_center(self):
|
||||
x = np.array([[0, 0, 0],
|
||||
@@ -101,5 +107,179 @@ 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, 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] == 2, \
|
||||
"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, 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)
|
||||
|
||||
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([[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)
|
||||
|
||||
def test_background_one_region_center(self):
|
||||
x = np.zeros((3, 3, 3), int)
|
||||
x[1, 1, 1] = 1
|
||||
|
||||
lb = np.ones_like(x) * BG
|
||||
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)
|
||||
|
||||
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()
|
||||
|
||||
Reference in New Issue
Block a user