8000 Merge pull request #3022 from lagru/maxima2d · scikit-image/scikit-image@912d9d2 · GitHub
[go: up one dir, main page]
Skip to content

Commit 912d9d2

Browse files
stefanvstefanv
authored
Merge pull request #3022 from lagru/maxima2d
Replace morphology.local_maxima with faster flood-fill based Cython version
2 parents e00b1a0 + 1b67372 commit 912d9d2

File tree

9 files changed

+1014
-286
lines changed

9 files changed

+1014
-286
lines changed

‎CONTRIBUTORS.txt

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -256,3 +256,6 @@
256256

257257
- David Volgyes
258258
Unsharp masking
259+
260+
- Lars Grüter
261+
Flood-fill based local maxima detection

‎doc/release/release_dev.rst

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -17,11 +17,15 @@ New Features
1717
------------
1818

1919
- unsharp mask filtering (#2772)
20+
- New options ``connectivity``, ``indices`` and ``allow_borders`` for
21+
``skimage.morphology.local_maxima`` and ``.local_minima``. #3022
2022

2123

2224
Improvements
2325
------------
2426

27+
- Performance of ``skimage.morphology.local_maxima`` and ``.local_minima`` was
28+
improved with a new Cython-based implementation. #3022
2529
- ``skivi`` is now using ``qtpy`` for Qt4/Qt5/PySide/PySide2 compatibility (a
2630
new optional dependency).
2731
- Performance is now monitored by

‎skimage/morphology/__init__.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,7 @@
1111
from .convex_hull import convex_hull_image, convex_hull_object
1212
from .greyreconstruct import reconstruction
1313
from .misc import remove_small_objects, remove_small_holes
14-
from .extrema import (h_minima, h_maxima, local_maxima, local_minima)
14+
from .extrema import h_minima, h_maxima, local_maxima, local_minima
1515

1616

1717
__all__ = ['binary_erosion',

‎skimage/morphology/_extrema_cy.pyx

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

‎skimage/morphology/_queue_with_history.pxi

Lines changed: 129 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,129 @@
1+
"""FIFO queue that preserves history until explicitly cleared.
2+
3+
.. warning::
4+
5+
One must define the type "QueueItem" before including this file. This makes
6+
it possible to store different types as QueueItems.
7+
8+
This queue can be operated like a class. The structure stores the state of the
9+
instance and the contained functions act as instance methods. The important
10+
distinction compared to normal queues is that popping an element doesn't remove
11+
it internally. Thus unless the queue's internal buffer is explicitly cleared
12+
(see `queue_clear`) popped items can be restored at anytime using
13+
`queue_restore`.
14+
15+
A possible application of this special functionality might be when actions
16+
performed on popped items must be undone at a later stage.
17+
18+
Example
19+
-------
20+
::
21+
22+
cdef QueueWithHistory queue
23+
cdef QueueItem item
24+
25+
queue_init(&queue, 5)
26+
try:
27+
for i in range(10):
28+
item = <QueueItem>i
29+
queue_push(&queue, &item)
30+
31+
while queue_pop(&queue, &item):
32+
print(item) # Prints 0 to 9
33+
34+
queue_restore(&queue)
35+
36+
while queue_pop(&queue, &item):
37+
print(item) # Prints 0 to 9 again
38+
finally:
39+
queue_exit(&queue)
40+
"""
41+
42+
from libc.stdlib cimport malloc, realloc, free
43+
44+
45+
# Store state of queue
46+
cdef struct QueueWithHistory:
47+
QueueItem* _buffer_ptr
48+
Py_ssize_t _buffer_size # Maximal number of elements the buffer can store
49+
Py_ssize_t _index_valid # Index to most recently inserted item
50+
Py_ssize_t _index_consumed # Index to most recently consumed item
51+
52+
53+
cdef inline void queue_init(QueueWithHistory* self, Py_ssize_t buffer_size) nogil:
54+
"""Initialize the queue and its buffer size.
55+
56+
The size is defined as the number of queue items to fit into the initial
57+
buffer, thus its true memory size is `buffer_size * sizeof(QueueItem)`. Be
58+
sure to call `queue_exit` after calling this function to free the allocated
59+
memory!
60+
"""
61+
self._buffer_ptr = <QueueItem*>malloc(buffer_size * sizeof(QueueItem))
62+
if not self._buffer_ptr:
63+
with gil:
64+
raise MemoryError("couldn't allocate buffer")
65+
self._buffer_size = buffer_size
66+
self._index_consumed = -1
67+
self._index_valid = -1
68+
69+
70+
cdef inline void queue_push(QueueWithHistory* self, QueueItem* item_ptr) nogil:
71+
"""Enqueue a new item."""
72+
self._index_valid += 1
73+
if self._buffer_size <= self._index_valid:
74+
_queue_grow_buffer(self)
75+
self._buffer_ptr[self._index_valid] = item_ptr[0]
76+
77+
78+
cdef inline unsigned char queue_pop(QueueWithHistory* self,
79+
QueueItem* item_ptr) nogil:
80+
"""If not empty pop an item and return 1 otherwise return 0.
81+
82+
The item is still preserved in the internal buffer and can be restored with
83+
`queue_restore`. To truly clear the internal buffer use `queue_clear`.
84+
"""
85+
if 0 <= self._index_consumed + 1 <= self._index_valid:
86+
self._index_consumed += 1
87+
item_ptr[0] = self._buffer_ptr[self._index_consumed]
88+
return 1
89+
return 0
90+
91+
92+
cdef inline void queue_restore(QueueWithHistory* self) nogil:
93+
"""Restore all consumed items to the queue.
94+
95+
The order of the restored queue is the same as previous one, meaning older
96+
items are popped first.
97+
"""
98+
self._index_consumed = -1
99+
100+
101+
cdef inline void queue_clear(QueueWithHistory* self) nogil:#
102+
"""Remove all consumable items.
103+
104+
After this the old items can't be restored with `queue_restore`.
105+
"""
106+
self._index_consumed = -1
107+
self._index_valid = -1
108+
109+
110+
cdef inline void queue_exit(QueueWithHistory* self) nogil:
111+
"""Free the buffer of the queue.
112+
113+
Don't use the queue after this command unless `queue_init` is called again.
114+
"""
115+
free(self._buffer_ptr)
116+
117+
118+
cdef inline void _queue_grow_buffer(QueueWithHistory* self) nogil:
119+
"""Double the memory used for the buffer."""
120+
cdef QueueItem* new_buffer
121+
self._buffer_size *= 2
122+
new_buffer_ptr = <QueueItem*>realloc(
123+
self._buffer_ptr,
124+
self._buffer_size * sizeof(QueueItem)
125+
)
126+
if not new_buffer_ptr:
127+
with gil:
128+
raise MemoryError("couldn't reallocate buffer")
129+
self._buffer_ptr = new_buffer_ptr

0 commit comments

Comments
 (0)
0