diff --git a/skimage/measure/_ccomp.pxd b/skimage/measure/_ccomp.pxd index 98ecf15a..b514e550 100644 --- a/skimage/measure/_ccomp.pxd +++ b/skimage/measure/_ccomp.pxd @@ -4,6 +4,6 @@ cimport numpy as cnp DTYPE = cnp.intp 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 DTYPE_t find_root(DTYPE_t *forest, DTYPE_t n) nogil +cdef void set_root(DTYPE_t *forest, DTYPE_t n, DTYPE_t root) nogil +cdef void join_trees(DTYPE_t *forest, DTYPE_t n, DTYPE_t m) nogil diff --git a/skimage/measure/_ccomp.pyx b/skimage/measure/_ccomp.pyx index a3b27593..90c921cc 100644 --- a/skimage/measure/_ccomp.pyx +++ b/skimage/measure/_ccomp.pyx @@ -234,7 +234,7 @@ cdef int ravel_index3D(int x, int y, int z, shape_info *shapeinfo): # 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): +cdef DTYPE_t find_root(DTYPE_t *forest, DTYPE_t n) nogil: """Find the root of node n. Given the example above, for any integer from 1 to 9, 1 is always returned """ @@ -244,7 +244,7 @@ cdef DTYPE_t find_root(DTYPE_t *forest, DTYPE_t n): return root -cdef inline void set_root(DTYPE_t *forest, DTYPE_t n, DTYPE_t root): +cdef inline void set_root(DTYPE_t *forest, DTYPE_t n, DTYPE_t root) nogil: """ 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. @@ -261,7 +261,7 @@ cdef inline void set_root(DTYPE_t *forest, DTYPE_t n, DTYPE_t root): forest[n] = root -cdef inline void join_trees(DTYPE_t *forest, DTYPE_t n, DTYPE_t m): +cdef inline void join_trees(DTYPE_t *forest, DTYPE_t n, DTYPE_t m) nogil: """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. @@ -416,7 +416,7 @@ def label(input, neighbors=None, background=None, return_num=False, [0 1 0] [0 0 1]] >>> from skimage.measure import label - >>> print(label(x, connectivity=1)) + >>> print(label(x, connectivity=1)) [[0 1 1] [2 3 1] [2 2 4]] diff --git a/skimage/segmentation/_felzenszwalb_cy.pyx b/skimage/segmentation/_felzenszwalb_cy.pyx index 37229225..0ab28e48 100644 --- a/skimage/segmentation/_felzenszwalb_cy.pyx +++ b/skimage/segmentation/_felzenszwalb_cy.pyx @@ -42,7 +42,9 @@ def _felzenszwalb_grey(image, double scale=1, sigma=0.8, if image.ndim != 2: raise ValueError("This algorithm works only on single-channel 2d" "images. Got image of shape %s" % str(image.shape)) + image = img_as_float(image) + # rescale scale to behave like in reference implementation scale = float(scale) / 255. image = scipy.ndimage.gaussian_filter(image, sigma=sigma) @@ -55,6 +57,7 @@ def _felzenszwalb_grey(image, double scale=1, sigma=0.8, cdef cnp.ndarray[cnp.float_t, ndim=1] costs = np.hstack([right_cost.ravel(), down_cost.ravel(), dright_cost.ravel(), uright_cost.ravel()]).astype(np.float) + # compute edges between pixels: height, width = image.shape[:2] cdef cnp.ndarray[cnp.intp_t, ndim=2] segments \ @@ -65,6 +68,7 @@ def _felzenszwalb_grey(image, double scale=1, sigma=0.8, uright_edges = np.c_[segments[:-1, 1:].ravel(), segments[1:, :-1].ravel()] cdef cnp.ndarray[cnp.intp_t, ndim=2] edges \ = np.vstack([right_edges, down_edges, dright_edges, uright_edges]) + # initialize data structures for segment size # and inner cost, then start greedy iteration over edges. edge_queue = np.argsort(costs) @@ -75,39 +79,43 @@ def _felzenszwalb_grey(image, double scale=1, sigma=0.8, cdef cnp.float_t *costs_p = costs.data cdef cnp.ndarray[cnp.intp_t, ndim=1] segment_size \ = np.ones(width * height, dtype=np.intp) + # inner cost of segments cdef cnp.ndarray[cnp.float_t, ndim=1] cint = np.zeros(width * height) cdef cnp.intp_t seg0, seg1, seg_new, e cdef float cost, inner_cost0, inner_cost1 - # set costs_p back one. we increase it before we use it - # since we might continue before that. - costs_p -= 1 - for e in range(costs.size): - seg0 = find_root(segments_p, edges_p[0]) - seg1 = find_root(segments_p, edges_p[1]) - edges_p += 2 - costs_p += 1 - if seg0 == seg1: - continue - inner_cost0 = cint[seg0] + scale / segment_size[seg0] - inner_cost1 = cint[seg1] + scale / segment_size[seg1] - if costs_p[0] < min(inner_cost0, inner_cost1): - # update size and cost - join_trees(segments_p, seg0, seg1) - seg_new = find_root(segments_p, seg0) - segment_size[seg_new] = segment_size[seg0] + segment_size[seg1] - cint[seg_new] = costs_p[0] + cdef Py_ssize_t num_costs = costs.size - # postprocessing to remove small segments - edges_p = edges.data - for e in range(costs.size): - seg0 = find_root(segments_p, edges_p[0]) - seg1 = find_root(segments_p, edges_p[1]) - edges_p += 2 - if seg0 == seg1: - continue - if segment_size[seg0] < min_size or segment_size[seg1] < min_size: - join_trees(segments_p, seg0, seg1) + with nogil: + # set costs_p back one. we increase it before we use it + # since we might continue before that. + costs_p -= 1 + for e in range(num_costs): + seg0 = find_root(segments_p, edges_p[0]) + seg1 = find_root(segments_p, edges_p[1]) + edges_p += 2 + costs_p += 1 + if seg0 == seg1: + continue + inner_cost0 = cint[seg0] + scale / segment_size[seg0] + inner_cost1 = cint[seg1] + scale / segment_size[seg1] + if costs_p[0] < min(inner_cost0, inner_cost1): + # update size and cost + join_trees(segments_p, seg0, seg1) + seg_new = find_root(segments_p, seg0) + segment_size[seg_new] = segment_size[seg0] + segment_size[seg1] + cint[seg_new] = costs_p[0] + + # postprocessing to remove small segments + edges_p = edges.data + for e in range(num_costs): + seg0 = find_root(segments_p, edges_p[0]) + seg1 = find_root(segments_p, edges_p[1]) + edges_p += 2 + if seg0 == seg1: + continue + if segment_size[seg0] < min_size or segment_size[seg1] < min_size: + join_trees(segments_p, seg0, seg1) # unravel the union find tree flat = segments.ravel() diff --git a/skimage/segmentation/tests/test_felzenszwalb.py b/skimage/segmentation/tests/test_felzenszwalb.py index 5324995d..c6f76a49 100644 --- a/skimage/segmentation/tests/test_felzenszwalb.py +++ b/skimage/segmentation/tests/test_felzenszwalb.py @@ -1,11 +1,11 @@ import numpy as np from numpy.testing import assert_equal, assert_array_equal -from skimage._shared.testing import assert_greater +from skimage._shared.testing import assert_greater, test_parallel from skimage.segmentation import felzenszwalb from skimage import data - +@test_parallel() def test_grey(): # very weak tests. This algorithm is pretty unstable. img = np.zeros((20, 21))