[Numpy-discussion] Example Usage of Neighborhood Iterator in Cython

[T J]

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').