From cd9f3bd92e013f58de5e860979c08222481646cf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sat, 17 Aug 2013 22:07:45 +0200 Subject: [PATCH 1/5] Use typed memoryviews in draw package --- skimage/draw/_draw.pyx | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/skimage/draw/_draw.pyx b/skimage/draw/_draw.pyx index e95a2366..17dee0f7 100644 --- a/skimage/draw/_draw.pyx +++ b/skimage/draw/_draw.pyx @@ -47,8 +47,6 @@ def line(Py_ssize_t y, Py_ssize_t x, Py_ssize_t y2, Py_ssize_t x2): """ - cdef cnp.ndarray[cnp.intp_t, ndim=1, mode="c"] rr, cc - cdef char steep = 0 cdef Py_ssize_t dx = abs(x2 - x) cdef Py_ssize_t dy = abs(y2 - y) @@ -69,8 +67,8 @@ def line(Py_ssize_t y, Py_ssize_t x, Py_ssize_t y2, Py_ssize_t x2): sx, sy = sy, sx d = (2 * dy) - dx - rr = np.zeros(int(dx) + 1, dtype=np.intp) - cc = np.zeros(int(dx) + 1, dtype=np.intp) + cdef Py_ssize_t[:] rr = np.zeros(int(dx) + 1, dtype=np.intp) + cdef Py_ssize_t[:] cc = np.zeros(int(dx) + 1, dtype=np.intp) for i in range(dx): if steep: @@ -149,8 +147,8 @@ def polygon(y, x, shape=None): # make contigous arrays for r, c coordinates cdef cnp.ndarray contiguous_rdata, contiguous_cdata - contiguous_rdata = np.ascontiguousarray(y, 'double') - contiguous_cdata = np.ascontiguousarray(x, 'double') + contiguous_rdata = np.ascontiguousarray(y, dtype=np.double) + contiguous_cdata = np.ascontiguousarray(x, dtype=np.double) cdef cnp.double_t* rptr = contiguous_rdata.data cdef cnp.double_t* cptr = contiguous_cdata.data From 3172f508d6f31247aada306b559b850658e2bb8b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sat, 17 Aug 2013 22:35:13 +0200 Subject: [PATCH 2/5] Use typed memoryviews in feature package --- skimage/feature/_template.pyx | 8 +++---- skimage/feature/_texture.pyx | 42 ++++++++++++++++------------------- skimage/feature/corner_cy.pyx | 21 +++++------------- skimage/feature/texture.py | 2 +- 4 files changed, 30 insertions(+), 43 deletions(-) diff --git a/skimage/feature/_template.pyx b/skimage/feature/_template.pyx index 03695959..855ece23 100644 --- a/skimage/feature/_template.pyx +++ b/skimage/feature/_template.pyx @@ -50,9 +50,9 @@ from skimage.transform import integral def match_template(cnp.ndarray[float, ndim=2, mode="c"] image, cnp.ndarray[float, ndim=2, mode="c"] template): - cdef cnp.ndarray[float, ndim=2, mode="c"] corr - cdef cnp.ndarray[float, ndim=2, mode="c"] image_sat - cdef cnp.ndarray[float, ndim=2, mode="c"] image_sqr_sat + cdef float[:, ::1] corr + cdef float[:, ::1] image_sat + cdef float[:, ::1] image_sqr_sat cdef float template_mean = np.mean(template) cdef float template_ssd cdef float inv_area @@ -94,4 +94,4 @@ def match_template(cnp.ndarray[float, ndim=2, mode="c"] image, den = sqrt((window_sqr_sum - window_mean_sqr) * template_ssd) corr[r, c] /= den - return corr + return np.asarray(corr) diff --git a/skimage/feature/_texture.pyx b/skimage/feature/_texture.pyx index f98ed4ca..0e17fd1d 100644 --- a/skimage/feature/_texture.pyx +++ b/skimage/feature/_texture.pyx @@ -8,15 +8,9 @@ from libc.math cimport sin, cos, abs from skimage._shared.interpolation cimport bilinear_interpolation -def _glcm_loop(cnp.ndarray[dtype=cnp.uint8_t, ndim=2, - negative_indices=False, mode='c'] image, - cnp.ndarray[dtype=cnp.float64_t, ndim=1, - negative_indices=False, mode='c'] distances, - cnp.ndarray[dtype=cnp.float64_t, ndim=1, - negative_indices=False, mode='c'] angles, - int levels, - cnp.ndarray[dtype=cnp.uint32_t, ndim=4, - negative_indices=False, mode='c'] out): +def _glcm_loop(cnp.uint8_t[:, ::1] image, double[:] distances, + double[:] angles, Py_ssize_t levels, + cnp.uint32_t[:, :, :, ::1] out): """Perform co-occurrence matrix accumulation. Parameters @@ -81,7 +75,7 @@ cdef inline int _bit_rotate_right(int value, int length): return (value >> 1) | ((value & 1) << (length - 1)) -def _local_binary_pattern(cnp.ndarray[double, ndim=2] image, +def _local_binary_pattern(double[:, ::1] image, int P, float R, char method='D'): """Gray scale and rotation invariant LBP (Local Binary Patterns). @@ -92,8 +86,8 @@ def _local_binary_pattern(cnp.ndarray[double, ndim=2] image, image : (N, M) double array Graylevel image. P : int - Number of circularly symmetric neighbour set points (quantization of the - angular space). + Number of circularly symmetric neighbour set points (quantization of + the angular space). R : float Radius of circle (spatial resolution of the operator). method : {'D', 'R', 'U', 'V'} @@ -111,19 +105,20 @@ def _local_binary_pattern(cnp.ndarray[double, ndim=2] image, """ # texture weights - cdef cnp.ndarray[int, ndim=1] weights = 2 ** np.arange(P, dtype=np.int32) + cdef int[:] weights = 2 ** np.arange(P, dtype=np.int32) # local position of texture elements - rp = - R * np.sin(2 * np.pi * np.arange(P, dtype=np.double) / P) - cp = R * np.cos(2 * np.pi * np.arange(P, dtype=np.double) / P) - cdef cnp.ndarray[double, ndim=2] coords = np.round(np.vstack([rp, cp]).T, 5) + rr = - R * np.sin(2 * np.pi * np.arange(P, dtype=np.double) / P) + cc = R * np.cos(2 * np.pi * np.arange(P, dtype=np.double) / P) + cdef double[:] rp = np.round(rr, 5) + cdef double[:] cp = np.round(cc, 5) - # pre allocate arrays for computation - cdef cnp.ndarray[double, ndim=1] texture = np.zeros(P, np.double) - cdef cnp.ndarray[char, ndim=1] signed_texture = np.zeros(P, np.int8) - cdef cnp.ndarray[int, ndim=1] rotation_chain = np.zeros(P, np.int32) + # pre-allocate arrays for computation + cdef double[:] texture = np.zeros(P, dtype=np.double) + cdef char[:] signed_texture = np.zeros(P, dtype=np.int8) + cdef int[:] rotation_chain = np.zeros(P, dtype=np.int32) output_shape = (image.shape[0], image.shape[1]) - cdef cnp.ndarray[double, ndim=2] output = np.zeros(output_shape, np.double) + cdef double[:, ::1] output = np.zeros(output_shape, dtype=np.double) cdef Py_ssize_t rows = image.shape[0] cdef Py_ssize_t cols = image.shape[1] @@ -133,8 +128,9 @@ def _local_binary_pattern(cnp.ndarray[double, ndim=2] image, for r in range(image.shape[0]): for c in range(image.shape[1]): for i in range(P): - texture[i] = bilinear_interpolation(image.data, - rows, cols, r + coords[i, 0], c + coords[i, 1], 'C', 0) + texture[i] = bilinear_interpolation(&image[0, 0], rows, cols, + r + rp[i], c + cp[i], + 'C', 0) # signed / thresholded texture for i in range(P): if texture[i] - image[r, c] >= 0: diff --git a/skimage/feature/corner_cy.pyx b/skimage/feature/corner_cy.pyx index d4bc5e5a..7d558e52 100644 --- a/skimage/feature/corner_cy.pyx +++ b/skimage/feature/corner_cy.pyx @@ -59,16 +59,8 @@ def corner_moravec(image, Py_ssize_t window_size=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 - - if image.ndim == 3: - cimage = rgb2grey(image) - cimage = np.ascontiguousarray(img_as_float(image)) - - out = np.zeros(image.shape, dtype=np.double) - - cdef double* image_data = cimage.data - cdef double* out_data = out.data + cdef double[:, ::1] cimage = np.ascontiguousarray(img_as_float(image)) + cdef double[:, ::1] out = np.zeros(image.shape, dtype=np.double) cdef double msum, min_msum cdef Py_ssize_t r, c, br, bc, mr, mc, a, b @@ -81,11 +73,10 @@ def corner_moravec(image, Py_ssize_t window_size=1): msum = 0 for mr in range(- window_size, window_size + 1): for mc in range(- window_size, window_size + 1): - a = (r + mr) * cols + c + mc - b = (br + mr) * cols + bc + mc - msum += (image_data[a] - image_data[b]) ** 2 + msum += (cimage[r + mr, c + mc] + - cimage[br + mr, bc + mc]) ** 2 min_msum = min(msum, min_msum) - out_data[r * cols + c] = min_msum + out[r, c] = min_msum - return out + return np.asarray(out) diff --git a/skimage/feature/texture.py b/skimage/feature/texture.py index 7655b82a..7549cfdd 100644 --- a/skimage/feature/texture.py +++ b/skimage/feature/texture.py @@ -271,6 +271,6 @@ def local_binary_pattern(image, P, R, method='default'): 'uniform': ord('U'), 'var': ord('V') } - image = np.array(image, dtype='double', copy=True) + image = np.ascontiguousarray(image, dtype=np.double) output = _local_binary_pattern(image, P, R, methods[method.lower()]) return output From f76ba9dff3d3645f380b894ddeb4d73b17f69401 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sat, 17 Aug 2013 22:47:03 +0200 Subject: [PATCH 3/5] Use typed memoryviews in filter package --- skimage/filter/_ctmf.pyx | 21 ++++-------- skimage/filter/_denoise_cy.pyx | 60 ++++++++++++++-------------------- 2 files changed, 32 insertions(+), 49 deletions(-) diff --git a/skimage/filter/_ctmf.pyx b/skimage/filter/_ctmf.pyx index e3e845e4..70a59aa3 100644 --- a/skimage/filter/_ctmf.pyx +++ b/skimage/filter/_ctmf.pyx @@ -731,13 +731,8 @@ cdef int c_median_filter(Py_ssize_t rows, return 0 -def median_filter(cnp.ndarray[dtype=cnp.uint8_t, ndim=2, - negative_indices=False, mode='c'] data, - cnp.ndarray[dtype=cnp.uint8_t, ndim=2, - negative_indices=False, mode='c'] mask, - cnp.ndarray[dtype=cnp.uint8_t, ndim=2, - negative_indices=False, mode='c'] output, - int radius, +def median_filter(cnp.uint8_t[:, ::1] data, cnp.uint8_t[:, ::1] mask, + cnp.uint8_t[:, ::1] output, int radius, cnp.int32_t percent): """Median filter with octagon shape and masking. @@ -773,12 +768,10 @@ def median_filter(cnp.ndarray[dtype=cnp.uint8_t, ndim=2, raise ValueError('Data shape (%d, %d) is not output shape (%d, %d)' % (data.shape[0], data.shape[1], output.shape[0], output.shape[1])) - if c_median_filter(data.shape[0], - data.shape[1], - data.strides[0], - data.strides[1], + if c_median_filter(data.shape[0], data.shape[1], + data.strides[0], data.strides[1], radius, percent, - data.data, - mask.data, - output.data): + &data[0, 0], + &mask[0, 0], + &output[0, 0]): raise MemoryError('Failed to allocate scratchpad memory') diff --git a/skimage/filter/_denoise_cy.pyx b/skimage/filter/_denoise_cy.pyx index 3459171c..0c4f2539 100644 --- a/skimage/filter/_denoise_cy.pyx +++ b/skimage/filter/_denoise_cy.pyx @@ -113,11 +113,8 @@ def denoise_bilateral(image, Py_ssize_t win_size=5, sigma_range=None, double max_value - cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] cimage - cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] out - - double* image_data - double* out_data + double[:, :, ::1] cimage + double[:, :, ::1] out double* color_lut double* range_lut @@ -143,8 +140,6 @@ def denoise_bilateral(image, Py_ssize_t win_size=5, sigma_range=None, cimage = np.ascontiguousarray(image) out = np.zeros((rows, cols, dims), dtype=np.double) - image_data = cimage.data - out_data = out.data color_lut = _compute_color_lut(bins, csigma_range, max_value) range_lut = _compute_range_lut(win_size, sigma_spatial) dist_scale = bins / dims / max_value @@ -159,11 +154,10 @@ def denoise_bilateral(image, Py_ssize_t win_size=5, sigma_range=None, for r in range(rows): for c in range(cols): - pixel_addr = r * cols * dims + c * dims total_weight = 0 for d in range(dims): total_values[d] = 0 - centres[d] = image_data[pixel_addr + d] + centres[d] = cimage[r, c, d] for wr in range(-window_ext, window_ext + 1): rr = wr + r kr = wr + window_ext @@ -175,7 +169,7 @@ def denoise_bilateral(image, Py_ssize_t win_size=5, sigma_range=None, # distance between centre stack and current position dist = 0 for d in range(dims): - value = get_pixel3d(image_data, rows, cols, dims, + value = get_pixel3d(&cimage[0, 0, 0], rows, cols, dims, rr, cc, d, cmode, cval) values[d] = value dist += (centres[d] - value)**2 @@ -189,7 +183,7 @@ def denoise_bilateral(image, Py_ssize_t win_size=5, sigma_range=None, total_values[d] += values[d] * weight total_weight += weight for d in range(dims): - out_data[pixel_addr + d] = total_values[d] / total_weight + out[r, c, d] = total_values[d] / total_weight free(color_lut) free(range_lut) @@ -197,11 +191,11 @@ def denoise_bilateral(image, Py_ssize_t win_size=5, sigma_range=None, free(centres) free(total_values) - return np.squeeze(out) + return np.squeeze(np.asarray(out)) def denoise_tv_bregman(image, double weight, int max_iter=100, double eps=1e-3, - isotropic=True): + char isotropic=True): """Perform total-variation denoising using split-Bregman optimization. Total-variation denoising (also know as total-variation regularization) @@ -258,21 +252,17 @@ def denoise_tv_bregman(image, double weight, int max_iter=100, double eps=1e-3, Py_ssize_t total = rows * cols * dims - shape_ext = (rows2, cols2, dims) + shape_ext = (rows2, cols2, dims) + u = np.zeros(shape_ext, dtype=np.double) - cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] cimage = \ - np.ascontiguousarray(image) - cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] u = \ - np.zeros(shape_ext, dtype=np.double) + cdef: + double[:, :, ::1] cimage = np.ascontiguousarray(image) + double[:, :, ::1] cu = u - cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] dx = \ - np.zeros(shape_ext, dtype=np.double) - cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] dy = \ - np.zeros(shape_ext, dtype=np.double) - cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] bx = \ - np.zeros(shape_ext, dtype=np.double) - cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] by = \ - np.zeros(shape_ext, dtype=np.double) + double[:, :, ::1] dx = np.zeros(shape_ext, dtype=np.double) + double[:, :, ::1] dy = np.zeros(shape_ext, dtype=np.double) + double[:, :, ::1] bx = np.zeros(shape_ext, dtype=np.double) + double[:, :, ::1] by = np.zeros(shape_ext, dtype=np.double) double ux, uy, uprev, unew, bxx, byy, dxx, dyy, s int i = 0 @@ -296,19 +286,19 @@ def denoise_tv_bregman(image, double weight, int max_iter=100, double eps=1e-3, for r in range(1, rows + 1): for c in range(1, cols + 1): - uprev = u[r, c, k] + uprev = cu[r, c, k] # forward derivatives - ux = u[r, c + 1, k] - uprev - uy = u[r + 1, c, k] - uprev + ux = cu[r, c + 1, k] - uprev + uy = cu[r + 1, c, k] - uprev # Gauss-Seidel method unew = ( lam * ( - + u[r + 1, c, k] - + u[r - 1, c, k] - + u[r, c + 1, k] - + u[r, c - 1, k] + + cu[r + 1, c, k] + + cu[r - 1, c, k] + + cu[r, c + 1, k] + + cu[r, c - 1, k] + dx[r, c - 1, k] - dx[r, c, k] @@ -321,7 +311,7 @@ def denoise_tv_bregman(image, double weight, int max_iter=100, double eps=1e-3, + by[r, c, k] ) + weight * cimage[r - 1, c - 1, k] ) / norm - u[r, c, k] = unew + cu[r, c, k] = unew # update root mean square error rmse += (unew - uprev)**2 @@ -360,4 +350,4 @@ def denoise_tv_bregman(image, double weight, int max_iter=100, double eps=1e-3, rmse = sqrt(rmse / total) i += 1 - return np.squeeze(u[1:-1, 1:-1]) + return np.squeeze(np.asarray(u[1:-1, 1:-1])) From 5b38bdac59a00d6a089cb934437caf4de09d7145 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sat, 17 Aug 2013 23:07:00 +0200 Subject: [PATCH 4/5] Fix missing conversion of typed memoryview to numpy array --- skimage/feature/_texture.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/feature/_texture.pyx b/skimage/feature/_texture.pyx index 0e17fd1d..9daa215d 100644 --- a/skimage/feature/_texture.pyx +++ b/skimage/feature/_texture.pyx @@ -177,4 +177,4 @@ def _local_binary_pattern(double[:, ::1] image, output[r, c] = lbp - return output + return np.asarray(output) From a8f0d46ab299dd7d30af32890e5d4c4508dd9b3a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sat, 17 Aug 2013 23:17:05 +0200 Subject: [PATCH 5/5] Fix missing conversion of typed memoryview to numpy array --- skimage/draw/_draw.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/draw/_draw.pyx b/skimage/draw/_draw.pyx index 17dee0f7..132ba1d4 100644 --- a/skimage/draw/_draw.pyx +++ b/skimage/draw/_draw.pyx @@ -86,7 +86,7 @@ def line(Py_ssize_t y, Py_ssize_t x, Py_ssize_t y2, Py_ssize_t x2): rr[dx] = y2 cc[dx] = x2 - return rr, cc + return np.asarray(rr), np.asarray(cc) def polygon(y, x, shape=None):