diff --git a/skimage/_shared/geometry.pxd b/skimage/_shared/geometry.pxd index afdc6b5b..3379318c 100644 --- a/skimage/_shared/geometry.pxd +++ b/skimage/_shared/geometry.pxd @@ -1,6 +1,6 @@ -cdef unsigned char point_in_polygon(int nr_verts, double *xp, double *yp, +cdef unsigned char point_in_polygon(Py_ssize_t nr_verts, double *xp, double *yp, double x, double y) -cdef void points_in_polygon(int nr_verts, double *xp, double *yp, - int nr_points, double *x, double *y, +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) diff --git a/skimage/_shared/geometry.pyx b/skimage/_shared/geometry.pyx index 3f4850b0..beb07e14 100644 --- a/skimage/_shared/geometry.pyx +++ b/skimage/_shared/geometry.pyx @@ -4,8 +4,8 @@ #cython: wraparound=False -cdef inline unsigned char point_in_polygon(int nr_verts, double *xp, double *yp, - double x, double y): +cdef inline unsigned char point_in_polygon(Py_ssize_t nr_verts, double *xp, + double *yp, double x, double y): """Test whether point lies inside a polygon. Parameters @@ -17,9 +17,9 @@ cdef inline unsigned char point_in_polygon(int nr_verts, double *xp, double *yp, x, y : double Coordinates of point. """ - cdef int i + cdef Py_ssize_t i cdef unsigned char c = 0 - cdef int j = nr_verts - 1 + cdef Py_ssize_t j = nr_verts - 1 for i in range(nr_verts): if ( (((yp[i] <= y) and (y < yp[j])) or @@ -31,8 +31,8 @@ cdef inline unsigned char point_in_polygon(int nr_verts, double *xp, double *yp, return c -cdef void points_in_polygon(int nr_verts, double *xp, double *yp, - int nr_points, double *x, double *y, +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): """Test whether points lie inside a polygon. @@ -49,6 +49,6 @@ cdef void points_in_polygon(int nr_verts, double *xp, double *yp, result : unsigned char array Test results for each point. """ - cdef int n + cdef Py_ssize_t n for n in range(nr_points): result[n] = point_in_polygon(nr_verts, xp, yp, x[n], y[n]) diff --git a/skimage/_shared/interpolation.pxd b/skimage/_shared/interpolation.pxd index 4e00dfb0..c5f32b6a 100644 --- a/skimage/_shared/interpolation.pxd +++ b/skimage/_shared/interpolation.pxd @@ -1,27 +1,27 @@ -cdef double nearest_neighbour_interpolation(double* image, int rows, - int cols, double r, +cdef double nearest_neighbour_interpolation(double* image, Py_ssize_t rows, + Py_ssize_t cols, double r, double c, char mode, double cval) -cdef double bilinear_interpolation(double* image, int rows, int cols, +cdef double bilinear_interpolation(double* image, Py_ssize_t rows, Py_ssize_t cols, double r, double c, char mode, double cval) cdef double quadratic_interpolation(double x, double[3] f) -cdef double biquadratic_interpolation(double* image, int rows, int cols, +cdef double biquadratic_interpolation(double* image, Py_ssize_t rows, Py_ssize_t cols, double r, double c, char mode, double cval) cdef double cubic_interpolation(double x, double[4] f) -cdef double bicubic_interpolation(double* image, int rows, int cols, +cdef double bicubic_interpolation(double* image, Py_ssize_t rows, Py_ssize_t cols, double r, double c, char mode, double cval) -cdef double get_pixel2d(double* image, int rows, int cols, int r, int c, - char mode, double cval) +cdef double get_pixel2d(double* image, Py_ssize_t rows, Py_ssize_t cols, Py_ssize_t r, + Py_ssize_t c, char mode, double cval) -cdef double get_pixel3d(double* image, int rows, int cols, int dims, int r, - int c, int d, char mode, double cval) +cdef double get_pixel3d(double* image, Py_ssize_t rows, Py_ssize_t cols, Py_ssize_t dims, + Py_ssize_t r, Py_ssize_t c, Py_ssize_t d, char mode, double cval) -cdef int coord_map(int dim, int coord, char mode) +cdef Py_ssize_t coord_map(Py_ssize_t dim, Py_ssize_t coord, char mode) diff --git a/skimage/_shared/interpolation.pyx b/skimage/_shared/interpolation.pyx index b34ba507..a8b96014 100644 --- a/skimage/_shared/interpolation.pyx +++ b/skimage/_shared/interpolation.pyx @@ -5,12 +5,12 @@ from libc.math cimport ceil, floor -cdef inline int round(double r): - return ((r + 0.5) if (r > 0.0) else (r - 0.5)) +cdef inline Py_ssize_t round(double r): + return ((r + 0.5) if (r > 0.0) else (r - 0.5)) -cdef inline double nearest_neighbour_interpolation(double* image, int rows, - int cols, double r, +cdef inline double nearest_neighbour_interpolation(double* image, Py_ssize_t rows, + Py_ssize_t cols, double r, double c, char mode, double cval): """Nearest neighbour interpolation at a given position in the image. @@ -35,13 +35,12 @@ cdef inline double nearest_neighbour_interpolation(double* image, int rows, """ - return get_pixel2d(image, rows, cols, round(r), round(c), - mode, cval) + return get_pixel2d(image, rows, cols, round(r), round(c), mode, cval) -cdef inline double bilinear_interpolation(double* image, int rows, int cols, - double r, double c, char mode, - double cval): +cdef inline double bilinear_interpolation(double* image, Py_ssize_t rows, + Py_ssize_t cols, double r, double c, + char mode, double cval): """Bilinear interpolation at a given position in the image. Parameters @@ -64,12 +63,12 @@ cdef inline double bilinear_interpolation(double* image, int rows, int cols, """ cdef double dr, dc - cdef int minr, minc, maxr, maxc + cdef Py_ssize_t minr, minc, maxr, maxc - minr = floor(r) - minc = floor(c) - maxr = ceil(r) - maxc = ceil(c) + minr = floor(r) + minc = floor(c) + maxr = ceil(r) + maxc = ceil(c) dr = r - minr dc = c - minc top = (1 - dc) * get_pixel2d(image, rows, cols, minr, minc, mode, cval) \ @@ -98,9 +97,9 @@ cdef inline double quadratic_interpolation(double x, double[3] f): return f[1] - 0.25 * (f[0] - f[2]) * x -cdef inline double biquadratic_interpolation(double* image, int rows, int cols, - double r, double c, char mode, - double cval): +cdef inline double biquadratic_interpolation(double* image, Py_ssize_t rows, + Py_ssize_t cols, double r, double c, + char mode, double cval): """Biquadratic interpolation at a given position in the image. Parameters @@ -123,8 +122,8 @@ cdef inline double biquadratic_interpolation(double* image, int rows, int cols, """ - cdef int r0 = round(r) - cdef int c0 = round(c) + cdef Py_ssize_t r0 = round(r) + cdef Py_ssize_t c0 = round(c) if r < 0: r0 -= 1 if c < 0: @@ -139,7 +138,7 @@ cdef inline double biquadratic_interpolation(double* image, int rows, int cols, cdef double fc[3], fr[3] - cdef int pr, pc + cdef Py_ssize_t pr, pc # row-wise cubic interpolation for pr in range(r0, r0 + 3): @@ -174,9 +173,9 @@ cdef inline double cubic_interpolation(double x, double[4] f): (3.0 * (f[1] - f[2]) + f[3] - f[0]))) -cdef inline double bicubic_interpolation(double* image, int rows, int cols, - double r, double c, char mode, - double cval): +cdef inline double bicubic_interpolation(double* image, Py_ssize_t rows, + Py_ssize_t cols, double r, double c, + char mode, double cval): """Bicubic interpolation at a given position in the image. Parameters @@ -199,8 +198,8 @@ cdef inline double bicubic_interpolation(double* image, int rows, int cols, """ - cdef int r0 = r - 1 - cdef int c0 = c - 1 + cdef Py_ssize_t r0 = r - 1 + cdef Py_ssize_t c0 = c - 1 if r < 0: r0 -= 1 if c < 0: @@ -211,7 +210,7 @@ cdef inline double bicubic_interpolation(double* image, int rows, int cols, cdef double fc[4], fr[4] - cdef int pr, pc + cdef Py_ssize_t pr, pc # row-wise cubic interpolation for pr in range(r0, r0 + 4): @@ -223,8 +222,8 @@ cdef inline double bicubic_interpolation(double* image, int rows, int cols, return cubic_interpolation(xr, fr) -cdef inline double get_pixel2d(double* image, int rows, int cols, int r, int c, - char mode, double cval): +cdef inline double get_pixel2d(double* image, Py_ssize_t rows, Py_ssize_t cols, + Py_ssize_t r, Py_ssize_t c, char mode, double cval): """Get a pixel from the image, taking wrapping mode into consideration. Parameters @@ -255,8 +254,9 @@ cdef inline double get_pixel2d(double* image, int rows, int cols, int r, int c, return image[coord_map(rows, r, mode) * cols + coord_map(cols, c, mode)] -cdef inline double get_pixel3d(double* image, int rows, int cols, int dims, int r, - int c, int d, char mode, double cval): +cdef inline double get_pixel3d(double* image, Py_ssize_t rows, Py_ssize_t cols, + Py_ssize_t dims, Py_ssize_t r, Py_ssize_t c, Py_ssize_t d, + char mode, double cval): """Get a pixel from the image, taking wrapping mode into consideration. Parameters @@ -289,7 +289,7 @@ cdef inline double get_pixel3d(double* image, int rows, int cols, int dims, int + d] -cdef inline int coord_map(int dim, int coord, char mode): +cdef inline Py_ssize_t coord_map(Py_ssize_t dim, Py_ssize_t coord, char mode): """ Wrap a coordinate, according to a given mode. @@ -308,20 +308,20 @@ cdef inline int coord_map(int dim, int coord, char mode): if mode == 'R': # reflect if coord < 0: # How many times times does the coordinate wrap? - if (-coord / dim) % 2 != 0: - return dim - (-coord % dim) + if (-coord / dim) % 2 != 0: + return dim - (-coord % dim) else: - return (-coord % dim) + return (-coord % dim) elif coord > dim: - if (coord / dim) % 2 != 0: - return (dim - (coord % dim)) + if (coord / dim) % 2 != 0: + return (dim - (coord % dim)) else: - return (coord % dim) + return (coord % dim) elif mode == 'W': # wrap if coord < 0: - return (dim - (-coord % dim)) + return (dim - (-coord % dim)) elif coord > dim: - return (coord % dim) + return (coord % dim) elif mode == 'N': # nearest if coord < 0: return 0 diff --git a/skimage/_shared/transform.pxd b/skimage/_shared/transform.pxd index 0edc22a4..ccb16ff3 100644 --- a/skimage/_shared/transform.pxd +++ b/skimage/_shared/transform.pxd @@ -2,4 +2,4 @@ cimport numpy as cnp cdef float integrate(cnp.ndarray[float, ndim=2, mode="c"] sat, - int r0, int c0, int r1, int c1) + Py_ssize_t r0, Py_ssize_t c0, Py_ssize_t r1, Py_ssize_t c1) diff --git a/skimage/_shared/transform.pyx b/skimage/_shared/transform.pyx index ba0efc71..8ce2ab67 100644 --- a/skimage/_shared/transform.pyx +++ b/skimage/_shared/transform.pyx @@ -6,7 +6,7 @@ cimport numpy as cnp cdef float integrate(cnp.ndarray[float, ndim=2, mode="c"] sat, - int r0, int c0, int r1, int c1): + Py_ssize_t r0, Py_ssize_t c0, Py_ssize_t r1, Py_ssize_t c1): """ Using a summed area table / integral image, calculate the sum over a given window. diff --git a/skimage/draw/_draw.pyx b/skimage/draw/_draw.pyx index e0c9231e..efd5fdc5 100644 --- a/skimage/draw/_draw.pyx +++ b/skimage/draw/_draw.pyx @@ -2,15 +2,15 @@ #cython: boundscheck=False #cython: nonecheck=False #cython: wraparound=False -import numpy as np -import math -from libc.math cimport sqrt -cimport numpy as np cimport cython +cimport numpy as np +from libc.math cimport sqrt +import math +import numpy as np from skimage._shared.geometry cimport point_in_polygon -def line(int y, int x, int y2, int x2): +def line(Py_ssize_t y, Py_ssize_t x, Py_ssize_t y2, Py_ssize_t x2): """Generate line pixel coordinates. Parameters @@ -29,12 +29,12 @@ def line(int y, int x, int y2, int x2): """ - cdef np.ndarray[np.int32_t, ndim=1, mode="c"] rr, cc + cdef np.ndarray[np.intp_t, ndim=1, mode="c"] rr, cc - cdef int steep = 0 - cdef int dx = abs(x2 - x) - cdef int dy = abs(y2 - y) - cdef int sx, sy, d, i + cdef char steep = 0 + cdef Py_ssize_t dx = abs(x2 - x) + cdef Py_ssize_t dy = abs(y2 - y) + cdef Py_ssize_t sx, sy, d, i if (x2 - x) > 0: sx = 1 @@ -51,8 +51,8 @@ def line(int y, int x, int y2, int x2): sx, sy = sy, sx d = (2 * dy) - dx - rr = np.zeros(int(dx) + 1, dtype=np.int32) - cc = np.zeros(int(dx) + 1, dtype=np.int32) + rr = np.zeros(int(dx) + 1, dtype=np.intp) + cc = np.zeros(int(dx) + 1, dtype=np.intp) for i in range(dx): if steep: @@ -96,18 +96,18 @@ def polygon(y, x, shape=None): """ - cdef int nr_verts = x.shape[0] - cdef int minr = max(0, y.min()) - cdef int maxr = math.ceil(y.max()) - cdef int minc = max(0, x.min()) - cdef int maxc = math.ceil(x.max()) + cdef Py_ssize_t nr_verts = x.shape[0] + cdef Py_ssize_t minr = int(max(0, y.min())) + cdef Py_ssize_t maxr = int(math.ceil(y.max())) + cdef Py_ssize_t minc = int(max(0, x.min())) + cdef Py_ssize_t maxc = int(math.ceil(x.max())) # make sure output coordinates do not exceed image size if shape is not None: maxr = min(shape[0] - 1, maxr) maxc = min(shape[1] - 1, maxc) - cdef int r, c + cdef Py_ssize_t r, c #: make contigous arrays for r, c coordinates cdef np.ndarray contiguous_rdata, contiguous_cdata @@ -148,17 +148,17 @@ def ellipse(double cy, double cx, double yradius, double xradius, shape=None): """ - cdef int minr = max(0, cy - yradius) - cdef int maxr = math.ceil(cy + yradius) - cdef int minc = max(0, cx - xradius) - cdef int maxc = math.ceil(cx + xradius) + cdef Py_ssize_t minr = int(max(0, cy - yradius)) + cdef Py_ssize_t maxr = int(math.ceil(cy + yradius)) + cdef Py_ssize_t minc = int(max(0, cx - xradius)) + cdef Py_ssize_t maxc = int(math.ceil(cx + xradius)) # make sure output coordinates do not exceed image size if shape is not None: maxr = min(shape[0] - 1, maxr) maxc = min(shape[1] - 1, maxc) - cdef int r, c + cdef Py_ssize_t r, c #: output coordinate arrays cdef list rr = list() @@ -195,7 +195,8 @@ def circle(double cy, double cx, double radius, shape=None): return ellipse(cy, cx, radius, radius, shape) -def circle_perimeter(int cy, int cx, int radius, method='bresenham'): +def circle_perimeter(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t radius, + method='bresenham'): """Generate circle perimeter coordinates. Parameters @@ -234,9 +235,9 @@ def circle_perimeter(int cy, int cx, int radius, method='bresenham'): cdef list rr = list() cdef list cc = list() - cdef int x = 0 - cdef int y = radius - cdef int d = 0 + cdef Py_ssize_t x = 0 + cdef Py_ssize_t y = radius + cdef Py_ssize_t d = 0 cdef char cmethod if method == 'bresenham': d = 3 - 2 * radius @@ -273,7 +274,8 @@ def circle_perimeter(int cy, int cx, int radius, method='bresenham'): return np.array(rr) + cy, np.array(cc) + cx -def ellipse_perimeter(int cy, int cx, int yradius, int xradius): +def ellipse_perimeter(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t yradius, + Py_ssize_t xradius): """Generate ellipse perimeter coordinates. Parameters @@ -302,8 +304,8 @@ def ellipse_perimeter(int cy, int cx, int yradius, int xradius): return np.array(cy), np.array(cx) # a and b are xradius an yradius compute 2a^2 and 2b^2 - cdef int twoasquared = 2 * xradius**2 - cdef int twobsquared = 2 * yradius**2 + cdef Py_ssize_t twoasquared = 2 * xradius**2 + cdef Py_ssize_t twobsquared = 2 * yradius**2 # Pixels cdef list px = list() @@ -311,14 +313,14 @@ def ellipse_perimeter(int cy, int cx, int yradius, int xradius): # First set of points: # start at the top - cdef int x = xradius - cdef int y = 0 + cdef Py_ssize_t x = xradius + cdef Py_ssize_t y = 0 - cdef int err = 0 - cdef int xstop = twobsquared * xradius - cdef int ystop = 0 - cdef int xchange = yradius * yradius * (1 - 2 * xradius) - cdef int ychange = xradius * xradius + cdef Py_ssize_t err = 0 + cdef Py_ssize_t xstop = twobsquared * xradius + cdef Py_ssize_t ystop = 0 + cdef Py_ssize_t xchange = yradius * yradius * (1 - 2 * xradius) + cdef Py_ssize_t ychange = xradius * xradius while xstop > ystop: px.extend([x, -x, -x, x]) diff --git a/skimage/feature/_daisy.py b/skimage/feature/_daisy.py index f83605b5..1bb8cbf5 100644 --- a/skimage/feature/_daisy.py +++ b/skimage/feature/_daisy.py @@ -180,7 +180,7 @@ def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8, color = (1, 0, 0) desc_y = i * step + radius desc_x = j * step + radius - coords = draw.circle_perimeter(desc_y, desc_x, sigmas[0]) + coords = draw.circle_perimeter(desc_y, desc_x, int(sigmas[0])) draw.set_color(descs_img, coords, color) max_bin = np.max(descs[i, j, :]) for o_num, o in enumerate(orientation_angles): @@ -188,8 +188,8 @@ def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8, bin_size = descs[i, j, o_num] / max_bin dy = sigmas[0] * bin_size * sin(o) dx = sigmas[0] * bin_size * cos(o) - coords = draw.line(desc_y, desc_x, desc_y + dy, - desc_x + dx) + coords = draw.line(desc_y, desc_x, int(desc_y + dy), + int(desc_x + dx)) draw.set_color(descs_img, coords, color) for r_num, r in enumerate(ring_radii): color_offset = float(1 + r_num) / rings @@ -199,7 +199,7 @@ def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8, hist_y = desc_y + int(round(r * sin(t))) hist_x = desc_x + int(round(r * cos(t))) coords = draw.circle_perimeter(hist_y, hist_x, - sigmas[r_num + 1]) + int(sigmas[r_num + 1])) draw.set_color(descs_img, coords, color) for o_num, o in enumerate(orientation_angles): # Draw histogram bins @@ -209,8 +209,9 @@ def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8, bin_size /= max_bin dy = sigmas[r_num + 1] * bin_size * sin(o) dx = sigmas[r_num + 1] * bin_size * cos(o) - coords = draw.line(hist_y, hist_x, hist_y + dy, - hist_x + dx) + coords = draw.line(hist_y, hist_x, + int(hist_y + dy), + int(hist_x + dx)) draw.set_color(descs_img, coords, color) return descs, descs_img else: diff --git a/skimage/feature/_hog.py b/skimage/feature/_hog.py index 8e8ea7f2..9fa018b7 100644 --- a/skimage/feature/_hog.py +++ b/skimage/feature/_hog.py @@ -142,8 +142,10 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8), centre = tuple([y * cy + cy // 2, x * cx + cx // 2]) dx = radius * cos(float(o) / orientations * np.pi) dy = radius * sin(float(o) / orientations * np.pi) - rr, cc = draw.bresenham(centre[0] - dx, centre[1] - dy, - centre[0] + dx, centre[1] + dy) + rr, cc = draw.bresenham(int(centre[0] - dx), + int(centre[1] - dy), + int(centre[0] + dx), + int(centre[1] + dy)) hog_image[rr, cc] += orientation_histogram[y, x, o] """ diff --git a/skimage/feature/_template.pyx b/skimage/feature/_template.pyx index 58d48524..842e1eda 100644 --- a/skimage/feature/_template.pyx +++ b/skimage/feature/_template.pyx @@ -42,12 +42,17 @@ from skimage._shared.transform cimport integrate @cython.boundscheck(False) def match_template(np.ndarray[float, ndim=2, mode="c"] image, np.ndarray[float, ndim=2, mode="c"] template): + cdef np.ndarray[float, ndim=2, mode="c"] corr cdef np.ndarray[float, ndim=2, mode="c"] image_sat cdef np.ndarray[float, ndim=2, mode="c"] image_sqr_sat cdef float template_mean = np.mean(template) cdef float template_ssd cdef float inv_area + cdef Py_ssize_t r, c, r_end, c_end + cdef Py_ssize_t template_rows = template.shape[0] + cdef Py_ssize_t template_cols = template.shape[1] + cdef float den, window_sqr_sum, window_mean_sqr, window_sum image_sat = integral.integral_image(image) image_sqr_sat = integral.integral_image(image**2) @@ -63,24 +68,24 @@ def match_template(np.ndarray[float, ndim=2, mode="c"] image, mode="valid"), dtype=np.float32) - cdef int i, j - cdef float den, window_sqr_sum, window_mean_sqr, window_sum, - # move window through convolution results, normalizing in the process - for i in range(corr.shape[0]): - for j in range(corr.shape[1]): - # subtract 1 because `i_end` and `j_end` are used for indexing into - # summed-area table, instead of slicing windows of the image. - i_end = i + template.shape[0] - 1 - j_end = j + template.shape[1] - 1 - window_sum = integrate(image_sat, i, j, i_end, j_end) + # move window through convolution results, normalizing in the process + for r in range(corr.shape[0]): + for c in range(corr.shape[1]): + # subtract 1 because `i_end` and `c_end` are used for indexing into + # summed-area table, instead of slicing windows of the image. + r_end = r + template_rows - 1 + c_end = c + template_cols - 1 + + window_sum = integrate(image_sat, r, c, r_end, c_end) window_mean_sqr = window_sum * window_sum * inv_area - window_sqr_sum = integrate(image_sqr_sat, i, j, i_end, j_end) + window_sqr_sum = integrate(image_sqr_sat, r, c, r_end, c_end) if window_sqr_sum <= window_mean_sqr: - corr[i, j] = 0 + corr[r, c] = 0 continue den = sqrt((window_sqr_sum - window_mean_sqr) * template_ssd) - corr[i, j] /= den + corr[r, c] /= den + return corr diff --git a/skimage/feature/_texture.pyx b/skimage/feature/_texture.pyx index 70a446bb..ffff38d7 100644 --- a/skimage/feature/_texture.pyx +++ b/skimage/feature/_texture.pyx @@ -16,8 +16,7 @@ def _glcm_loop(np.ndarray[dtype=np.uint8_t, ndim=2, negative_indices=False, mode='c'] angles, int levels, np.ndarray[dtype=np.uint32_t, ndim=4, - negative_indices=False, mode='c'] out - ): + negative_indices=False, mode='c'] out): """Perform co-occurrence matrix accumulation. Parameters @@ -37,23 +36,26 @@ def _glcm_loop(np.ndarray[dtype=np.uint8_t, ndim=2, the results of the GLCM computation. """ + cdef: - np.int32_t a_inx, d_idx - np.int32_t r, c, rows, cols, row, col - np.int32_t i, j + Py_ssize_t a_idx, d_idx, r, c, rows, cols, row, col + np.uint8_t i, j + np.float64_t angle, distance rows = image.shape[0] cols = image.shape[1] - for a_idx, angle in enumerate(angles): - for d_idx, distance in enumerate(distances): + 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] # compute the location of the offset pixel row = r + (sin(angle) * distance + 0.5) - col = c + (cos(angle) * distance + 0.5); + col = c + (cos(angle) * distance + 0.5) # make sure the offset is within bounds if row >= 0 and row < rows and \ @@ -123,11 +125,11 @@ def _local_binary_pattern(np.ndarray[double, ndim=2] image, output_shape = (image.shape[0], image.shape[1]) cdef np.ndarray[double, ndim=2] output = np.zeros(output_shape, np.double) - cdef int rows = image.shape[0] - cdef int cols = image.shape[1] + cdef Py_ssize_t rows = image.shape[0] + cdef Py_ssize_t cols = image.shape[1] cdef double lbp - cdef int r, c, changes, i + cdef Py_ssize_t r, c, changes, i for r in range(image.shape[0]): for c in range(image.shape[1]): for i in range(P): diff --git a/skimage/feature/corner_cy.pyx b/skimage/feature/corner_cy.pyx index 2f001ea4..407e6642 100644 --- a/skimage/feature/corner_cy.pyx +++ b/skimage/feature/corner_cy.pyx @@ -10,7 +10,7 @@ from skimage.color import rgb2grey from skimage.util import img_as_float -def corner_moravec(image, int window_size=1): +def corner_moravec(image, Py_ssize_t window_size=1): """Compute Moravec corner measure response image. This is one of the simplest corner detectors and is comparatively fast but @@ -56,8 +56,8 @@ def corner_moravec(image, int window_size=1): [ 0., 0., 0., 0., 0., 0., 0.]]) """ - cdef int rows = image.shape[0] - cdef int cols = image.shape[1] + cdef Py_ssize_t rows = image.shape[0] + cdef Py_ssize_t cols = image.shape[1] cdef cnp.ndarray[dtype=cnp.double_t, ndim=2, mode='c'] cimage, out @@ -71,7 +71,7 @@ def corner_moravec(image, int window_size=1): cdef double* out_data = out.data cdef double msum, min_msum - cdef int r, c, br, bc, mr, mc, a, b + 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 diff --git a/skimage/feature/tests/test_hog.py b/skimage/feature/tests/test_hog.py index 6f2d4cdf..cfef2da0 100644 --- a/skimage/feature/tests/test_hog.py +++ b/skimage/feature/tests/test_hog.py @@ -102,7 +102,7 @@ def test_hog_orientations_circle(): width = height = 100 image = np.zeros((height, width)) - rr, cc = draw.circle(height/2, width/2, width/3) + rr, cc = draw.circle(int(height / 2), int(width / 2), int(width / 3)) image[rr, cc] = 100 image = ndimage.gaussian_filter(image, 2) diff --git a/skimage/filter/_ctmf.pyx b/skimage/filter/_ctmf.pyx index 5a81cb98..65feb9e4 100644 --- a/skimage/filter/_ctmf.pyx +++ b/skimage/filter/_ctmf.pyx @@ -84,9 +84,9 @@ cdef struct PixelCount: # relative offsets from the octagon center # cdef struct SCoord: - np.int32_t stride # add the stride to the memory location - np.int32_t x - np.int32_t y + Py_ssize_t stride # add the stride to the memory location + Py_ssize_t x + Py_ssize_t y cdef struct Histograms: void *memory # pointer to the allocated memory @@ -95,16 +95,16 @@ cdef struct Histograms: np.uint8_t *data # pointer to the image data np.uint8_t *mask # pointer to the image mask np.uint8_t *output # pointer to the output array - np.int32_t column_count # number of columns represented by this + Py_ssize_t column_count # number of columns represented by this # structure - np.int32_t stripe_length # number of columns including "radius" before + Py_ssize_t stripe_length # number of columns including "radius" before # and after - np.int32_t row_count # number of rows available in image - np.int32_t current_column # the column being processed - np.int32_t current_row # the row being processed - np.int32_t current_stride # offset in data and mask to current location - np.int32_t radius # the "radius" of the octagon - np.int32_t a_2 # 1/2 of the length of a side of the octagon + Py_ssize_t row_count # number of rows available in image + Py_ssize_t current_column # the column being processed + Py_ssize_t current_row # the row being processed + Py_ssize_t current_stride # offset in data and mask to current location + Py_ssize_t radius # the "radius" of the octagon + Py_ssize_t a_2 # 1/2 of the length of a side of the octagon # # # The strides are the offsets in the array to the points that need to @@ -138,50 +138,50 @@ cdef struct Histograms: SCoord last_bottom_left # (-) trailing edge bottom - 1 col SCoord bottom_left # (+) left side of octagon's bottom - 1 col - np.int32_t row_stride # stride between one row and the next - np.int32_t col_stride # stride between one column and the next + Py_ssize_t row_stride # stride between one row and the next + Py_ssize_t col_stride # stride between one column and the next # The accumulator holds the running histogram # HistogramPiece accumulator # # The running count of pixels in the accumulator # - np.uint32_t accumulator_count + Py_ssize_t accumulator_count # # The percent of pixels within the octagon whose value is # less than or equal to the median-filtered value (e.g. for # median, this is 50, for lower quartile it's 25) # - np.int32_t percent + Py_ssize_t percent # # last_update_column keeps track of the column # of the last update # to the fine histogram accumulator. Short-term, the median # stays in one coarse block so only one fine histogram might # need to be updated # - np.int32_t last_update_column[16] + Py_ssize_t last_update_column[16] ############################################################################ # # allocate_histograms - allocates the Histograms structure for the run # ############################################################################ -cdef Histograms *allocate_histograms(np.int32_t rows, - np.int32_t columns, - np.int32_t row_stride, - np.int32_t col_stride, - np.int32_t radius, - np.int32_t percent, +cdef Histograms *allocate_histograms(Py_ssize_t rows, + Py_ssize_t columns, + Py_ssize_t row_stride, + Py_ssize_t col_stride, + Py_ssize_t radius, + Py_ssize_t percent, np.uint8_t *data, np.uint8_t *mask, np.uint8_t *output): cdef: - unsigned int adjusted_stripe_length = columns + 2*radius + 1 - unsigned int memory_size + Py_ssize_t adjusted_stripe_length = columns + 2*radius + 1 + Py_ssize_t memory_size void *ptr Histograms *ph size_t roundoff - int a + Py_ssize_t a SCoord *psc memory_size = (adjusted_stripe_length * @@ -199,7 +199,7 @@ cdef Histograms *allocate_histograms(np.int32_t rows, # # Align histogram memory to a 32-byte boundary # - roundoff = ptr + roundoff = ptr roundoff += 31 roundoff -= roundoff % 32 ptr = roundoff @@ -232,7 +232,7 @@ cdef Histograms *allocate_histograms(np.int32_t rows, # a_2 is the offset from the center to each of the octagon # corners # - a = ( radius * 2.0 / 2.414213) + a = (radius * 2.0 / 2.414213) a_2 = a / 2 if a_2 == 0: a_2 = 1 @@ -326,19 +326,18 @@ cdef void set_stride(Histograms *ph, SCoord *psc): # a column that is "radius" to the left. # ############################################################################ -cdef inline np.int32_t tl_br_colidx(Histograms *ph, np.int32_t colidx): - return (colidx + 3 * ph.radius + ph.current_row) % \ - ph.stripe_length +cdef inline Py_ssize_t tl_br_colidx(Histograms *ph, Py_ssize_t colidx): + return (colidx + 3*ph.radius + ph.current_row) % ph.stripe_length -cdef inline np.int32_t tr_bl_colidx(Histograms *ph, np.int32_t colidx): - return (colidx + 3 * ph.radius + ph.row_count - ph.current_row) % \ - ph.stripe_length +cdef inline Py_ssize_t tr_bl_colidx(Histograms *ph, Py_ssize_t colidx): + return (colidx + 3*ph.radius + ph.row_count-ph.current_row) % \ + ph.stripe_length -cdef inline np.int32_t leading_edge_colidx(Histograms *ph, np.int32_t colidx): - return (colidx + 5 * ph.radius) % ph.stripe_length +cdef inline Py_ssize_t leading_edge_colidx(Histograms *ph, Py_ssize_t colidx): + return (colidx + 5*ph.radius) % ph.stripe_length -cdef inline np.int32_t trailing_edge_colidx(Histograms *ph, np.int32_t colidx): - return (colidx + 3 * ph.radius - 1) % ph.stripe_length +cdef inline Py_ssize_t trailing_edge_colidx(Histograms *ph, Py_ssize_t colidx): + return (colidx + 3*ph.radius - 1) % ph.stripe_length ############################################################################ # @@ -349,9 +348,8 @@ cdef inline np.int32_t trailing_edge_colidx(Histograms *ph, np.int32_t colidx): # colidx - the index of the column to add # ############################################################################ -cdef inline void accumulate_coarse_histogram(Histograms *ph, np.int32_t colidx): - cdef: - int offset +cdef inline void accumulate_coarse_histogram(Histograms *ph, Py_ssize_t colidx): + cdef Py_ssize_t offset offset = tr_bl_colidx(ph, colidx) if ph.pixel_count[offset].top_right > 0: @@ -372,9 +370,8 @@ cdef inline void accumulate_coarse_histogram(Histograms *ph, np.int32_t colidx): # for a given column # ############################################################################ -cdef inline void deaccumulate_coarse_histogram(Histograms *ph, np.int32_t colidx): - cdef: - int offset +cdef inline void deaccumulate_coarse_histogram(Histograms *ph, Py_ssize_t colidx): + cdef Py_ssize_t offset # # The trailing diagonals don't appear until here # @@ -403,11 +400,11 @@ cdef inline void deaccumulate_coarse_histogram(Histograms *ph, np.int32_t colidx # ############################################################################ cdef inline void accumulate_fine_histogram(Histograms *ph, - np.int32_t colidx, - np.uint32_t fineidx): + Py_ssize_t colidx, + Py_ssize_t fineidx): cdef: - int fineoffset = fineidx * 16 - int offset + Py_ssize_t fineoffset = fineidx * 16 + Py_ssize_t offset offset = tr_bl_colidx(ph, colidx) add16(ph.accumulator.fine + fineoffset, @@ -427,11 +424,11 @@ cdef inline void accumulate_fine_histogram(Histograms *ph, # ############################################################################ cdef inline void deaccumulate_fine_histogram(Histograms *ph, - np.int32_t colidx, - np.uint32_t fineidx): + Py_ssize_t colidx, + Py_ssize_t fineidx): cdef: - int fineoffset = fineidx * 16 - int offset + Py_ssize_t fineoffset = fineidx * 16 + Py_ssize_t offset # # The trailing diagonals don't appear until here @@ -459,10 +456,7 @@ cdef inline void deaccumulate_fine_histogram(Histograms *ph, ############################################################################ cdef inline void accumulate(Histograms *ph): - cdef: - int i - int j - np.int32_t accumulator + cdef np.int32_t accumulator accumulate_coarse_histogram(ph, ph.current_column) deaccumulate_coarse_histogram(ph, ph.current_column) @@ -486,11 +480,11 @@ cdef inline void accumulate(Histograms *ph): # to choose remains to be done. ############################################################################ -cdef inline void update_fine(Histograms *ph, int fineidx): +cdef inline void update_fine(Histograms *ph, Py_ssize_t fineidx): cdef: - int first_update_column = ph.last_update_column[fineidx]+1 - int update_limit = ph.current_column+1 - int i + Py_ssize_t first_update_column = ph.last_update_column[fineidx]+1 + Py_ssize_t update_limit = ph.current_column+1 + Py_ssize_t i for i in range(first_update_column, update_limit): accumulate_fine_histogram(ph, i, fineidx) @@ -515,15 +509,15 @@ cdef inline void update_histogram(Histograms *ph, SCoord *last_coord, SCoord *coord): cdef: - np.int32_t current_column = ph.current_column - np.int32_t current_row = ph.current_row - np.int32_t current_stride = ph.current_stride - np.int32_t column_count = ph.column_count - np.int32_t row_count = ph.row_count + Py_ssize_t current_column = ph.current_column + Py_ssize_t current_row = ph.current_row + Py_ssize_t current_stride = ph.current_stride + Py_ssize_t column_count = ph.column_count + Py_ssize_t row_count = ph.row_count np.uint8_t value - np.int32_t stride - np.int32_t x - np.int32_t y + Py_ssize_t stride + Py_ssize_t x + Py_ssize_t y x = last_coord.x + current_column y = last_coord.y + current_row @@ -556,21 +550,21 @@ cdef inline void update_histogram(Histograms *ph, ############################################################################ cdef inline void update_current_location(Histograms *ph): cdef: - np.int32_t current_column = ph.current_column - np.int32_t radius = ph.radius - np.int32_t top_left_off = tl_br_colidx(ph, current_column) - np.int32_t top_right_off = tr_bl_colidx(ph, current_column) - np.int32_t bottom_left_off = tr_bl_colidx(ph, current_column) - np.int32_t bottom_right_off = tl_br_colidx(ph, current_column) - np.int32_t leading_edge_off = leading_edge_colidx(ph, current_column) + Py_ssize_t current_column = ph.current_column + Py_ssize_t radius = ph.radius + Py_ssize_t top_left_off = tl_br_colidx(ph, current_column) + Py_ssize_t top_right_off = tr_bl_colidx(ph, current_column) + Py_ssize_t bottom_left_off = tr_bl_colidx(ph, current_column) + Py_ssize_t bottom_right_off = tl_br_colidx(ph, current_column) + Py_ssize_t leading_edge_off = leading_edge_colidx(ph, current_column) np.int32_t *coarse_histogram np.int32_t *fine_histogram - np.int32_t last_xoff - np.int32_t last_yoff - np.int32_t last_stride - np.int32_t xoff - np.int32_t yoff - np.int32_t stride + Py_ssize_t last_xoff + Py_ssize_t last_yoff + Py_ssize_t last_stride + Py_ssize_t xoff + Py_ssize_t yoff + Py_ssize_t stride update_histogram(ph, &ph.histogram[top_left_off].top_left, &ph.pixel_count[top_left_off].top_left, @@ -605,16 +599,18 @@ cdef inline void update_current_location(Histograms *ph): cdef inline np.uint8_t find_median(Histograms *ph): cdef: - np.uint32_t pixels_below # of pixels below the median - int i - int j - int k + Py_ssize_t pixels_below # of pixels below the median + Py_ssize_t i + Py_ssize_t j + Py_ssize_t k np.uint32_t accumulator if ph.accumulator_count == 0: return 0 - pixels_below = ((ph.accumulator_count * ph.percent + 50) - / 100) # +50 for roundoff + + # +50 for roundoff + pixels_below = (ph.accumulator_count * ph.percent + 50) / 100 + if pixels_below > 0: pixels_below -= 1 @@ -626,10 +622,10 @@ cdef inline np.uint8_t find_median(Histograms *ph): accumulator -= ph.accumulator.coarse[i] update_fine(ph, i) - for j in range(i * 16, (i + 1) * 16): + for j in range(i*16, (i + 1)*16): accumulator += ph.accumulator.fine[j] if accumulator > pixels_below: - return j + return j return 0 @@ -648,23 +644,25 @@ cdef inline np.uint8_t find_median(Histograms *ph): # output - array to be filled with filtered pixels # ############################################################################ -cdef int c_median_filter(np.int32_t rows, - np.int32_t columns, - np.int32_t row_stride, - np.int32_t col_stride, - np.int32_t radius, - np.int32_t percent, +cdef int c_median_filter(Py_ssize_t rows, + Py_ssize_t columns, + Py_ssize_t row_stride, + Py_ssize_t col_stride, + Py_ssize_t radius, + Py_ssize_t percent, np.uint8_t *data, np.uint8_t *mask, np.uint8_t *output): cdef: Histograms *ph Histogram *phistogram - int row, col, i - np.int32_t top_left_off - np.int32_t top_right_off - np.int32_t bottom_left_off - np.int32_t bottom_right_off + Py_ssize_t row + Py_ssize_t col + Py_ssize_t i + Py_ssize_t top_left_off + Py_ssize_t top_right_off + Py_ssize_t bottom_left_off + Py_ssize_t bottom_right_off ph = allocate_histograms(rows, columns, row_stride, col_stride, radius, percent, data, mask, output) diff --git a/skimage/filter/_denoise_cy.pyx b/skimage/filter/_denoise_cy.pyx index 8ebf5c73..5a85e84a 100644 --- a/skimage/filter/_denoise_cy.pyx +++ b/skimage/filter/_denoise_cy.pyx @@ -17,7 +17,7 @@ cdef inline double _gaussian_weight(double sigma, double value): return exp(-0.5 * (value / sigma)**2) -cdef double* _compute_color_lut(int bins, double sigma, double max_value): +cdef double* _compute_color_lut(Py_ssize_t bins, double sigma, double max_value): cdef: double* color_lut = malloc(bins * sizeof(double)) @@ -29,7 +29,7 @@ cdef double* _compute_color_lut(int bins, double sigma, double max_value): return color_lut -cdef double* _compute_range_lut(int win_size, double sigma): +cdef double* _compute_range_lut(Py_ssize_t win_size, double sigma): cdef: double* range_lut = malloc(win_size**2 * sizeof(double)) @@ -45,9 +45,9 @@ cdef double* _compute_range_lut(int win_size, double sigma): return range_lut -def denoise_bilateral(image, int win_size=5, sigma_range=None, - double sigma_spatial=1, int bins=10000, mode='constant', - double cval=0): +def denoise_bilateral(image, Py_ssize_t win_size=5, sigma_range=None, + double sigma_spatial=1, Py_ssize_t bins=10000, + mode='constant', double cval=0): """Denoise image using bilateral filter. This is an edge-preserving and noise reducing denoising filter. It averages diff --git a/skimage/graph/heap.pxd b/skimage/graph/heap.pxd index 51b4f79c..e31aa150 100644 --- a/skimage/graph/heap.pxd +++ b/skimage/graph/heap.pxd @@ -1,7 +1,7 @@ """ This is the definition file for heap.pyx. It contains the definitions of the heap classes, such that other cython modules can "cimport heap" and thus use the -C versions of pop(), push(), and value_of(): pop_fast(), push_fast() and +C versions of pop(), push(), and value_of(): pop_fast(), push_fast() and value_of_fast() """ @@ -14,16 +14,16 @@ ctypedef unsigned char LEVELS_T cdef class BinaryHeap: cdef readonly INDEX_T count - cdef readonly LEVELS_T levels, min_levels + cdef readonly LEVELS_T levels, min_levels cdef VALUE_T *_values 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 INDEX_T push_fast(self, VALUE_T value, REFERENCE_T reference) cdef VALUE_T pop_fast(self) @@ -32,8 +32,7 @@ cdef class FastUpdateBinaryHeap(BinaryHeap): cdef INDEX_T *_crossref cdef BOOL_T _invalid_ref cdef BOOL_T _pushed - + cdef VALUE_T value_of_fast(self, REFERENCE_T reference) - cdef INDEX_T push_if_lower_fast(self, VALUE_T value, + cdef INDEX_T push_if_lower_fast(self, VALUE_T value, REFERENCE_T reference) - \ No newline at end of file diff --git a/skimage/measure/_find_contours.pyx b/skimage/measure/_find_contours.pyx index 4f1b3cee..91ba318c 100644 --- a/skimage/measure/_find_contours.pyx +++ b/skimage/measure/_find_contours.pyx @@ -1,10 +1,7 @@ -# -*- python -*- # cython: cdivision=True - import numpy as np cimport numpy as np -np.import_array() cdef inline double _get_fraction(double from_value, double to_value, double level): @@ -14,7 +11,7 @@ cdef inline double _get_fraction(double from_value, double to_value, def iterate_and_store(np.ndarray[double, ndim=2] array, - double level, int vertex_connect_high): + double level, Py_ssize_t vertex_connect_high): """Iterate across the given array in a marching-squares fashion, looking for segments that cross 'level'. If such a segment is found, its coordinates are added to a growing list of segments, @@ -27,7 +24,7 @@ def iterate_and_store(np.ndarray[double, ndim=2] array, raise ValueError("Input array must be at least 2x2.") cdef list arc_list = [] - cdef int n + cdef Py_ssize_t n # The plan is to iterate a 2x2 square across the input array. This means # that the upper-left corner of the square needs to iterate across a @@ -39,17 +36,17 @@ def iterate_and_store(np.ndarray[double, ndim=2] array, # index varies the fastest). # Current coords start at 0,0. - cdef int[2] coords + cdef Py_ssize_t[2] coords coords[0] = 0 coords[1] = 0 # Calculate the number of iterations we'll need - cdef int num_square_steps = (array.shape[0] - 1) * (array.shape[1] - 1) + cdef Py_ssize_t num_square_steps = (array.shape[0] - 1) * (array.shape[1] - 1) cdef unsigned char square_case = 0 cdef tuple top, bottom, left, right cdef double ul, ur, ll, lr - cdef int r0, r1, c0, c1 + cdef Py_ssize_t r0, r1, c0, c1 for n in range(num_square_steps): # There are sixteen different possible square types, diagramed below. diff --git a/skimage/measure/_moments.pyx b/skimage/measure/_moments.pyx index f84e14dd..b94a7b89 100644 --- a/skimage/measure/_moments.pyx +++ b/skimage/measure/_moments.pyx @@ -7,7 +7,7 @@ cimport numpy as np def central_moments(np.ndarray[np.double_t, ndim=2] array, double cr, double cc, int order): - cdef int p, q, r, c + cdef Py_ssize_t p, q, r, c cdef np.ndarray[np.double_t, ndim=2] mu mu = np.zeros((order + 1, order + 1), 'double') for p in range(order + 1): @@ -18,7 +18,7 @@ def central_moments(np.ndarray[np.double_t, ndim=2] array, double cr, double cc, return mu def normalized_moments(np.ndarray[np.double_t, ndim=2] mu, int order): - cdef int p, q + cdef Py_ssize_t p, q cdef np.ndarray[np.double_t, ndim=2] nu nu = np.zeros((order + 1, order + 1), 'double') for p in range(order + 1): diff --git a/skimage/morphology/_convex_hull.pyx b/skimage/morphology/_convex_hull.pyx index dc426900..9bec825b 100644 --- a/skimage/morphology/_convex_hull.pyx +++ b/skimage/morphology/_convex_hull.pyx @@ -13,47 +13,44 @@ def possible_hull(np.ndarray[dtype=np.uint8_t, ndim=2, mode="c"] img): Returns ------- - coords : ndarray (N, 2) + coords : ndarray (cols, 2) The ``(row, column)`` coordinates of all pixels that possibly belong to the convex hull. """ - cdef int i, j, k - cdef unsigned int M, N - - M = img.shape[0] - N = img.shape[1] + cdef Py_ssize_t r, c + cdef Py_ssize_t rows = img.shape[0] + cdef Py_ssize_t cols = img.shape[1] - # Output: M storage slots for left boundary pixels - # N storage slots for top boundary pixels - # M storage slots for right boundary pixels - # N storage slots for bottom boundary pixels - cdef np.ndarray[dtype=np.int_t, ndim=2] nonzero = \ - np.ones((2 * (M + N), 2), dtype=np.int) - nonzero *= -1 + # Output: rows storage slots for left boundary pixels + # cols storage slots for top boundary pixels + # rows storage slots for right boundary pixels + # cols storage slots for bottom boundary pixels + cdef np.ndarray[dtype=np.intp_t, ndim=2] nonzero = \ + np.ones((2 * (rows + cols), 2), dtype=np.int) + nonzero *= -1 - k = 0 - for i in range(M): - for j in range(N): - if img[i, j] != 0: + for r in range(rows): + for c in range(cols): + if img[r, c] != 0: # Left check - if nonzero[i, 1] == -1: - nonzero[i, 0] = i - nonzero[i, 1] = j + if nonzero[r, 1] == -1: + nonzero[r, 0] = r + nonzero[r, 1] = c # Right check - elif nonzero[M + N + i, 1] < j: - nonzero[M + N + i, 0] = i - nonzero[M + N + i, 1] = j + elif nonzero[rows + cols + r, 1] < c: + nonzero[rows + cols + r, 0] = r + nonzero[rows + cols + r, 1] = c # Top check - if nonzero[M + j, 1] == -1: - nonzero[M + j, 0] = i - nonzero[M + j, 1] = j + if nonzero[rows + c, 1] == -1: + nonzero[rows + c, 0] = r + nonzero[rows + c, 1] = c # Bottom check - elif nonzero[2 * M + N + j, 0] < i: - nonzero[2 * M + N + j, 0] = i - nonzero[2 * M + N + j, 1] = j - + elif nonzero[2 * rows + cols + c, 0] < r: + nonzero[2 * rows + cols + c, 0] = r + nonzero[2 * rows + cols + c, 1] = c + return nonzero[nonzero[:, 0] != -1] diff --git a/skimage/morphology/_pnpoly.pyx b/skimage/morphology/_pnpoly.pyx index 7deb6a50..f6cf386f 100644 --- a/skimage/morphology/_pnpoly.pyx +++ b/skimage/morphology/_pnpoly.pyx @@ -31,18 +31,19 @@ def grid_points_inside_poly(shape, verts): vx = verts[:, 0].astype(np.double) vy = verts[:, 1].astype(np.double) - cdef int V = vx.shape[0] + cdef Py_ssize_t V = vx.shape[0] - cdef int M = shape[0] - cdef int N = shape[1] - cdef int m, n + cdef Py_ssize_t M = shape[0] + cdef Py_ssize_t N = shape[1] + cdef Py_ssize_t m, n cdef np.ndarray[dtype=np.uint8_t, ndim=2, mode="c"] out = \ np.zeros((M, N), dtype=np.uint8) for m in range(M): for n in range(N): - out[m, n] = point_in_polygon(V, vx.data, vy.data, m, n) + out[m, n] = point_in_polygon(V, vx.data, vy.data, + m, n) return out.view(bool) @@ -76,7 +77,7 @@ def points_inside_poly(points, verts): vy = verts[:, 1].astype(np.double) cdef np.ndarray[np.uint8_t, ndim=1] out = \ - np.zeros(x.shape[0], dtype=np.uint8) + np.zeros(x.shape[0], dtype=np.uint8) points_in_polygon(vx.shape[0], vx.data, vy.data, x.shape[0], x.data, y.data, diff --git a/skimage/morphology/_skeletonize.py b/skimage/morphology/_skeletonize.py index 04a65da6..b48beb86 100644 --- a/skimage/morphology/_skeletonize.py +++ b/skimage/morphology/_skeletonize.py @@ -277,8 +277,8 @@ def medial_axis(image, mask=None, return_distance=False): i, j = np.mgrid[0:image.shape[0], 0:image.shape[1]] result = masked_image.copy() distance = distance[result] - i = np.ascontiguousarray(i[result], np.int32) - j = np.ascontiguousarray(j[result], np.int32) + i = np.ascontiguousarray(i[result], np.intp) + j = np.ascontiguousarray(j[result], np.intp) result = np.ascontiguousarray(result, np.uint8) # Determine the order in which pixels are processed. diff --git a/skimage/morphology/_skeletonize_cy.pyx b/skimage/morphology/_skeletonize_cy.pyx index ff5fcdf2..9bef86b9 100644 --- a/skimage/morphology/_skeletonize_cy.pyx +++ b/skimage/morphology/_skeletonize_cy.pyx @@ -15,11 +15,11 @@ cimport cython @cython.boundscheck(False) -def _skeletonize_loop(np.ndarray[dtype=np.uint8_t, ndim=2, +def _skeletonize_loop(np.ndarray[dtype=np.uint8_t, ndim=2, negative_indices=False, mode='c'] result, - np.ndarray[dtype=np.int32_t, ndim=1, + np.ndarray[dtype=np.intp_t, ndim=1, negative_indices=False, mode='c'] i, - np.ndarray[dtype=np.int32_t, ndim=1, + np.ndarray[dtype=np.intp_t, ndim=1, negative_indices=False, mode='c'] j, np.ndarray[dtype=np.int32_t, ndim=1, negative_indices=False, mode='c'] order, @@ -37,13 +37,13 @@ def _skeletonize_loop(np.ndarray[dtype=np.uint8_t, ndim=2, i, j : ndarrays The coordinates of each foreground pixel in the image - + order : ndarray The index of each pixel, in the order of processing (order[0] is the first pixel to process, etc.) - + table : ndarray - The 512-element lookup table of values after transformation + The 512-element lookup table of values after transformation (whether to keep or not each configuration in a binary 3x3 array) Notes @@ -55,15 +55,15 @@ def _skeletonize_loop(np.ndarray[dtype=np.uint8_t, ndim=2, the quench-line of the brushfire will be evaluated later than a point closer to the edge. - Note that the neighbourhood of a pixel may evolve before the loop - arrives at this pixel. This is why it is possible to compute the + Note that the neighbourhood of a pixel may evolve before the loop + arrives at this pixel. This is why it is possible to compute the skeleton in only one pass, thanks to an adapted ordering of the pixels. """ cdef: np.int32_t accumulator - np.int32_t index, order_index - np.int32_t ii, jj + Py_ssize_t index, order_index + Py_ssize_t ii, jj for index in range(order.shape[0]): accumulator = 16 @@ -110,21 +110,21 @@ def _table_lookup_index(np.ndarray[dtype=np.uint8_t, ndim=2, 256 128 64 32 16 8 4 2 1 - + but this runs about twice as fast because of inlining and the hardwired kernel. """ cdef: - np.ndarray[dtype=np.int32_t, ndim=2, + np.ndarray[dtype=np.int32_t, ndim=2, negative_indices=False, mode='c'] indexer np.int32_t *p_indexer np.uint8_t *p_image - np.int32_t i_stride - np.int32_t i_shape - np.int32_t j_shape - np.int32_t i - np.int32_t j - np.int32_t offset + Py_ssize_t i_stride + Py_ssize_t i_shape + Py_ssize_t j_shape + Py_ssize_t i + Py_ssize_t j + Py_ssize_t offset i_shape = image.shape[0] j_shape = image.shape[1] diff --git a/skimage/morphology/_watershed.pyx b/skimage/morphology/_watershed.pyx index c86d8744..122f0262 100644 --- a/skimage/morphology/_watershed.pyx +++ b/skimage/morphology/_watershed.pyx @@ -9,39 +9,33 @@ All rights reserved. Original author: Lee Kamentsky """ - -cdef extern from "numpy/arrayobject.h": - cdef void import_array() -import_array() - import numpy as np cimport numpy as np cimport cython -DTYPE_INT32 = np.int32 + ctypedef np.int32_t DTYPE_INT32_t DTYPE_BOOL = np.bool ctypedef np.int8_t DTYPE_BOOL_t + include "heap_watershed.pxi" + @cython.boundscheck(False) -def watershed(np.ndarray[DTYPE_INT32_t, ndim=1, negative_indices=False, - mode='c'] image, - np.ndarray[DTYPE_INT32_t, ndim=2, negative_indices=False, - mode='c'] pq, - DTYPE_INT32_t age, - np.ndarray[DTYPE_INT32_t, ndim=2, negative_indices=False, - mode='c'] structure, - DTYPE_INT32_t ndim, - np.ndarray[DTYPE_BOOL_t, ndim=1, negative_indices=False, - mode='c'] mask, - np.ndarray[DTYPE_INT32_t, ndim=1, negative_indices=False, - mode='c'] image_shape, - np.ndarray[DTYPE_INT32_t, ndim=1, negative_indices=False, - mode='c'] output): +def watershed(np.ndarray[DTYPE_INT32_t, ndim=1, negative_indices=False, + mode='c'] image, + np.ndarray[DTYPE_INT32_t, ndim=2, negative_indices=False, + mode='c'] pq, + Py_ssize_t age, + np.ndarray[DTYPE_INT32_t, ndim=2, negative_indices=False, + mode='c'] structure, + np.ndarray[DTYPE_BOOL_t, ndim=1, negative_indices=False, + mode='c'] mask, + np.ndarray[DTYPE_INT32_t, ndim=1, negative_indices=False, + mode='c'] output): """Do heavy lifting of watershed algorithm - + Parameters ---------- @@ -58,20 +52,17 @@ def watershed(np.ndarray[DTYPE_INT32_t, ndim=1, negative_indices=False, in a flattened array. The remaining elements are the offsets from the point to its neighbor in the various dimensions - ndim - # of dimensions in the image mask - numpy boolean (char) array indicating which pixels to consider and which to ignore. Also flattened. - image_shape - the dimensions of the image, for boundary checking, - a numpy array of np.int32 output - put the image labels in here """ cdef Heapitem elem cdef Heapitem new_elem - cdef DTYPE_INT32_t nneighbors = structure.shape[0] - cdef DTYPE_INT32_t i = 0 - cdef DTYPE_INT32_t index = 0 - cdef DTYPE_INT32_t old_index = 0 - cdef DTYPE_INT32_t max_index = image.shape[0] + cdef Py_ssize_t nneighbors = structure.shape[0] + cdef Py_ssize_t i = 0 + cdef Py_ssize_t index = 0 + cdef Py_ssize_t old_index = 0 + cdef Py_ssize_t max_index = image.shape[0] cdef Heap *hp = heap_from_numpy2() diff --git a/skimage/morphology/ccomp.pxd b/skimage/morphology/ccomp.pxd index 0b431832..b50703e9 100644 --- a/skimage/morphology/ccomp.pxd +++ b/skimage/morphology/ccomp.pxd @@ -1,10 +1,10 @@ """Export fast union find in Cython""" cimport numpy as np -DTYPE = np.int -ctypedef np.int_t DTYPE_t +DTYPE = np.intp +ctypedef np.intp_t DTYPE_t -cdef DTYPE_t find_root(np.int_t *forest, np.int_t n) -cdef set_root(np.int_t *forest, np.int_t n, np.int_t root) -cdef join_trees(np.int_t *forest, np.int_t n, np.int_t m) -cdef link_bg(np.int_t *forest, np.int_t n, np.int_t *background_node) +cdef DTYPE_t find_root(DTYPE_t *forest, DTYPE_t n) +cdef set_root(DTYPE_t *forest, DTYPE_t n, DTYPE_t root) +cdef join_trees(DTYPE_t *forest, DTYPE_t n, DTYPE_t m) +cdef link_bg(DTYPE_t *forest, DTYPE_t n, DTYPE_t *background_node) diff --git a/skimage/morphology/ccomp.pyx b/skimage/morphology/ccomp.pyx index 6a4fb1f2..dc85196b 100644 --- a/skimage/morphology/ccomp.pyx +++ b/skimage/morphology/ccomp.pyx @@ -23,23 +23,25 @@ See also: # Tree operations implemented by an array as described in Wu et al. # The term "forest" is used to indicate an array that stores one or more trees -DTYPE = np.int +DTYPE = np.intp -cdef DTYPE_t find_root(np.int_t *forest, np.int_t n): + +cdef DTYPE_t find_root(DTYPE_t *forest, DTYPE_t n): """Find the root of node n. """ - cdef np.int_t root = n + cdef DTYPE_t root = n while (forest[root] < root): root = forest[root] return root -cdef set_root(np.int_t *forest, np.int_t n, np.int_t root): + +cdef set_root(DTYPE_t *forest, DTYPE_t n, DTYPE_t root): """ Set all nodes on a path to point to new_root. """ - cdef np.int_t j + cdef DTYPE_t j while (forest[n] < n): j = forest[n] forest[n] = root @@ -48,12 +50,12 @@ cdef set_root(np.int_t *forest, np.int_t n, np.int_t root): forest[n] = root -cdef join_trees(np.int_t *forest, np.int_t n, np.int_t m): +cdef join_trees(DTYPE_t *forest, DTYPE_t n, DTYPE_t m): """Join two trees containing nodes n and m. """ - cdef np.int_t root = find_root(forest, n) - cdef np.int_t root_m + cdef DTYPE_t root = find_root(forest, n) + cdef DTYPE_t root_m if (n != m): root_m = find_root(forest, m) @@ -64,7 +66,8 @@ cdef join_trees(np.int_t *forest, np.int_t n, np.int_t m): set_root(forest, n, root) set_root(forest, m, root) -cdef link_bg(np.int_t *forest, np.int_t n, np.int_t *background_node): + +cdef link_bg(DTYPE_t *forest, DTYPE_t n, DTYPE_t *background_node): """ Link a node to the background node. @@ -76,7 +79,7 @@ cdef link_bg(np.int_t *forest, np.int_t n, np.int_t *background_node): # Connected components search as described in Fiorio et al. -def label(input, np.int_t neighbors=8, np.int_t background=-1): +def label(input, DTYPE_t neighbors=8, DTYPE_t background=-1): """Label connected regions of an integer array. Two pixels are connected when they are neighbors and have the same value. @@ -134,8 +137,8 @@ def label(input, np.int_t neighbors=8, np.int_t background=-1): [-1 -1 -1]] """ - cdef np.int_t rows = input.shape[0] - cdef np.int_t cols = input.shape[1] + cdef DTYPE_t rows = input.shape[0] + cdef DTYPE_t cols = input.shape[1] cdef np.ndarray[DTYPE_t, ndim=2] data = np.array(input, copy=True, dtype=DTYPE) @@ -143,12 +146,12 @@ def label(input, np.int_t neighbors=8, np.int_t background=-1): forest = np.arange(data.size, dtype=DTYPE).reshape((rows, cols)) - cdef np.int_t *forest_p = forest.data - cdef np.int_t *data_p = data.data + cdef DTYPE_t *forest_p = forest.data + cdef DTYPE_t *data_p = data.data - cdef np.int_t i, j + cdef DTYPE_t i, j - cdef np.int_t background_node = -999 + cdef DTYPE_t background_node = -999 if neighbors != 4 and neighbors != 8: raise ValueError('Neighbors must be either 4 or 8.') @@ -197,7 +200,7 @@ def label(input, np.int_t neighbors=8, np.int_t background=-1): # Label output - cdef np.int_t ctr = 0 + cdef DTYPE_t ctr = 0 for i in range(rows): for j in range(cols): if (i*cols + j) == background_node: diff --git a/skimage/morphology/cmorph.pyx b/skimage/morphology/cmorph.pyx index 9b8b3a27..a09a39a3 100644 --- a/skimage/morphology/cmorph.pyx +++ b/skimage/morphology/cmorph.pyx @@ -13,13 +13,13 @@ def dilate(np.ndarray[np.uint8_t, ndim=2] image, np.ndarray[np.uint8_t, ndim=2] out=None, char shift_x=0, char shift_y=0): - cdef int rows = image.shape[0] - cdef int cols = image.shape[1] - cdef int srows = selem.shape[0] - cdef int scols = selem.shape[1] + cdef Py_ssize_t rows = image.shape[0] + cdef Py_ssize_t cols = image.shape[1] + cdef Py_ssize_t srows = selem.shape[0] + cdef Py_ssize_t scols = selem.shape[1] - cdef int centre_r = int(selem.shape[0] / 2) - shift_y - cdef int centre_c = int(selem.shape[1] / 2) - shift_x + cdef Py_ssize_t centre_r = int(selem.shape[0] / 2) - shift_y + cdef Py_ssize_t centre_c = int(selem.shape[1] / 2) - shift_x image = np.ascontiguousarray(image) if out is None: @@ -30,11 +30,11 @@ def dilate(np.ndarray[np.uint8_t, ndim=2] image, cdef np.uint8_t* out_data = out.data cdef np.uint8_t* image_data = image.data - cdef int r, c, rr, cc, s, value, local_max + cdef Py_ssize_t r, c, rr, cc, s, value, local_max - cdef int selem_num = np.sum(selem != 0) - cdef int* sr = malloc(selem_num * sizeof(int)) - cdef int* sc = malloc(selem_num * sizeof(int)) + cdef Py_ssize_t selem_num = np.sum(selem != 0) + cdef Py_ssize_t* sr = malloc(selem_num * sizeof(Py_ssize_t)) + cdef Py_ssize_t* sc = malloc(selem_num * sizeof(Py_ssize_t)) s = 0 for r in range(srows): @@ -68,13 +68,13 @@ def erode(np.ndarray[np.uint8_t, ndim=2] image, np.ndarray[np.uint8_t, ndim=2] out=None, char shift_x=0, char shift_y=0): - cdef int rows = image.shape[0] - cdef int cols = image.shape[1] - cdef int srows = selem.shape[0] - cdef int scols = selem.shape[1] + cdef Py_ssize_t rows = image.shape[0] + cdef Py_ssize_t cols = image.shape[1] + cdef Py_ssize_t srows = selem.shape[0] + cdef Py_ssize_t scols = selem.shape[1] - cdef int centre_r = int(selem.shape[0] / 2) - shift_y - cdef int centre_c = int(selem.shape[1] / 2) - shift_x + cdef Py_ssize_t centre_r = int(selem.shape[0] / 2) - shift_y + cdef Py_ssize_t centre_c = int(selem.shape[1] / 2) - shift_x image = np.ascontiguousarray(image) if out is None: @@ -87,9 +87,9 @@ def erode(np.ndarray[np.uint8_t, ndim=2] image, cdef int r, c, rr, cc, s, value, local_min - cdef int selem_num = np.sum(selem != 0) - cdef int* sr = malloc(selem_num * sizeof(int)) - cdef int* sc = malloc(selem_num * sizeof(int)) + cdef Py_ssize_t selem_num = np.sum(selem != 0) + cdef Py_ssize_t* sr = malloc(selem_num * sizeof(Py_ssize_t)) + cdef Py_ssize_t* sc = malloc(selem_num * sizeof(Py_ssize_t)) s = 0 for r in range(srows): diff --git a/skimage/morphology/heap_general.pxi b/skimage/morphology/heap_general.pxi index a113b98e..67b21fa6 100644 --- a/skimage/morphology/heap_general.pxi +++ b/skimage/morphology/heap_general.pxi @@ -10,21 +10,18 @@ All rights reserved. Original author: Lee Kamentsky """ -cdef extern from "stdlib.h": - ctypedef unsigned long size_t - void free(void *ptr) - void *malloc(size_t size) - void *realloc(void *ptr, size_t size) +from libc.stdlib cimport free, malloc, realloc + cdef struct Heap: - unsigned int items - unsigned int space + Py_ssize_t items + Py_ssize_t space Heapitem *data Heapitem **ptrs cdef inline Heap *heap_from_numpy2(): - cdef unsigned int k - cdef Heap *heap + cdef Py_ssize_t k + cdef Heap *heap heap = malloc(sizeof (Heap)) heap.items = 0 heap.space = 1000 @@ -39,7 +36,7 @@ cdef inline void heap_done(Heap *heap): free(heap.ptrs) free(heap) -cdef inline void swap(unsigned int a, unsigned int b, Heap *h): +cdef inline void swap(Py_ssize_t a, Py_ssize_t b, Heap *h): h.ptrs[a], h.ptrs[b] = h.ptrs[b], h.ptrs[a] @@ -47,13 +44,13 @@ cdef inline void swap(unsigned int a, unsigned int b, Heap *h): # heappop - inlined # # pop an element off the heap, maintaining heap invariant -# +# # Note: heap ordering is the same as python heapq, i.e., smallest first. ###################################################### -cdef inline void heappop(Heap *heap, - Heapitem *dest): - cdef unsigned int i, smallest, l, r # heap indices - +cdef inline void heappop(Heap *heap, Heapitem *dest): + + cdef Py_ssize_t i, smallest, l, r # heap indices + # # Start by copying the first element to the destination # @@ -76,10 +73,10 @@ cdef inline void heappop(Heap *heap, smallest = i while True: # loop invariant here: smallest == i - + # find smallest of (i, l, r), and swap it to i's position if necessary - l = i*2+1 #__left(i) - r = i*2+2 #__right(i) + l = i * 2 + 1 #__left(i) + r = i * 2 + 2 #__right(i) if l < heap.items: if smaller(heap.ptrs[l], heap.ptrs[i]): smallest = l @@ -88,13 +85,14 @@ cdef inline void heappop(Heap *heap, else: # this is unnecessary, but trims 0.04 out of 0.85 seconds... break - # the element at i is smaller than either of its children, heap invariant restored. + # the element at i is smaller than either of its children, heap + # invariant restored. if smallest == i: break # swap swap(i, smallest, heap) i = smallest - + ################################################## # heappush - inlined # @@ -102,34 +100,36 @@ cdef inline void heappop(Heap *heap, # # Note: heap ordering is the same as python heapq, i.e., smallest first. ################################################## -cdef inline void heappush(Heap *heap, - Heapitem *new_elem): - cdef unsigned int child = heap.items - cdef unsigned int parent - cdef unsigned int k - cdef Heapitem *new_data +cdef inline void heappush(Heap *heap, Heapitem *new_elem): - # grow if necessary - if heap.items == heap.space: + cdef Py_ssize_t child = heap.items + cdef Py_ssize_t parent + cdef Py_ssize_t k + cdef Heapitem *new_data + + # grow if necessary + if heap.items == heap.space: heap.space = heap.space * 2 - new_data = realloc( heap.data, (heap.space * sizeof(Heapitem))) - heap.ptrs = realloc( heap.ptrs, (heap.space * sizeof(Heapitem *))) + new_data = realloc(heap.data, + (heap.space * sizeof(Heapitem))) + heap.ptrs = realloc(heap.ptrs, + (heap.space * sizeof(Heapitem *))) for k in range(heap.items): heap.ptrs[k] = new_data + (heap.ptrs[k] - heap.data) for k in range(heap.items, heap.space): heap.ptrs[k] = new_data + k heap.data = new_data - # insert new data at child - heap.ptrs[child][0] = new_elem[0] - heap.items += 1 + # insert new data at child + heap.ptrs[child][0] = new_elem[0] + heap.items += 1 - # restore heap invariant, all parents <= children - while child>0: - parent = (child + 1) / 2 - 1 # __parent(i) - - if smaller(heap.ptrs[child], heap.ptrs[parent]): - swap(parent, child, heap) - child = parent - else: - break + # restore heap invariant, all parents <= children + while child > 0: + parent = (child + 1) / 2 - 1 # __parent(i) + + if smaller(heap.ptrs[child], heap.ptrs[parent]): + swap(parent, child, heap) + child = parent + else: + break diff --git a/skimage/morphology/heap_watershed.pxi b/skimage/morphology/heap_watershed.pxi index ea66da26..573d32bb 100644 --- a/skimage/morphology/heap_watershed.pxi +++ b/skimage/morphology/heap_watershed.pxi @@ -13,14 +13,17 @@ import numpy as np cimport numpy as np cimport cython + cdef struct Heapitem: np.int32_t value np.int32_t age - np.int32_t index + Py_ssize_t index + cdef inline int smaller(Heapitem *a, Heapitem *b): if a.value <> b.value: - return a.value < b.value + return a.value < b.value return a.age < b.age + include "heap_general.pxi" diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py index b6f2a4fb..fbd63281 100644 --- a/skimage/morphology/watershed.py +++ b/skimage/morphology/watershed.py @@ -214,9 +214,7 @@ def watershed(image, markers, connectivity=None, offset=None, mask=None): c_mask = c_mask.astype(np.int8).flatten() _watershed.watershed(c_image.flatten(), pq, age, c, - c_image.ndim, c_mask, - np.array(c_image.shape, np.int32), c_output) c_output = c_output.reshape(c_image.shape)[[slice(1, -1, None)] * image.ndim] diff --git a/skimage/segmentation/_felzenszwalb_cy.pyx b/skimage/segmentation/_felzenszwalb_cy.pyx index efae45a4..5a8ac012 100644 --- a/skimage/segmentation/_felzenszwalb_cy.pyx +++ b/skimage/segmentation/_felzenszwalb_cy.pyx @@ -10,7 +10,7 @@ from ..util import img_as_float @cython.boundscheck(False) @cython.wraparound(False) @cython.cdivision(True) -def _felzenszwalb_grey(image, double scale=1, sigma=0.8, int min_size=20): +def _felzenszwalb_grey(image, double scale=1, sigma=0.8, Py_ssize_t min_size=20): """Felzenszwalb's efficient graph based segmentation for a single channel. Produces an oversegmentation of a 2d image using a fast, minimum spanning @@ -54,23 +54,23 @@ def _felzenszwalb_grey(image, double scale=1, sigma=0.8, int min_size=20): uright_cost.ravel()]).astype(np.float) # compute edges between pixels: height, width = image.shape[:2] - cdef np.ndarray[np.int_t, ndim=2] segments \ - = np.arange(width * height, dtype=np.int).reshape(height, width) + cdef np.ndarray[np.intp_t, ndim=2] segments \ + = np.arange(width * height, dtype=np.intp).reshape(height, width) right_edges = np.c_[segments[1:, :].ravel(), segments[:-1, :].ravel()] down_edges = np.c_[segments[:, 1:].ravel(), segments[:, :-1].ravel()] dright_edges = np.c_[segments[1:, 1:].ravel(), segments[:-1, :-1].ravel()] uright_edges = np.c_[segments[:-1, 1:].ravel(), segments[1:, :-1].ravel()] - cdef np.ndarray[np.int_t, ndim=2] edges \ + cdef np.ndarray[np.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) edges = np.ascontiguousarray(edges[edge_queue]) costs = np.ascontiguousarray(costs[edge_queue]) - cdef np.int_t *segments_p = segments.data - cdef np.int_t *edges_p = edges.data + cdef np.intp_t *segments_p = segments.data + cdef np.intp_t *edges_p = edges.data cdef np.float_t *costs_p = costs.data - cdef np.ndarray[np.int_t, ndim=1] segment_size \ + cdef np.ndarray[np.intp_t, ndim=1] segment_size \ = np.ones(width * height, dtype=np.int) # inner cost of segments cdef np.ndarray[np.float_t, ndim=1] cint = np.zeros(width * height) @@ -96,7 +96,7 @@ def _felzenszwalb_grey(image, double scale=1, sigma=0.8, int min_size=20): cint[seg_new] = costs_p[0] # postprocessing to remove small segments - edges_p = edges.data + 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]) diff --git a/skimage/segmentation/_quickshift.pyx b/skimage/segmentation/_quickshift.pyx index b465eb08..61e4d0e0 100644 --- a/skimage/segmentation/_quickshift.pyx +++ b/skimage/segmentation/_quickshift.pyx @@ -85,18 +85,19 @@ def quickshift(image, ratio=1., float kernel_size=5, max_dist=10, raise ValueError("Sigma should be >= 1") cdef int w = int(3 * kernel_size) - cdef int height = image_c.shape[0] - cdef int width = image_c.shape[1] - cdef int channels = image_c.shape[2] + cdef Py_ssize_t height = image_c.shape[0] + cdef Py_ssize_t width = image_c.shape[1] + cdef Py_ssize_t channels = image_c.shape[2] cdef double current_density, closest, dist - cdef int r, c, r_, c_, channel + cdef Py_ssize_t r, c, r_, c_, channel, r_min, c_min cdef np.float_t* image_p = image_c.data cdef np.float_t* current_pixel_p = image_p cdef np.ndarray[dtype=np.float_t, ndim=2] densities \ = np.zeros((height, width)) + # compute densities for r in range(height): for c in range(width): @@ -120,6 +121,7 @@ def quickshift(image, ratio=1., float kernel_size=5, max_dist=10, = np.arange(width * height).reshape(height, width) cdef np.ndarray[dtype=np.float_t, ndim=2] dist_parent \ = np.zeros((height, width)) + # find nearest node with higher density current_pixel_p = image_p for r in range(height): diff --git a/skimage/segmentation/_slic.pyx b/skimage/segmentation/_slic.pyx index dfb37a7c..e91ed21d 100644 --- a/skimage/segmentation/_slic.pyx +++ b/skimage/segmentation/_slic.pyx @@ -1,3 +1,4 @@ +#cython: boundscheck=False import numpy as np cimport numpy as np from time import time @@ -7,7 +8,7 @@ from ..color import rgb2lab, gray2rgb def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1, - convert2lab=True): + convert2lab=True): """Segments image using k-means clustering in Color-(x,y) space. Parameters @@ -62,10 +63,10 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1, image = rgb2lab(image) # initialize on grid: - cdef int height, width + cdef Py_ssize_t height, width height, width = image.shape[:2] # approximate grid size for desired n_segments - cdef int step = np.ceil(np.sqrt(height * width / n_segments)) + cdef Py_ssize_t step = int(np.ceil(np.sqrt(height * width / n_segments))) grid_y, grid_x = np.mgrid[:height, :width] means_y = grid_y[::step, ::step] means_x = grid_x[::step, ::step] @@ -81,11 +82,11 @@ def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1, ratio = (ratio / float(step)) ** 2 cdef np.ndarray[dtype=np.float_t, ndim=3] image_yx \ = np.dstack([grid_y, grid_x, image / ratio]).copy("C") - cdef int i, k, x, y, x_min, x_max, y_min, y_max, changes + cdef Py_ssize_t i, k, x, y, x_min, x_max, y_min, y_max, changes cdef double dist_mean - cdef np.ndarray[dtype=np.int_t, ndim=2] nearest_mean \ - = np.zeros((height, width), dtype=np.int) + cdef np.ndarray[dtype=np.intp_t, ndim=2] nearest_mean \ + = np.zeros((height, width), dtype=np.intp) cdef np.ndarray[dtype=np.float_t, ndim=2] distance \ = np.empty((height, width)) cdef np.float_t* image_p = image_yx.data diff --git a/skimage/segmentation/tests/test_slic.py b/skimage/segmentation/tests/test_slic.py index 59088fcd..89dee59b 100644 --- a/skimage/segmentation/tests/test_slic.py +++ b/skimage/segmentation/tests/test_slic.py @@ -13,8 +13,8 @@ def test_color(): img[img > 1] = 1 img[img < 0] = 0 seg = slic(img, sigma=0, n_segments=4) - # we expect 4 segments: - print(seg) + + # we expect 4 segments assert_equal(len(np.unique(seg)), 4) assert_array_equal(seg[:10, :10], 0) assert_array_equal(seg[10:, :10], 2) @@ -31,7 +31,7 @@ def test_gray(): img[img > 1] = 1 img[img < 0] = 0 seg = slic(img, sigma=0, n_segments=4, ratio=50.0) - print(seg) + assert_equal(len(np.unique(seg)), 4) assert_array_equal(seg[:10, :10], 0) assert_array_equal(seg[10:, :10], 2) diff --git a/skimage/transform/_hough_transform.pyx b/skimage/transform/_hough_transform.pyx index ef1cf700..662b3b88 100644 --- a/skimage/transform/_hough_transform.pyx +++ b/skimage/transform/_hough_transform.pyx @@ -1,23 +1,22 @@ +#cython: cdivision=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False cimport cython import numpy as np cimport numpy as np -from random import randint -from libc.math cimport abs, fabs, sqrt, ceil, floor +from libc.math cimport abs, fabs, sqrt, ceil from libc.stdlib cimport rand -np.import_array() - - cdef double PI_2 = 1.5707963267948966 cdef double NEG_PI_2 = -PI_2 -cdef inline int round(double r): - return ((r + 0.5) if (r > 0.0) else (r - 0.5)) +cdef inline Py_ssize_t round(double r): + return ((r + 0.5) if (r > 0.0) else (r - 0.5)) -@cython.boundscheck(False) def _hough(np.ndarray img, np.ndarray[ndim=1, dtype=np.double_t] theta=None): if img.ndim != 2: @@ -36,21 +35,20 @@ def _hough(np.ndarray img, np.ndarray[ndim=1, dtype=np.double_t] theta=None): # compute the bins and allocate the accumulator array cdef np.ndarray[ndim=2, dtype=np.uint64_t] accum cdef np.ndarray[ndim=1, dtype=np.double_t] bins - cdef int max_distance, offset + cdef Py_ssize_t max_distance, offset - max_distance = 2 * ceil((sqrt(img.shape[0] * img.shape[0] + - img.shape[1] * img.shape[1]))) + max_distance = 2 * ceil(sqrt(img.shape[0] * img.shape[0] + + img.shape[1] * img.shape[1])) accum = np.zeros((max_distance, theta.shape[0]), dtype=np.uint64) bins = np.linspace(-max_distance / 2.0, max_distance / 2.0, max_distance) offset = max_distance / 2 # compute the nonzero indexes cdef np.ndarray[ndim=1, dtype=np.npy_intp] x_idxs, y_idxs - y_idxs, x_idxs = np.PyArray_Nonzero(img) - + y_idxs, x_idxs = np.nonzero(img) # finally, run the transform - cdef int nidxs, nthetas, i, j, x, y, accum_idx + 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): @@ -61,67 +59,75 @@ def _hough(np.ndarray img, np.ndarray[ndim=1, dtype=np.double_t] theta=None): accum[accum_idx, j] += 1 return accum, theta, bins -import math -@cython.cdivision(True) -@cython.boundscheck(False) -def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \ - int line_gap, np.ndarray[ndim=1, dtype=np.double_t] theta=None): +def _probabilistic_hough(np.ndarray img, int value_threshold, + int line_length, int line_gap, + np.ndarray[ndim=1, dtype=np.double_t] theta=None): + if img.ndim != 2: raise ValueError('The input image must be 2D.') - # compute the array of angles and their sine and cosine - cdef np.ndarray[ndim=1, dtype=np.double_t] ctheta - cdef np.ndarray[ndim=1, dtype=np.double_t] stheta - # calculate thetas if none specified + if theta is None: - theta = np.linspace(math.pi/2, -math.pi/2, 180) - theta = math.pi/2-np.arange(180)/180.0* math.pi - ctheta = np.cos(theta) - stheta = np.sin(theta) - cdef int height = img.shape[0] - cdef int width = img.shape[1] + theta = PI_2 - np.arange(180) / 180.0 * 2 * PI_2 + + cdef Py_ssize_t height = img.shape[0] + cdef Py_ssize_t width = img.shape[1] + # compute the bins and allocate the accumulator array cdef np.ndarray[ndim=2, dtype=np.int64_t] accum - cdef np.ndarray[ndim=2, dtype=np.uint8_t] mask = np.zeros((height, width), dtype=np.uint8) - cdef np.ndarray[ndim=2, dtype=np.int32_t] line_end = np.zeros((2, 2), dtype=np.int32) - cdef int max_distance, offset, num_indexes, index + cdef np.ndarray[ndim=1, dtype=np.double_t] ctheta, stheta + cdef np.ndarray[ndim=2, dtype=np.uint8_t] mask = \ + np.zeros((height, width), dtype=np.uint8) + cdef np.ndarray[ndim=2, dtype=np.int32_t] line_end = \ + np.zeros((2, 2), dtype=np.int32) + cdef Py_ssize_t max_distance, offset, num_indexes, index cdef double a, b - cdef int nidxs, nthetas, i, j, x, y, px, py, accum_idx, value, max_value, max_theta + cdef Py_ssize_t nidxs, i, j, x, y, px, py, accum_idx + cdef int value, max_value, max_theta cdef int shift = 16 # maximum line number cutoff - cdef int lines_max = 2 ** 15 - cdef int xflag, x0, y0, dx0, dy0, dx, dy, gap, x1, y1, good_line, count + cdef Py_ssize_t lines_max = 2 ** 15 + cdef Py_ssize_t xflag, x0, y0, dx0, dy0, dx, dy, gap, x1, y1, \ + good_line, count + cdef list lines = list() + max_distance = 2 * ceil((sqrt(img.shape[0] * img.shape[0] + img.shape[1] * img.shape[1]))) accum = np.zeros((max_distance, theta.shape[0]), dtype=np.int64) offset = max_distance / 2 - # find the nonzero indexes - cdef np.ndarray[ndim=1, dtype=np.npy_intp] x_idxs, y_idxs - y_idxs, x_idxs = np.nonzero(img) - num_indexes = y_idxs.shape[0] # x and y are the same shape nthetas = theta.shape[0] - points = [] - for i in range(num_indexes): - points.append((x_idxs[i], y_idxs[i])) - lines = [] - # create mask of all non-zero indexes - for i in range(num_indexes): - mask[y_idxs[i], x_idxs[i]] = 1 + + # compute sine and cosine of angles + ctheta = np.cos(theta) + stheta = np.sin(theta) + + # find the nonzero indexes + y_idxs, x_idxs = np.nonzero(img) + points = list(zip(x_idxs, y_idxs)) + # mask all non-zero indexes + mask[y_idxs, x_idxs] = 1 + while 1: - # select random non-zero point + + # quit if no remaining points count = len(points) if count == 0: break - index = rand() % (count) + + # select random non-zero point + index = rand() % count x = points[index][0] y = points[index][1] del points[index] + # if previously eliminated, skip if not mask[y, x]: continue + value = 0 - max_value = value_threshold-1 + max_value = value_threshold - 1 max_theta = -1 + # apply hough transform on point for j in range(nthetas): accum_idx = round((ctheta[j] * x + stheta[j] * y)) + offset @@ -132,7 +138,9 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \ max_theta = j if max_value < value_threshold: continue - # from the random point walk in opposite directions and find line beginning and end + + # from the random point walk in opposite directions and find line + # beginning and end a = -stheta[max_theta] b = ctheta[max_theta] x0 = x @@ -188,6 +196,7 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \ # confirm line length is sufficient good_line = abs(line_end[1, 1] - line_end[0, 1]) >= line_length or \ abs(line_end[1, 0] - line_end[0, 0]) >= line_length + # pass 2: walk the line again and reset accumulator and mask for k in range(2): px = x0 @@ -207,7 +216,8 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \ # if non-zero point found, continue the line if mask[y1, x1]: if good_line: - accum_idx = round((ctheta[j] * x1 + stheta[j] * y1)) + offset + accum_idx = round((ctheta[j] * x1 \ + + stheta[j] * y1)) + offset accum[accum_idx, max_theta] -= 1 mask[y1, x1] = 0 # exit when the point is the line end @@ -218,9 +228,9 @@ def _probabilistic_hough(np.ndarray img, int value_threshold, int line_length, \ # add line to the result if good_line: - lines.append(((line_end[0, 0], line_end[0, 1]), (line_end[1, 0], line_end[1, 1]))) + lines.append(((line_end[0, 0], line_end[0, 1]), + (line_end[1, 0], line_end[1, 1]))) if len(lines) > lines_max: return lines + return lines - - diff --git a/skimage/transform/_warps_cy.pyx b/skimage/transform/_warps_cy.pyx index ce400ed6..951eb861 100644 --- a/skimage/transform/_warps_cy.pyx +++ b/skimage/transform/_warps_cy.pyx @@ -93,7 +93,7 @@ def _warp_fast(np.ndarray image, np.ndarray H, output_shape=None, int order=1, "`constant`, `nearest`, `wrap` or `reflect`.") cdef char mode_c = ord(mode[0].upper()) - cdef int out_r, out_c + cdef Py_ssize_t out_r, out_c if output_shape is None: out_r = img.shape[0] out_c = img.shape[1] @@ -104,12 +104,12 @@ def _warp_fast(np.ndarray image, np.ndarray H, output_shape=None, int order=1, cdef np.ndarray[dtype=np.double_t, ndim=2] out = \ np.zeros((out_r, out_c), dtype=np.double) - cdef int tfr, tfc + cdef Py_ssize_t tfr, tfc cdef double r, c - cdef int rows = img.shape[0] - cdef int cols = img.shape[1] + cdef Py_ssize_t rows = img.shape[0] + cdef Py_ssize_t cols = img.shape[1] - cdef double (*interp_func)(double*, int, int, double, double, + cdef double (*interp_func)(double*, Py_ssize_t, Py_ssize_t, double, double, char, double) if order == 0: interp_func = nearest_neighbour_interpolation