Merge pull request #424 from ahojnnes/ssize_t

RF: Globally change all index variables to ssize_t.
This commit is contained in:
Stefan van der Walt
2013-02-24 01:27:45 -08:00
35 changed files with 521 additions and 509 deletions
+3 -3
View File
@@ -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)
+7 -7
View File
@@ -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])
+10 -10
View File
@@ -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)
+39 -39
View File
@@ -5,12 +5,12 @@
from libc.math cimport ceil, floor
cdef inline int round(double r):
return <int>((r + 0.5) if (r > 0.0) else (r - 0.5))
cdef inline Py_ssize_t round(double r):
return <Py_ssize_t>((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, <int>round(r), <int>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 = <int>floor(r)
minc = <int>floor(c)
maxr = <int>ceil(r)
maxc = <int>ceil(c)
minr = <Py_ssize_t>floor(r)
minc = <Py_ssize_t>floor(c)
maxr = <Py_ssize_t>ceil(r)
maxc = <Py_ssize_t>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 = <int>round(r)
cdef int c0 = <int>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 = <int>r - 1
cdef int c0 = <int>c - 1
cdef Py_ssize_t r0 = <Py_ssize_t>r - 1
cdef Py_ssize_t c0 = <Py_ssize_t>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 <int>(-coord / dim) % 2 != 0:
return dim - <int>(-coord % dim)
if <Py_ssize_t>(-coord / dim) % 2 != 0:
return dim - <Py_ssize_t>(-coord % dim)
else:
return <int>(-coord % dim)
return <Py_ssize_t>(-coord % dim)
elif coord > dim:
if <int>(coord / dim) % 2 != 0:
return <int>(dim - (coord % dim))
if <Py_ssize_t>(coord / dim) % 2 != 0:
return <Py_ssize_t>(dim - (coord % dim))
else:
return <int>(coord % dim)
return <Py_ssize_t>(coord % dim)
elif mode == 'W': # wrap
if coord < 0:
return <int>(dim - (-coord % dim))
return <Py_ssize_t>(dim - (-coord % dim))
elif coord > dim:
return <int>(coord % dim)
return <Py_ssize_t>(coord % dim)
elif mode == 'N': # nearest
if coord < 0:
return 0
+1 -1
View File
@@ -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)
+1 -1
View File
@@ -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.
+39 -37
View File
@@ -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 = <int>max(0, y.min())
cdef int maxr = <int>math.ceil(y.max())
cdef int minc = <int>max(0, x.min())
cdef int maxc = <int>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 = <int>max(0, cy - yradius)
cdef int maxr = <int>math.ceil(cy + yradius)
cdef int minc = <int>max(0, cx - xradius)
cdef int maxc = <int>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])
+7 -6
View File
@@ -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:
+4 -2
View File
@@ -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]
"""
+18 -13
View File
@@ -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
+13 -11
View File
@@ -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 + <int>(sin(angle) * distance + 0.5)
col = c + <int>(cos(angle) * distance + 0.5);
col = c + <int>(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):
+4 -4
View File
@@ -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 = <double*>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
+1 -1
View File
@@ -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)
+97 -99
View File
@@ -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 = <size_t> ptr
roundoff = <Py_ssize_t>ptr
roundoff += 31
roundoff -= roundoff % 32
ptr = <void *> 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 = <int> (<np.float64_t> radius * 2.0 / 2.414213)
a = <Py_ssize_t>(<np.float64_t>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 <np.int32_t> (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 <np.int32_t> (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 <np.int32_t> (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 <np.int32_t> (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 = <np.uint32_t> ((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 <np.uint8_t> j
return <np.uint8_t>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)
+5 -5
View File
@@ -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 = <double*>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 = <double*>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
+6 -7
View File
@@ -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)
+5 -8
View File
@@ -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.
+2 -2
View File
@@ -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):
+27 -30
View File
@@ -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]
+7 -6
View File
@@ -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, <double*>vx.data, <double*>vy.data, m, n)
out[m, n] = point_in_polygon(V, <double*>vx.data, <double*>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], <double*>vx.data, <double*>vy.data,
x.shape[0], <double*>x.data, <double*>y.data,
+2 -2
View File
@@ -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.
+18 -18
View File
@@ -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]
+20 -29
View File
@@ -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 *> heap_from_numpy2()
+6 -6
View File
@@ -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)
+20 -17
View File
@@ -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 = <np.int_t*>forest.data
cdef np.int_t *data_p = <np.int_t*>data.data
cdef DTYPE_t *forest_p = <DTYPE_t*>forest.data
cdef DTYPE_t *data_p = <DTYPE_t*>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:
+19 -19
View File
@@ -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 = <np.uint8_t*>out.data
cdef np.uint8_t* image_data = <np.uint8_t*>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 = <int*>malloc(selem_num * sizeof(int))
cdef int* sc = <int*>malloc(selem_num * sizeof(int))
cdef Py_ssize_t selem_num = np.sum(selem != 0)
cdef Py_ssize_t* sr = <Py_ssize_t*>malloc(selem_num * sizeof(Py_ssize_t))
cdef Py_ssize_t* sc = <Py_ssize_t*>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 = <int*>malloc(selem_num * sizeof(int))
cdef int* sc = <int*>malloc(selem_num * sizeof(int))
cdef Py_ssize_t selem_num = np.sum(selem != 0)
cdef Py_ssize_t* sr = <Py_ssize_t*>malloc(selem_num * sizeof(Py_ssize_t))
cdef Py_ssize_t* sc = <Py_ssize_t*>malloc(selem_num * sizeof(Py_ssize_t))
s = 0
for r in range(srows):
+42 -42
View File
@@ -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 = <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 = <Heapitem *> realloc(<void *> heap.data, <size_t> (heap.space * sizeof(Heapitem)))
heap.ptrs = <Heapitem **> realloc(<void *> heap.ptrs, <size_t> (heap.space * sizeof(Heapitem *)))
new_data = <Heapitem*>realloc(<void*>heap.data,
<Py_ssize_t>(heap.space * sizeof(Heapitem)))
heap.ptrs = <Heapitem**>realloc(<void*>heap.ptrs,
<Py_ssize_t>(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
+5 -2
View File
@@ -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"
-2
View File
@@ -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]
+8 -8
View File
@@ -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 = <np.int_t*>segments.data
cdef np.int_t *edges_p = <np.int_t*>edges.data
cdef np.intp_t *segments_p = <np.intp_t*>segments.data
cdef np.intp_t *edges_p = <np.intp_t*>edges.data
cdef np.float_t *costs_p = <np.float_t*>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 = <np.int_t*>edges.data
edges_p = <np.intp_t*>edges.data
for e in range(costs.size):
seg0 = find_root(segments_p, edges_p[0])
seg1 = find_root(segments_p, edges_p[1])
+6 -4
View File
@@ -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 = <np.float_t*> 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):
+7 -6
View File
@@ -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 = <np.float_t*> image_yx.data
+3 -3
View File
@@ -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)
+64 -54
View File
@@ -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 <int>((r + 0.5) if (r > 0.0) else (r - 0.5))
cdef inline Py_ssize_t round(double r):
return <Py_ssize_t>((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 * <int>ceil((sqrt(img.shape[0] * img.shape[0] +
img.shape[1] * img.shape[1])))
max_distance = 2 * <Py_ssize_t>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 * <int>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 = <int>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 = <int>round((ctheta[j] * x1 + stheta[j] * y1)) + offset
accum_idx = <int>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
+5 -5
View File
@@ -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