1#cython: cdivision=True
2#cython: boundscheck=False
3#cython: nonecheck=False
4#cython: wraparound=False
5
6"""Cython code used in extrema.py."""
7
8cimport numpy as cnp
9cnp.import_array()
10
11
12# Must be defined to use QueueWithHistory
13ctypedef Py_ssize_t QueueItem
14
15include "_queue_with_history.pxi"
16
17
18ctypedef fused dtype_t:
19    cnp.uint8_t
20    cnp.uint16_t
21    cnp.uint32_t
22    cnp.uint64_t
23    cnp.int8_t
24    cnp.int16_t
25    cnp.int32_t
26    cnp.int64_t
27    cnp.float32_t
28    cnp.float64_t
29
30
31# Definition of flag values used for `flags` in _local_maxima & _fill_plateau
32cdef:
33    # First or last index in a dimension
34    unsigned char BORDER_INDEX = 3
35    # Potentially part of a maximum
36    unsigned char CANDIDATE = 2
37    # Index was queued (flood-fill) and might still be part of maximum OR
38    # when evaluation is complete this flag value marks definite local maxima
39    unsigned char QUEUED_CANDIDATE = 1
40    # None of the above is true
41    unsigned char NOT_MAXIMUM = 0
42
43
44def _local_maxima(dtype_t[::1] image not None,
45                  unsigned char[::1] flags,
46                  Py_ssize_t[::1] neighbor_offsets not None):
47    """Detect local maxima in n-dimensional array.
48
49    Inner function to `local_maxima` that detects all local maxima (including
50    plateaus) in the image. The result is stored inplace inside `flags` with
51    the value of "QUEUED_CANDIDATE".
52
53    Parameters
54    ----------
55    image : ndarray, one-dimensional
56        The raveled view of a n-dimensional array.
57    flags : ndarray
58        An array of flags that is used to store the state of each pixel during
59        evaluation and is MODIFIED INPLACE. Initially, pixels that border the
60        image edge must be marked as "BORDER_INDEX" while all other pixels
61        should be marked with "NOT_MAXIMUM".
62    neighbor_offsets : ndarray
63        A one-dimensional array that contains the offsets to find the
64        connected neighbors for any index in `image`.
65    """
66    cdef:
67        QueueWithHistory queue
68        unsigned char last_dim_in_neighbors
69
70    # This needs the GIL
71    last_dim_in_neighbors = -1 in neighbor_offsets and 1 in neighbor_offsets
72    with nogil:
73        if last_dim_in_neighbors:
74            # If adjacent pixels in the last dimension are part of the
75            # neighborhood, the number of candidates which have to be evaluated
76            # in the second algorithmic step `_fill_plateau` can be reduced
77            # ahead of time with this function
78            _mark_candidates_in_last_dimension(image, flags)
79        else:
80            # Otherwise simply mark all pixels (except border) as candidates
81            _mark_candidates_all(flags)
82
83        # Initialize a buffer used to queue positions while evaluating each
84        # potential maximum (flagged with 2)
85        queue_init(&queue, 64)
86        try:
87            for i in range(image.shape[0]):
88                if flags[i] == CANDIDATE:
89                    # Index is potentially part of a maximum:
90                    # Find all samples which are part of the plateau and fill
91                    # with 0 or 1 depending on whether it's a true maximum
92                    _fill_plateau(image, flags, neighbor_offsets, &queue, i)
93        finally:
94            # Ensure that memory is released again
95            queue_exit(&queue)
96
97
98cdef inline void _mark_candidates_in_last_dimension(
99        dtype_t[::1] image, unsigned char[::1] flags) nogil:
100    """Mark local maxima in last dimension.
101
102    This function considers only the last dimension of the image and marks
103    pixels with the "CANDIDATE" flag if it is a local maximum.
104
105    Parameters
106    ----------
107    image :
108        The raveled view of a n-dimensional array.
109    flags :
110        An array of flags that is used to store the state of each pixel during
111        evaluation.
112
113    Notes
114    -----
115    By evaluating this necessary but not sufficient condition first, usually a
116    significant amount of pixels can be rejected without having to evaluate the
117    entire neighborhood of their plateau. This can reduces the number of
118    candidates that need to be evaluated with the more expensive flood-fill
119    performed in `_fill_plateaus`.
120
121    However this is only possible if the adjacent pixels in the last dimension
122    are part of the defined neighborhood (see argument `neighbor_offsets`
123    in `_local_maxima`).
124    """
125    cdef Py_ssize_t i, i_ahead
126    i = 1
127    while i < image.shape[0]:
128        if image[i - 1] < image[i] and flags[i] != BORDER_INDEX:
129            # Potential maximum (in last dimension) is found, find
130            # other edge of current plateau or "edge of dimension"
131            i_ahead = i + 1
132            while (
133                image[i] == image[i_ahead] and
134                flags[i_ahead] != BORDER_INDEX
135            ):
136                i_ahead += 1
137            if image[i] > image[i_ahead]:
138                # Found local maximum (in one dimension), mark all
139                # parts of the plateau as potential maximum
140                flags[i:i_ahead] = CANDIDATE
141            i = i_ahead
142        else:
143            i += 1
144
145
146cdef inline void _mark_candidates_all(unsigned char[::1] flags) nogil:
147    """Mark all pixels as potential maxima, exclude border pixels.
148
149    This function marks pixels with the "CANDIDATE" flag if they aren't the
150    first or last index in any dimension (not flagged with "BORDER_INDEX").
151
152    Parameters
153    ----------
154    flags :
155        An array of flags that is used to store the state of each pixel during
156        evaluation.
157    """
158    cdef Py_ssize_t i = 1
159    while i < flags.shape[0]:
160        if flags[i] != BORDER_INDEX:
161            flags[i] = CANDIDATE
162        i += 1
163
164
165cdef inline void _fill_plateau(
166        dtype_t[::1] image, unsigned char[::1] flags,
167        Py_ssize_t[::1] neighbor_offsets, QueueWithHistory* queue_ptr,
168        Py_ssize_t start_index) nogil:
169    """Fill with 1 if plateau is local maximum else with 0.
170
171    Parameters
172    ----------
173    image :
174        The raveled view of a n-dimensional array.
175    flags :
176        An array of flags that is used to store the state of each pixel during
177        evaluation.
178    neighbor_offsets :
179        A one-dimensional array that contains the offsets to find the
180        connected neighbors for any index in `image`.
181    queue_ptr :
182        Pointer to initialized queue.
183    start_index :
184        Start position for the flood-fill.
185    """
186    cdef:
187        dtype_t height
188        unsigned char is_maximum
189        QueueItem current_index, neighbor
190
191    height = image[start_index]  # Height of the evaluated plateau
192    is_maximum = 1 # Boolean flag, initially assume true
193
194    # Queue start position after clearing the buffer
195    # which might have been used already
196    queue_clear(queue_ptr)
197    queue_push(queue_ptr, &start_index)
198    flags[start_index] = QUEUED_CANDIDATE
199
200    # Break loop if all queued positions were evaluated
201    while queue_pop(queue_ptr, &current_index):
202        # Look at all neighboring samples
203        for i in range(neighbor_offsets.shape[0]):
204            neighbor = current_index + neighbor_offsets[i]
205
206            if image[neighbor] == height:
207                # Value is part of plateau
208                if flags[neighbor] == BORDER_INDEX:
209                    # Plateau touches border and can't be maximum
210                    is_maximum = 0
211                elif flags[neighbor] != QUEUED_CANDIDATE:
212                    # Index wasn't queued already, do so now
213                    queue_push(queue_ptr, &neighbor)
214                    flags[neighbor] = QUEUED_CANDIDATE
215
216            elif image[neighbor] > height:
217                # Current plateau can't be maximum because it borders a
218                # larger one
219                is_maximum = 0
220
221    if not is_maximum:
222        queue_restore(queue_ptr)
223        # Initial guess was wrong -> flag as NOT_MAXIMUM
224        while queue_pop(queue_ptr, &neighbor):
225            flags[neighbor] = NOT_MAXIMUM
226