scikit-image

1
cimport numpy as cnp
2
cimport cython
3
cnp.import_array()
4

5

6
ctypedef fused np_ints:
7
    cnp.int32_t
8
    cnp.int64_t
9

10
ctypedef fused np_uints:
11
    cnp.uint32_t
12
    cnp.uint64_t
13

14

15
@cython.boundscheck(False)
16
@cython.nonecheck(False)
17
@cython.wraparound(False)
18
def reconstruction_loop(cnp.ndarray[dtype=np_uints, ndim=1,
19
                                    negative_indices=False, mode='c'] aranks,
20
                        cnp.ndarray[dtype=np_ints, ndim=1,
21
                                    negative_indices=False, mode='c'] aprev,
22
                        cnp.ndarray[dtype=np_ints, ndim=1,
23
                                    negative_indices=False, mode='c'] anext,
24
                        cnp.ndarray[dtype=np_ints, ndim=1,
25
                                    negative_indices=False, mode='c'] astrides,
26
                        Py_ssize_t current_idx,
27
                        Py_ssize_t image_stride):
28
    """The inner loop for reconstruction.
29

30
    This algorithm uses the rank-order of pixels. If low intensity pixels have
31
    a low rank and high intensity pixels have a high rank, then this loop
32
    performs reconstruction by dilation. If this ranking is reversed, the
33
    result is reconstruction by erosion.
34

35
    For each pixel in the seed image, check its neighbors. If its neighbor's
36
    rank is below that of the current pixel, replace the neighbor's rank with
37
    the rank of the current pixel. This dilation is limited by the mask, i.e.
38
    the rank at each pixel cannot exceed the mask as that pixel.
39

40
    Parameters
41
    ----------
42
    aranks : array
43
        The rank order of the flattened seed and mask images.
44
    aprev, anext: arrays
45
        Indices of previous and next pixels in rank sorted order.
46
    astrides : array
47
        Strides to neighbors of the current pixel.
48
    current_idx : int
49
        Index of highest-ranked pixel used as starting point in loop.
50
    image_stride : int
51
        Stride between seed image and mask image in `aranks`.
52
    """
53
    cdef unsigned int neighbor_rank, current_rank, mask_rank
54
    cdef int i, neighbor_idx, current_link, nprev, nnext
55
    cdef int nstrides = astrides.shape[0]
56
    cdef np_uints *ranks = <np_uints *>(aranks.data)
57
    cdef np_ints *prev = <np_ints *>(aprev.data)
58
    cdef np_ints *next = <np_ints *>(anext.data)
59
    cdef np_ints *strides = <np_ints *>(astrides.data)
60

61
    with nogil:
62
        while current_idx != -1:
63
            if current_idx < image_stride:
64
                current_rank = ranks[current_idx]
65
                if current_rank == 0:
66
                    break
67
                for i in range(nstrides):
68
                    neighbor_idx = current_idx + strides[i]
69
                    neighbor_rank = ranks[neighbor_idx]
70
                    # Only propagate neighbors ranked below the current rank
71
                    if neighbor_rank < current_rank:
72
                        mask_rank = ranks[neighbor_idx + image_stride]
73
                        # Only propagate neighbors ranked below the mask rank
74
                        if neighbor_rank < mask_rank:
75
                            # Raise the neighbor to the mask rank if
76
                            # the mask ranked below the current rank
77
                            if mask_rank < current_rank:
78
                                current_link = neighbor_idx + image_stride
79
                                ranks[neighbor_idx] = mask_rank
80
                            else:
81
                                current_link = current_idx
82
                                ranks[neighbor_idx] = current_rank
83
                            # unlink the neighbor
84
                            nprev = prev[neighbor_idx]
85
                            nnext = next[neighbor_idx]
86
                            next[nprev] = nnext
87
                            if nnext != -1:
88
                                prev[nnext] = nprev
89
                            # link to the neighbor after the current link
90
                            nnext = next[current_link]
91
                            next[neighbor_idx] = nnext
92
                            prev[neighbor_idx] = current_link
93
                            if nnext >= 0:
94
                                prev[nnext] = neighbor_idx
95
                                next[current_link] = neighbor_idx
96
            current_idx = next[current_idx]
97