MAINT: Rewrote scipy.ndimage.find_object function in cython

.gitignore

@@ -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

scipy/ndimage/measurements.py

@@ -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',
@@ -264,10 +265,13 @@ def find_objects(input, max_label=0):
     if numpy.iscomplexobj(input):
         raise TypeError('Complex type not supported')
 
+    if input.dtype == np.bool:
+        input = input.astype(int)
+
     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):

scipy/ndimage/setup.py

@@ -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

scipy/ndimage/src/_ni_measure.pyx

@@ -0,0 +1,173 @@
+######################################################################
+# 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 PyArrayObject:
+        int nd
+
+    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, np.intp_t max_label, 
+               int *regions, int rank) nogil
+
+def get_funcs(np.ndarray[data_t] input):
+    return (<Py_intptr_t> findObjectsPoint[data_t])
+
+######################################################################
+# Update Regions According to Input Data Type
+######################################################################
+
+cdef inline findObjectsPoint(data_t *data, np.flatiter _iti, np.intp_t max_label, 
+            int *regions, int rank):
+    cdef int kk =0
+    cdef np.intp_t cc
+
+    # only integer or boolean values are allowed, since s_index is being used in indexing
+    cdef np.intp_t s_index = <np.intp_t> (<data_t *> data)[0] - 1
+
+    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 = (<PyArrayIterObject *>_iti).coordinates[kk]
+                    regions[s_index + kk] = cc
+                    regions[s_index + kk + rank] = cc + 1
+
+            else:
+                for kk in range(rank):
+                    cc = (<PyArrayIterObject *>_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_raw, np.intp_t max_label):
+    cdef:
+        funcs = get_funcs(input_raw.take([0]))
+
+        int ii, rank, size_regions
+        int start, jj, idx, end
+        int *regions = NULL
+
+        np.flatiter _iti
+        PyArrayIterObject *iti
+
+        func_p findObjectsPoint = <func_p> <void *> <Py_intptr_t> funcs
+
+        int flags = np.NPY_CONTIGUOUS | np.NPY_NOTSWAPPED | np.NPY_ALIGNED
+
+    input = np.PyArray_FROM_OF(input_raw, flags)
+
+    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 = <int *> PyDataMem_NEW(size_regions * sizeof(int))
+        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, 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)
+        raise
+
+    return result

scipy/ndimage/src/nd_image.c

@@ -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;
@@ -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,