Example Usage of Neighborhood Iterator in Cython
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').
Hi, On Mon, Oct 17, 2011 at 9:19 PM, T J <tjhnson@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
On Mon, Oct 17, 2011 at 12:45 PM, eat <e.antero.tammi@gmail.com> wrote:
Just wondering what are the main benefits, of your approach, comparing to simple:
As I hinted, my goal was not to construct a "practical" example, but rather, to demonstrate how to use the neighborhood iterator in Cython. Roll and mod are quite nice. :) Now imagine working with higher dimensional arrays with more exotic neighborhoods (like the letter X).
Just in time! I was just working on a cythonic replacement to ndimage.generic_filter (well, I took a a short two years break in the middle). thank you very much, Nadav. ________________________________________ From: numpy-discussion-bounces@scipy.org [numpy-discussion-bounces@scipy.org] On Behalf Of T J [tjhnson@gmail.com] Sent: 17 October 2011 23:16 To: Discussion of Numerical Python Subject: Re: [Numpy-discussion] Example Usage of Neighborhood Iterator in Cython On Mon, Oct 17, 2011 at 12:45 PM, eat <e.antero.tammi@gmail.com> wrote:
Just wondering what are the main benefits, of your approach, comparing to simple:
As I hinted, my goal was not to construct a "practical" example, but rather, to demonstrate how to use the neighborhood iterator in Cython. Roll and mod are quite nice. :) Now imagine working with higher dimensional arrays with more exotic neighborhoods (like the letter X). _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
participants (3)
-
eat -
Nadav Horesh -
T J