Skip to content

Commit

Permalink
MAINT: Rewriting of find_objects in cython
Browse files Browse the repository at this point in the history
  • Loading branch information
bewithaman committed Jul 2, 2015
1 parent 9f936e4 commit a83b747
Show file tree
Hide file tree
Showing 7 changed files with 215 additions and 186 deletions.
1 change: 1 addition & 0 deletions .gitignore
View file
Expand Up @@ -155,6 +155,7 @@ scipy/linalg/fblas.pyf
scipy/linalg/flapack.pyf
scipy/linalg/src/id_dist/src/*_subr_*.f
scipy/ndimage/src/_ni_label.c
scipy/ndimage/src/_ni_measure.c
scipy/optimize/cobyla/_cobylamodule.c
scipy/optimize/lbfgsb/_lbfgsbmodule.c
scipy/optimize/minpack2/minpack2module.c
Expand Down
3 changes: 2 additions & 1 deletion scipy/ndimage/measurements.py
View file
Expand Up @@ -35,6 +35,7 @@
from . import _ni_support
from . import _ni_label
from . import _nd_image
from . import _ni_measure
from . import morphology

__all__ = ['label', 'find_objects', 'labeled_comprehension', 'sum', 'mean',
Expand Down Expand Up @@ -267,7 +268,7 @@ def find_objects(input, max_label=0):
if max_label < 1:
max_label = input.max()

return _nd_image.find_objects(input, max_label)
return _ni_measure._findObjects(input, max_label)


def labeled_comprehension(input, labels, index, func, out_dtype, default, pass_positions=False):
Expand Down
5 changes: 5 additions & 0 deletions scipy/ndimage/setup.py
View file
Expand Up @@ -23,6 +23,11 @@ def configuration(parent_package='', top_path=None):
include_dirs=['src']+[get_include()],
)

config.add_extension("_ni_measure",
sources=["src/_ni_measure.c",],
include_dirs=['src']+[get_include()],
)

config.add_data_dir('tests')

return config
Expand Down
207 changes: 207 additions & 0 deletions scipy/ndimage/src/_ni_measure.pyx
View file
@@ -0,0 +1,207 @@
######################################################################
# Cython version of scipy.ndimage.measurement.find_objects() function.
######################################################################

cimport cython
from cython cimport sizeof
import numpy as np
cimport numpy as np

np.import_array()

cdef extern from *:
ctypedef int Py_intptr_t

ctypedef void (*PyArray_CopySwapFunc)(void *, void *, int, void *)

cdef extern from "numpy/arrayobject.h" nogil:
ctypedef struct PyArrayIterObject:
np.intp_t *coordinates
char *dataptr
np.npy_bool contiguous

ctypedef struct PyArray_ArrFuncs:
PyArray_CopySwapFunc *copyswap

ctypedef struct PyArray_Descr:
PyArray_ArrFuncs *f

ctypedef struct PyArrayObject:
PyArray_Descr *descr
int nd

PyArray_Descr *PyArray_DESCR(PyArrayObject* arr)

void *PyDataMem_NEW(size_t)
void PyDataMem_FREE(void *)


######################################################################
# Use of Cython's type templates for type specialization
######################################################################

# Only integer values are allowed.
ctypedef fused data_t:
np.int8_t
np.int16_t
np.int32_t
np.int64_t
np.uint8_t
np.uint16_t
np.uint32_t
np.uint64_t

#####################################################################
# Function Specializers and asociate function for using fused type
#####################################################################

ctypedef void (*func_p)(void *data, np.flatiter iti, PyArrayIterObject *iti,
np.ndarray input, np.intp_t max_label, np.intp_t* regions,
int rank) nogil

def get_funcs(np.ndarray[data_t] input):
return (<Py_intptr_t> findObjectsPoint[data_t])


######################################################################
# Dereferencing and Dealing with Misalligned pointers
######################################################################

ctypedef data_t (* func2_p)(data_t *, np.flatiter, np.ndarray)

cdef data_t get_from_iter(data_t *data, np.flatiter iter, np.ndarray arr):
return (<data_t *>np.PyArray_ITER_DATA(iter))[0]

cdef data_t get_misaligned_from_iter(data_t *data, np.flatiter iter, np.ndarray arr):

cdef data_t ret = 0
cdef PyArray_CopySwapFunc copyswap = <PyArray_CopySwapFunc> <void *> PyArray_DESCR(<PyArrayObject*> arr).f.copyswap

copyswap(&ret, np.PyArray_ITER_DATA(iter), 1,<void *> arr)

return ret

######################################################################
# Update Regions According to Input Data Type
######################################################################

cdef inline findObjectsPoint(data_t *data, np.flatiter _iti, PyArrayIterObject *iti,
np.ndarray input, np.intp_t max_label, np.intp_t* regions,
int rank):
cdef int kk =0
cdef np.intp_t cc

cdef func2_p deref_p
if np.PyArray_ISBYTESWAPPED(input) == True:
deref_p = get_misaligned_from_iter

else:
deref_p = get_from_iter

# only integer or boolean values are allowed, since s_index is being used in indexing
# cdef np.uintp_t s_index = deref_p(data, _iti, input) - 1
cdef np.intp_t s_index = <np.intp_t>((<data_t *>np.PyArray_ITER_DATA(_iti))[0])

if s_index >=0 and s_index < max_label:
if rank > 0:
s_index *= 2 * rank
if regions[s_index] < 0:
for kk in range(rank):
cc = iti.coordinates[kk]
regions[s_index + kk] = cc
regions[s_index + kk + rank] = cc + 1

else:
for kk in range(rank):
cc = iti.coordinates[kk]
if cc < regions[s_index + kk]:
regions[s_index + kk] = cc
if cc +1 > regions[s_index + kk + rank]:
regions[s_index + kk + rank] = cc + 1
else:
regions[s_index] = 1

return 1

######################################################################
# Implementaion of find_Objects function:-
######################################################################

cpdef _findObjects(np.ndarray input, np.intp_t max_label):
cdef:
funcs = get_funcs(input.take([0]))

int ii, rank, size_regions
int start, jj, idx, end
np.intp_t *regions = NULL

np.flatiter _iti
PyArrayIterObject *iti

func_p findObjectsPoint = <func_p> <void *> <Py_intptr_t> funcs

rank = input.ndim

# Array Iterator defining and Initialization:
_iti = np.PyArray_IterNew(input)
iti = <PyArrayIterObject *> _iti

# If iterator is contiguous, PyArray_ITER_NEXT will treat it as 1D Array
iti.contiguous = 0

# Declaring output array
size_regions = 0

if max_label < 0:
max_label = 0
elif max_label > 0:
if rank > 0:
size_regions = 2 * max_label * rank

else:
size_regions = max_label

regions = <np.intp_t *> PyDataMem_NEW(size_regions * sizeof(np.intp_t))
if regions == NULL:
raise MemoryError()

else:
regions = NULL

try:
with nogil:
if rank > 0:
for jj in range(size_regions):
regions[jj] = -1

# Iteration over array:
while np.PyArray_ITER_NOTDONE(_iti):
findObjectsPoint(np.PyArray_ITER_DATA(_iti), _iti, iti, input,
max_label, regions, rank)
np.PyArray_ITER_NEXT(_iti)

result = []

for ii in range(max_label):
if rank > 0:
idx = 2 * rank * ii

else:
idx = ii

if regions[idx] >= 0:
slc = ()
for jj in range(rank):
start = regions[idx + jj]
end = regions[idx + jj + rank]

slc += (slice(start, end),)
result.append(slc)

else:
result.append(None)
except:
# clean up and re-raise
PyDataMem_FREE(regions)

return result
91 changes: 0 additions & 91 deletions scipy/ndimage/src/nd_image.c
View file
Expand Up @@ -629,95 +629,6 @@ static PyObject *Py_ZoomShift(PyObject *obj, PyObject *args)
return PyErr_Occurred() ? NULL : Py_BuildValue("");
}

static PyObject *Py_FindObjects(PyObject *obj, PyObject *args)
{
PyArrayObject *input = NULL;
PyObject *result = NULL, *tuple = NULL, *start = NULL, *end = NULL;
PyObject *slc = NULL;
int jj;
npy_intp max_label;
npy_intp ii, *regions = NULL;

if (!PyArg_ParseTuple(args, "O&n",
NI_ObjectToInputArray, &input, &max_label))
goto exit;

if (max_label < 0)
max_label = 0;
if (max_label > 0) {
if (input->nd > 0) {
regions = (npy_intp*)malloc(2 * max_label * input->nd *
sizeof(npy_intp));
} else {
regions = (npy_intp*)malloc(max_label * sizeof(npy_intp));
}
if (!regions) {
PyErr_NoMemory();
goto exit;
}
}

if (!NI_FindObjects(input, max_label, regions))
goto exit;

result = PyList_New(max_label);
if (!result) {
PyErr_NoMemory();
goto exit;
}

for(ii = 0; ii < max_label; ii++) {
npy_intp idx = input->nd > 0 ? 2 * input->nd * ii : ii;
if (regions[idx] >= 0) {
PyObject *tuple = PyTuple_New(input->nd);
if (!tuple) {
PyErr_NoMemory();
goto exit;
}
for(jj = 0; jj < input->nd; jj++) {
start = PyLong_FromSsize_t(regions[idx + jj]);
end = PyLong_FromSsize_t(regions[idx + jj + input->nd]);
if (!start || !end) {
PyErr_NoMemory();
goto exit;
}
slc = PySlice_New(start, end, NULL);
if (!slc) {
PyErr_NoMemory();
goto exit;
}
Py_XDECREF(start);
Py_XDECREF(end);
start = end = NULL;
PyTuple_SetItem(tuple, jj, slc);
slc = NULL;
}
PyList_SetItem(result, ii, tuple);
tuple = NULL;
} else {
Py_INCREF(Py_None);
PyList_SetItem(result, ii, Py_None);
}
}

Py_INCREF(result);

exit:
Py_XDECREF(input);
Py_XDECREF(result);
Py_XDECREF(tuple);
Py_XDECREF(start);
Py_XDECREF(end);
Py_XDECREF(slc);
free(regions);
if (PyErr_Occurred()) {
Py_XDECREF(result);
return NULL;
} else {
return result;
}
}

static PyObject *Py_WatershedIFT(PyObject *obj, PyObject *args)
{
PyArrayObject *input = NULL, *output = NULL, *markers = NULL;
Expand Down Expand Up @@ -917,8 +828,6 @@ static PyMethodDef methods[] = {
METH_VARARGS, NULL},
{"zoom_shift", (PyCFunction)Py_ZoomShift,
METH_VARARGS, NULL},
{"find_objects", (PyCFunction)Py_FindObjects,
METH_VARARGS, NULL},
{"watershed_ift", (PyCFunction)Py_WatershedIFT,
METH_VARARGS, NULL},
{"distance_transform_bf", (PyCFunction)Py_DistanceTransformBruteForce,
Expand Down

0 comments on commit a83b747

Please sign in to comment.