Merge pull request #697 from ahojnnes/memoryviews

Memoryview refactoring: Part 2.
This commit is contained in:
Stefan van der Walt
2013-08-20 23:56:25 -07:00
13 changed files with 95 additions and 110 deletions
+1 -1
View File
@@ -152,4 +152,4 @@ def _test_random(shape):
if __name__ == "__main__":
run_module_suite()
np.testing.run_module_suite()
+2 -2
View File
@@ -1,5 +1,5 @@
import numpy as np
from numpy.testing import *
from numpy.testing import assert_equal, assert_array_equal
import skimage.graph.spath as spath
@@ -33,4 +33,4 @@ def test_non_square():
if __name__ == "__main__":
run_module_suite()
np.testing.run_module_suite()
+2 -4
View File
@@ -4,8 +4,6 @@
#cython: wraparound=False
import numpy as np
cimport numpy as cnp
cdef inline double _get_fraction(double from_value, double to_value,
double level):
@@ -14,7 +12,7 @@ cdef inline double _get_fraction(double from_value, double to_value,
return ((level - from_value) / (to_value - from_value))
def iterate_and_store(cnp.ndarray[double, ndim=2] array,
def iterate_and_store(double[:, :] array,
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
@@ -46,7 +44,7 @@ def iterate_and_store(cnp.ndarray[double, ndim=2] array,
# Calculate the number of iterations we'll need
cdef Py_ssize_t num_square_steps = (array.shape[0] - 1) \
* (array.shape[1] - 1)
* (array.shape[1] - 1)
cdef unsigned char square_case = 0
cdef tuple top, bottom, left, right
+26 -14
View File
@@ -7,7 +7,7 @@ import numpy as np
cimport numpy as cnp
def possible_hull(cnp.ndarray[dtype=cnp.uint8_t, ndim=2, mode="c"] img):
def possible_hull(cnp.uint8_t[:, ::1] img):
"""Return positions of pixels that possibly belong to the convex hull.
Parameters
@@ -30,31 +30,43 @@ def possible_hull(cnp.ndarray[dtype=cnp.uint8_t, ndim=2, mode="c"] img):
# cols storage slots for top boundary pixels
# rows storage slots for right boundary pixels
# cols storage slots for bottom boundary pixels
cdef cnp.ndarray[dtype=cnp.intp_t, ndim=2] nonzero = \
np.ones((2 * (rows + cols), 2), dtype=np.intp)
nonzero *= -1
coords = np.ones((2 * (rows + cols), 2), dtype=np.intp)
coords *= -1
cdef Py_ssize_t[:, ::1] nonzero = coords
cdef Py_ssize_t rows_cols = rows + cols
cdef Py_ssize_t rows_2_cols = 2 * rows + cols
cdef Py_ssize_t rows_cols_r, rows_c
for r in range(rows):
rows_cols_r = rows_cols + r
for c in range(cols):
if img[r, c] != 0:
rows_c = rows + c
rows_2_cols_c = rows_2_cols + c
# Left check
if nonzero[r, 1] == -1:
nonzero[r, 0] = r
nonzero[r, 1] = c
# Right check
elif nonzero[rows + cols + r, 1] < c:
nonzero[rows + cols + r, 0] = r
nonzero[rows + cols + r, 1] = c
elif nonzero[rows_cols_r, 1] < c:
nonzero[rows_cols_r, 0] = r
nonzero[rows_cols_r, 1] = c
# Top check
if nonzero[rows + c, 1] == -1:
nonzero[rows + c, 0] = r
nonzero[rows + c, 1] = c
if nonzero[rows_c, 1] == -1:
nonzero[rows_c, 0] = r
nonzero[rows_c, 1] = c
# Bottom check
elif nonzero[2 * rows + cols + c, 0] < r:
nonzero[2 * rows + cols + c, 0] = r
nonzero[2 * rows + cols + c, 1] = c
elif nonzero[rows_2_cols_c, 0] < r:
nonzero[rows_2_cols_c, 0] = r
nonzero[rows_2_cols_c, 1] = c
return nonzero[nonzero[:, 0] != -1]
return coords[coords[:, 0] != -1]
+2 -2
View File
@@ -21,8 +21,8 @@ def reconstruction_loop(cnp.ndarray[dtype=cnp.uint32_t, ndim=1,
negative_indices=False, mode='c'] anext,
cnp.ndarray[dtype=cnp.int32_t, ndim=1,
negative_indices=False, mode='c'] astrides,
int current_idx,
int image_stride):
Py_ssize_t current_idx,
Py_ssize_t image_stride):
"""The inner loop for reconstruction.
This algorithm uses the rank-order of pixels. If low intensity pixels have
+9 -10
View File
@@ -29,7 +29,7 @@ def grid_points_inside_poly(shape, verts):
True where the grid falls inside the polygon.
"""
cdef cnp.ndarray[cnp.double_t, ndim=1, mode="c"] vx, vy
cdef double[:] vx, vy
verts = np.asarray(verts)
vx = verts[:, 0].astype(np.double)
@@ -45,8 +45,7 @@ def grid_points_inside_poly(shape, verts):
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, &vx[0], &vy[0], m, n)
return out.view(bool)
@@ -57,18 +56,18 @@ def points_inside_poly(points, verts):
Parameters
----------
points : (N, 2) array
Input points, ``(x, y)``.
Input points, ``(x, y)``.
verts : (M, 2) array
Vertices of the polygon, sorted either clockwise or anti-clockwise.
The first point may (but does not need to be) duplicated.
Vertices of the polygon, sorted either clockwise or anti-clockwise.
The first point may (but does not need to be) duplicated.
Returns
-------
mask : (N,) array of bool
True if corresponding point is inside the polygon.
True if corresponding point is inside the polygon.
"""
cdef cnp.ndarray[cnp.double_t, ndim=1, mode="c"] x, y, vx, vy
cdef double[:] x, y, vx, vy
points = np.asarray(points)
verts = np.asarray(verts)
@@ -82,8 +81,8 @@ def points_inside_poly(points, verts):
cdef cnp.ndarray[cnp.uint8_t, ndim=1] out = \
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,
points_in_polygon(vx.shape[0], &vx[0], &vy[0],
x.shape[0], &x[0], &y[0],
<unsigned char*>out.data)
return out.astype(bool)
+17 -24
View File
@@ -19,16 +19,9 @@ import numpy as np
cimport numpy as cnp
def _skeletonize_loop(cnp.ndarray[dtype=cnp.uint8_t, ndim=2,
negative_indices=False, mode='c'] result,
cnp.ndarray[dtype=cnp.intp_t, ndim=1,
negative_indices=False, mode='c'] i,
cnp.ndarray[dtype=cnp.intp_t, ndim=1,
negative_indices=False, mode='c'] j,
cnp.ndarray[dtype=cnp.int32_t, ndim=1,
negative_indices=False, mode='c'] order,
cnp.ndarray[dtype=cnp.uint8_t, ndim=1,
negative_indices=False, mode='c'] table):
def _skeletonize_loop(cnp.uint8_t[:, ::1] result,
Py_ssize_t[:] i, Py_ssize_t[:] j,
cnp.int32_t[:] order, cnp.uint8_t[:] table):
"""
Inner loop of skeletonize function
@@ -65,9 +58,11 @@ def _skeletonize_loop(cnp.ndarray[dtype=cnp.uint8_t, ndim=2,
pixels.
"""
cdef:
cnp.int32_t accumulator
Py_ssize_t accumulator
Py_ssize_t index, order_index
Py_ssize_t ii, jj
Py_ssize_t rows = result.shape[0]
Py_ssize_t cols = result.shape[1]
for index in range(order.shape[0]):
accumulator = 16
@@ -80,26 +75,25 @@ def _skeletonize_loop(cnp.ndarray[dtype=cnp.uint8_t, ndim=2,
accumulator += 1
if result[ii - 1, jj]:
accumulator += 2
if jj < result.shape[1] - 1 and result[ii - 1, jj + 1]:
if jj < cols - 1 and result[ii - 1, jj + 1]:
accumulator += 4
if jj > 0 and result[ii, jj - 1]:
accumulator += 8
if jj < result.shape[1] - 1 and result[ii, jj + 1]:
if jj < cols - 1 and result[ii, jj + 1]:
accumulator += 32
if ii < result.shape[0]-1:
if ii < rows - 1:
if jj > 0 and result[ii + 1, jj - 1]:
accumulator += 64
if result[ii + 1, jj]:
accumulator += 128
if jj < result.shape[1] - 1 and result[ii + 1, jj + 1]:
if jj < cols - 1 and result[ii + 1, jj + 1]:
accumulator += 256
# Assign the value of table corresponding to the configuration
result[ii, jj] = table[accumulator]
def _table_lookup_index(cnp.ndarray[dtype=cnp.uint8_t, ndim=2,
negative_indices=False, mode='c'] image):
def _table_lookup_index(cnp.uint8_t[:, ::1] image):
"""
Return an index into a table per pixel of a binary image
@@ -120,9 +114,8 @@ def _table_lookup_index(cnp.ndarray[dtype=cnp.uint8_t, ndim=2,
hardwired kernel.
"""
cdef:
cnp.ndarray[dtype=cnp.int32_t, ndim=2,
negative_indices=False, mode='c'] indexer
cnp.int32_t *p_indexer
Py_ssize_t[:, ::1] indexer
Py_ssize_t *p_indexer
cnp.uint8_t *p_image
Py_ssize_t i_stride
Py_ssize_t i_shape
@@ -133,9 +126,9 @@ def _table_lookup_index(cnp.ndarray[dtype=cnp.uint8_t, ndim=2,
i_shape = image.shape[0]
j_shape = image.shape[1]
indexer = np.zeros((i_shape, j_shape), np.int32)
p_indexer = <cnp.int32_t *>indexer.data
p_image = <cnp.uint8_t *>image.data
indexer = np.zeros((i_shape, j_shape), dtype=np.intp)
p_indexer = &indexer[0, 0]
p_image = &image[0, 0]
i_stride = image.strides[0]
assert i_shape >= 3 and j_shape >= 3, \
"Please use the slow method for arrays < 3x3"
@@ -214,4 +207,4 @@ def _table_lookup_index(cnp.ndarray[dtype=cnp.uint8_t, ndim=2,
indexer[i - 1, j_shape - 1] += 128
indexer[i, j_shape - 1] += 16
indexer[i + 1, j_shape - 1] += 2
return indexer
return np.asarray(indexer)
+5 -10
View File
@@ -23,17 +23,12 @@ 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,
def watershed(DTYPE_INT32_t[:] image,
DTYPE_INT32_t[:, ::1] 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):
DTYPE_INT32_t[:, ::1] structure,
DTYPE_BOOL_t[:] mask,
DTYPE_INT32_t[:] output):
"""Do heavy lifting of watershed algorithm
Parameters
+14 -20
View File
@@ -8,9 +8,9 @@ cimport numpy as np
from libc.stdlib cimport malloc, free
def _dilate(np.ndarray[np.uint8_t, ndim=2] image,
np.ndarray[np.uint8_t, ndim=2] selem,
np.ndarray[np.uint8_t, ndim=2] out=None,
def _dilate(np.uint8_t[:, :] image,
np.uint8_t[:, :] selem,
np.uint8_t[:, :] out=None,
char shift_x=0, char shift_y=0):
"""Return greyscale morphological dilation of an image.
@@ -52,12 +52,9 @@ def _dilate(np.ndarray[np.uint8_t, ndim=2] image,
else:
out = np.ascontiguousarray(out)
cdef np.uint8_t* out_data = <np.uint8_t*>out.data
cdef np.uint8_t* image_data = <np.uint8_t*>image.data
cdef Py_ssize_t r, c, rr, cc, s, value, local_max
cdef Py_ssize_t selem_num = np.sum(selem != 0)
cdef Py_ssize_t selem_num = np.sum(np.asarray(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))
@@ -76,21 +73,21 @@ def _dilate(np.ndarray[np.uint8_t, ndim=2] image,
rr = r + sr[s]
cc = c + sc[s]
if 0 <= rr < rows and 0 <= cc < cols:
value = image_data[rr * cols + cc]
value = image[rr, cc]
if value > local_max:
local_max = value
out_data[r * cols + c] = local_max
out[r, c] = local_max
free(sr)
free(sc)
return out
return np.asarray(out)
def _erode(np.ndarray[np.uint8_t, ndim=2] image,
np.ndarray[np.uint8_t, ndim=2] selem,
np.ndarray[np.uint8_t, ndim=2] out=None,
def _erode(np.uint8_t[:, :] image,
np.uint8_t[:, :] selem,
np.uint8_t[:, :] out=None,
char shift_x=0, char shift_y=0):
"""Return greyscale morphological erosion of an image.
@@ -131,12 +128,9 @@ def _erode(np.ndarray[np.uint8_t, ndim=2] image,
else:
out = np.ascontiguousarray(out)
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_min
cdef Py_ssize_t selem_num = np.sum(selem != 0)
cdef Py_ssize_t selem_num = np.sum(np.asarray(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))
@@ -155,13 +149,13 @@ def _erode(np.ndarray[np.uint8_t, ndim=2] image,
rr = r + sr[s]
cc = c + sc[s]
if 0 <= rr < rows and 0 <= cc < cols:
value = image_data[rr * cols + cc]
value = image[rr, cc]
if value < local_min:
local_min = value
out_data[r * cols + c] = local_min
out[r, c] = local_min
free(sr)
free(sc)
return out
return np.asarray(out)
+1 -1
View File
@@ -159,7 +159,7 @@ def reconstruction(seed, mask, method='dilation', selem=None, offset=None):
# Create a list of strides across the array to get the neighbors within
# a flattened array
value_stride = np.array(images.strides[1:]) / images.dtype.itemsize
image_stride = images.strides[0] / images.dtype.itemsize
image_stride = images.strides[0] // images.dtype.itemsize
selem_mgrid = np.mgrid[[slice(-o, d - o)
for d, o in zip(selem.shape, offset)]]
selem_offsets = selem_mgrid[:, selem].transpose()
+2 -1
View File
@@ -12,7 +12,8 @@ from skimage.morphology.ccomp cimport find_root, join_trees
from ..util import img_as_float
def _felzenszwalb_grey(image, double scale=1, sigma=0.8, Py_ssize_t 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
+8 -12
View File
@@ -2,15 +2,11 @@
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
import collections as coll
import numpy as np
from time import time
from scipy import ndimage
cimport numpy as cnp
from ..util import img_as_float, regular_grid
from ..color import rgb2lab, gray2rgb
from skimage.util import img_as_float, regular_grid
from skimage.color import rgb2lab, gray2rgb
def _slic_cython(double[:, :, :, ::1] image_zyx,
@@ -19,18 +15,18 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
double[:, ::1] means,
Py_ssize_t max_iter, Py_ssize_t n_segments):
"""Helper function for SLIC segmentation.
Parameters
----------
image_zyx : 4D np.ndarray of double, shape (Z, Y, X, 6)
image_zyx : 4D array of double, shape (Z, Y, X, 6)
The image with embedded coordinates, that is, `image_zyx[i, j, k]` is
`array([i, j, k, r, g, b])` or `array([i, j, k, L, a, b])`, depending
on the colorspace.
nearest_mean : 3D np.ndarray of int, shape (Z, Y, X)
nearest_mean : 3D array of int, shape (Z, Y, X)
The (initially empty) label field.
distance : 3D np.ndarray of double, shape (Z, Y, X)
distance : 3D array of double, shape (Z, Y, X)
The (initially infinity) array of distances to the nearest centroid.
means : 2D np.ndarray of double, shape (n_segments, 6)
means : 2D array of double, shape (n_segments, 6)
The centroids obtained by SLIC.
max_iter : int
The maximum number of k-means iterations.
@@ -39,7 +35,7 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
Returns
-------
nearest_mean : 3D np.ndarray of int, shape (Z, Y, X)
nearest_mean : 3D array of int, shape (Z, Y, X)
The label field/superpixels found by SLIC.
"""
+6 -9
View File
@@ -83,10 +83,8 @@ def _warp_fast(cnp.ndarray image, cnp.ndarray H, output_shape=None,
"""
cdef cnp.ndarray[dtype=cnp.double_t, ndim=2, mode="c"] img = \
np.ascontiguousarray(image, dtype=np.double)
cdef cnp.ndarray[dtype=cnp.double_t, ndim=2, mode="c"] M = \
np.ascontiguousarray(H)
cdef double[:, ::1] img = np.ascontiguousarray(image, dtype=np.double)
cdef double[:, ::1] M = np.ascontiguousarray(H)
if mode not in ('constant', 'wrap', 'reflect', 'nearest'):
raise ValueError("Invalid mode specified. Please use "
@@ -101,8 +99,7 @@ def _warp_fast(cnp.ndarray image, cnp.ndarray H, output_shape=None,
out_r = output_shape[0]
out_c = output_shape[1]
cdef cnp.ndarray[dtype=cnp.double_t, ndim=2] out = \
np.zeros((out_r, out_c), dtype=np.double)
cdef double[:, ::1] out = np.zeros((out_r, out_c), dtype=np.double)
cdef Py_ssize_t tfr, tfc
cdef double r, c
@@ -122,8 +119,8 @@ def _warp_fast(cnp.ndarray image, cnp.ndarray H, output_shape=None,
for tfr in range(out_r):
for tfc in range(out_c):
_matrix_transform(tfc, tfr, <double*>M.data, &c, &r)
out[tfr, tfc] = interp_func(<double*>img.data, rows, cols, r, c,
_matrix_transform(tfc, tfr, &M[0, 0], &c, &r)
out[tfr, tfc] = interp_func(&img[0, 0], rows, cols, r, c,
mode_c, cval)
return out
return np.asarray(out)