Refactor and fix peak_local_max

doc/release/release_dev.rst

@@ -73,7 +73,7 @@ Improvements
 - ``measure.label`` has been accelerated for boolean input images, by using
   ``scipy.ndimage``'s implementation for this case (#4945).
 - ``util.apply_parallel`` now works with multichannel data (#4927).
-
+- ``skimage.feature.peak_local_max`` supports now any Minkowski distance.
 
 
 API Changes
@@ -93,6 +93,8 @@ API Changes
   descriptive name.
 - The ``level`` parameter of ``measure.find_contours`` is now a keyword
   argument, with a default value set to (max(image) - min(image)) / 2.
+- ``p_norm`` argument was added to ``skimage.feature.peak_local_max``
+  to add support for Minkowski distances.
 
 
 Bugfixes
@@ -122,6 +124,12 @@ Bugfixes
   NaNs are found in the computation (as a result of NaNs in input images).
   Before this fix, an incorrect value could be returned where the input images
   had NaNs (#4886).
+- ``min_distance`` is now enforced for ``skimage.feature.peak_local_max``
+  (#2592).
+- Peak detection in labels is fixed in ``skimage.feature.peak_local_max``
+  (#4756).
+- Input ``labels`` argument renumbering in ``skimage.feature.peak_local_max``
+  is avoided (#5047).
 
 
 Deprecations

skimage/_shared/coord.py

@@ -0,0 +1,51 @@
+import numpy as np
+from scipy.spatial import cKDTree, distance
+
+
+def ensure_spacing(coord, spacing=1, p_norm=np.inf):
+    """Returns a subset of coord where a minimum spacing is guaranteed.
+
+    Parameters
+    ----------
+    coord : ndarray
+        The coordinates of the considered points.
+    spacing : float
+        the maximum allowed spacing between the points.
+    p_norm : float
+        Which Minkowski p-norm to use. Should be in the range [1, inf].
+        A finite large p may cause a ValueError if overflow can occur.
+        ``inf`` corresponds to the Chebyshev distance and 2 to the
+        Euclidean distance.
+
+    Returns
+    -------
+    output : ndarray
+        A subset of coord where a minimum spacing is guaranteed.
+
+    """
+
+    output = coord
+    if len(coord):
+        # Use KDtree to find the peaks that are too close to each other
+        tree = cKDTree(coord)
+
+        indices = tree.query_ball_point(coord, r=spacing, p=p_norm)
+        rejected_peaks_indices = set()
+        for idx, candidates in enumerate(indices):
+            if idx not in rejected_peaks_indices:
+                # keep current point and the points at exactly spacing from it
+                candidates.remove(idx)
+                dist = distance.cdist([coord[idx]],
+                                      coord[candidates],
+                                      distance.minkowski,
+                                      p=p_norm).reshape(-1)
+                candidates = [c for c, d in zip(candidates, dist)
+                              if d < spacing]
+
+                # candidates.remove(keep)
+                rejected_peaks_indices.update(candidates)
+
+        # Remove the peaks that are too close to each other
+        output = np.delete(coord, tuple(rejected_peaks_indices), axis=0)
+
+    return output

skimage/_shared/tests/test_coord.py

@@ -0,0 +1,43 @@
+import pytest
+import numpy as np
+from scipy.spatial.distance import pdist, minkowski
+from skimage._shared.coord import ensure_spacing
+
+
+@pytest.mark.parametrize("p", [1, 2, np.inf])
+def test_ensure_spacing_trivial(p):
+    # --- Empty input
+    assert ensure_spacing([], p_norm=p) == []
+
+    # --- Verified spacing
+    coord = np.random.randn(100, 2)
+
+    # --- 0 spacing
+    assert np.array_equal(coord, ensure_spacing(coord, spacing=0, p_norm=p))
+
+    # Spacing is chosen to be half the minimum distance
+    spacing = pdist(coord, minkowski, p=p).min() * 0.5
+
+    out = ensure_spacing(coord, spacing=spacing, p_norm=p)
+
+    assert np.array_equal(coord, out)
+
+
+@pytest.mark.parametrize("ndim", [1, 2, 3, 4, 5])
+def test_ensure_spacing_nD(ndim):
+    coord = np.ones((5, ndim))
+
+    expected = np.ones((1, ndim))
+
+    assert np.array_equal(ensure_spacing(coord), expected)
+
+
+@pytest.mark.parametrize("p", [1, 2, np.inf])
+def test_ensure_spacing_p_norm(p):
+    coord = np.random.randn(100, 2)
+
+    # --- Consider the average distance btween the point as spacing
+    spacing = np.median(pdist(coord, minkowski, p))
+    out = ensure_spacing(coord, spacing=spacing, p_norm=p)
+
+    assert pdist(out, minkowski, p).min() > spacing