From 0d4c76948bbe1bdfca432076f30fb9aeea80db95 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 19 May 2015 18:22:32 -0700 Subject: [PATCH 01/20] Do not acquire GIL for warp_fast --- skimage/_shared/interpolation.pxd | 21 +++++++++++---------- skimage/transform/_warps_cy.pyx | 15 ++++++++------- 2 files changed, 19 insertions(+), 17 deletions(-) diff --git a/skimage/_shared/interpolation.pxd b/skimage/_shared/interpolation.pxd index 088be1ef..3b4ac538 100644 --- a/skimage/_shared/interpolation.pxd +++ b/skimage/_shared/interpolation.pxd @@ -5,7 +5,7 @@ from libc.math cimport ceil, floor -cdef inline Py_ssize_t round(double r): +cdef inline Py_ssize_t round(double r) nogil: return ((r + 0.5) if (r > 0.0) else (r - 0.5)) @@ -13,7 +13,7 @@ cdef inline double nearest_neighbour_interpolation(double* image, Py_ssize_t rows, Py_ssize_t cols, double r, double c, char mode, - double cval): + double cval) nogil: """Nearest neighbour interpolation at a given position in the image. Parameters @@ -41,7 +41,7 @@ cdef inline double nearest_neighbour_interpolation(double* image, cdef inline double bilinear_interpolation(double* image, Py_ssize_t rows, Py_ssize_t cols, double r, double c, - char mode, double cval): + char mode, double cval) nogil: """Bilinear interpolation at a given position in the image. Parameters @@ -80,7 +80,7 @@ cdef inline double bilinear_interpolation(double* image, Py_ssize_t rows, return (1 - dr) * top + dr * bottom -cdef inline double quadratic_interpolation(double x, double[3] f): +cdef inline double quadratic_interpolation(double x, double[3] f) nogil: """WARNING: Do not use, not implemented correctly. Quadratic interpolation. @@ -105,7 +105,8 @@ cdef inline double quadratic_interpolation(double x, double[3] f): cdef inline double biquadratic_interpolation(double* image, Py_ssize_t rows, Py_ssize_t cols, double r, - double c, char mode, double cval): + double c, char mode, + double cval) nogil: """WARNING: Do not use, not implemented correctly. Biquadratic interpolation at a given position in the image. @@ -152,7 +153,7 @@ cdef inline double biquadratic_interpolation(double* image, Py_ssize_t rows, return quadratic_interpolation(xr, fr) -cdef inline double cubic_interpolation(double x, double[4] f): +cdef inline double cubic_interpolation(double x, double[4] f) nogil: """Cubic interpolation. Parameters @@ -177,7 +178,7 @@ cdef inline double cubic_interpolation(double x, double[4] f): cdef inline double bicubic_interpolation(double* image, Py_ssize_t rows, Py_ssize_t cols, double r, double c, - char mode, double cval): + char mode, double cval) nogil: """Bicubic interpolation at a given position in the image. Interpolation using Catmull-Rom splines, based on the bicubic convolution @@ -236,7 +237,7 @@ cdef inline double bicubic_interpolation(double* image, Py_ssize_t rows, cdef inline double get_pixel2d(double* image, Py_ssize_t rows, Py_ssize_t cols, long r, long c, char mode, - double cval): + double cval) nogil: """Get a pixel from the image, taking wrapping mode into consideration. Parameters @@ -269,7 +270,7 @@ cdef inline double get_pixel2d(double* image, Py_ssize_t rows, Py_ssize_t cols, cdef inline double get_pixel3d(double* image, Py_ssize_t rows, Py_ssize_t cols, Py_ssize_t dims, long r, long c, - long d, char mode, double cval): + long d, char mode, double cval) nogil: """Get a pixel from the image, taking wrapping mode into consideration. Parameters @@ -302,7 +303,7 @@ cdef inline double get_pixel3d(double* image, Py_ssize_t rows, Py_ssize_t cols, + coord_map(dims, d, mode)] -cdef inline Py_ssize_t coord_map(Py_ssize_t dim, long coord, char mode): +cdef inline Py_ssize_t coord_map(Py_ssize_t dim, long coord, char mode) nogil: """Wrap a coordinate, according to a given mode. Parameters diff --git a/skimage/transform/_warps_cy.pyx b/skimage/transform/_warps_cy.pyx index 27b6c95a..aa2f3b6e 100644 --- a/skimage/transform/_warps_cy.pyx +++ b/skimage/transform/_warps_cy.pyx @@ -11,7 +11,7 @@ from .._shared.interpolation cimport (nearest_neighbour_interpolation, cdef inline void _matrix_transform(double x, double y, double* H, double *x_, - double *y_): + double *y_) nogil: """Apply a homography to a coordinate. Parameters @@ -102,7 +102,7 @@ def _warp_fast(cnp.ndarray image, cnp.ndarray H, output_shape=None, cdef Py_ssize_t cols = img.shape[1] cdef double (*interp_func)(double*, Py_ssize_t, Py_ssize_t, double, double, - char, double) + char, double) nogil if order == 0: interp_func = nearest_neighbour_interpolation elif order == 1: @@ -112,10 +112,11 @@ def _warp_fast(cnp.ndarray image, cnp.ndarray H, output_shape=None, elif order == 3: interp_func = bicubic_interpolation - for tfr in range(out_r): - for tfc in range(out_c): - _matrix_transform(tfc, tfr, &M[0, 0], &c, &r) - out[tfr, tfc] = interp_func(&img[0, 0], rows, cols, r, c, - mode_c, cval) + with nogil: + for tfr in range(out_r): + for tfc in range(out_c): + _matrix_transform(tfc, tfr, &M[0, 0], &c, &r) + out[tfr, tfc] = interp_func(&img[0, 0], rows, cols, r, c, + mode_c, cval) return np.asarray(out) From 07f2e4b93f06e52689cdb9a73cbbf9556e93a0cd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 19 May 2015 18:25:01 -0700 Subject: [PATCH 02/20] Test parallel execution of warp_fast --- skimage/transform/tests/test_warps.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/skimage/transform/tests/test_warps.py b/skimage/transform/tests/test_warps.py index afa78da3..4de59359 100644 --- a/skimage/transform/tests/test_warps.py +++ b/skimage/transform/tests/test_warps.py @@ -11,6 +11,7 @@ from skimage.transform import (warp, warp_coords, rotate, resize, rescale, from skimage import transform as tf, data, img_as_float from skimage.color import rgb2gray from skimage._shared._warnings import expected_warnings +from skimage._shared.testing import test_parallel np.random.seed(0) @@ -41,6 +42,7 @@ def test_warp_callable(): assert_almost_equal(outx, refx) +@test_parallel() def test_warp_matrix(): x = np.zeros((5, 5), dtype=np.double) x[2, 2] = 1 From 7a8afbddf3e4f3d68dbc63c4c359027746d141bc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 19 May 2015 18:37:36 -0700 Subject: [PATCH 03/20] Do not acquire GIL for iradon_sart --- skimage/transform/_radon_transform.pyx | 158 +++++++++--------- .../transform/tests/test_radon_transform.py | 2 + 2 files changed, 84 insertions(+), 76 deletions(-) diff --git a/skimage/transform/_radon_transform.pyx b/skimage/transform/_radon_transform.pyx index 91b943b9..67ddffed 100644 --- a/skimage/transform/_radon_transform.pyx +++ b/skimage/transform/_radon_transform.pyx @@ -40,48 +40,51 @@ cpdef bilinear_ray_sum(cnp.double_t[:, :] image, cnp.double_t theta, # s0 is the half-length of the ray's path in the reconstruction circle cdef cnp.double_t s0 s0 = sqrt(radius**2 - t**2) if radius**2 >= t**2 else 0. - cdef Py_ssize_t Ns = 2 * ( ceil(2 * s0)) # number of steps - # along the ray + cdef Py_ssize_t Ns = 2 * (ceil(2 * s0)) # number of steps + # along the ray cdef cnp.double_t ray_sum = 0. cdef cnp.double_t weight_norm = 0. cdef cnp.double_t ds, dx, dy, x0, y0, x, y, di, dj, cdef cnp.double_t index_i, index_j, weight cdef Py_ssize_t k, i, j - if Ns > 0: - # step length between samples - ds = 2 * s0 / Ns - dx = -ds * cos(theta) - dy = -ds * sin(theta) - # point of entry of the ray into the reconstruction circle - x0 = s0 * cos(theta) - t * sin(theta) - y0 = s0 * sin(theta) + t * cos(theta) - for k in range(Ns+1): - x = x0 + k * dx - y = y0 + k * dy - index_i = x + rotation_center - index_j = y + rotation_center - i = floor(index_i) - j = floor(index_j) - di = index_i - floor(index_i) - dj = index_j - floor(index_j) - # Use linear interpolation between values - # Where values fall outside the array, assume zero - if i > 0 and j > 0: - weight = (1. - di) * (1. - dj) * ds - ray_sum += weight * image[i, j] - weight_norm += weight**2 - if i > 0 and j < image.shape[1] - 1: - weight = (1. - di) * dj * ds - ray_sum += weight * image[i, j+1] - weight_norm += weight**2 - if i < image.shape[0] - 1 and j > 0: - weight = di * (1 - dj) * ds - ray_sum += weight * image[i+1, j] - weight_norm += weight**2 - if i < image.shape[0] - 1 and j < image.shape[1] - 1: - weight = di * dj * ds - ray_sum += weight * image[i+1, j+1] - weight_norm += weight**2 + + with nogil: + if Ns > 0: + # step length between samples + ds = 2 * s0 / Ns + dx = -ds * cos(theta) + dy = -ds * sin(theta) + # point of entry of the ray into the reconstruction circle + x0 = s0 * cos(theta) - t * sin(theta) + y0 = s0 * sin(theta) + t * cos(theta) + for k in range(Ns + 1): + x = x0 + k * dx + y = y0 + k * dy + index_i = x + rotation_center + index_j = y + rotation_center + i = floor(index_i) + j = floor(index_j) + di = index_i - floor(index_i) + dj = index_j - floor(index_j) + # Use linear interpolation between values + # Where values fall outside the array, assume zero + if i > 0 and j > 0: + weight = (1. - di) * (1. - dj) * ds + ray_sum += weight * image[i, j] + weight_norm += weight**2 + if i > 0 and j < image.shape[1] - 1: + weight = (1. - di) * dj * ds + ray_sum += weight * image[i, j+1] + weight_norm += weight**2 + if i < image.shape[0] - 1 and j > 0: + weight = di * (1 - dj) * ds + ray_sum += weight * image[i+1, j] + weight_norm += weight**2 + if i < image.shape[0] - 1 and j < image.shape[1] - 1: + weight = di * dj * ds + ray_sum += weight * image[i+1, j+1] + weight_norm += weight**2 + return ray_sum, weight_norm @@ -89,26 +92,25 @@ cpdef bilinear_ray_update(cnp.double_t[:, :] image, cnp.double_t[:, :] image_update, cnp.double_t theta, cnp.double_t ray_position, cnp.double_t projected_value): - """ - Compute the update along a ray using bilinear interpolation. + """Compute the update along a ray using bilinear interpolation. Parameters ---------- image : 2D array, dtype=float - Current reconstruction estimate + Current reconstruction estimate. image_update : 2D array, dtype=float Array of same shape as ``image``. Updates will be added to this array. theta : float - Angle of the projection + Angle of the projection. ray_position : float - Position of the ray within the projection + Position of the ray within the projection. projected_value : float - Projected value (from the sinogram) + Projected value (from the sinogram). Returns ------- deviation : - Deviation before updating the image + Deviation before updating the image. """ cdef cnp.double_t ray_sum, weight_norm, deviation ray_sum, weight_norm = bilinear_ray_sum(image, theta, ray_position) @@ -125,43 +127,47 @@ cpdef bilinear_ray_update(cnp.double_t[:, :] image, # s0 is the half-length of the ray's path in the reconstruction circle cdef cnp.double_t s0 s0 = sqrt(radius*radius - t*t) if radius**2 >= t**2 else 0. - cdef Py_ssize_t Ns = 2 * ( ceil(2 * s0)) - cdef cnp.double_t hamming_beta = 0.46164 # beta for equiripple Hamming window + cdef Py_ssize_t Ns = 2 * (ceil(2 * s0)) + # beta for equiripple Hamming window + cdef cnp.double_t hamming_beta = 0.46164 cdef cnp.double_t ds, dx, dy, x0, y0, x, y, di, dj, index_i, index_j cdef cnp.double_t hamming_window cdef Py_ssize_t k, i, j - if Ns > 0: - # Step length between samples - ds = 2 * s0 / Ns - dx = -ds * cos(theta) - dy = -ds * sin(theta) - # Point of entry of the ray into the reconstruction circle - x0 = s0 * cos(theta) - t * sin(theta) - y0 = s0 * sin(theta) + t * cos(theta) - for k in range(Ns+1): - x = x0 + k * dx - y = y0 + k * dy - index_i = x + rotation_center - index_j = y + rotation_center - i = floor(index_i) - j = floor(index_j) - di = index_i - floor(index_i) - dj = index_j - floor(index_j) - hamming_window = ((1 - hamming_beta) - - hamming_beta * cos(2 * M_PI * k / (Ns - 1))) - if i > 0 and j > 0: - image_update[i, j] += (deviation * (1. - di) * (1. - dj) - * ds * hamming_window) - if i > 0 and j < image.shape[1] - 1: - image_update[i, j+1] += (deviation * (1. - di) * dj - * ds * hamming_window) - if i < image.shape[0] - 1 and j > 0: - image_update[i+1, j] += (deviation * di * (1 - dj) - * ds * hamming_window) - if i < image.shape[0] - 1 and j < image.shape[1] - 1: - image_update[i+1, j+1] += (deviation * di * dj + + with nogil: + if Ns > 0: + # Step length between samples + ds = 2 * s0 / Ns + dx = -ds * cos(theta) + dy = -ds * sin(theta) + # Point of entry of the ray into the reconstruction circle + x0 = s0 * cos(theta) - t * sin(theta) + y0 = s0 * sin(theta) + t * cos(theta) + for k in range(Ns + 1): + x = x0 + k * dx + y = y0 + k * dy + index_i = x + rotation_center + index_j = y + rotation_center + i = floor(index_i) + j = floor(index_j) + di = index_i - floor(index_i) + dj = index_j - floor(index_j) + hamming_window = ((1 - hamming_beta) + - hamming_beta * cos(2 * M_PI * k / (Ns - 1))) + if i > 0 and j > 0: + image_update[i, j] += (deviation * (1. - di) * (1. - dj) * ds * hamming_window) + if i > 0 and j < image.shape[1] - 1: + image_update[i, j+1] += (deviation * (1. - di) * dj + * ds * hamming_window) + if i < image.shape[0] - 1 and j > 0: + image_update[i+1, j] += (deviation * di * (1 - dj) + * ds * hamming_window) + if i < image.shape[0] - 1 and j < image.shape[1] - 1: + image_update[i+1, j+1] += (deviation * di * dj + * ds * hamming_window) + return deviation diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py index 9e51ad0a..6df8b7f8 100644 --- a/skimage/transform/tests/test_radon_transform.py +++ b/skimage/transform/tests/test_radon_transform.py @@ -8,6 +8,7 @@ import os.path from skimage.transform import radon, iradon, iradon_sart, rescale from skimage.io import imread from skimage import data_dir +from skimage._shared.testing import test_parallel PHANTOM = imread(os.path.join(data_dir, "phantom.png"), @@ -310,6 +311,7 @@ def test_order_angles_golden_ratio(): assert len(indices) == len(set(indices)) +@test_parallel() def test_iradon_sart(): debug = False From 83fe07097f4d503184244a203476ed6a3ef8de2e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 19 May 2015 18:42:53 -0700 Subject: [PATCH 04/20] Do not acquire GIL for hough_circle --- skimage/transform/_hough_transform.pyx | 32 ++++++++++--------- .../transform/tests/test_hough_transform.py | 2 ++ 2 files changed, 19 insertions(+), 15 deletions(-) diff --git a/skimage/transform/_hough_transform.pyx b/skimage/transform/_hough_transform.pyx index e51085cb..601ac0ea 100644 --- a/skimage/transform/_hough_transform.pyx +++ b/skimage/transform/_hough_transform.pyx @@ -79,22 +79,24 @@ def _hough_circle(cnp.ndarray img, num_circle_pixels = circle_x.size - if normalize: - incr = 1.0 / num_circle_pixels - else: - incr = 1 + with nogil: - # For each non zero pixel - for p in range(num_pixels): - # Plug the circle at (px, py), - # its coordinates are (tx, ty) - for c in range(num_circle_pixels): - tx = circle_x[c] + x[p] - ty = circle_y[c] + y[p] - if offset: - acc[i, tx, ty] += incr - elif 0 <= tx < xmax and 0 <= ty < ymax: - acc[i, tx, ty] += incr + if normalize: + incr = 1.0 / num_circle_pixels + else: + incr = 1 + + # For each non zero pixel + for p in range(num_pixels): + # Plug the circle at (px, py), + # its coordinates are (tx, ty) + for c in range(num_circle_pixels): + tx = circle_x[c] + x[p] + ty = circle_y[c] + y[p] + if offset: + acc[i, tx, ty] += incr + elif 0 <= tx < xmax and 0 <= ty < ymax: + acc[i, tx, ty] += incr return acc diff --git a/skimage/transform/tests/test_hough_transform.py b/skimage/transform/tests/test_hough_transform.py index 884d7957..ea9eac42 100644 --- a/skimage/transform/tests/test_hough_transform.py +++ b/skimage/transform/tests/test_hough_transform.py @@ -4,6 +4,7 @@ from numpy.testing import assert_almost_equal, assert_equal import skimage.transform as tf from skimage.draw import line, circle_perimeter, ellipse_perimeter from skimage._shared._warnings import expected_warnings +from skimage._shared.testing import test_parallel def append_desc(func, description): @@ -129,6 +130,7 @@ def test_hough_line_peaks_num(): min_angle=0, num_peaks=1)[0]) == 1 +@test_parallel() def test_hough_circle(): # Prepare picture img = np.zeros((120, 100), dtype=int) From 7337a3e14b09f54985eb9b7878eabf95de95e4fb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 19 May 2015 18:57:20 -0700 Subject: [PATCH 05/20] Do not acquire GIL for hough_line --- skimage/transform/_hough_transform.pyx | 19 +++++++++++-------- .../transform/tests/test_hough_transform.py | 1 + 2 files changed, 12 insertions(+), 8 deletions(-) diff --git a/skimage/transform/_hough_transform.pyx b/skimage/transform/_hough_transform.pyx index 601ac0ea..6ce69d8e 100644 --- a/skimage/transform/_hough_transform.pyx +++ b/skimage/transform/_hough_transform.pyx @@ -308,14 +308,17 @@ def hough_line(cnp.ndarray img, # finally, run the transform cdef Py_ssize_t nidxs, nthetas, i, j, x, y, accum_idx - nidxs = y_idxs.shape[0] # x and y are the same shape - nthetas = theta.shape[0] - for i in range(nidxs): - x = x_idxs[i] - y = y_idxs[i] - for j in range(nthetas): - accum_idx = round((ctheta[j] * x + stheta[j] * y)) + offset - accum[accum_idx, j] += 1 + + with nogil: + nidxs = y_idxs.shape[0] # x and y are the same shape + nthetas = theta.shape[0] + for i in range(nidxs): + x = x_idxs[i] + y = y_idxs[i] + for j in range(nthetas): + accum_idx = round((ctheta[j] * x + stheta[j] * y)) + offset + accum[accum_idx, j] += 1 + return accum, theta, bins diff --git a/skimage/transform/tests/test_hough_transform.py b/skimage/transform/tests/test_hough_transform.py index ea9eac42..e7030494 100644 --- a/skimage/transform/tests/test_hough_transform.py +++ b/skimage/transform/tests/test_hough_transform.py @@ -16,6 +16,7 @@ def append_desc(func, description): return func +@test_parallel() def test_hough_line(): # Generate a test image img = np.zeros((100, 150), dtype=int) From 785176b4c2dac3ddd2a7d5d7ee9781af21d686a3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 19 May 2015 22:40:19 -0700 Subject: [PATCH 06/20] Do not acquire GIL for slic --- skimage/segmentation/_slic.pyx | 231 ++++++++++++------------ skimage/segmentation/tests/test_slic.py | 2 + 2 files changed, 118 insertions(+), 115 deletions(-) diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx index 0935550c..19fc1627 100644 --- a/skimage/segmentation/_slic.pyx +++ b/skimage/segmentation/_slic.pyx @@ -101,88 +101,89 @@ def _slic_cython(double[:, :, :, ::1] image_zyx, # The reference implementation (Achanta et al.) calls this invxywt cdef double spatial_weight = float(1) / (step ** 2) - for i in range(max_iter): - change = 0 - distance[:, :, :] = DBL_MAX + with nogil: + for i in range(max_iter): + change = 0 + distance[:, :, :] = DBL_MAX - # assign pixels to segments - for k in range(n_segments): + # assign pixels to segments + for k in range(n_segments): - # segment coordinate centers - cz = segments[k, 0] - cy = segments[k, 1] - cx = segments[k, 2] + # segment coordinate centers + cz = segments[k, 0] + cy = segments[k, 1] + cx = segments[k, 2] - # compute windows - z_min = max(cz - 2 * step_z, 0) - z_max = min(cz + 2 * step_z + 1, depth) - y_min = max(cy - 2 * step_y, 0) - y_max = min(cy + 2 * step_y + 1, height) - x_min = max(cx - 2 * step_x, 0) - x_max = min(cx + 2 * step_x + 1, width) + # compute windows + z_min = max(cz - 2 * step_z, 0) + z_max = min(cz + 2 * step_z + 1, depth) + y_min = max(cy - 2 * step_y, 0) + y_max = min(cy + 2 * step_y + 1, height) + x_min = max(cx - 2 * step_x, 0) + x_max = min(cx + 2 * step_x + 1, width) - for z in range(z_min, z_max): - dz = (sz * (cz - z)) ** 2 - for y in range(y_min, y_max): - dy = (sy * (cy - y)) ** 2 - for x in range(x_min, x_max): - dist_center = (dz + dy + (sx * (cx - x)) ** 2) * spatial_weight - dist_color = 0 - for c in range(3, n_features): - dist_color += (image_zyx[z, y, x, c - 3] - - segments[k, c]) ** 2 - if slic_zero: - dist_center += dist_color / max_dist_color[k] - else: - dist_center += dist_color + for z in range(z_min, z_max): + dz = (sz * (cz - z)) ** 2 + for y in range(y_min, y_max): + dy = (sy * (cy - y)) ** 2 + for x in range(x_min, x_max): + dist_center = (dz + dy + (sx * (cx - x)) ** 2) * spatial_weight + dist_color = 0 + for c in range(3, n_features): + dist_color += (image_zyx[z, y, x, c - 3] + - segments[k, c]) ** 2 + if slic_zero: + dist_center += dist_color / max_dist_color[k] + else: + dist_center += dist_color - if distance[z, y, x] > dist_center: - nearest_segments[z, y, x] = k - distance[z, y, x] = dist_center - change = 1 + if distance[z, y, x] > dist_center: + nearest_segments[z, y, x] = k + distance[z, y, x] = dist_center + change = 1 - # stop if no pixel changed its segment - if change == 0: - break + # stop if no pixel changed its segment + if change == 0: + break - # recompute segment centers + # recompute segment centers - # sum features for all segments - n_segment_elems[:] = 0 - segments[:, :] = 0 - for z in range(depth): - for y in range(height): - for x in range(width): - k = nearest_segments[z, y, x] - n_segment_elems[k] += 1 - segments[k, 0] += z - segments[k, 1] += y - segments[k, 2] += x - for c in range(3, n_features): - segments[k, c] += image_zyx[z, y, x, c - 3] - - # divide by number of elements per segment to obtain mean - for k in range(n_segments): - for c in range(n_features): - segments[k, c] /= n_segment_elems[k] - - # If in SLICO mode, update the color distance maxima - if slic_zero: + # sum features for all segments + n_segment_elems[:] = 0 + segments[:, :] = 0 for z in range(depth): for y in range(height): for x in range(width): - k = nearest_segments[z, y, x] - dist_color = 0 - + n_segment_elems[k] += 1 + segments[k, 0] += z + segments[k, 1] += y + segments[k, 2] += x for c in range(3, n_features): - dist_color += (image_zyx[z, y, x, c - 3] - - segments[k, c]) ** 2 + segments[k, c] += image_zyx[z, y, x, c - 3] - # The reference implementation seems to only change - # the color if it increases from previous iteration - if max_dist_color[k] < dist_color: - max_dist_color[k] = dist_color + # divide by number of elements per segment to obtain mean + for k in range(n_segments): + for c in range(n_features): + segments[k, c] /= n_segment_elems[k] + + # If in SLICO mode, update the color distance maxima + if slic_zero: + for z in range(depth): + for y in range(height): + for x in range(width): + + k = nearest_segments[z, y, x] + dist_color = 0 + + for c in range(3, n_features): + dist_color += (image_zyx[z, y, x, c - 3] - + segments[k, c]) ** 2 + + # The reference implementation seems to only change + # the color if it increases from previous iteration + if max_dist_color[k] < dist_color: + max_dist_color[k] = dist_color return np.asarray(nearest_segments) @@ -237,54 +238,54 @@ def _enforce_label_connectivity_cython(Py_ssize_t[:, :, ::1] segments, cdef Py_ssize_t[:, ::1] coord_list = np.zeros((max_size, 3), dtype=np.intp) - # loop through all image - for z in range(depth): - for y in range(height): - for x in range(width): - if connected_segments[z, y, x] >= 0: - continue - # find the component size - adjacent = 0 - label = segments[z, y, x] - connected_segments[z, y, x] = current_new_label - current_segment_size = 1 - bfs_visited = 0 - coord_list[bfs_visited, 0] = z - coord_list[bfs_visited, 1] = y - coord_list[bfs_visited, 2] = x + with nogil: + for z in range(depth): + for y in range(height): + for x in range(width): + if connected_segments[z, y, x] >= 0: + continue + # find the component size + adjacent = 0 + label = segments[z, y, x] + connected_segments[z, y, x] = current_new_label + current_segment_size = 1 + bfs_visited = 0 + coord_list[bfs_visited, 0] = z + coord_list[bfs_visited, 1] = y + coord_list[bfs_visited, 2] = x - #perform a breadth first search to find - # the size of the connected component - while bfs_visited < current_segment_size < max_size: - for i in range(6): - zz = coord_list[bfs_visited, 0] + ddz[i] - yy = coord_list[bfs_visited, 1] + ddy[i] - xx = coord_list[bfs_visited, 2] + ddx[i] - if (0 <= xx < width and - 0 <= yy < height and - 0 <= zz < depth): - if (segments[zz, yy, xx] == label and - connected_segments[zz, yy, xx] == -1): - connected_segments[zz, yy, xx] = \ - current_new_label - coord_list[current_segment_size, 0] = zz - coord_list[current_segment_size, 1] = yy - coord_list[current_segment_size, 2] = xx - current_segment_size += 1 - if current_segment_size >= max_size: - break - elif (connected_segments[zz, yy, xx] >= 0 and - connected_segments[zz, yy, xx] != current_new_label): - adjacent = connected_segments[zz, yy, xx] - bfs_visited += 1 + #perform a breadth first search to find + # the size of the connected component + while bfs_visited < current_segment_size < max_size: + for i in range(6): + zz = coord_list[bfs_visited, 0] + ddz[i] + yy = coord_list[bfs_visited, 1] + ddy[i] + xx = coord_list[bfs_visited, 2] + ddx[i] + if (0 <= xx < width and + 0 <= yy < height and + 0 <= zz < depth): + if (segments[zz, yy, xx] == label and + connected_segments[zz, yy, xx] == -1): + connected_segments[zz, yy, xx] = \ + current_new_label + coord_list[current_segment_size, 0] = zz + coord_list[current_segment_size, 1] = yy + coord_list[current_segment_size, 2] = xx + current_segment_size += 1 + if current_segment_size >= max_size: + break + elif (connected_segments[zz, yy, xx] >= 0 and + connected_segments[zz, yy, xx] != current_new_label): + adjacent = connected_segments[zz, yy, xx] + bfs_visited += 1 - # change to an adjacent one, like in the original paper - if current_segment_size < min_size: - for i in range(current_segment_size): - connected_segments[coord_list[i, 0], - coord_list[i, 1], - coord_list[i, 2]] = adjacent - else: - current_new_label += 1 + # change to an adjacent one, like in the original paper + if current_segment_size < min_size: + for i in range(current_segment_size): + connected_segments[coord_list[i, 0], + coord_list[i, 1], + coord_list[i, 2]] = adjacent + else: + current_new_label += 1 return np.asarray(connected_segments) diff --git a/skimage/segmentation/tests/test_slic.py b/skimage/segmentation/tests/test_slic.py index 2dd4cde7..73629103 100644 --- a/skimage/segmentation/tests/test_slic.py +++ b/skimage/segmentation/tests/test_slic.py @@ -2,8 +2,10 @@ import itertools as it import numpy as np from numpy.testing import assert_equal, assert_raises from skimage.segmentation import slic +from skimage._shared.testing import test_parallel +@test_parallel() def test_color_2d(): rnd = np.random.RandomState(0) img = np.zeros((20, 21, 3)) From 2d0b4dc37a59e698268d73313a54540b76e377be Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 19 May 2015 22:43:46 -0700 Subject: [PATCH 07/20] Do not acquire GIL for quickshift --- skimage/segmentation/_quickshift.pyx | 75 ++++++++++--------- skimage/segmentation/tests/test_quickshift.py | 3 +- 2 files changed, 41 insertions(+), 37 deletions(-) diff --git a/skimage/segmentation/_quickshift.pyx b/skimage/segmentation/_quickshift.pyx index b43775b6..fef74a9e 100644 --- a/skimage/segmentation/_quickshift.pyx +++ b/skimage/segmentation/_quickshift.pyx @@ -8,6 +8,7 @@ from itertools import product cimport numpy as cnp from libc.math cimport exp, sqrt +from libc.float cimport DBL_MAX from ..util import img_as_float from ..color import rgb2lab @@ -99,19 +100,20 @@ def quickshift(image, ratio=1., float kernel_size=5, max_dist=10, = np.zeros((height, width)) # compute densities - for r in range(height): - for c in range(width): - r_min, r_max = max(r - w, 0), min(r + w + 1, height) - c_min, c_max = max(c - w, 0), min(c + w + 1, width) - for r_ in range(r_min, r_max): - for c_ in range(c_min, c_max): - dist = 0 - for channel in range(channels): - dist += (current_pixel_p[channel] - - image_c[r_, c_, channel])**2 - dist += (r - r_)**2 + (c - c_)**2 - densities[r, c] += exp(-dist / (2 * kernel_size**2)) - current_pixel_p += channels + with nogil: + for r in range(height): + for c in range(width): + r_min, r_max = max(r - w, 0), min(r + w + 1, height) + c_min, c_max = max(c - w, 0), min(c + w + 1, width) + for r_ in range(r_min, r_max): + for c_ in range(c_min, c_max): + dist = 0 + for channel in range(channels): + dist += (current_pixel_p[channel] - + image_c[r_, c_, channel])**2 + dist += (r - r_)**2 + (c - c_)**2 + densities[r, c] += exp(-dist / (2 * kernel_size**2)) + current_pixel_p += channels # this will break ties that otherwise would give us headache densities += random_state.normal(scale=0.00001, size=(height, width)) @@ -123,29 +125,30 @@ def quickshift(image, ratio=1., float kernel_size=5, max_dist=10, = np.zeros((height, width)) # find nearest node with higher density - current_pixel_p = image_p - for r in range(height): - for c in range(width): - current_density = densities[r, c] - closest = np.inf - r_min, r_max = max(r - w, 0), min(r + w + 1, height) - c_min, c_max = max(c - w, 0), min(c + w + 1, width) - for r_ in range(r_min, r_max): - for c_ in range(c_min, c_max): - if densities[r_, c_] > current_density: - dist = 0 - # We compute the distances twice since otherwise - # we get crazy memory overhead - # (width * height * windowsize**2) - for channel in range(channels): - dist += (current_pixel_p[channel] - - image_c[r_, c_, channel])**2 - dist += (r - r_)**2 + (c - c_)**2 - if dist < closest: - closest = dist - parent[r, c] = r_ * width + c_ - dist_parent[r, c] = sqrt(closest) - current_pixel_p += channels + with nogil: + current_pixel_p = image_p + for r in range(height): + for c in range(width): + current_density = densities[r, c] + closest = DBL_MAX + r_min, r_max = max(r - w, 0), min(r + w + 1, height) + c_min, c_max = max(c - w, 0), min(c + w + 1, width) + for r_ in range(r_min, r_max): + for c_ in range(c_min, c_max): + if densities[r_, c_] > current_density: + dist = 0 + # We compute the distances twice since otherwise + # we get crazy memory overhead + # (width * height * windowsize**2) + for channel in range(channels): + dist += (current_pixel_p[channel] - + image_c[r_, c_, channel])**2 + dist += (r - r_)**2 + (c - c_)**2 + if dist < closest: + closest = dist + parent[r, c] = r_ * width + c_ + dist_parent[r, c] = sqrt(closest) + current_pixel_p += channels dist_parent_flat = dist_parent.ravel() flat = parent.ravel() diff --git a/skimage/segmentation/tests/test_quickshift.py b/skimage/segmentation/tests/test_quickshift.py index a940f859..21dbe8ea 100644 --- a/skimage/segmentation/tests/test_quickshift.py +++ b/skimage/segmentation/tests/test_quickshift.py @@ -1,10 +1,11 @@ import numpy as np from numpy.testing import assert_equal, assert_array_equal from nose.tools import assert_true -from skimage._shared.testing import assert_greater +from skimage._shared.testing import assert_greater, test_parallel from skimage.segmentation import quickshift +@test_parallel() def test_grey(): rnd = np.random.RandomState(0) img = np.zeros((20, 21)) From 0758fc9cd6a2d77157988545ca68c4b7d2950492 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 19 May 2015 22:52:29 -0700 Subject: [PATCH 08/20] Do not acquire GIL for felzenswalb --- skimage/measure/_ccomp.pxd | 6 +- skimage/measure/_ccomp.pyx | 8 +-- skimage/segmentation/_felzenszwalb_cy.pyx | 64 +++++++++++-------- .../segmentation/tests/test_felzenszwalb.py | 4 +- 4 files changed, 45 insertions(+), 37 deletions(-) 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)) From 8e6542c86ccba0d1eae86d44ff43fff71ee22679 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 19 May 2015 23:03:46 -0700 Subject: [PATCH 09/20] Do not acquire GIL for heap --- skimage/graph/heap.pxd | 14 +++++++------- skimage/graph/heap.pyx | 22 ++++++++++++---------- skimage/graph/tests/test_heap.py | 2 ++ 3 files changed, 21 insertions(+), 17 deletions(-) diff --git a/skimage/graph/heap.pxd b/skimage/graph/heap.pxd index e31aa150..71b8f2b0 100644 --- a/skimage/graph/heap.pxd +++ b/skimage/graph/heap.pxd @@ -19,13 +19,13 @@ cdef class BinaryHeap: cdef REFERENCE_T *_references cdef REFERENCE_T _popped_ref - cdef void _add_or_remove_level(self, LEVELS_T add_or_remove) - cdef void _update(self) - cdef void _update_one(self, INDEX_T i) - cdef void _remove(self, INDEX_T i) + cdef void _add_or_remove_level(self, LEVELS_T add_or_remove) nogil + cdef void _update(self) nogil + cdef void _update_one(self, INDEX_T i) nogil + cdef void _remove(self, INDEX_T i) nogil - cdef INDEX_T push_fast(self, VALUE_T value, REFERENCE_T reference) - cdef VALUE_T pop_fast(self) + cdef INDEX_T push_fast(self, VALUE_T value, REFERENCE_T reference) nogil + cdef VALUE_T pop_fast(self) nogil cdef class FastUpdateBinaryHeap(BinaryHeap): cdef readonly REFERENCE_T max_reference @@ -35,4 +35,4 @@ cdef class FastUpdateBinaryHeap(BinaryHeap): cdef VALUE_T value_of_fast(self, REFERENCE_T reference) cdef INDEX_T push_if_lower_fast(self, VALUE_T value, - REFERENCE_T reference) + REFERENCE_T reference) nogil diff --git a/skimage/graph/heap.pyx b/skimage/graph/heap.pyx index a4368595..4be48e46 100644 --- a/skimage/graph/heap.pyx +++ b/skimage/graph/heap.pyx @@ -43,7 +43,8 @@ cdef extern from "pyport.h": cdef VALUE_T inf = Py_HUGE_VAL # this is handy -cdef inline INDEX_T index_min(INDEX_T a, INDEX_T b): return a if a <= b else b +cdef inline INDEX_T index_min(INDEX_T a, INDEX_T b) nogil: + return a if a <= b else b cdef class BinaryHeap: @@ -185,7 +186,7 @@ cdef class BinaryHeap: ## C Maintanance methods - cdef void _add_or_remove_level(self, LEVELS_T add_or_remove): + cdef void _add_or_remove_level(self, LEVELS_T add_or_remove) nogil: # init indexing ints cdef INDEX_T i, i1, i2, n @@ -228,7 +229,7 @@ cdef class BinaryHeap: self._update() - cdef void _update(self): + cdef void _update(self) nogil: """Update the full tree from the bottom up. This should be done after resizing. """ @@ -251,7 +252,7 @@ cdef class BinaryHeap: values[ii] = values[i+1] - cdef void _update_one(self, INDEX_T i): + cdef void _update_one(self, INDEX_T i) nogil: """Update the tree for one value.""" # shorter name for values @@ -279,7 +280,7 @@ cdef class BinaryHeap: i = ii-1 - cdef void _remove(self, INDEX_T i1): + cdef void _remove(self, INDEX_T i1) nogil: """Remove a value from the heap. By index.""" cdef LEVELS_T levels = self.levels @@ -314,7 +315,7 @@ cdef class BinaryHeap: ## C Public methods - cdef INDEX_T push_fast(self, VALUE_T value, REFERENCE_T reference): + cdef INDEX_T push_fast(self, VALUE_T value, REFERENCE_T reference) nogil: """The c-method for fast pushing. Returns the index relative to the start of the last level in the heap.""" @@ -338,7 +339,7 @@ cdef class BinaryHeap: return count - cdef VALUE_T pop_fast(self): + cdef VALUE_T pop_fast(self) nogil: """The c-method for fast popping. Returns the minimum value. The reference is put in self._popped_ref""" @@ -543,7 +544,7 @@ cdef class FastUpdateBinaryHeap(BinaryHeap): self._crossref[i] = -1 - cdef void _remove(self, INDEX_T i1): + cdef void _remove(self, INDEX_T i1) nogil: """ Remove a value from the heap. By index. """ cdef LEVELS_T levels = self.levels cdef INDEX_T count = self.count @@ -581,7 +582,7 @@ cdef class FastUpdateBinaryHeap(BinaryHeap): self._update_one(i2) - cdef INDEX_T push_fast(self, VALUE_T value, REFERENCE_T reference): + cdef INDEX_T push_fast(self, VALUE_T value, REFERENCE_T reference) nogil: """The c method for fast pushing. If the reference is already present, will update its value, otherwise @@ -611,7 +612,8 @@ cdef class FastUpdateBinaryHeap(BinaryHeap): self._crossref[reference] = ir return ir - cdef INDEX_T push_if_lower_fast(self, VALUE_T value, REFERENCE_T reference): + cdef INDEX_T push_if_lower_fast(self, VALUE_T value, + REFERENCE_T reference) nogil: """If the reference is already present, will update its value ONLY if the new value is lower than the old one. If the reference is not present, this append it. If a value was appended, self._pushed is diff --git a/skimage/graph/tests/test_heap.py b/skimage/graph/tests/test_heap.py index cf7e3d0b..856b8e30 100644 --- a/skimage/graph/tests/test_heap.py +++ b/skimage/graph/tests/test_heap.py @@ -1,8 +1,10 @@ import time import random import skimage.graph.heap as heap +from skimage._shared.testing import test_parallel +@test_parallel() def test_heap(): _test_heap(100000, True) _test_heap(100000, False) From a7c0b11edac93527789113dbe067984b67e57e8d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 20 May 2015 07:10:51 -0700 Subject: [PATCH 10/20] Do not acquire GIL for median_filter --- skimage/filters/_ctmf.pyx | 52 +++++++++++++++++++++++---------------- 1 file changed, 31 insertions(+), 21 deletions(-) diff --git a/skimage/filters/_ctmf.pyx b/skimage/filters/_ctmf.pyx index 70a59aa3..106a0495 100644 --- a/skimage/filters/_ctmf.pyx +++ b/skimage/filters/_ctmf.pyx @@ -21,8 +21,8 @@ from libc.string cimport memset cdef extern from "../_shared/vectorized_ops.h": - void add16(cnp.uint16_t *dest, cnp.uint16_t *src) - void sub16(cnp.uint16_t *dest, cnp.uint16_t *src) + void add16(cnp.uint16_t *dest, cnp.uint16_t *src) nogil + void sub16(cnp.uint16_t *dest, cnp.uint16_t *src) nogil ############################################################################## @@ -174,14 +174,14 @@ cdef Histograms *allocate_histograms(Py_ssize_t rows, Py_ssize_t percent, cnp.uint8_t *data, cnp.uint8_t *mask, - cnp.uint8_t *output): + cnp.uint8_t *output) nogil: cdef: Py_ssize_t adjusted_stripe_length = columns + 2*radius + 1 Py_ssize_t memory_size void *ptr Histograms *ph Py_ssize_t roundoff - Py_ssize_t a + Py_ssize_t a, a_2 SCoord *psc memory_size = (adjusted_stripe_length * @@ -292,7 +292,7 @@ cdef Histograms *allocate_histograms(Py_ssize_t rows, # free_histograms - frees the Histograms structure # ############################################################################ -cdef void free_histograms(Histograms *ph): +cdef void free_histograms(Histograms *ph) nogil: free(ph.memory) ############################################################################ @@ -301,7 +301,7 @@ cdef void free_histograms(Histograms *ph): # ############################################################################ -cdef void set_stride(Histograms *ph, SCoord *psc): +cdef void set_stride(Histograms *ph, SCoord *psc) nogil: psc.stride = psc.x * ph.col_stride + psc.y * ph.row_stride ############################################################################ @@ -326,17 +326,23 @@ cdef void set_stride(Histograms *ph, SCoord *psc): # a column that is "radius" to the left. # ############################################################################ -cdef inline Py_ssize_t tl_br_colidx(Histograms *ph, Py_ssize_t colidx): +@cython.cdivision(True) +cdef inline Py_ssize_t tl_br_colidx(Histograms *ph, Py_ssize_t colidx) nogil: return (colidx + 3*ph.radius + ph.current_row) % ph.stripe_length -cdef inline Py_ssize_t tr_bl_colidx(Histograms *ph, Py_ssize_t colidx): +@cython.cdivision(True) +cdef inline Py_ssize_t tr_bl_colidx(Histograms *ph, Py_ssize_t colidx) nogil: return (colidx + 3*ph.radius + ph.row_count-ph.current_row) % \ ph.stripe_length -cdef inline Py_ssize_t leading_edge_colidx(Histograms *ph, Py_ssize_t colidx): +@cython.cdivision(True) +cdef inline Py_ssize_t leading_edge_colidx(Histograms *ph, + Py_ssize_t colidx) nogil: return (colidx + 5*ph.radius) % ph.stripe_length -cdef inline Py_ssize_t trailing_edge_colidx(Histograms *ph, Py_ssize_t colidx): +@cython.cdivision(True) +cdef inline Py_ssize_t trailing_edge_colidx(Histograms *ph, + Py_ssize_t colidx) nogil: return (colidx + 3*ph.radius - 1) % ph.stripe_length ############################################################################ @@ -348,7 +354,8 @@ cdef inline Py_ssize_t trailing_edge_colidx(Histograms *ph, Py_ssize_t colidx): # colidx - the index of the column to add # ############################################################################ -cdef inline void accumulate_coarse_histogram(Histograms *ph, Py_ssize_t colidx): +cdef inline void accumulate_coarse_histogram(Histograms *ph, + Py_ssize_t colidx) nogil: cdef Py_ssize_t offset offset = tr_bl_colidx(ph, colidx) @@ -370,7 +377,8 @@ cdef inline void accumulate_coarse_histogram(Histograms *ph, Py_ssize_t colidx): # for a given column # ############################################################################ -cdef inline void deaccumulate_coarse_histogram(Histograms *ph, Py_ssize_t colidx): +cdef inline void deaccumulate_coarse_histogram(Histograms *ph, + Py_ssize_t colidx) nogil: cdef Py_ssize_t offset # # The trailing diagonals don't appear until here @@ -401,7 +409,7 @@ cdef inline void deaccumulate_coarse_histogram(Histograms *ph, Py_ssize_t colidx ############################################################################ cdef inline void accumulate_fine_histogram(Histograms *ph, Py_ssize_t colidx, - Py_ssize_t fineidx): + Py_ssize_t fineidx) nogil: cdef: Py_ssize_t fineoffset = fineidx * 16 Py_ssize_t offset @@ -425,7 +433,7 @@ cdef inline void accumulate_fine_histogram(Histograms *ph, ############################################################################ cdef inline void deaccumulate_fine_histogram(Histograms *ph, Py_ssize_t colidx, - Py_ssize_t fineidx): + Py_ssize_t fineidx) nogil: cdef: Py_ssize_t fineoffset = fineidx * 16 Py_ssize_t offset @@ -455,7 +463,7 @@ cdef inline void deaccumulate_fine_histogram(Histograms *ph, # ############################################################################ -cdef inline void accumulate(Histograms *ph): +cdef inline void accumulate(Histograms *ph) nogil: cdef cnp.int32_t accumulator accumulate_coarse_histogram(ph, ph.current_column) deaccumulate_coarse_histogram(ph, ph.current_column) @@ -480,7 +488,7 @@ cdef inline void accumulate(Histograms *ph): # to choose remains to be done. ############################################################################ -cdef inline void update_fine(Histograms *ph, Py_ssize_t fineidx): +cdef inline void update_fine(Histograms *ph, Py_ssize_t fineidx) nogil: cdef: Py_ssize_t first_update_column = ph.last_update_column[fineidx]+1 Py_ssize_t update_limit = ph.current_column+1 @@ -507,7 +515,7 @@ cdef inline void update_histogram(Histograms *ph, HistogramPiece *hist_piece, pixel_count_t *pixel_count, SCoord *last_coord, - SCoord *coord): + SCoord *coord) nogil: cdef: Py_ssize_t current_column = ph.current_column Py_ssize_t current_row = ph.current_row @@ -548,7 +556,7 @@ cdef inline void update_histogram(Histograms *ph, # update_current_location - update the histograms at the current location # ############################################################################ -cdef inline void update_current_location(Histograms *ph): +cdef inline void update_current_location(Histograms *ph) nogil: cdef: Py_ssize_t current_column = ph.current_column Py_ssize_t radius = ph.radius @@ -597,7 +605,7 @@ cdef inline void update_current_location(Histograms *ph): # ############################################################################ -cdef inline cnp.uint8_t find_median(Histograms *ph): +cdef inline cnp.uint8_t find_median(Histograms *ph) nogil: cdef: Py_ssize_t pixels_below # of pixels below the median Py_ssize_t i @@ -652,12 +660,13 @@ cdef int c_median_filter(Py_ssize_t rows, Py_ssize_t percent, cnp.uint8_t *data, cnp.uint8_t *mask, - cnp.uint8_t *output): + cnp.uint8_t *output) nogil: cdef: Histograms *ph Histogram *phistogram Py_ssize_t row Py_ssize_t col + Py_ssize_t end_col Py_ssize_t i Py_ssize_t top_left_off Py_ssize_t top_right_off @@ -706,7 +715,8 @@ cdef int c_median_filter(Py_ssize_t rows, # Update locations and coarse accumulator for the octagon # for points before 0 # - for col in range(-radius, 0 if row >= 0 else columns+radius): + end_col = 0 if row >= 0 else columns + radius + for col in range(-radius, end_col): ph.current_column = col ph.current_stride = row * row_stride + col * col_stride update_current_location(ph) From aa38f2485b73dfa84ea4cc90d1329c607e9cd16d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 20 May 2015 07:26:24 -0700 Subject: [PATCH 11/20] Do not acquire GIL for rank filters --- skimage/filters/rank/bilateral_cy.pyx | 6 +- skimage/filters/rank/core_cy.pxd | 6 +- skimage/filters/rank/core_cy.pyx | 212 ++++++++++++------------ skimage/filters/rank/generic_cy.pyx | 38 ++--- skimage/filters/rank/percentile_cy.pyx | 18 +- skimage/filters/rank/tests/test_rank.py | 5 +- 6 files changed, 145 insertions(+), 140 deletions(-) diff --git a/skimage/filters/rank/bilateral_cy.pyx b/skimage/filters/rank/bilateral_cy.pyx index a008e682..df37f6c0 100644 --- a/skimage/filters/rank/bilateral_cy.pyx +++ b/skimage/filters/rank/bilateral_cy.pyx @@ -14,7 +14,7 @@ cdef inline void _kernel_mean(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t i cdef Py_ssize_t bilat_pop = 0 @@ -38,7 +38,7 @@ cdef inline void _kernel_pop(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t i cdef Py_ssize_t bilat_pop = 0 @@ -57,7 +57,7 @@ cdef inline void _kernel_sum(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t i cdef Py_ssize_t bilat_pop = 0 diff --git a/skimage/filters/rank/core_cy.pxd b/skimage/filters/rank/core_cy.pxd index 795f0819..a0b07654 100644 --- a/skimage/filters/rank/core_cy.pxd +++ b/skimage/filters/rank/core_cy.pxd @@ -11,13 +11,13 @@ ctypedef fused dtype_t_out: double_t -cdef dtype_t _max(dtype_t a, dtype_t b) -cdef dtype_t _min(dtype_t a, dtype_t b) +cdef dtype_t _max(dtype_t a, dtype_t b) nogil +cdef dtype_t _min(dtype_t a, dtype_t b) nogil cdef void _core(void kernel(dtype_t_out*, Py_ssize_t, Py_ssize_t*, double, dtype_t, Py_ssize_t, Py_ssize_t, double, - double, Py_ssize_t, Py_ssize_t), + double, Py_ssize_t, Py_ssize_t) nogil, dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, diff --git a/skimage/filters/rank/core_cy.pyx b/skimage/filters/rank/core_cy.pyx index 822c7251..e7464063 100644 --- a/skimage/filters/rank/core_cy.pyx +++ b/skimage/filters/rank/core_cy.pyx @@ -9,29 +9,29 @@ cimport numpy as cnp from libc.stdlib cimport malloc, free -cdef inline dtype_t _max(dtype_t a, dtype_t b): +cdef inline dtype_t _max(dtype_t a, dtype_t b) nogil: return a if a >= b else b -cdef inline dtype_t _min(dtype_t a, dtype_t b): +cdef inline dtype_t _min(dtype_t a, dtype_t b) nogil: return a if a <= b else b cdef inline void histogram_increment(Py_ssize_t* histo, double* pop, - dtype_t value): + dtype_t value) nogil: histo[value] += 1 pop[0] += 1 cdef inline void histogram_decrement(Py_ssize_t* histo, double* pop, - dtype_t value): + dtype_t value) nogil: histo[value] -= 1 pop[0] -= 1 cdef inline char is_in_mask(Py_ssize_t rows, Py_ssize_t cols, Py_ssize_t r, Py_ssize_t c, - char* mask): + char* mask) nogil: """Check whether given coordinate is within image and mask is true.""" if r < 0 or r > rows - 1 or c < 0 or c > cols - 1: return 0 @@ -44,7 +44,7 @@ cdef inline char is_in_mask(Py_ssize_t rows, Py_ssize_t cols, cdef void _core(void kernel(dtype_t_out*, Py_ssize_t, Py_ssize_t*, double, dtype_t, Py_ssize_t, Py_ssize_t, double, - double, Py_ssize_t, Py_ssize_t), + double, Py_ssize_t, Py_ssize_t) nogil, dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, @@ -123,123 +123,125 @@ cdef void _core(void kernel(dtype_t_out*, Py_ssize_t, Py_ssize_t*, double, t = np.vstack((np.zeros((1, selem.shape[1])), selem)) cdef unsigned char[:, :] t_n = (np.diff(t, axis=0) > 0).view(np.uint8) - for r in range(srows): - for c in range(scols): - if t_e[r, c]: - se_e_r[num_se_e] = r - centre_r - se_e_c[num_se_e] = c - centre_c - num_se_e += 1 - if t_w[r, c]: - se_w_r[num_se_w] = r - centre_r - se_w_c[num_se_w] = c - centre_c - num_se_w += 1 - if t_n[r, c]: - se_n_r[num_se_n] = r - centre_r - se_n_c[num_se_n] = c - centre_c - num_se_n += 1 - if t_s[r, c]: - se_s_r[num_se_s] = r - centre_r - se_s_c[num_se_s] = c - centre_c - num_se_s += 1 + with nogil: - for r in range(srows): - for c in range(scols): - rr = r - centre_r - cc = c - centre_c - if selem[r, c]: + for r in range(srows): + for c in range(scols): + if t_e[r, c]: + se_e_r[num_se_e] = r - centre_r + se_e_c[num_se_e] = c - centre_c + num_se_e += 1 + if t_w[r, c]: + se_w_r[num_se_w] = r - centre_r + se_w_c[num_se_w] = c - centre_c + num_se_w += 1 + if t_n[r, c]: + se_n_r[num_se_n] = r - centre_r + se_n_c[num_se_n] = c - centre_c + num_se_n += 1 + if t_s[r, c]: + se_s_r[num_se_s] = r - centre_r + se_s_c[num_se_s] = c - centre_c + num_se_s += 1 + + for r in range(srows): + for c in range(scols): + rr = r - centre_r + cc = c - centre_c + if selem[r, c]: + if is_in_mask(rows, cols, rr, cc, mask_data): + histogram_increment(histo, &pop, image[rr, cc]) + + r = 0 + c = 0 + kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin, mid_bin, + p0, p1, s0, s1) + + # main loop + r = 0 + for even_row in range(0, rows, 2): + + # ---> west to east + for c in range(1, cols): + for s in range(num_se_e): + rr = r + se_e_r[s] + cc = c + se_e_c[s] + if is_in_mask(rows, cols, rr, cc, mask_data): + histogram_increment(histo, &pop, image[rr, cc]) + + for s in range(num_se_w): + rr = r + se_w_r[s] + cc = c + se_w_c[s] - 1 + if is_in_mask(rows, cols, rr, cc, mask_data): + histogram_decrement(histo, &pop, image[rr, cc]) + + kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin, + mid_bin, p0, p1, s0, s1) + + r += 1 # pass to the next row + if r >= rows: + break + + # ---> north to south + for s in range(num_se_s): + rr = r + se_s_r[s] + cc = c + se_s_c[s] if is_in_mask(rows, cols, rr, cc, mask_data): histogram_increment(histo, &pop, image[rr, cc]) - r = 0 - c = 0 - kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin, mid_bin, - p0, p1, s0, s1) - - # main loop - r = 0 - for even_row in range(0, rows, 2): - - # ---> west to east - for c in range(1, cols): - for s in range(num_se_e): - rr = r + se_e_r[s] - cc = c + se_e_c[s] - if is_in_mask(rows, cols, rr, cc, mask_data): - histogram_increment(histo, &pop, image[rr, cc]) - - for s in range(num_se_w): - rr = r + se_w_r[s] - cc = c + se_w_c[s] - 1 + for s in range(num_se_n): + rr = r + se_n_r[s] - 1 + cc = c + se_n_c[s] if is_in_mask(rows, cols, rr, cc, mask_data): histogram_decrement(histo, &pop, image[rr, cc]) kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin, mid_bin, p0, p1, s0, s1) - r += 1 # pass to the next row - if r >= rows: - break + # ---> east to west + for c in range(cols - 2, -1, -1): + for s in range(num_se_w): + rr = r + se_w_r[s] + cc = c + se_w_c[s] + if is_in_mask(rows, cols, rr, cc, mask_data): + histogram_increment(histo, &pop, image[rr, cc]) - # ---> north to south - for s in range(num_se_s): - rr = r + se_s_r[s] - cc = c + se_s_c[s] - if is_in_mask(rows, cols, rr, cc, mask_data): - histogram_increment(histo, &pop, image[rr, cc]) + for s in range(num_se_e): + rr = r + se_e_r[s] + cc = c + se_e_c[s] + 1 + if is_in_mask(rows, cols, rr, cc, mask_data): + histogram_decrement(histo, &pop, image[rr, cc]) - for s in range(num_se_n): - rr = r + se_n_r[s] - 1 - cc = c + se_n_c[s] - if is_in_mask(rows, cols, rr, cc, mask_data): - histogram_decrement(histo, &pop, image[rr, cc]) + kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin, + mid_bin, p0, p1, s0, s1) - kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin, - mid_bin, p0, p1, s0, s1) + r += 1 # pass to the next row + if r >= rows: + break - # ---> east to west - for c in range(cols - 2, -1, -1): - for s in range(num_se_w): - rr = r + se_w_r[s] - cc = c + se_w_c[s] + # ---> north to south + for s in range(num_se_s): + rr = r + se_s_r[s] + cc = c + se_s_c[s] if is_in_mask(rows, cols, rr, cc, mask_data): histogram_increment(histo, &pop, image[rr, cc]) - for s in range(num_se_e): - rr = r + se_e_r[s] - cc = c + se_e_c[s] + 1 + for s in range(num_se_n): + rr = r + se_n_r[s] - 1 + cc = c + se_n_c[s] if is_in_mask(rows, cols, rr, cc, mask_data): histogram_decrement(histo, &pop, image[rr, cc]) - kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin, - mid_bin, p0, p1, s0, s1) + kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], + max_bin, mid_bin, p0, p1, s0, s1) - r += 1 # pass to the next row - if r >= rows: - break - - # ---> north to south - for s in range(num_se_s): - rr = r + se_s_r[s] - cc = c + se_s_c[s] - if is_in_mask(rows, cols, rr, cc, mask_data): - histogram_increment(histo, &pop, image[rr, cc]) - - for s in range(num_se_n): - rr = r + se_n_r[s] - 1 - cc = c + se_n_c[s] - if is_in_mask(rows, cols, rr, cc, mask_data): - histogram_decrement(histo, &pop, image[rr, cc]) - - kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], - max_bin, mid_bin, p0, p1, s0, s1) - - # release memory allocated by malloc - free(se_e_r) - free(se_e_c) - free(se_w_r) - free(se_w_c) - free(se_n_r) - free(se_n_c) - free(se_s_r) - free(se_s_c) - free(histo) + # release memory allocated by malloc + free(se_e_r) + free(se_e_c) + free(se_w_r) + free(se_w_c) + free(se_n_r) + free(se_n_c) + free(se_s_r) + free(se_s_c) + free(histo) diff --git a/skimage/filters/rank/generic_cy.pyx b/skimage/filters/rank/generic_cy.pyx index 3b28005a..1356f369 100644 --- a/skimage/filters/rank/generic_cy.pyx +++ b/skimage/filters/rank/generic_cy.pyx @@ -14,7 +14,7 @@ cdef inline void _kernel_autolevel(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t i, imin, imax, delta @@ -41,7 +41,7 @@ cdef inline void _kernel_bottomhat(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t i @@ -59,7 +59,7 @@ cdef inline void _kernel_equalize(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t i cdef Py_ssize_t sum = 0 @@ -79,7 +79,7 @@ cdef inline void _kernel_gradient(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t i, imin, imax @@ -102,7 +102,7 @@ cdef inline void _kernel_maximum(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t i @@ -120,7 +120,7 @@ cdef inline void _kernel_mean(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t i cdef Py_ssize_t mean = 0 @@ -138,7 +138,7 @@ cdef inline void _kernel_subtract_mean(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t i cdef Py_ssize_t mean = 0 @@ -156,7 +156,7 @@ cdef inline void _kernel_median(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t i cdef double sum = pop / 2.0 @@ -177,7 +177,7 @@ cdef inline void _kernel_minimum(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t i @@ -195,7 +195,7 @@ cdef inline void _kernel_modal(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t hmax = 0, imax = 0 @@ -217,7 +217,7 @@ cdef inline void _kernel_enhance_contrast(dtype_t_out* out, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, Py_ssize_t s0, - Py_ssize_t s1): + Py_ssize_t s1) nogil: cdef Py_ssize_t i, imin, imax @@ -243,7 +243,7 @@ cdef inline void _kernel_pop(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: out[0] = pop @@ -253,7 +253,7 @@ cdef inline void _kernel_sum(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t i cdef Py_ssize_t sum = 0 @@ -271,7 +271,7 @@ cdef inline void _kernel_threshold(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t i cdef Py_ssize_t mean = 0 @@ -289,7 +289,7 @@ cdef inline void _kernel_tophat(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t i @@ -307,7 +307,7 @@ cdef inline void _kernel_noise_filter(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t i cdef Py_ssize_t min_i @@ -334,7 +334,7 @@ cdef inline void _kernel_entropy(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t i cdef double e, p @@ -354,7 +354,7 @@ cdef inline void _kernel_otsu(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t i cdef Py_ssize_t max_i cdef double P, mu1, mu2, q1, new_q1, sigma_b, max_sigma_b @@ -394,7 +394,7 @@ cdef inline void _kernel_win_hist(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t i cdef Py_ssize_t max_i cdef double scale diff --git a/skimage/filters/rank/percentile_cy.pyx b/skimage/filters/rank/percentile_cy.pyx index 4d90555f..b230b660 100644 --- a/skimage/filters/rank/percentile_cy.pyx +++ b/skimage/filters/rank/percentile_cy.pyx @@ -12,7 +12,7 @@ cdef inline void _kernel_autolevel(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t i, imin, imax, sum, delta @@ -46,7 +46,7 @@ cdef inline void _kernel_gradient(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t i, imin, imax, sum, delta @@ -75,7 +75,7 @@ cdef inline void _kernel_mean(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t i, sum, mean, n @@ -101,7 +101,7 @@ cdef inline void _kernel_sum(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t i, sum, sum_g, n @@ -128,7 +128,7 @@ cdef inline void _kernel_subtract_mean(dtype_t_out* out, Py_ssize_t odepth, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, Py_ssize_t s0, - Py_ssize_t s1): + Py_ssize_t s1) nogil: cdef Py_ssize_t i, sum, mean, n @@ -156,7 +156,7 @@ cdef inline void _kernel_enhance_contrast(dtype_t_out* out, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, Py_ssize_t s0, - Py_ssize_t s1): + Py_ssize_t s1) nogil: cdef Py_ssize_t i, imin, imax, sum, delta @@ -191,7 +191,7 @@ cdef inline void _kernel_percentile(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t i cdef Py_ssize_t sum = 0 @@ -216,7 +216,7 @@ cdef inline void _kernel_pop(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef Py_ssize_t i, sum, n @@ -237,7 +237,7 @@ cdef inline void _kernel_threshold(dtype_t_out* out, Py_ssize_t odepth, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, double p0, double p1, - Py_ssize_t s0, Py_ssize_t s1): + Py_ssize_t s0, Py_ssize_t s1) nogil: cdef int i cdef Py_ssize_t sum = 0 diff --git a/skimage/filters/rank/tests/test_rank.py b/skimage/filters/rank/tests/test_rank.py index 48628432..021f3d36 100644 --- a/skimage/filters/rank/tests/test_rank.py +++ b/skimage/filters/rank/tests/test_rank.py @@ -8,6 +8,7 @@ from skimage import data, util, morphology from skimage.morphology import grey, disk from skimage.filters import rank from skimage._shared._warnings import expected_warnings +from skimage._shared.testing import test_parallel def test_all(): @@ -15,7 +16,9 @@ def test_all(): check_all() +@test_parallel() def check_all(): + np.random.seed(0) image = np.random.rand(25, 25) selem = morphology.disk(1) refs = np.load(os.path.join(skimage.data_dir, "rank_filter_tests.npz")) @@ -489,7 +492,7 @@ def test_entropy(): assert(np.max(rank.entropy(data, selem)) == 12) # make sure output is of dtype double - with expected_warnings(['Bitdepth of 11']): + with expected_warnings(['Bitdepth of 11']): out = rank.entropy(data, np.ones((16, 16), dtype=np.uint8)) assert out.dtype == np.double From c52c7f4bfb1133456b225e25a4ada94ad25e1e91 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 20 May 2015 07:38:47 -0700 Subject: [PATCH 12/20] Do not acquire GIL for ORB --- skimage/feature/orb_cy.pyx | 35 ++++++++++++++++--------------- skimage/feature/tests/test_orb.py | 2 ++ 2 files changed, 20 insertions(+), 17 deletions(-) diff --git a/skimage/feature/orb_cy.pyx b/skimage/feature/orb_cy.pyx index 49c4e650..b3140aa5 100644 --- a/skimage/feature/orb_cy.pyx +++ b/skimage/feature/orb_cy.pyx @@ -30,27 +30,28 @@ def _orb_loop(double[:, ::1] image, Py_ssize_t[:, ::1] keypoints, cdef signed char[:, ::1] cpos0 = POS0 cdef signed char[:, ::1] cpos1 = POS1 - for i in range(descriptors.shape[0]): + with nogil: + for i in range(descriptors.shape[0]): - angle = orientations[i] - sin_a = sin(angle) - cos_a = cos(angle) + angle = orientations[i] + sin_a = sin(angle) + cos_a = cos(angle) - kr = keypoints[i, 0] - kc = keypoints[i, 1] + kr = keypoints[i, 0] + kc = keypoints[i, 1] - for j in range(descriptors.shape[1]): - pr0 = cpos0[j, 0] - pc0 = cpos0[j, 1] - pr1 = cpos1[j, 0] - pc1 = cpos1[j, 1] + for j in range(descriptors.shape[1]): + pr0 = cpos0[j, 0] + pc0 = cpos0[j, 1] + pr1 = cpos1[j, 0] + pc1 = cpos1[j, 1] - spr0 = round(sin_a * pr0 + cos_a * pc0) - spc0 = round(cos_a * pr0 - sin_a * pc0) - spr1 = round(sin_a * pr1 + cos_a * pc1) - spc1 = round(cos_a * pr1 - sin_a * pc1) + spr0 = round(sin_a * pr0 + cos_a * pc0) + spc0 = round(cos_a * pr0 - sin_a * pc0) + spr1 = round(sin_a * pr1 + cos_a * pc1) + spc1 = round(cos_a * pr1 - sin_a * pc1) - if image[kr + spr0, kc + spc0] < image[kr + spr1, kc + spc1]: - descriptors[i, j] = True + if image[kr + spr0, kc + spc0] < image[kr + spr1, kc + spc1]: + descriptors[i, j] = True return np.asarray(descriptors) diff --git a/skimage/feature/tests/test_orb.py b/skimage/feature/tests/test_orb.py index 5bed2bf3..8895d857 100644 --- a/skimage/feature/tests/test_orb.py +++ b/skimage/feature/tests/test_orb.py @@ -3,11 +3,13 @@ from numpy.testing import assert_equal, assert_almost_equal, run_module_suite from skimage.feature import ORB from skimage import data from skimage.color import rgb2gray +from skimage._shared.testing import test_parallel img = rgb2gray(data.lena()) +@test_parallel() def test_keypoints_orb_desired_no_of_keypoints(): detector_extractor = ORB(n_keypoints=10, fast_n=12, fast_threshold=0.20) detector_extractor.detect(img) From a489b85b22d077a6e78fbf099faee00a77ca54a2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 20 May 2015 07:45:04 -0700 Subject: [PATCH 13/20] Do not acquire GIL for corner detectors --- skimage/feature/corner_cy.pyx | 137 ++++++++++++++------------- skimage/feature/tests/test_corner.py | 4 + 2 files changed, 75 insertions(+), 66 deletions(-) diff --git a/skimage/feature/corner_cy.pyx b/skimage/feature/corner_cy.pyx index cef104e0..4f5c32c6 100644 --- a/skimage/feature/corner_cy.pyx +++ b/skimage/feature/corner_cy.pyx @@ -5,7 +5,7 @@ import numpy as np cimport numpy as cnp from libc.float cimport DBL_MAX -from libc.math cimport atan2 +from libc.math cimport atan2, fabs from ..util import img_as_float, pad from ..color import rgb2grey @@ -67,27 +67,30 @@ def corner_moravec(image, Py_ssize_t window_size=1): cdef double msum, min_msum cdef Py_ssize_t r, c, br, bc, mr, mc, a, b - for r in range(2 * window_size, rows - 2 * window_size): - for c in range(2 * window_size, cols - 2 * window_size): - min_msum = DBL_MAX - for br in range(r - window_size, r + window_size + 1): - for bc in range(c - window_size, c + window_size + 1): - if br != r and bc != c: - msum = 0 - for mr in range(- window_size, window_size + 1): - for mc in range(- window_size, window_size + 1): - msum += (cimage[r + mr, c + mc] - - cimage[br + mr, bc + mc]) ** 2 - min_msum = min(msum, min_msum) - out[r, c] = min_msum + with nogil: + for r in range(2 * window_size, rows - 2 * window_size): + for c in range(2 * window_size, cols - 2 * window_size): + min_msum = DBL_MAX + for br in range(r - window_size, r + window_size + 1): + for bc in range(c - window_size, c + window_size + 1): + if br != r and bc != c: + msum = 0 + for mr in range(- window_size, window_size + 1): + for mc in range(- window_size, window_size + 1): + msum += (cimage[r + mr, c + mc] + - cimage[br + mr, bc + mc]) ** 2 + min_msum = min(msum, min_msum) + + out[r, c] = min_msum return np.asarray(out) cdef inline double _corner_fast_response(double curr_pixel, double* circle_intensities, - signed char* bins, signed char state, char n): + signed char* bins, signed char state, + char n) nogil: cdef char consecutive_count = 0 cdef double curr_response cdef Py_ssize_t l, m @@ -97,7 +100,7 @@ cdef inline double _corner_fast_response(double curr_pixel, if consecutive_count == n: curr_response = 0 for m in range(16): - curr_response += abs(circle_intensities[m] - curr_pixel) + curr_response += fabs(circle_intensities[m] - curr_pixel) return curr_response else: consecutive_count = 0 @@ -124,49 +127,50 @@ def _corner_fast(double[:, ::1] image, signed char n, double threshold): cdef double curr_response - for i in range(3, rows - 3): - for j in range(3, cols - 3): + with nogil: + for i in range(3, rows - 3): + for j in range(3, cols - 3): - curr_pixel = image[i, j] - lower_threshold = curr_pixel - threshold - upper_threshold = curr_pixel + threshold + curr_pixel = image[i, j] + lower_threshold = curr_pixel - threshold + upper_threshold = curr_pixel + threshold - for k in range(16): - circle_intensities[k] = image[i + rp[k], j + cp[k]] - if circle_intensities[k] > upper_threshold: - # Brighter pixel - bins[k] = 'b' - elif circle_intensities[k] < lower_threshold: - # Darker pixel - bins[k] = 'd' - else: - # Similar pixel - bins[k] = 's' + for k in range(16): + circle_intensities[k] = image[i + rp[k], j + cp[k]] + if circle_intensities[k] > upper_threshold: + # Brighter pixel + bins[k] = 'b' + elif circle_intensities[k] < lower_threshold: + # Darker pixel + bins[k] = 'd' + else: + # Similar pixel + bins[k] = 's' - # High speed test for n >= 12 - if n >= 12: - speed_sum_b = 0 - speed_sum_d = 0 - for k in range(0, 16, 4): - if bins[k] == 'b': - speed_sum_b += 1 - elif bins[k] == 'd': - speed_sum_d += 1 - if speed_sum_d < 3 and speed_sum_b < 3: - continue + # High speed test for n >= 12 + if n >= 12: + speed_sum_b = 0 + speed_sum_d = 0 + for k in range(0, 16, 4): + if bins[k] == 'b': + speed_sum_b += 1 + elif bins[k] == 'd': + speed_sum_d += 1 + if speed_sum_d < 3 and speed_sum_b < 3: + continue - # Test for bright pixels - curr_response = \ - _corner_fast_response(curr_pixel, circle_intensities, - bins, 'b', n) - - # Test for dark pixels - if curr_response == 0: + # Test for bright pixels curr_response = \ _corner_fast_response(curr_pixel, circle_intensities, - bins, 'd', n) + bins, 'b', n) - corner_response[i, j] = curr_response + # Test for dark pixels + if curr_response == 0: + curr_response = \ + _corner_fast_response(curr_pixel, circle_intensities, + bins, 'd', n) + + corner_response[i, j] = curr_response return np.asarray(corner_response) @@ -254,22 +258,23 @@ def corner_orientations(image, Py_ssize_t[:, :] corners, mask): cdef double curr_pixel cdef double m01, m10, m01_tmp - for i in range(corners.shape[0]): - r0 = corners[i, 0] - c0 = corners[i, 1] + with nogil: + for i in range(corners.shape[0]): + r0 = corners[i, 0] + c0 = corners[i, 1] - m01 = 0 - m10 = 0 + m01 = 0 + m10 = 0 - for r in range(mrows): - m01_tmp = 0 - for c in range(mcols): - if cmask[r, c]: - curr_pixel = cimage[r0 + r, c0 + c] - m10 += curr_pixel * (c - mcols2) - m01_tmp += curr_pixel - m01 += m01_tmp * (r - mrows2) + for r in range(mrows): + m01_tmp = 0 + for c in range(mcols): + if cmask[r, c]: + curr_pixel = cimage[r0 + r, c0 + c] + m10 += curr_pixel * (c - mcols2) + m01_tmp += curr_pixel + m01 += m01_tmp * (r - mrows2) - orientations[i] = atan2(m01, m10) + orientations[i] = atan2(m01, m10) return np.asarray(orientations) diff --git a/skimage/feature/tests/test_corner.py b/skimage/feature/tests/test_corner.py index a9f5bb91..aa7e1409 100644 --- a/skimage/feature/tests/test_corner.py +++ b/skimage/feature/tests/test_corner.py @@ -6,6 +6,7 @@ from skimage import data from skimage import img_as_float from skimage.color import rgb2gray from skimage.morphology import octagon +from skimage._shared.testing import test_parallel from skimage.feature import (corner_moravec, corner_harris, corner_shi_tomasi, corner_subpix, peak_local_max, corner_peaks, @@ -99,6 +100,7 @@ def test_hessian_matrix_det(): assert_almost_equal(det, 0, decimal = 3) +@test_parallel() def test_square_image(): im = np.zeros((50, 50)).astype(float) im[:25, :25] = 1. @@ -280,6 +282,7 @@ def test_corner_fast_image_unsupported_error(): assert_raises(ValueError, corner_fast, img) +@test_parallel() def test_corner_fast_lena(): img = rgb2gray(data.astronaut()) expected = np.array([[101, 198], @@ -335,6 +338,7 @@ def test_corner_orientations_even_shape_error(): np.asarray([[7, 7]]), np.ones((4, 4))) +@test_parallel() def test_corner_orientations_lena(): img = rgb2gray(data.lena()) corners = corner_peaks(corner_fast(img, 11, 0.35)) From bca458bc1ce41d97177877faf157a8ad1ab59ecc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 20 May 2015 09:04:58 -0700 Subject: [PATCH 14/20] Do not acquire GIL for CENSURE --- skimage/feature/censure_cy.pyx | 86 ++++++++++++++------------- skimage/feature/tests/test_censure.py | 2 + 2 files changed, 46 insertions(+), 42 deletions(-) diff --git a/skimage/feature/censure_cy.pyx b/skimage/feature/censure_cy.pyx index 1c352fde..f8c712f8 100644 --- a/skimage/feature/censure_cy.pyx +++ b/skimage/feature/censure_cy.pyx @@ -18,55 +18,57 @@ def _censure_dob_loop(Py_ssize_t n, cdef Py_ssize_t n2 = 2 * n cdef double total_weight = inner_weight + outer_weight - # top-left pixel - inner = (integral_img[n2 + n, n2 + n] - + integral_img[n2 - n - 1, n2 - n - 1] - - integral_img[n2 + n, n2 - n - 1] - - integral_img[n2 - n - 1, n2 + n]) + with nogil: - outer = integral_img[2 * n2, 2 * n2] + # top-left pixel + inner = (integral_img[n2 + n, n2 + n] + + integral_img[n2 - n - 1, n2 - n - 1] + - integral_img[n2 + n, n2 - n - 1] + - integral_img[n2 - n - 1, n2 + n]) - filtered_image[n2, n2] = (outer_weight * outer - - total_weight * inner) + outer = integral_img[2 * n2, 2 * n2] - # left column - for i in range(n2 + 1, integral_img.shape[0] - n2): - inner = (integral_img[i + n, n2 + n] - + integral_img[i - n - 1, n2 - n - 1] - - integral_img[i + n, n2 - n - 1] - - integral_img[i - n - 1, n2 + n]) + filtered_image[n2, n2] = (outer_weight * outer + - total_weight * inner) - outer = (integral_img[i + n2, 2 * n2] - - integral_img[i - n2 - 1, 2 * n2]) + # left column + for i in range(n2 + 1, integral_img.shape[0] - n2): + inner = (integral_img[i + n, n2 + n] + + integral_img[i - n - 1, n2 - n - 1] + - integral_img[i + n, n2 - n - 1] + - integral_img[i - n - 1, n2 + n]) - filtered_image[i, n2] = (outer_weight * outer - - total_weight * inner) + outer = (integral_img[i + n2, 2 * n2] + - integral_img[i - n2 - 1, 2 * n2]) - # top row - for j in range(n2 + 1, integral_img.shape[1] - n2): - inner = (integral_img[n2 + n, j + n] - + integral_img[n2 - n - 1, j - n - 1] - - integral_img[n2 + n, j - n - 1] - - integral_img[n2 - n - 1, j + n]) + filtered_image[i, n2] = (outer_weight * outer + - total_weight * inner) - outer = (integral_img[2 * n2, j + n2] - - integral_img[2 * n2, j - n2 - 1]) - - filtered_image[n2, j] = (outer_weight * outer - - total_weight * inner) - - # remaining block - for i in range(n2 + 1, integral_img.shape[0] - n2): + # top row for j in range(n2 + 1, integral_img.shape[1] - n2): - inner = (integral_img[i + n, j + n] - + integral_img[i - n - 1, j - n - 1] - - integral_img[i + n, j - n - 1] - - integral_img[i - n - 1, j + n]) + inner = (integral_img[n2 + n, j + n] + + integral_img[n2 - n - 1, j - n - 1] + - integral_img[n2 + n, j - n - 1] + - integral_img[n2 - n - 1, j + n]) - outer = (integral_img[i + n2, j + n2] - + integral_img[i - n2 - 1, j - n2 - 1] - - integral_img[i + n2, j - n2 - 1] - - integral_img[i - n2 - 1, j + n2]) + outer = (integral_img[2 * n2, j + n2] + - integral_img[2 * n2, j - n2 - 1]) - filtered_image[i, j] = (outer_weight * outer - - total_weight * inner) + filtered_image[n2, j] = (outer_weight * outer + - total_weight * inner) + + # remaining block + for i in range(n2 + 1, integral_img.shape[0] - n2): + for j in range(n2 + 1, integral_img.shape[1] - n2): + inner = (integral_img[i + n, j + n] + + integral_img[i - n - 1, j - n - 1] + - integral_img[i + n, j - n - 1] + - integral_img[i - n - 1, j + n]) + + outer = (integral_img[i + n2, j + n2] + + integral_img[i - n2 - 1, j - n2 - 1] + - integral_img[i + n2, j - n2 - 1] + - integral_img[i - n2 - 1, j + n2]) + + filtered_image[i, j] = (outer_weight * outer + - total_weight * inner) diff --git a/skimage/feature/tests/test_censure.py b/skimage/feature/tests/test_censure.py index 95bfcda3..ecd8a90f 100644 --- a/skimage/feature/tests/test_censure.py +++ b/skimage/feature/tests/test_censure.py @@ -2,6 +2,7 @@ import numpy as np from numpy.testing import assert_array_equal, assert_raises from skimage.data import moon from skimage.feature import CENSURE +from skimage._shared.testing import test_parallel img = moon() @@ -53,6 +54,7 @@ def test_keypoints_censure_moon_image_dob(): assert_array_equal(expected_scales, detector.scales) +@test_parallel() def test_keypoints_censure_moon_image_octagon(): """Verify the actual Censure keypoints and their corresponding scale with the expected values for Octagon filter.""" From 9cd7ca5881d92c51c5e67ed15f7294c0d29b56ab Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 20 May 2015 17:39:46 -0700 Subject: [PATCH 15/20] Do not acquire GIL for BRIEF --- skimage/feature/brief.py | 11 ++++++----- skimage/feature/brief_cy.pyx | 21 +++++++++++---------- skimage/feature/tests/test_brief.py | 1 + 3 files changed, 18 insertions(+), 15 deletions(-) diff --git a/skimage/feature/brief.py b/skimage/feature/brief.py index 783cdbd0..3a53f89f 100644 --- a/skimage/feature/brief.py +++ b/skimage/feature/brief.py @@ -140,7 +140,8 @@ class BRIEF(DescriptorExtractor): """ assert_nD(image, 2) - np.random.seed(self.sample_seed) + random = np.random.RandomState() + random.seed(self.sample_seed) image = _prepare_grayscale_input_2D(image) @@ -151,7 +152,7 @@ class BRIEF(DescriptorExtractor): desc_size = self.descriptor_size patch_size = self.patch_size if self.mode == 'normal': - samples = (patch_size / 5.0) * np.random.randn(desc_size * 8) + samples = (patch_size / 5.0) * random.randn(desc_size * 8) samples = np.array(samples, dtype=np.int32) samples = samples[(samples < (patch_size // 2)) & (samples > - (patch_size - 2) // 2)] @@ -159,9 +160,9 @@ class BRIEF(DescriptorExtractor): pos1 = samples[:desc_size * 2].reshape(desc_size, 2) pos2 = samples[desc_size * 2:desc_size * 4].reshape(desc_size, 2) elif self.mode == 'uniform': - samples = np.random.randint(-(patch_size - 2) // 2, - (patch_size // 2) + 1, - (desc_size * 2, 2)) + samples = random.randint(-(patch_size - 2) // 2, + (patch_size // 2) + 1, + (desc_size * 2, 2)) samples = np.array(samples, dtype=np.int32) pos1, pos2 = np.split(samples, 2) diff --git a/skimage/feature/brief_cy.pyx b/skimage/feature/brief_cy.pyx index 8cd1afa7..1773894d 100644 --- a/skimage/feature/brief_cy.pyx +++ b/skimage/feature/brief_cy.pyx @@ -12,13 +12,14 @@ def _brief_loop(double[:, ::1] image, unsigned char[:, ::1] descriptors, cdef Py_ssize_t k, d, kr, kc, pr0, pr1, pc0, pc1 - for p in range(pos0.shape[0]): - pr0 = pos0[p, 0] - pc0 = pos0[p, 1] - pr1 = pos1[p, 0] - pc1 = pos1[p, 1] - for k in range(keypoints.shape[0]): - kr = keypoints[k, 0] - kc = keypoints[k, 1] - if image[kr + pr0, kc + pc0] < image[kr + pr1, kc + pc1]: - descriptors[k, p] = True + with nogil: + for p in range(pos0.shape[0]): + pr0 = pos0[p, 0] + pc0 = pos0[p, 1] + pr1 = pos1[p, 0] + pc1 = pos1[p, 1] + for k in range(keypoints.shape[0]): + kr = keypoints[k, 0] + kc = keypoints[k, 1] + if image[kr + pr0, kc + pc0] < image[kr + pr1, kc + pc1]: + descriptors[k, p] = True diff --git a/skimage/feature/tests/test_brief.py b/skimage/feature/tests/test_brief.py index 554301b2..e70a1148 100644 --- a/skimage/feature/tests/test_brief.py +++ b/skimage/feature/tests/test_brief.py @@ -4,6 +4,7 @@ from skimage import data from skimage import transform as tf from skimage.color import rgb2gray from skimage.feature import BRIEF, corner_peaks, corner_harris +from skimage._shared.testing import test_parallel def test_color_image_unsupported_error(): From 4a9534548956c80692f745b807041f138be08e97 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 20 May 2015 17:51:41 -0700 Subject: [PATCH 16/20] Do not acquire GIL for texture functions --- skimage/feature/_texture.pyx | 279 +++++++++++++------------- skimage/feature/tests/test_texture.py | 3 + 2 files changed, 146 insertions(+), 136 deletions(-) diff --git a/skimage/feature/_texture.pyx b/skimage/feature/_texture.pyx index 5eb91486..d454a812 100644 --- a/skimage/feature/_texture.pyx +++ b/skimage/feature/_texture.pyx @@ -4,7 +4,7 @@ #cython: wraparound=False import numpy as np cimport numpy as cnp -from libc.math cimport sin, cos, abs +from libc.math cimport sin, cos, abs, NAN from .._shared.interpolation cimport bilinear_interpolation, round @@ -36,32 +36,33 @@ def _glcm_loop(cnp.uint8_t[:, ::1] image, double[:] distances, cnp.uint8_t i, j cnp.float64_t angle, distance - rows = image.shape[0] - cols = image.shape[1] + with nogil: + rows = image.shape[0] + cols = image.shape[1] - for a_idx in range(len(angles)): - angle = angles[a_idx] - for d_idx in range(len(distances)): - distance = distances[d_idx] - for r in range(rows): - for c in range(cols): - i = image[r, c] + for a_idx in range(angles.shape[0]): + angle = angles[a_idx] + for d_idx in range(distances.shape[0]): + distance = distances[d_idx] + for r in range(rows): + for c in range(cols): + i = image[r, c] - # compute the location of the offset pixel - row = r + round(sin(angle) * distance) - col = c + round(cos(angle) * distance) + # compute the location of the offset pixel + row = r + round(sin(angle) * distance) + col = c + round(cos(angle) * distance) - # make sure the offset is within bounds - if row >= 0 and row < rows and \ - col >= 0 and col < cols: - j = image[row, col] + # make sure the offset is within bounds + if row >= 0 and row < rows and \ + col >= 0 and col < cols: + j = image[row, col] - if i >= 0 and i < levels and \ - j >= 0 and j < levels: - out[i, j, d_idx, a_idx] += 1 + if i >= 0 and i < levels and \ + j >= 0 and j < levels: + out[i, j, d_idx, a_idx] += 1 -cdef inline int _bit_rotate_right(int value, int length): +cdef inline int _bit_rotate_right(int value, int length) nogil: """Cyclic bit shift to the right. Parameters @@ -132,124 +133,130 @@ def _local_binary_pattern(double[:, ::1] image, # To compute the variance features cdef double sum_, var_, texture_i - for r in range(image.shape[0]): - for c in range(image.shape[1]): - for i in range(P): - texture[i] = bilinear_interpolation(&image[0, 0], rows, cols, - r + rp[i], c + cp[i], - 'C', 0) - # signed / thresholded texture - for i in range(P): - if texture[i] - image[r, c] >= 0: - signed_texture[i] = 1 - else: - signed_texture[i] = 0 - - lbp = 0 - - # if method == 'var': - if method == 'V': - # Compute the variance without passing from numpy. - # Following the LBP paper, we're taking a biased estimate - # of the variance (ddof=0) - sum_ = 0.0 - var_ = 0.0 + with nogil: + for r in range(image.shape[0]): + for c in range(image.shape[1]): for i in range(P): - texture_i = texture[i] - sum_ += texture_i - var_ += texture_i * texture_i - var_ = (var_ - (sum_ * sum_) / P) / P - if var_ != 0: - lbp = var_ - else: - lbp = np.nan - # if method == 'uniform': - elif method == 'U' or method == 'N': - # determine number of 0 - 1 changes - changes = 0 - for i in range(P - 1): - changes += abs(signed_texture[i] - signed_texture[i + 1]) - if method == 'N': - # Uniform local binary patterns are defined as patterns - # with at most 2 value changes (from 0 to 1 or from 1 to - # 0). Uniform patterns can be caraterized by their number - # `n_ones` of 1. The possible values for `n_ones` range - # from 0 to P. - # Here is an example for P = 4: - # n_ones=0: 0000 - # n_ones=1: 0001, 1000, 0100, 0010 - # n_ones=2: 0011, 1001, 1100, 0110 - # n_ones=3: 0111, 1011, 1101, 1110 - # n_ones=4: 1111 - # - # For a pattern of size P there are 2 constant patterns - # corresponding to n_ones=0 and n_ones=P. For each other - # value of `n_ones` , i.e n_ones=[1..P-1], there are P - # possible patterns which are related to each other through - # circular permutations. The total number of uniform - # patterns is thus (2 + P * (P - 1)). - # Given any pattern (uniform or not) we must be able to - # associate a unique code: - # 1. Constant patterns patterns (with n_ones=0 and - # n_ones=P) and non uniform patterns are given fixed - # code values. - # 2. Other uniform patterns are indexed considering the - # value of n_ones, and an index called 'rot_index' - # reprenting the number of circular right shifts - # required to obtain the pattern starting from a - # reference position (corresponding to all zeros stacked - # on the right). This number of rotations (or circular - # right shifts) 'rot_index' is efficiently computed by - # considering the positions of the first 1 and the first - # 0 found in the pattern. - - if changes <= 2: - # We have a uniform pattern - n_ones = 0 # determines the number of ones - first_one = -1 # position was the first one - first_zero = -1 # position of the first zero - for i in range(P): - if signed_texture[i]: - n_ones += 1 - if first_one == -1: - first_one = i - else: - if first_zero == -1: - first_zero = i - if n_ones == 0: - lbp = 0 - elif n_ones == P: - lbp = P * (P - 1) + 1 - else: - if first_one == 0: - rot_index = n_ones - first_zero - else: - rot_index = P - first_one - lbp = 1 + (n_ones - 1) * P + rot_index - else: # changes > 2 - lbp = P * (P - 1) + 2 - else: # method != 'N' - if changes <= 2: - for i in range(P): - lbp += signed_texture[i] + texture[i] = bilinear_interpolation(&image[0, 0], rows, cols, + r + rp[i], c + cp[i], + 'C', 0) + # signed / thresholded texture + for i in range(P): + if texture[i] - image[r, c] >= 0: + signed_texture[i] = 1 else: - lbp = P + 1 - else: - # method == 'default' - for i in range(P): - lbp += signed_texture[i] * weights[i] + signed_texture[i] = 0 - # method == 'ror' - if method == 'R': - # shift LBP P times to the right and get minimum value - rotation_chain[0] = lbp - for i in range(1, P): - rotation_chain[i] = \ - _bit_rotate_right(rotation_chain[i - 1], P) - lbp = rotation_chain[0] - for i in range(1, P): - lbp = min(lbp, rotation_chain[i]) + lbp = 0 - output[r, c] = lbp + # if method == 'var': + if method == 'V': + # Compute the variance without passing from numpy. + # Following the LBP paper, we're taking a biased estimate + # of the variance (ddof=0) + sum_ = 0.0 + var_ = 0.0 + for i in range(P): + texture_i = texture[i] + sum_ += texture_i + var_ += texture_i * texture_i + var_ = (var_ - (sum_ * sum_) / P) / P + if var_ != 0: + lbp = var_ + else: + lbp = NAN + # if method == 'uniform': + elif method == 'U' or method == 'N': + # determine number of 0 - 1 changes + changes = 0 + for i in range(P - 1): + changes += (signed_texture[i] + - signed_texture[i + 1]) != 0 + if method == 'N': + # Uniform local binary patterns are defined as patterns + # with at most 2 value changes (from 0 to 1 or from 1 to + # 0). Uniform patterns can be characterized by their + # number `n_ones` of 1. The possible values for + # `n_ones` range from 0 to P. + # + # Here is an example for P = 4: + # n_ones=0: 0000 + # n_ones=1: 0001, 1000, 0100, 0010 + # n_ones=2: 0011, 1001, 1100, 0110 + # n_ones=3: 0111, 1011, 1101, 1110 + # n_ones=4: 1111 + # + # For a pattern of size P there are 2 constant patterns + # corresponding to n_ones=0 and n_ones=P. For each other + # value of `n_ones` , i.e n_ones=[1..P-1], there are P + # possible patterns which are related to each other + # through circular permutations. The total number of + # uniform patterns is thus (2 + P * (P - 1)). + + # Given any pattern (uniform or not) we must be able to + # associate a unique code: + # + # 1. Constant patterns patterns (with n_ones=0 and + # n_ones=P) and non uniform patterns are given fixed + # code values. + # + # 2. Other uniform patterns are indexed considering the + # value of n_ones, and an index called 'rot_index' + # reprenting the number of circular right shifts + # required to obtain the pattern starting from a + # reference position (corresponding to all zeros stacked + # on the right). This number of rotations (or circular + # right shifts) 'rot_index' is efficiently computed by + # considering the positions of the first 1 and the first + # 0 found in the pattern. + + if changes <= 2: + # We have a uniform pattern + n_ones = 0 # determines the number of ones + first_one = -1 # position was the first one + first_zero = -1 # position of the first zero + for i in range(P): + if signed_texture[i]: + n_ones += 1 + if first_one == -1: + first_one = i + else: + if first_zero == -1: + first_zero = i + if n_ones == 0: + lbp = 0 + elif n_ones == P: + lbp = P * (P - 1) + 1 + else: + if first_one == 0: + rot_index = n_ones - first_zero + else: + rot_index = P - first_one + lbp = 1 + (n_ones - 1) * P + rot_index + else: # changes > 2 + lbp = P * (P - 1) + 2 + else: # method != 'N' + if changes <= 2: + for i in range(P): + lbp += signed_texture[i] + else: + lbp = P + 1 + else: + # method == 'default' + for i in range(P): + lbp += signed_texture[i] * weights[i] + + # method == 'ror' + if method == 'R': + # shift LBP P times to the right and get minimum value + rotation_chain[0] = lbp + for i in range(1, P): + rotation_chain[i] = \ + _bit_rotate_right(rotation_chain[i - 1], P) + lbp = rotation_chain[0] + for i in range(1, P): + lbp = min(lbp, rotation_chain[i]) + + output[r, c] = lbp return np.asarray(output) diff --git a/skimage/feature/tests/test_texture.py b/skimage/feature/tests/test_texture.py index 788ccbc6..10444aa8 100644 --- a/skimage/feature/tests/test_texture.py +++ b/skimage/feature/tests/test_texture.py @@ -1,5 +1,6 @@ import numpy as np from skimage.feature import greycomatrix, greycoprops, local_binary_pattern +from skimage._shared.testing import test_parallel class TestGLCM(): @@ -10,6 +11,7 @@ class TestGLCM(): [0, 2, 2, 2], [2, 2, 3, 3]], dtype=np.uint8) + @test_parallel() def test_output_angles(self): result = greycomatrix(self.image, [1], [0, np.pi / 4, np.pi / 2, 3 * np.pi / 4], 4) assert result.shape == (4, 4, 1, 4) @@ -162,6 +164,7 @@ class TestLBP(): [ 0, 255, 30, 34, 255, 24], [146, 241, 255, 0, 189, 126]], dtype='double') + @test_parallel() def test_default(self): lbp = local_binary_pattern(self.image, 8, 1, 'default') ref = np.array([[ 0, 251, 0, 255, 96, 255], From a2ddba21a56bf9f09bfc2800eb9a20aac2f31968 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 20 May 2015 20:22:02 -0700 Subject: [PATCH 17/20] Do not acquire GIL for hessian_matrix_det --- skimage/feature/_hessian_det_appx.pyx | 44 ++++++++++++++------------- skimage/feature/tests/test_corner.py | 1 + 2 files changed, 24 insertions(+), 21 deletions(-) diff --git a/skimage/feature/_hessian_det_appx.pyx b/skimage/feature/_hessian_det_appx.pyx index 0c54cba7..3c8cec4d 100644 --- a/skimage/feature/_hessian_det_appx.pyx +++ b/skimage/feature/_hessian_det_appx.pyx @@ -6,7 +6,8 @@ import numpy as np cimport numpy as cnp -cdef inline Py_ssize_t _clip(Py_ssize_t x, Py_ssize_t low, Py_ssize_t high): +cdef inline Py_ssize_t _clip(Py_ssize_t x, Py_ssize_t low, + Py_ssize_t high) nogil: """Clips coordinate between high and low. This method was created so that `hessian_det_appx` does not have to make @@ -36,7 +37,7 @@ cdef inline Py_ssize_t _clip(Py_ssize_t x, Py_ssize_t low, Py_ssize_t high): cdef inline cnp.double_t _integ( cnp.double_t[:, ::1] img, Py_ssize_t r, Py_ssize_t c, - Py_ssize_t rl, Py_ssize_t cl): + Py_ssize_t rl, Py_ssize_t cl) nogil: """Integrate over the integral image in the given window This method was created so that `hessian_det_appx` does not have to make @@ -123,31 +124,32 @@ def _hessian_matrix_det(cnp.double_t[:, ::1] img, double sigma): cdef float dxx, dyy, dxy - if not size % 2: - size += 1 + with nogil: + if size % 2 == 0: + size += 1 - for r in range(height): - for c in range(width): - tl = _integ(img, r - s3, c - s3, s3, s3) # top left - br = _integ(img, r + 1, c + 1, s3, s3) # bottom right - bl = _integ(img, r - s3, c + 1, s3, s3) # bottom left - tr = _integ(img, r + 1, c - s3, s3, s3) # top right + for r in range(height): + for c in range(width): + tl = _integ(img, r - s3, c - s3, s3, s3) # top left + br = _integ(img, r + 1, c + 1, s3, s3) # bottom right + bl = _integ(img, r - s3, c + 1, s3, s3) # bottom left + tr = _integ(img, r + 1, c - s3, s3, s3) # top right - dxy = bl + tr - tl - br - dxy = -dxy * w_i + dxy = bl + tr - tl - br + dxy = -dxy * w_i - mid = _integ(img, r - s3 + 1, c - s2, 2 * s3 - 1, w) # middle box - side = _integ(img, r - s3 + 1, c - s3 / 2, 2 * s3 - 1, s3) # sides + mid = _integ(img, r - s3 + 1, c - s2, 2 * s3 - 1, w) # middle box + side = _integ(img, r - s3 + 1, c - s3 / 2, 2 * s3 - 1, s3) # sides - dxx = mid - 3 * side - dxx = -dxx * w_i + dxx = mid - 3 * side + dxx = -dxx * w_i - mid = _integ(img, r - s2, c - s3 + 1, w, 2 * s3 - 1) - side = _integ(img, r - s3 / 2, c - s3 + 1, s3, 2 * s3 - 1) + mid = _integ(img, r - s2, c - s3 + 1, w, 2 * s3 - 1) + side = _integ(img, r - s3 / 2, c - s3 + 1, s3, 2 * s3 - 1) - dyy = mid - 3 * side - dyy = -dyy * w_i + dyy = mid - 3 * side + dyy = -dyy * w_i - out[r, c] = (dxx * dyy - 0.81 * (dxy * dxy)) + out[r, c] = (dxx * dyy - 0.81 * (dxy * dxy)) return out diff --git a/skimage/feature/tests/test_corner.py b/skimage/feature/tests/test_corner.py index aa7e1409..d6c3418f 100644 --- a/skimage/feature/tests/test_corner.py +++ b/skimage/feature/tests/test_corner.py @@ -93,6 +93,7 @@ def test_hessian_matrix_eigvals(): [0, 0, 0, 0, 0]])) +@test_parallel() def test_hessian_matrix_det(): image = np.zeros((5, 5)) image[2, 2] = 1 From fcdffd7e3d41790ab97c05e93d42ead74c6e6f55 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 20 May 2015 20:26:57 -0700 Subject: [PATCH 18/20] Do not acquire GIL for draw functions --- skimage/draw/_draw.pyx | 66 +++++++++++++++++---------------- skimage/draw/tests/test_draw.py | 2 + 2 files changed, 36 insertions(+), 32 deletions(-) diff --git a/skimage/draw/_draw.pyx b/skimage/draw/_draw.pyx index d99c10ae..0ac540df 100644 --- a/skimage/draw/_draw.pyx +++ b/skimage/draw/_draw.pyx @@ -6,7 +6,7 @@ import math import numpy as np cimport numpy as cnp -from libc.math cimport sqrt, sin, cos, floor, ceil +from libc.math cimport sqrt, sin, cos, floor, ceil, fabs from .._shared.geometry cimport point_in_polygon @@ -83,39 +83,41 @@ def line(Py_ssize_t y, Py_ssize_t x, Py_ssize_t y2, Py_ssize_t x2): cdef Py_ssize_t dy = abs(y2 - y) cdef Py_ssize_t sx, sy, d, i - if (x2 - x) > 0: - sx = 1 - else: - sx = -1 - if (y2 - y) > 0: - sy = 1 - else: - sy = -1 - if dy > dx: - steep = 1 - x, y = y, x - dx, dy = dy, dx - sx, sy = sy, sx - d = (2 * dy) - dx + with nogil: + if (x2 - x) > 0: + sx = 1 + else: + sx = -1 + if (y2 - y) > 0: + sy = 1 + else: + sy = -1 + if dy > dx: + steep = 1 + x, y = y, x + dx, dy = dy, dx + sx, sy = sy, sx + d = (2 * dy) - dx cdef Py_ssize_t[::1] rr = np.zeros(int(dx) + 1, dtype=np.intp) cdef Py_ssize_t[::1] cc = np.zeros(int(dx) + 1, dtype=np.intp) - for i in range(dx): - if steep: - rr[i] = x - cc[i] = y - else: - rr[i] = y - cc[i] = x - while d >= 0: - y = y + sy - d = d - (2 * dx) - x = x + sx - d = d + (2 * dy) + with nogil: + for i in range(dx): + if steep: + rr[i] = x + cc[i] = y + else: + rr[i] = y + cc[i] = x + while d >= 0: + y = y + sy + d = d - (2 * dx) + x = x + sx + d = d + (2 * dy) - rr[dx] = y2 - cc[dx] = x2 + rr[dx] = y2 + cc[dx] = x2 return np.asarray(rr), np.asarray(cc) @@ -187,7 +189,7 @@ def line_aa(Py_ssize_t y1, Py_ssize_t x1, Py_ssize_t y2, Py_ssize_t x2): while True: cc.append(x) rr.append(y) - val.append(1. * abs(err - dx + dy) / (ed)) + val.append(1. * fabs(err - dx + dy) / (ed)) e = err if 2 * e >= -dx: if x == x2: @@ -195,7 +197,7 @@ def line_aa(Py_ssize_t y1, Py_ssize_t x1, Py_ssize_t y2, Py_ssize_t x2): if e + dy < ed: cc.append(x) rr.append(y + sign_y) - val.append(1. * abs(e + dy) / (ed)) + val.append(1. * fabs(e + dy) / (ed)) err -= dy x += sign_x if 2 * e <= dy: @@ -204,7 +206,7 @@ def line_aa(Py_ssize_t y1, Py_ssize_t x1, Py_ssize_t y2, Py_ssize_t x2): if dx - e < ed: cc.append(x) rr.append(y) - val.append(abs(dx - e) / (ed)) + val.append(fabs(dx - e) / (ed)) err += dx y += sign_y diff --git a/skimage/draw/tests/test_draw.py b/skimage/draw/tests/test_draw.py index fd759b4e..c3e49664 100644 --- a/skimage/draw/tests/test_draw.py +++ b/skimage/draw/tests/test_draw.py @@ -1,5 +1,6 @@ from numpy.testing import assert_array_equal, assert_equal import numpy as np +from skimage._shared.testing import test_parallel from skimage.draw import (set_color, line, line_aa, polygon, circle, circle_perimeter, circle_perimeter_aa, @@ -19,6 +20,7 @@ def test_set_color(): assert_array_equal(img, img_) +@test_parallel() def test_line_horizontal(): img = np.zeros((10, 10)) From 6e500b173377db76c77c27af7496e3619992b4a9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 20 May 2015 20:29:19 -0700 Subject: [PATCH 19/20] Do not acquire GIL for shared functions --- skimage/_shared/geometry.pxd | 4 ++-- skimage/_shared/geometry.pyx | 5 +++-- skimage/_shared/transform.pxd | 2 +- skimage/_shared/transform.pyx | 2 +- 4 files changed, 7 insertions(+), 6 deletions(-) diff --git a/skimage/_shared/geometry.pxd b/skimage/_shared/geometry.pxd index 3379318c..4524fa37 100644 --- a/skimage/_shared/geometry.pxd +++ b/skimage/_shared/geometry.pxd @@ -1,6 +1,6 @@ cdef unsigned char point_in_polygon(Py_ssize_t nr_verts, double *xp, double *yp, - double x, double y) + double x, double y) nogil cdef void points_in_polygon(Py_ssize_t nr_verts, double *xp, double *yp, Py_ssize_t nr_points, double *x, double *y, - unsigned char *result) + unsigned char *result) nogil diff --git a/skimage/_shared/geometry.pyx b/skimage/_shared/geometry.pyx index beb07e14..ea6695a9 100644 --- a/skimage/_shared/geometry.pyx +++ b/skimage/_shared/geometry.pyx @@ -5,7 +5,8 @@ cdef inline unsigned char point_in_polygon(Py_ssize_t nr_verts, double *xp, - double *yp, double x, double y): + double *yp, double x, + double y) nogil: """Test whether point lies inside a polygon. Parameters @@ -33,7 +34,7 @@ cdef inline unsigned char point_in_polygon(Py_ssize_t nr_verts, double *xp, cdef void points_in_polygon(Py_ssize_t nr_verts, double *xp, double *yp, Py_ssize_t nr_points, double *x, double *y, - unsigned char *result): + unsigned char *result) nogil: """Test whether points lie inside a polygon. Parameters diff --git a/skimage/_shared/transform.pxd b/skimage/_shared/transform.pxd index 4978e843..b48fe92e 100644 --- a/skimage/_shared/transform.pxd +++ b/skimage/_shared/transform.pxd @@ -2,4 +2,4 @@ cimport numpy as cnp cdef float integrate(float[:, ::1] sat, Py_ssize_t r0, Py_ssize_t c0, - Py_ssize_t r1, Py_ssize_t c1) + Py_ssize_t r1, Py_ssize_t c1) nogil diff --git a/skimage/_shared/transform.pyx b/skimage/_shared/transform.pyx index d77ec583..cbb40144 100644 --- a/skimage/_shared/transform.pyx +++ b/skimage/_shared/transform.pyx @@ -6,7 +6,7 @@ cimport numpy as cnp cdef float integrate(float[:, ::1] sat, Py_ssize_t r0, Py_ssize_t c0, - Py_ssize_t r1, Py_ssize_t c1): + Py_ssize_t r1, Py_ssize_t c1) nogil: """ Using a summed area table / integral image, calculate the sum over a given window. From 990892fcdec52229149d81781a7366243cf16253 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 21 May 2015 08:27:00 -0700 Subject: [PATCH 20/20] Replace libc.math.NAN with numpy.math.NAN for MSVC --- skimage/feature/_texture.pyx | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/skimage/feature/_texture.pyx b/skimage/feature/_texture.pyx index d454a812..d716e690 100644 --- a/skimage/feature/_texture.pyx +++ b/skimage/feature/_texture.pyx @@ -4,10 +4,14 @@ #cython: wraparound=False import numpy as np cimport numpy as cnp -from libc.math cimport sin, cos, abs, NAN +from libc.math cimport sin, cos, abs from .._shared.interpolation cimport bilinear_interpolation, round +cdef extern from "numpy/npy_math.h": + double NAN "NPY_NAN" + + def _glcm_loop(cnp.uint8_t[:, ::1] image, double[:] distances, double[:] angles, Py_ssize_t levels, cnp.uint32_t[:, :, :, ::1] out):