Hi, On Mon, Oct 17, 2011 at 9:19 PM, T J <tjhn...@gmail.com> wrote:
> 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'). > Just wondering what are the main benefits, of your approach, comparing to simple: In []: a= arange(5) In []: n= 10 In []: b= arange(n)[:, None] In []: mod(a+ roll(b, 1), n) Out[]: array([[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]]) In []: mod(a+ roll(b, 2), n) Out[]: array([[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]]) Regards, eat > _______________________________________________ > NumPy-Discussion mailing list > NumPy-Discussion@scipy.org > http://mail.scipy.org/mailman/listinfo/numpy-discussion >
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
