[Numpy-discussion] Example Usage of Neighborhood Iterator in Cython
T J
tjhnson@gmail....
Mon Oct 17 13:19:10 CDT 2011
I recently put together a Cython example which uses the neighborhood
iterator. It was trickier than I thought it would be, so I thought to
share it with the community. The function takes a 1-dimensional array
and returns a 2-dimensional array of neighborhoods in the original
area. This is somewhat similar to the functionality provided by
segment_axis (http://projects.scipy.org/numpy/ticket/901), but I
believe this slightly different in that neighborhood can extend to the
left of the current index (assuming circular boundaries). Keep in
mind that this is just an example, and normal usage probably is not
concerned with creating a new array.
External link: http://codepad.org/czRIzXQl
--------------
import numpy as np
cimport numpy as np
cdef extern from "numpy/arrayobject.h":
ctypedef extern class numpy.flatiter [object PyArrayIterObject]:
cdef int nd_m1
cdef np.npy_intp index, size
cdef np.ndarray ao
cdef char *dataptr
# This isn't exposed to the Python API.
# So we can't use the same approach we used to define flatiter
ctypedef struct PyArrayNeighborhoodIterObject:
int nd_m1
np.npy_intp index, size
np.PyArrayObject *ao # note the change from np.ndarray
char *dataptr
object PyArray_NeighborhoodIterNew(flatiter it, np.npy_intp* bounds,
int mode, np.ndarray fill_value)
int PyArrayNeighborhoodIter_Next(PyArrayNeighborhoodIterObject *it)
int PyArrayNeighborhoodIter_Reset(PyArrayNeighborhoodIterObject *it)
object PyArray_IterNew(object arr)
void PyArray_ITER_NEXT(flatiter it)
np.npy_intp PyArray_SIZE(np.ndarray arr)
cdef enum:
NPY_NEIGHBORHOOD_ITER_ZERO_PADDING,
NPY_NEIGHBORHOOD_ITER_ONE_PADDING,
NPY_NEIGHBORHOOD_ITER_CONSTANT_PADDING,
NPY_NEIGHBORHOOD_ITER_CIRCULAR_PADDING,
NPY_NEIGHBORHOOD_ITER_MIRROR_PADDING
np.import_array()
def windowed(np.ndarray[np.int_t, ndim=1] arr, bounds):
cdef flatiter iterx = <flatiter>PyArray_IterNew(<object>arr)
cdef np.npy_intp size = PyArray_SIZE(arr)
cdef np.npy_intp* boundsPtr = [bounds[0], bounds[1]]
cdef int hoodSize = bounds[1] - bounds[0] + 1
# Create the Python object and keep a reference to it
cdef object niterx_ = PyArray_NeighborhoodIterNew(iterx,
boundsPtr, NPY_NEIGHBORHOOD_ITER_CIRCULAR_PADDING, None)
cdef PyArrayNeighborhoodIterObject *niterx = \
<PyArrayNeighborhoodIterObject *>niterx_
cdef int i,j
cdef np.ndarray[np.int_t, ndim=2] hoods
hoods = np.empty((arr.shape[0], hoodSize), dtype=np.int)
for i in range(iterx.size):
for j in range(niterx.size):
hoods[i,j] = (niterx.dataptr)[0]
PyArrayNeighborhoodIter_Next(niterx)
PyArray_ITER_NEXT(iterx)
PyArrayNeighborhoodIter_Reset(niterx)
return hoods
def test():
x = np.arange(10)
print x
print
print windowed(x, [-1, 3])
print
print windowed(x, [-2, 2])
----------
If you run test(), this is what you should see:
[0 1 2 3 4 5 6 7 8 9]
[[9 0 1 2 3]
[0 1 2 3 4]
[1 2 3 4 5]
[2 3 4 5 6]
[3 4 5 6 7]
[4 5 6 7 8]
[5 6 7 8 9]
[6 7 8 9 0]
[7 8 9 0 1]
[8 9 0 1 2]]
[[8 9 0 1 2]
[9 0 1 2 3]
[0 1 2 3 4]
[1 2 3 4 5]
[2 3 4 5 6]
[3 4 5 6 7]
[4 5 6 7 8]
[5 6 7 8 9]
[6 7 8 9 0]
[7 8 9 0 1]]
windowed(x, [0, 2]) is almost like segment_axis(x, 3, 2, end='wrap').
More information about the NumPy-Discussion
mailing list