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, ¤t_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