scikit-image · stefanv · Jun 19, 2018 · Apr 22, 2018 · Apr 23, 2018 · Apr 23, 2018
diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt
@@ -256,3 +256,6 @@
 
 - David Volgyes
   Unsharp masking
+
+- Lars Grüter
+  Flood-fill based local maxima detection
diff --git a/doc/release/release_dev.rst b/doc/release/release_dev.rst
@@ -17,11 +17,15 @@ New Features
 ------------
 
 - unsharp mask filtering (#2772)
+- New options ``connectivity``, ``indices`` and ``allow_borders`` for
+  ``skimage.morphology.local_maxima`` and ``.local_minima``. #3022
 
 
 Improvements
 ------------
 
+- Performance of ``skimage.morphology.local_maxima`` and ``.local_minima`` was
+  improved with a new Cython-based implementation. #3022
 - ``skivi`` is now using ``qtpy`` for Qt4/Qt5/PySide/PySide2 compatibility (a
   new optional dependency).
 - Performance is now monitored by

diff --git a/skimage/morphology/__init__.py b/skimage/morphology/__init__.py
@@ -11,7 +11,7 @@
 from .convex_hull import convex_hull_image, convex_hull_object
 from .greyreconstruct import reconstruction
 from .misc import remove_small_objects, remove_small_holes
-from .extrema import (h_minima, h_maxima, local_maxima, local_minima)
+from .extrema import h_minima, h_maxima, local_maxima, local_minima
 
 
 __all__ = ['binary_erosion',

diff --git a/skimage/morphology/_extrema_cy.pyx b/skimage/morphology/_extrema_cy.pyx
@@ -0,0 +1,224 @@
+#cython: cdivision=True
+#cython: boundscheck=False
+#cython: nonecheck=False
+#cython: wraparound=False
+
+"""Cython code used in extrema.py."""
+
+cimport numpy as cnp
+
+
+# Must be defined to use QueueWithHistory
+ctypedef Py_ssize_t QueueItem
+
+include "_queue_with_history.pxi"
+
+
+ctypedef fused dtype_t:
+    cnp.uint8_t
+    cnp.uint16_t
+    cnp.uint32_t
+    cnp.uint64_t
+    cnp.int8_t
+    cnp.int16_t
+    cnp.int32_t
+    cnp.int64_t
+    cnp.float32_t
+    cnp.float64_t
+
+
+# Definition of flag values used for `flags` in _local_maxima & _fill_plateau
+cdef:
+    # First or last index in a dimension
+    unsigned char BORDER_INDEX = 3
+    # Potentially part of a maximum
+    unsigned char CANDIDATE = 2
+    # Index was queued (flood-fill) and might still be part of maximum OR
+    # when evaluation is complete this flag value marks definite local maxima
+    unsigned char QUEUED_CANDIDATE = 1
+    # None of the above is true
+    unsigned char NOT_MAXIMUM = 0
+
+
+def _local_maxima(dtype_t[::1] image not None,
+                  unsigned char[::1] flags,
+                  Py_ssize_t[::1] neighbor_offsets not None):
+    """Detect local maxima in n-dimensional array.
+
+    Inner function to `local_maxima` that detects all local maxima (including
+    plateaus) in the image. The result is stored inplace inside `flags` with
+    the value of "QUEUED_CANDIDATE".
+
+    Parameters
+    ----------
+    image : ndarray, one-dimensional
+        The raveled view of a n-dimensional array.
+    flags : ndarray
+        An array of flags that is used to store the state of each pixel during
+        evaluation and is MODIFIED INPLACE. Initially, pixels that border the
+        image edge must be marked as "BORDER_INDEX" while all other pixels
+        should be marked with "NOT_MAXIMUM".
+    neighbor_offsets : ndarray
+        A one-dimensional array that contains the offsets to find the
+        connected neighbors for any index in `image`.
+    """
+    cdef:
+        QueueWithHistory queue
+        unsigned char last_dim_in_neighbors
+
+    # This needs the GIL
+    last_dim_in_neighbors = -1 in neighbor_offsets and 1 in neighbor_offsets
+    with nogil:
+        if last_dim_in_neighbors:
+            # If adjacent pixels in the last dimension are part of the
+            # neighborhood, the number of candidates which have to be evaluated
+            # in the second algorithmic step `_fill_plateau` can be reduced
+            # ahead of time with this function
+            _mark_candidates_in_last_dimension(image, flags)
+        else:
+            # Otherwise simply mark all pixels (except border) as candidates
+            _mark_candidates_all(flags)
+
+        # Initialize a buffer used to queue positions while evaluating each
+        # potential maximum (flagged with 2)
+        queue_init(&queue, 64)
+        try:
+            for i in range(image.shape[0]):
+                if flags[i] == CANDIDATE:
+                    # Index is potentially part of a maximum:
+                    # Find all samples which are part of the plateau and fill
+                    # with 0 or 1 depending on whether it's a true maximum
+                    _fill_plateau(image, flags, neighbor_offsets, &queue, i)
+        finally:
+            # Ensure that memory is released again
+            queue_exit(&queue)
+
+
+cdef inline void _mark_candidates_in_last_dimension(
+        dtype_t[::1] image, unsigned char[::1] flags) nogil:
+    """Mark local maxima in last dimension.
+
+    This function considers only the last dimension of the image and marks 
+    pixels with the "CANDIDATE" flag if it is a local maximum. 
+
+    Parameters
+    ----------
+    image :
        The raveled view of a n-dimensional array.
+    flags :
+        An array of flags that is used to store the state of each pixel during
+        evaluation.
+
+    Notes
+    -----
+    By evaluating this necessary but not sufficient condition first, usually a
+    significant amount of pixels can be rejected without having to evaluate the
+    entire neighborhood of their plateau. This can reduces the number of 
+    candidates that need to be evaluated with the more expensive flood-fill 
+    performed in `_fill_plateaus`.
+
+    However this is only possible if the adjacent pixels in the last dimension
+    are part of the defined neighborhood (see argument `neighbor_offsets`
+    in `_local_maxima`).
+    """
+    cdef Py_ssize_t i, i_ahead
+    i = 1
+    while i < image.shape[0]:
+        if image[i - 1] < image[i] and flags[i] != BORDER_INDEX:
+            # Potential maximum (in last dimension) is found, find
+            # other edge of current plateau or "edge of dimension"
+            i_ahead = i + 1
+            while (
+                image[i] == image[i_ahead] and
+                flags[i_ahead] != BORDER_INDEX
+            ):
+                i_ahead += 1
+            if image[i] > image[i_ahead]:
+                # Found local maximum (in one dimension), mark all
+                # parts of the plateau as potential maximum
+                flags[i:i_ahead] = CANDIDATE
+            i = i_ahead
+        else:
+            i += 1
+
+
+cdef inline void _mark_candidates_all(unsigned char[::1] flags) nogil:
+    """Mark all pixels as potential maxima, exclude border pixels.
+
+    This function marks pixels with the "CANDIDATE" flag if they aren't the
+    first or last index in any dimension (not flagged with "BORDER_INDEX"). 
+
+    Parameters
+    ----------
+    flags :
+        An array of flags that is used to store the state of each pixel during
+        evaluation.
+    """
+    cdef Py_ssize_t i = 1
+    while i < flags.shape[0]:
+        if flags[i] != BORDER_INDEX:
+            flags[i] = CANDIDATE
+        i += 1
+
+
+cdef inline void _fill_plateau(
+        dtype_t[::1] image, unsigned char[::1] flags,
+        Py_ssize_t[::1] neighbor_offsets, QueueWithHistory* queue_ptr,
+        Py_ssize_t start_index) nogil:
+    """Fill with 1 if plateau is local maximum else with 0.
+
+    Parameters
+    ----------
+    image :
+        The raveled view of a n-dimensional array.
+    flags :
+        An array of flags that is used to store the state of each pixel during
+        evaluation.
+    neighbor_offsets :
+        A one-dimensional array that contains the offsets to find the
+        connected neighbors for any index in `image`.
+    queue_ptr :
+        Pointer to initialized queue.
+    start_index :
+        Start position for the flood-fill.
+    """
+    cdef:
+        dtype_t height
+        unsigned char is_maximum
+        QueueItem current_index, neighbor
+
+    height = image[start_index]  # Height of the evaluated plateau
+    is_maximum = 1 # Boolean flag, initially assume true
+
+    # Queue start position after clearing the buffer
+    # which might have been used already
+    queue_clear(queue_ptr)
+    queue_push(queue_ptr, &start_index)
+    flags[start_index] = QUEUED_CANDIDATE
+
+    # Break loop if all queued positions were evaluated
+    while queue_pop(queue_ptr, &current_index):
+        # Look at all neighboring samples
+        for i in range(neighbor_offsets.shape[0]):
+            neighbor = current_index + neighbor_offsets[i]
+
+            if image[neighbor] == height:
+                # Value is part of plateau
+                if flags[neighbor] == BORDER_INDEX:
+                    # Plateau touches border and can't be maximum
+                    is_maximum = 0
+                elif flags[neighbor] != QUEUED_CANDIDATE:
+                    # Index wasn't queued already, do so now
+                    queue_push(queue_ptr, &neighbor)
+                    flags[neighbor] = QUEUED_CANDIDATE
+
+            elif image[neighbor] > height:
+                # Current plateau can't be maximum because it borders a
+                # larger one
+                is_maximum = 0
+
+    if not is_maximum:
+        queue_restore(queue_ptr)
+        # Initial guess was wrong -> flag as NOT_MAXIMUM
+        while queue_pop(queue_ptr, &neighbor):
+            flags[neighbor] = NOT_MAXIMUM
diff --git a/skimage/morphology/_queue_with_history.pxi b/skimage/morphology/_queue_with_history.pxi
@@ -0,0 +1,129 @@
+"""FIFO queue that preserves history until explicitly cleared.
+
+.. warning::
+
+    One must define the type "QueueItem" before including this file. This makes
+    it possible to store different types as QueueItems.
+
+This queue can be operated like a class. The structure stores the state of the
+instance and the contained functions act as instance methods. The important
+distinction compared to normal queues is that popping an element doesn't remove
+it internally. Thus unless the queue's internal buffer is explicitly cleared
+(see `queue_clear`) popped items can be restored at anytime using
+`queue_restore`.
+
+A possible application of this special functionality might be when actions
+performed on popped items must be undone at a later stage.
+
+Example
+-------
+::
+
+    cdef QueueWithHistory queue
+    cdef QueueItem item
+
+    queue_init(&queue, 5)
+    try:
+        for i in range(10):
+            item = <QueueItem>i
+            queue_push(&queue, &item)
+
+        while queue_pop(&queue, &item):
+            print(item)  # Prints 0 to 9
+
+        queue_restore(&queue)
+
+        while queue_pop(&queue, &item):
+            print(item)  # Prints 0 to 9 again
+    finally:
+        queue_exit(&queue)
+"""
+
+from libc.stdlib cimport malloc, realloc, free
+
+
+# Store state of queue
+cdef struct QueueWithHistory:
+    QueueItem* _buffer_ptr
+    Py_ssize_t _buffer_size  # Maximal number of elements the buffer can store
+    Py_ssize_t _index_valid  # Index to most recently inserted item
+    Py_ssize_t _index_consumed  # Index to most recently consumed item
+
+
+cdef inline void queue_init(QueueWithHistory* self, Py_ssize_t buffer_size) nogil:
+    """Initialize the queue and its buffer size.
+
+    The size is defined as the number of queue items to fit into the initial 
+    buffer, thus its true memory size is `buffer_size * sizeof(QueueItem)`. Be
+    sure to call `queue_exit` after calling this function to free the allocated
+    memory!
+    """
+    self._buffer_ptr = <QueueItem*>malloc(buffer_size * sizeof(QueueItem))
+    if not self._buffer_ptr:
+        with gil:
+            raise MemoryError("couldn't allocate buffer")
+    self._buffer_size = buffer_size
+    self._index_consumed = -1
+    self._index_valid = -1
+
+
+cdef inline void queue_push(QueueWithHistory* self, QueueItem* item_ptr) nogil:
+    """Enqueue a new item."""
+    self._index_valid += 1
+    if self._buffer_size <= self._index_valid:
+        _queue_grow_buffer(self)
+    self._buffer_ptr[self._index_valid] = item_ptr[0]
+
+
+cdef inline unsigned char queue_pop(QueueWithHistory* self,
+                                    QueueItem* item_ptr) nogil:
+    """If not empty pop an item and return 1 otherwise return 0.
+
+    The item is still preserved in the internal buffer and can be restored with
+    `queue_restore`. To truly clear the internal buffer use `queue_clear`.
+    """
+    if 0 <= self._index_consumed + 1 <= self._index_valid:
+        self._index_consumed += 1
+        item_ptr[0] = self._buffer_ptr[self._index_consumed]
+        return 1
+    return 0
+
+
+cdef inline void queue_restore(QueueWithHistory* self) nogil:
+    """Restore all consumed items to the queue.
+
+    The order of the restored queue is the same as previous one, meaning older
+    items are popped first.
+    """
+    self._index_consumed = -1
+
+
+cdef inline void queue_clear(QueueWithHistory* self) nogil:#
+    """Remove all consumable items.
+
+    After this the old items can't be restored with `queue_restore`.
+    """
+    self._index_consumed = -1
+    self._index_valid = -1
+
+
+cdef inline void queue_exit(QueueWithHistory* self) nogil:
+    """Free the buffer of the queue.
+
+    Don't use the queue after this command unless `queue_init` is called again.
+    """
+    free(self._buffer_ptr)
+
+
+cdef inline void _queue_grow_buffer(QueueWithHistory* self) nogil:
+    """Double the memory used for the buffer."""
+    cdef QueueItem* new_buffer
+    self._buffer_size *= 2
+    new_buffer_ptr = <QueueItem*>realloc(
+        self._buffer_ptr,
+        self._buffer_size * sizeof(QueueItem)
+    )
+    if not new_buffer_ptr:
+        with gil:
+            raise MemoryError("couldn't reallocate buffer")
+    self._buffer_ptr = new_buffer_ptr