Do not acquire GIL for felzenswalb

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

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]]

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 = <cnp.float_t*>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 = <cnp.intp_t*>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 = <cnp.intp_t*>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()

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))