Add function to remove near objects in nd-image

benchmarks/benchmark_morphology.py

@@ -7,7 +7,7 @@
 from numpy.lib import NumpyVersion as Version
 
 import skimage
-from skimage import data, morphology, util
+from skimage import data, morphology, util, color
 
 
 class Skeletonize3d(object):
@@ -157,3 +157,42 @@ def time_erosion(
         self, shape, footprint, radius, *args
     ):
         morphology.erosion(self.image, self.footprint)
+
+
+class RemoveNearObjects(object):
+
+    param_names = ["minimal_distance"]
+    params = [5, 100]
+
+    def setup(self, *args):
+        image = data.hubble_deep_field()
+        image = color.rgb2gray(image)
+        self.objects = image > 0.18  # Chosen with threshold_li
+
+    def time_remove_near_objects(self, minimal_distance):
+        morphology.remove_near_objects(
+            self.objects, minimal_distance=minimal_distance
+        )
+
+    def peakmem_reference(self, *args):
+        """Provide reference for memory measurement with empty benchmark.
+
+        Peakmem benchmarks measure the maximum amount of RAM used by a
+        function. However, this maximum also includes the memory used
+        during the setup routine (as of asv 0.2.1; see [1]_).
+        Measuring an empty peakmem function might allow us to disambiguate
+        between the memory used by setup and the memory used by target (see
+        other ``peakmem_`` functions below).
+
+        References
+        ----------
+        .. [1]: https://asv.readthedocs.io/en/stable/writing_benchmarks.html#peak-memory
+        """
+        pass
+
+    def peakmem_remove_near_objects(self, minimal_distance):
+        morphology.remove_near_objects(
+            self.objects,
+            minimal_distance=minimal_distance,
+            priority=self.objects,
+        )

doc/examples/features_detection/plot_remove_objects.py

@@ -0,0 +1,45 @@
+"""
+================
+Removing objects
+================
+
+scikit-image supports several ways to remove objects inside N-dimensional
+images. In this context "objects" (and "holes") are defined as groups of
+connected samples that are distinct from the background. A binary image can
+contain several objects that are not connected to each other.
+
+The code snippet below demonstrates two ways to remove objects inside an image
+either by removing objects
+
+- based on the number of samples that make up each object
+- or based on their distance of each object to each other.
+"""
+
+import matplotlib.pyplot as plt
+from skimage import data, morphology, color, filters
+
+# Extract foreground by thresholding an image taken by the Hubble Telescope
+image = color.rgb2gray(data.hubble_deep_field())
+foreground = image > filters.threshold_li(image)
+
+# Separate objects into regions larger and smaller than 100 pixels
+large_objects = morphology.remove_small_objects(foreground, min_size=100)
+small_objects = foreground ^ large_objects
+
+# Remove objects until remaining ones are at least 100 pixels apart,
+# smaller ones are removed in favor of larger ones by default
+spaced_objects = morphology.remove_near_objects(
+    foreground, minimal_distance=100
+)
+
+# Plot the results
+fig, ax = plt.subplots(2, 2, figsize=(10, 10))
+ax[0, 0].set_title("original")
+ax[0, 0].imshow(foreground)
+ax[0, 1].set_title("large objects")
+ax[0, 1].imshow(large_objects)
+ax[1, 1].set_title("small objects")
+ax[1, 1].imshow(small_objects)
+ax[1, 0].set_title("spaced objects (nearest removed)")
+ax[1, 0].imshow(spaced_objects)
+plt.show()

doc/release/release_dev.rst

@@ -15,6 +15,8 @@ https://scikit-image.org
 
 New Features
 ------------
+- Added ``skimage.morphology.remove_near_objects`` to remove objects until a
+  minimal distance between remaining ones is ensured.
 
 
 

skimage/morphology/__init__.py

@@ -8,7 +8,8 @@
 from ._skeletonize import skeletonize, medial_axis, thin, skeletonize_3d
 from .convex_hull import convex_hull_image, convex_hull_object
 from .grayreconstruct import reconstruction
