MAINT: Added watersheft function in _ni_measure.pyx Note: The functio…

…n is yet not complete.

MAINT: Improvement in watersheft func

Worked on Watershed function

scipy/ndimage/src/_ni_measure.pyx

@@ -1,3 +1,5 @@
+## NOTE: the function watersheft is yetnot complete. It will take me around day 
+## to make it functional
 ######################################################################
 # Cython version of scipy.ndimage.measurement.find_objects() function.
 ######################################################################
@@ -35,7 +37,6 @@ cdef extern from "numpy/arrayobject.h" nogil:
     void *PyDataMem_NEW(size_t)
     void PyDataMem_FREE(void *)
 
-
 ######################################################################
 # Use of Cython's type templates for type specialization
 ######################################################################
@@ -205,3 +206,377 @@ cpdef _findObjects(np.ndarray input, np.intp_t max_label):
         raise
 
     return result
+
+
+## NOTE: the function watersheft is yetnot complete. It will take me around day 
+## to make it functional
+
+##############################################################################
+##############################################################################
+#   Py_WatershedIFT function in cython
+##############################################################################
+##############################################################################
+
+
+DEF WS_MAXDIM = 7
+#############################################################################
+# Fused type delarations
+#############################################################################
+ctypedef fused data_watershed:
+    np.int8_t
+    np.int16_t
+
+ctypedef fused data_markers:
+    np.int8_t
+    np.int16_t
+    np.int32_t
+    np.int64_t
+    np.uint8_t
+    np.uint16_t
+    np.uint32_t
+    np.uint64_t
+
+
+ctypedef fused data_output:
+    np.int8_t
+    np.int16_t
+    np.int32_t
+    np.int64_t
+    np.uint8_t
+    np.uint16_t
+    np.uint32_t
+    np.uint64_t
+
+def get_funcp_watershed(np.ndarray[data_watershed] input, np.ndarray[data_markers] markers,
+                       np.ndarray[data_output] output):
+    return (<Py_intptr_t> get_value[data_watershed],
+            <Py_intptr_t> case_windex1[data_output],
+            <Py_intptr_t> markers_to_output[data_markers, data_output],
+            <Py_intptr_t> case_windex2[data_output])
+
+################################################################################
+# Pointerdeclaration
+#################################################################################
+
+ctypedef np.intp_t (*funcp_watershed)(void *data, np.flatiter x,
+                    PyArrayIterObject *y, np.ndarray input) nogil
+
+ctypedef np.intp_t (*funcp_markers_to_output)(void * data_m, void * data_o, np.flatiter _mi, 
+                    np.flatiter _li)
+
+ctypedef np.intp_t (*funcp_windex1)(void *data, np.intp_t v_index, np.intp_t p_index, np.intp_t *strides, int i_contiguous, 
+                            np.intp_t p_idx, np.intp_t v_idx, np.flatiter _ii, np.intp_t *vval, np.intp_t *pval, np.ndarray input)
+
+ctypedef np.intp_t (*funcp_windex2) (void *data, np.intp_t v_index, np.intp_t p_index, np.intp_t *strides, np.ndarray output, np.ndarray input, 
+    np.intp_t idx, int o_contiguous, int *label)
+
+################################################################################
+# Asociate function Declaration
+#################################################################################
+
+cdef void get_index(data_output *data, np.intp_t index, np.intp_t *strides, np.ndarray array, 
+                    np.intp_t rank, np.intp_t *out, np.intp_t contiguous):
+    cdef int qq
+    cdef np.intp_t cc, idx
+
+    do:
+        if contiguous:
+            out[0] = index * sizeof(data_output)
+        else:
+            idx = index
+            out[0] = 0
+            for qq in range(rank):
+                cc = idx / strides[qq]
+                idx -= cc * strides[qq]
+                out[0] += array.strides[qq] * cc
+    while False 
+
+cdef np.intp_t get_value(data_watershed *data, np.flatiter _iti, PyArrayObject *iti, 
+                         np.ndarray array):
+    #use byteswapped technique
+    return <np.intp_t>((<data_watershed *> np.PyArray_ITER_DATA(_iti))[0])
+
+cdef np.intp_t markers_to_output(data_markers *data_m, data_output *data_p, np.flatiter _mi, 
+                                np.flatiter _li):
+    #use byteswapping technique
+    cdef np.intp_t temp = 4
+    temp = (<data_markers *> np.PyArray_ITER_DATA(_mi))[0]
+    (<data_output *> np.PyArray_ITER_DATA(_mi))[0] = < data_output > temp
+    return temp
+
+cdef np.intp_t case_windex1(data_output *data, np.intp_t v_index, np.intp_t p_index, np.intp_t *strides, int i_contiguous, 
+                            np.intp_t p_idx, np.intp_t v_idx, np.flatiter _ii, np.intp_t *vval, np.intp_t *pval, np.ndarray input):
+
+    # Another Function defination
+    vval[0] = (<data_output *>(np.PyArray_ITER_DATA(_ii) + v_idx))[0]
+    pval[0] = (<data_output *>(np.PyArray_ITER_DATA(_ii) + p_idx))[0]
+
+cdef np.intp_t case_windex2(data_output *data, np.intp_t v_index, np.intp_t p_index, np.intp_t *strides, np.ndarray output, np.ndarray input, 
+    np.intp_t idx, int o_contiguous, int *label):
+    # Another function definations
+    get_index(data, v_index, strides, output, input.ndim, idx, o_contiguous)
+    label[0] = (<data_output *> data)[0]
+    #same function call
+    get_index(data, p_index, strides, output.strides, out, )
+    (<data_output *> data)[0] = label[0]
+
+    # cdef int del = 0
+
+    return 1
+
+#############################################################################
+# Basic function
+#############################################################################
+
+cdef struct WatershedElement:
+    np.uintp_t index
+    np.uint16_t cost
+    void *next, *prev
+    np.uint8_t done
+
+cpdef int watershed_ift(np.ndarray input, np.ndarray markers, np.ndarray structure, 
+                        np.ndarray output):
+    cdef:
+        funcs = get_funcp_watershed(input.take([0]), markers.take([0]), output.take([0]))
+        int ll, jj, hh, kk, i_contiguous, o_contiguous, label
+        np.intp_t size, maxval, nneigh, ssize, ival
+        np.intp_t strides[WS_MAXDIM], coordinates[WS_MAXDIM]
+        np.intp_t *nstrides = NULL
+        bint *ps = NULL
+        np.flatiter _ii, _li, _mi
+        PyArrayIterObject *ii, *li, *mi
+        WatershedElement *temp = NULL, **first = NULL, **last = NULL
+
+
+    # if input.ndim > WS_MAXDIM:
+    #     Raise RuntimeError("Too many dimensions")
+
+    ssize = structure.ndim
+    size = input.ndim
+
+    temp = <WatershedElement *> PyDataMem_NEW(size * sizeof(WatershedElement))
+    # Error condition
+
+    # Iterator inititalization
+    _ii = np.PyArray_IterNew(input)
+    _li = np.PyArray_IterNew(output)
+    _mi = np.PyArray_IterNew(markers)
+
+    ii = <PyArrayIterObject *> _ii
+    li = <PyArrayIterObject *> _li
+    mi = <PyArrayIterObject *> _mi
+
+    cdef funcp_watershed get_value = <funcp_watershed> <void *> <Py_intptr_t> funcs[0]
+    cdef funcp_windex1 case_windex1 = <funcp_windex1> <void *> <Py_intptr_t> funcs[1]
+    cdef funcp_markers_to_output markers_to_output = <funcp_markers_to_output> <void *> <Py_intptr_t> funcs[2]
+    cdef funcp_windex2 case_windex2 = <funcp_windex2> <void *> <Py_intptr_t> funcs[3]
+
+    for jj in range(size):
+        # Need value in function in ival from pi using fused_type
+        ival = get_value(np.PyArray_ITER_DATA(_ii), _ii, ii, input)
+
+        temp[jj].index = jj
+        temp[jj].done = 0
+        if ival > maxval:
+            maxval = ival
+
+        np.PyArray_ITER_NEXT(_ii)
+
+    # Allocate and initialize the storage for the queue
+    first = <WatershedElement ** >  PyDataMem_NEW((maxval + 1) * sizeof(WatershedElement *))
+    last =  <WatershedElement ** > PyDataMem_NEW((maxval + 1) * sizeof(WatershedElement *))
+    # error in allocations
+
+
+    for hh in range(maxval):
+        first[hh] = last[hh] = NULL
+
+    for ll in range(input.ndim):
+        coordinates[ll] = 0
+
+    for jj in range(size):
+        label = markers_to_output(np.PyArray_ITER_DATA(_mi), np.PyArray_ITER_DATA(_li), _mi, _li)
+        np.PyArray_ITER_NEXT(_mi)
+        np.PyArray_ITER_NEXT(_li)
+        if label != 0:
+            temp[jj].cost = 0
+            if first[0] is NULL:
+                first[0] = &(temp[jj])
+                # beware here.. could get erreors
+                first[0].prev = NULL
+                first[0].next = NULL
+                last[0] = first[0]
+
+            else:
+                if label > 0:
+                    # object markers are enqueued at the beginning, so they
+                    # are processed first.
+                    temp[jj].next = first[0]
+                    temp[jj].prev = NULL
+                    first[0].prev = &(temp[jj])
+                    first[0] = &(temp[jj])
+
+                else:
+                    # background markers are enqueued at the end, so they are
+                    # processed after the object markers. 
+                    temp[jj].next = NULL
+                    temp[jj].prev = last[0]
+                    last[0].next = &(temp[jj])
+                    last[0] = &(temp[jj])
+
+        else:
+            # This node is not a marker 
+            temp[jj].cost = maxval + 1
+            temp[jj].next = NULL
+            temp[jj].prev = NULL
+
+        ll = input.ndim - 1
+        while ll >=0:
+            if coordinates[ll] < input.dimensions[ll] - 1:
+                coordinates[ll] += 1
+                break
+
+            else:
+                coordinates[ll] = 0
+            ll -= 1
+
+    nneigh = 0
+    for kk in range(ssize):
+        if ps[kk] and kk != ssize/2:
+            nneigh += 1
+
+    nstrides = <np.intp_t *> PyDataMem_NEW(nneigh * sizeof(np.intp_t))
+    #Error in allocation
+
+    strides[input.ndim -1] = 1
+
+    for ll in range(input.ndim -1) :
+        strides[ll] = input.dimensions[ll + 1] * strides[ll + 1]
+
+    for ll in range(input.ndim):
+        coordinates[ll] = -1
+
+    for kk in range(nneigh):
+        nstrides[kk] = 0
+
+    jj = 0
+    cdef int offset = 0
+
+    for kk in range(ssize):
+        if ps[kk] is not 0:
+            for ll in range(input.ndim):
+                offset += coordinates[ll] * strides[ll]
+            if offset is not 0:
+                nstrides[jj] += offset
+                jj += 1
+
+
+        ll = input.ndim -1
+        while ll >= 0:
+            if coordinates[ll] < 1:
+                coordinates[ll] += 1
+
+            else:
+                coordinates[ll] = -1
+
+
+    # Propogation Phase
+    cdef:
+        WatershedElement *v, *p, *prev, *next
+        np.intp_t v_index, p_index, idx, cc
+        int qq, outside
+        np.intp_t max, pval, vval, wvp, pcost, p_idx, v_idx
+
+    for jj in range(maxval +1):
+        while first[jj] is not NULL:
+            v = first[jj]
+            first[jj] = <WatershedElement *>first[jj].next
+            if first[jj] is not NULL:
+                first[jj].prev = NULL
+
+            v.prev = NULL
+            v.next = NULL
+
+            v.done = 1
+
+            for hh in range(nneigh):
+                v_index = v.index
+                p_index = v.index
+                outside = 0
+                p_index += nstrides[hh]
+                # Check if the neighbour is within the extent of the array
+                idx = p_index
+                for qq in range(input.ndim):
+                    cc = idx / strides[qq]
+                    if cc < 0 or cc >=input.dimensions[qq]:
+                        outside = 1
+                        break
+
+                if outside is not 0:
+                    p = &(temp[p_index])
+                    # If the neighbour is not processed Yet
+                    if p.done is 0:
+
+                        case_windex1(np.PyArray_ITER_DATA(_ii), v_index, p_index, strides,
+                                     i_contiguous, p_idx, v_idx, _ii, &vval, &pval, input)
+
+                        # Calculate Cost
+                        wvp = pval - vval
+                        if wvp < 0:
+                            wvp = -wvp
+                        # Find the maximum of this cost and the current element cost
+                        pcost = p.cost
+                        if v.cost > wvp:
+                            max = v.cost
+
+                        else:
+                            max = wvp
+
+                        if max < pcost:
+                            # If this maximum is less than the neighbors cost,
+                            # adapt the cost and the label of the neighbor: 
+                            p.cost = max
+                            case_windex2(np.PyArray_ITER_DATA(_li), v_index, p_index, strides, output, input, idx,
+                                        o_contiguous, &label)
+
+                            # If the neighbor is in a queue, remove it:
+                            if p.next or p.prev is not NULL:
+                                prev = <WatershedElement *> p.prev
+                                next = <WatershedElement *> p.next
+                                if first[pcost] == p:
+                                    first[pcost] = next
+
+                                if last[pcost] == p:
+                                    last[pcost] = p
+
+                                if prev is not NULL:
+                                    prev.next = next
+
+                                if next is not NULL:
+                                    next.prev = prev
+
+                            # Insert the neighbor in the appropiate queue:
+                            if label < 0:
+                                p.prev = last[max]
+                                p.next = NULL
+                                if last[max] is not NULL:
+                                    last[max].next = p
+                                last[max] = p
+                                if first[max] is NULL:
+                                    first[max] = p
+
+                            else:
+                                p.next = first[max]
+                                p.prev = NULL
+                                if first[max] is not NULL:
+                                    first[max].prev = p
+                                first[max] = p
+                                if last[max] is NULL:
+                                    last[max] = p
+
+    PyDataMem_FREE(temp)
+    PyDataMem_FREE(first)
+    PyDataMem_FREE(last)
+    PyDataMem_FREE(nstrides)
+
+    return 1