-from .misc import remove_small_objects, remove_small_holes
+from .misc import (remove_small_objects, remove_small_holes,
+                   remove_near_objects)
 from .extrema import h_minima, h_maxima, local_maxima, local_minima
 from ._flood_fill import flood, flood_fill
 from .max_tree import (max_tree, area_opening, area_closing,
@@ -44,6 +45,7 @@
            'convex_hull_object',
            'reconstruction',
            'remove_small_objects',
+           'remove_near_objects',
            'remove_small_holes',
            'h_minima',
            'h_maxima',

skimage/morphology/_near_objects_cy.pyx

@@ -0,0 +1,156 @@
+#cython: cdivision=True
+#cython: boundscheck=False
+#cython: nonecheck=False
+#cython: wraparound=False
+
+
+"""Cython code used in `remove_near_objects` function."""
+
+
+import numpy as np
+cimport numpy as cnp
+
+from _shared.fused_numerics cimport np_real_numeric
+
+# Must be defined to use QueueWithHistory
+ctypedef Py_ssize_t QueueItem
+
+
+include "_queue_with_history.pxi"
+
+
+def _remove_near_objects(
+    np_real_numeric[::1] image not None,
+    Py_ssize_t[::1] labels not None,
+    Py_ssize_t[::1] raveled_indices not None,
+    Py_ssize_t[::1] neighbor_offsets not None,
+    kdtree,
+    cnp.float64_t p_norm,
+    cnp.float64_t minimal_distance,
+    tuple shape,
+):
+    """Remove objects until a minimal distance is ensured.
+
+    Iterates over all objects (connected pixels that are True) inside an image
+    and removes neighboring objects until all remaining ones are at least a
+    minimal distance from each other.
+
+    Parameters
+    ----------
+    image :
+        The raveled view of a n-dimensional array. which is modified inplace.
+    labels :
+        An array with labels for each object in `image` matching it in shape.
+    raveled_indices :
+        Indices into `image` and `labels` that determines the iteration order
+        and thus which objects take precedence.
+    neighbor_offsets :
+        A one-dimensional array that contains the offsets to find the
+        connected neighbors for any index in `image`.
+    kdtree : scipy.spatial.cKDTree
+        A KDTree containing the coordinates of all objects in `image`.
+    minimal_distance :
+        The minimal allowed distance between objects.
+    p_norm :
+        Which Minkowski p-norm to use to calculate the distance between
+        objects. Defaults to 2 which corresponds to the Euclidean distance
+        while 1 corresponds to the Manatten distance.
+    shape :
+        The shape of the unraveled `image`.
+
+    Notes
+    -----
+    This function and its partner function :func:`~._remove_object` can deal
+    with objects where `labels` is 0 inside objects as long as its enclosing
+    surface points (in the sense of the neighborhood) are labeled.
+    This significantly improves the performance by reducing number of queries
+    to the KDTree and its size.
+    This effect grows with the size to surface ratio of all evaluated objects.
+    """
+    cdef:
+        Py_ssize_t i, j, index_i, index_j
+        list in_range
+        QueueWithHistory queue
+
+    queue_init(&queue, 64)
+    try:
+        for i in range(raveled_indices.shape[0]):
+            index_i = raveled_indices[i]
+
+            # Skip if point is part of a removed object
+            if image[index_i] == 0:
+                continue
+
+            in_range = kdtree.query_ball_point(
+                np.unravel_index(index_i, shape),
+                r=minimal_distance,
+                p=p_norm,
+            )
+
+            # Remove objects in `in_range` that don't share the same label id
+            for j in in_range:
+                index_j = raveled_indices[j]
+                if (
+                    image[index_j] != 0
+                    and labels[index_i] != labels[index_j]
+                ):
+                    _remove_object(
+                        image=image,
+                        labels=labels,
+                        start_index=index_j,
+                        neighbor_offsets=neighbor_offsets,
+                        queue_ptr=&queue,
+                    )
+    finally:
+        queue_exit(&queue)
+
+
+cdef inline void _remove_object(
+    np_real_numeric[::1] image,
+    Py_ssize_t[::1] labels,
+    Py_ssize_t start_index,
+    Py_ssize_t[::1] neighbor_offsets,
+    QueueWithHistory* queue_ptr,
+):
+    """Remove single connected object.
+
+    Performs a flood-fill on the object with the value 0. Samples with a label
+    id == 0 and an image value != 0 are considered to be inside the evaluated
+    object. 
+
+    Parameters
+    ----------
+    image :
+        The raveled view of a n-dimensional array. which is modified inplace.
+    labels :
+        An array with labels for each object in `image` matching it in shape.
+    start_index :
+        Start position for the flood-fill.
+    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.
+    """
+    cdef Py_ssize_t i, point, neighbor, max_index
+    cdef cnp.uint32_t label
+
+    max_index = image.shape[0]
+    queue_clear(queue_ptr)
+    queue_push(queue_ptr, &start_index)
+    image[start_index] = 0
+    label = labels[start_index]
+
+    while queue_pop(queue_ptr, &point):
+        for i in range(neighbor_offsets.shape[0]):
+            neighbor = point + neighbor_offsets[i]
+            # Bounds checking because image wasn't padded to signal the edge
+            if not 0 <= neighbor < max_index:
+                continue
+
+            # The algorithm might cross the image edge when two objects are
+            # neighbors in the raveled view -> check that the label id is
+            # either the same (object's surface) or 0 (inside object).
+            if image[neighbor] != 0 and labels[neighbor] in (0, label):
+                queue_push(queue_ptr, &neighbor)
+                image[neighbor] = 0