diff --git a/MIGRATION_GUIDE.md b/MIGRATION_GUIDE.md
index 72b1286..6c12a90 100644
--- a/MIGRATION_GUIDE.md
+++ b/MIGRATION_GUIDE.md
@@ -2047,6 +2047,53 @@ Important differences:
   axis-0 coordinate `-1` and `0` on all other axes. The first row of
   `indices` will then contain `-1` everywhere.
 
+### Non-Maximum Distance Suppression
+
+`nifty.filters.nonMaximumDistanceSuppression` filters a set of candidate
+points using a distance map: each point's suppression radius is the distance
+value at its own location, and from every group of points that fall within
+one another's radius only the one with the largest distance value is kept.
+`bioimage-cpp` exposes the same algorithm as
+`bic.distance.non_maximum_distance_suppression`.
+
+nifty:
+
+```python
+from nifty.filters import nonMaximumDistanceSuppression
+
+# distanceMap: float32 array; points: uint64 array of shape (N, ndim)
+kept = nonMaximumDistanceSuppression(distanceMap, points)
+```
+
+bioimage-cpp:
+
+```python
+import bioimage_cpp as bic
+
+kept = bic.distance.non_maximum_distance_suppression(distance_map, points)
+```
+
+Name mapping:
+
+| nifty name | bioimage-cpp name |
+| --- | --- |
+| `nifty.filters.nonMaximumDistanceSuppression` | `non_maximum_distance_suppression` |
+
+Important differences:
+
+- Snake_case naming, consistent with the rest of `bic.distance`.
+- `points` may be `int64`, `uint64`, `int32`, or `uint32`; the returned array
+  has shape `(K, ndim)` and preserves the input `points` dtype (nifty always
+  returned `uint64`). Output rows are the retained points in ascending
+  input-index order.
+- `distance_map` is coerced to C-contiguous `float32` if needed. The
+  per-point radius is dynamic (the distance value at each point), matching
+  nifty; there is no fixed-radius mode.
+- The algorithm is otherwise identical to nifty, including its float
+  arithmetic, so results match element-for-element. It uses an O(N²)
+  pairwise distance matrix; threshold the distance map first to keep the
+  candidate count modest.
+
 ## I/O and Build Dependencies
 
 `bioimage-cpp` intentionally does not replace nifty or affogato I/O helpers.
diff --git a/development/distance/check_non_maximum_distance_suppression.py b/development/distance/check_non_maximum_distance_suppression.py
new file mode 100644
index 0000000..0bfc98d
--- /dev/null
+++ b/development/distance/check_non_maximum_distance_suppression.py
@@ -0,0 +1,140 @@
+"""Cross-check bioimage-cpp's non_maximum_distance_suppression against nifty.
+
+Builds random binary masks, computes their Euclidean distance transform, picks
+candidate points by thresholding the distance map, and compares
+``bic.distance.non_maximum_distance_suppression`` against
+``nifty.filters.nonMaximumDistanceSuppression`` for 2D and 3D inputs. Reports
+both correctness (set + row order) and per-call runtime.
+
+Not part of the pytest suite; requires nifty and scipy.
+"""
+
+from __future__ import annotations
+
+import argparse
+import sys
+from statistics import median
+from time import perf_counter
+
+import numpy as np
+
+import bioimage_cpp as bic
+
+try:
+    from nifty.filters import nonMaximumDistanceSuppression
+except ImportError as error:  # pragma: no cover - dev script
+    sys.stderr.write(f"nifty not installed: {error}\n")
+    sys.exit(1)
+
+try:
+    from scipy.ndimage import distance_transform_edt
+except ImportError as error:  # pragma: no cover - dev script
+    sys.stderr.write(f"scipy not installed: {error}\n")
+    sys.exit(1)
+
+
+CASES = [
+    # (name, shape, foreground_fraction, threshold)
+    ("2d_small", (60, 60), 0.85, 2.0),
+    ("2d_large", (256, 256), 0.9, 3.0),
+    ("3d_small", (25, 25, 25), 0.85, 2.0),
+    ("3d_large", (40, 40, 40), 0.9, 3.0),
+]
+
+
+def time_call(fn, repeats):
+    timings = []
+    result = None
+    for _ in range(repeats):
+        start = perf_counter()
+        result = fn()
+        timings.append(perf_counter() - start)
+    return median(timings), result
+
+
+def run_case(name, shape, fg_fraction, threshold, n_trials, repeats, rng):
+    rows = []
+    for trial in range(n_trials):
+        mask = rng.random(shape) < fg_fraction
+        dm = distance_transform_edt(mask).astype(np.float32)
+        coords = np.argwhere(dm > threshold).astype(np.uint64)
+        if len(coords) == 0:
+            continue
+
+        ref_s, ref = time_call(
+            lambda: nonMaximumDistanceSuppression(dm, coords), repeats
+        )
+        ours_s, ours = time_call(
+            lambda: bic.distance.non_maximum_distance_suppression(dm, coords), repeats
+        )
+
+        exact = ref.shape == ours.shape and np.array_equal(ref, ours)
+        same_set = {tuple(r) for r in ref.tolist()} == {tuple(r) for r in ours.tolist()}
+        rows.append(
+            {
+                "case": name,
+                "trial": trial,
+                "n_points": len(coords),
+                "n_ref": len(ref),
+                "n_ours": len(ours),
+                "set_ok": same_set,
+                "order_ok": exact,
+                "ref_s": ref_s,
+                "ours_s": ours_s,
+            }
+        )
+    return rows
+
+
+def main():
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--trials", type=int, default=5)
+    parser.add_argument("--repeats", type=int, default=3)
+    parser.add_argument("--seed", type=int, default=0)
+    args = parser.parse_args()
+
+    rng = np.random.default_rng(args.seed)
+
+    all_rows = []
+    for name, shape, fg, thr in CASES:
+        all_rows.extend(
+            run_case(name, shape, fg, thr, args.trials, args.repeats, rng)
+        )
+
+    header = (
+        f"{'case':>10} {'trial':>5} {'n_pts':>7} {'n_ref':>6} {'n_ours':>6}"
+        f" {'set':>5} {'order':>6} {'nifty_ms':>9} {'bic_ms':>9} {'speedup':>8}"
+    )
+    print(header)
+    print("-" * len(header))
+    all_ok = True
+    speedups = []
+    for r in all_rows:
+        speedup = r["ref_s"] / r["ours_s"] if r["ours_s"] > 0 else float("nan")
+        speedups.append(speedup)
+        print(
+            f"{r['case']:>10} {r['trial']:>5d} {r['n_points']:>7d}"
+            f" {r['n_ref']:>6d} {r['n_ours']:>6d}"
+            f" {'OK' if r['set_ok'] else 'FAIL':>5}"
+            f" {'OK' if r['order_ok'] else 'FAIL':>6}"
+            f" {r['ref_s'] * 1e3:>9.3f} {r['ours_s'] * 1e3:>9.3f}"
+            f" {speedup:>7.2f}x"
+        )
+        all_ok = all_ok and r["set_ok"] and r["order_ok"]
+
+    finite = [s for s in speedups if np.isfinite(s)]
+    if finite:
+        geo_mean = float(np.exp(np.mean(np.log(finite))))
+        print(
+            f"\nSpeedup (bic vs nifty): min {min(finite):.2f}x, "
+            f"max {max(finite):.2f}x, geo-mean {geo_mean:.2f}x"
+        )
+
+    if not all_ok:
+        print("\nFAIL: output mismatch vs nifty", file=sys.stderr)
+        sys.exit(1)
+    print("All cases match nifty (set and row order).")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/include/bioimage_cpp/non_maximum_distance_suppression.hxx b/include/bioimage_cpp/non_maximum_distance_suppression.hxx
new file mode 100644
index 0000000..1239e73
--- /dev/null
+++ b/include/bioimage_cpp/non_maximum_distance_suppression.hxx
@@ -0,0 +1,143 @@
+#pragma once
+
+#include "bioimage_cpp/array_view.hxx"
+#include "bioimage_cpp/detail/grid.hxx"
+
+#include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <cstdint>
+#include <limits>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
+namespace bioimage_cpp::distance {
+
+namespace detail {
+
+template <class PointT>
+inline std::ptrdiff_t point_to_flat(
+    const PointT *coord_row,
+    const std::vector<std::ptrdiff_t> &strides,
+    std::ptrdiff_t ndim
+) {
+    std::ptrdiff_t flat = 0;
+    for (std::ptrdiff_t d = 0; d < ndim; ++d) {
+        flat += static_cast<std::ptrdiff_t>(coord_row[d]) *
+                strides[static_cast<std::size_t>(d)];
+    }
+    return flat;
+}
+
+} // namespace detail
+
+// Non-maximum suppression of candidate points by a distance map.
+//
+// For each input point p_i, let d_i = distance_map[p_i]. Among all input
+// points (including i itself) within Euclidean distance d_i of p_i, the one
+// with the largest distance_map value is selected. The unique set of selected
+// indices is returned in ascending order via `kept_indices`.
+//
+// Matches `nifty.filters.nonMaximumDistanceSuppression`, including its float
+// arithmetic: coordinate differences and their squared sum accumulate in
+// float, the Euclidean distance is `float(sqrt(sum))`, and the neighborhood
+// test compares that distance directly against d_i. Replicating this exactly
+// (rather than comparing squared distances) keeps boundary ties identical to
+// nifty.
+//
+// Complexity: O(N^2) time and O(N^2) memory for the symmetric distance matrix.
+template <class PointT>
+inline void non_maximum_distance_suppression(
+    const ConstArrayView<float> &distance_map,
+    const ConstArrayView<PointT> &points,
+    std::vector<std::size_t> &kept_indices
+) {
+    if (distance_map.ndim() < 1) {
+        throw std::invalid_argument(
+            "distance_map must have ndim >= 1, got ndim=0"
+        );
+    }
+    if (points.ndim() != 2) {
+        throw std::invalid_argument(
+            "points must have ndim == 2, got ndim=" + std::to_string(points.ndim())
+        );
+    }
+    const auto n_points = points.shape[0];
+    const auto coord_ndim = points.shape[1];
+    if (coord_ndim != distance_map.ndim()) {
+        throw std::invalid_argument(
+            "points second axis must match distance_map ndim, got points.shape[1]=" +
+            std::to_string(coord_ndim) + ", distance_map.ndim()=" +
+            std::to_string(distance_map.ndim())
+        );
+    }
+
+    kept_indices.clear();
+    if (n_points == 0) {
+        return;
+    }
+
+    const auto strides = bioimage_cpp::detail::c_order_strides(distance_map.shape);
+    const auto n = static_cast<std::size_t>(n_points);
+
+    // Precompute flat index and distance value at each point.
+    std::vector<float> point_dist(n);
+    for (std::size_t i = 0; i < n; ++i) {
+        const auto *row =
+            points.data + static_cast<std::ptrdiff_t>(i) * coord_ndim;
+        const auto flat = detail::point_to_flat(row, strides, coord_ndim);
+        point_dist[i] = distance_map.data[flat];
+    }
+
+    // Pairwise Euclidean distance matrix (symmetric, N x N). The float
+    // accumulation and float(sqrt(...)) match nifty bit-for-bit so the
+    // neighborhood test below produces identical boundary decisions.
+    std::vector<float> pd(n * n, 0.0f);
+    for (std::size_t i = 0; i < n; ++i) {
+        const auto *row_i =
+            points.data + static_cast<std::ptrdiff_t>(i) * coord_ndim;
+        for (std::size_t j = i + 1; j < n; ++j) {
+            const auto *row_j =
+                points.data + static_cast<std::ptrdiff_t>(j) * coord_ndim;
+            float sum_sq = 0.0f;
+            for (std::ptrdiff_t d = 0; d < coord_ndim; ++d) {
+                const float diff = static_cast<float>(row_i[d]) -
+                                   static_cast<float>(row_j[d]);
+                sum_sq += diff * diff;
+            }
+            const auto val = static_cast<float>(std::sqrt(static_cast<double>(sum_sq)));
+            pd[i * n + j] = val;
+            pd[j * n + i] = val;
+        }
+    }
+
+    // For each point, scan all neighbors within its dynamic radius and keep
+    // the one with the largest distance_map value. Strict `>` ensures the
+    // first index encountered wins ties, matching nifty's behavior.
+    std::vector<std::size_t> bests;
+    bests.reserve(n);
+    for (std::size_t i = 0; i < n; ++i) {
+        const float d_i = point_dist[i];
+        float best_val = -std::numeric_limits<float>::infinity();
+        std::size_t best_idx = i;
+        const auto *pd_row = pd.data() + i * n;
+        for (std::size_t j = 0; j < n; ++j) {
+            if (pd_row[j] > d_i) {
+                continue;
+            }
+            const float dj = point_dist[j];
+            if (dj > best_val) {
+                best_val = dj;
+                best_idx = j;
+            }
+        }
+        bests.push_back(best_idx);
+    }
+
+    std::sort(bests.begin(), bests.end());
+    bests.erase(std::unique(bests.begin(), bests.end()), bests.end());
+    kept_indices = std::move(bests);
+}
+
+} // namespace bioimage_cpp::distance
diff --git a/src/bindings/distance.cxx b/src/bindings/distance.cxx
index 9ad9975..b05bdb6 100644
--- a/src/bindings/distance.cxx
+++ b/src/bindings/distance.cxx
@@ -2,6 +2,7 @@
 
 #include "bioimage_cpp/array_view.hxx"
 #include "bioimage_cpp/distance_transform.hxx"
+#include "bioimage_cpp/non_maximum_distance_suppression.hxx"
 
 #include <nanobind/nanobind.h>
 #include <nanobind/ndarray.h>
@@ -13,6 +14,7 @@
 #include <optional>
 #include <stdexcept>
 #include <string>
+#include <type_traits>
 #include <vector>
 
 namespace nb = nanobind;
@@ -187,6 +189,82 @@ nb::tuple distance_transform_uint8(
     return nb::make_tuple(distances_result, indices_result, vectors_result);
 }
 
+template <class PointT>
+nb::ndarray<nb::numpy, PointT, nb::c_contig> non_maximum_distance_suppression_impl(
+    nb::ndarray<nb::numpy, const float, nb::c_contig> distance_map,
+    nb::ndarray<nb::numpy, const PointT, nb::c_contig> points,
+    const std::size_t n_threads
+) {
+    (void)n_threads;  // Reserved for future parallelization; single-threaded.
+
+    if (distance_map.ndim() == 0) {
+        throw std::invalid_argument("distance_map must have ndim >= 1, got ndim=0");
+    }
+    if (points.ndim() != 2) {
+        throw std::invalid_argument(
+            "points must have ndim == 2, got ndim=" + std::to_string(points.ndim())
+        );
+    }
+    const auto coord_ndim = static_cast<std::size_t>(points.shape(1));
+    if (coord_ndim != distance_map.ndim()) {
+        throw std::invalid_argument(
+            "points.shape[1] must match distance_map ndim, got points.shape[1]=" +
+            std::to_string(coord_ndim) + ", distance_map.ndim()=" +
+            std::to_string(distance_map.ndim())
+        );
+    }
+
+    const auto map_shape = ndarray_shape(distance_map);
+    const auto n_points = static_cast<std::size_t>(points.shape(0));
+
+    // Bounds-check every coordinate before dropping the GIL.
+    const PointT *points_data = points.data();
+    for (std::size_t i = 0; i < n_points; ++i) {
+        for (std::size_t d = 0; d < coord_ndim; ++d) {
+            const PointT coord = points_data[i * coord_ndim + d];
+            if constexpr (std::is_signed_v<PointT>) {
+                if (coord < 0) {
+                    throw std::invalid_argument(
+                        "points coordinate out of bounds: points[" + std::to_string(i) +
+                        ", " + std::to_string(d) + "]=" + std::to_string(coord) +
+                        " is negative"
+                    );
+                }
+            }
+            if (static_cast<std::ptrdiff_t>(coord) >= map_shape[d]) {
+                throw std::invalid_argument(
+                    "points coordinate out of bounds: points[" + std::to_string(i) +
+                    ", " + std::to_string(d) + "]=" + std::to_string(coord) +
+                    " >= distance_map.shape[" + std::to_string(d) + "]=" +
+                    std::to_string(map_shape[d])
+                );
+            }
+        }
+    }
+
+    ConstArrayView<float> map_view{distance_map.data(), map_shape, {}};
+    ConstArrayView<PointT> points_view{points_data, ndarray_shape(points), {}};
+
+    std::vector<std::size_t> kept_indices;
+    {
+        nb::gil_scoped_release release;
+        distance::non_maximum_distance_suppression(map_view, points_view, kept_indices);
+    }
+
+    const std::size_t n_kept = kept_indices.size();
+    std::vector<std::size_t> out_shape{n_kept, coord_ndim};
+    auto output =
+        make_array<PointT, nb::ndarray<nb::numpy, PointT, nb::c_contig>>(out_shape);
+    PointT *out_data = output.data();
+    for (std::size_t k = 0; k < n_kept; ++k) {
+        const std::size_t i = kept_indices[k];
+        for (std::size_t d = 0; d < coord_ndim; ++d) {
+            out_data[k * coord_ndim + d] = points_data[i * coord_ndim + d];
+        }
+    }
+    return output;
+}
+
 } // namespace
 
 void bind_distance(nb::module_ &m) {
@@ -206,6 +284,33 @@ void bind_distance(nb::module_ &m) {
         "combination of (distances, indices, vectors) in a single separable F&H\n"
         "sweep. Pre-allocated output buffers are written into directly."
     );
+
+    const char *nms_doc =
+        "Non-maximum distance suppression of candidate points by a float32\n"
+        "distance map. For each point p_i, keeps the point with the largest\n"
+        "distance value within Euclidean distance distance_map[p_i] of p_i.\n"
+        "Returns the unique selected points (shape (K, ndim)) in ascending\n"
+        "input-index order. O(N^2) time and memory.";
+    m.def(
+        "_non_maximum_distance_suppression_int64",
+        &non_maximum_distance_suppression_impl<std::int64_t>,
+        nb::arg("distance_map"), nb::arg("points"), nb::arg("n_threads"), nms_doc
+    );
+    m.def(
+        "_non_maximum_distance_suppression_uint64",
+        &non_maximum_distance_suppression_impl<std::uint64_t>,
+        nb::arg("distance_map"), nb::arg("points"), nb::arg("n_threads"), nms_doc
+    );
+    m.def(
+        "_non_maximum_distance_suppression_int32",
+        &non_maximum_distance_suppression_impl<std::int32_t>,
+        nb::arg("distance_map"), nb::arg("points"), nb::arg("n_threads"), nms_doc
+    );
+    m.def(
+        "_non_maximum_distance_suppression_uint32",
+        &non_maximum_distance_suppression_impl<std::uint32_t>,
+        nb::arg("distance_map"), nb::arg("points"), nb::arg("n_threads"), nms_doc
+    );
 }
 
 } // namespace bioimage_cpp::bindings
diff --git a/src/bioimage_cpp/distance/__init__.py b/src/bioimage_cpp/distance/__init__.py
index ed7f020..ac0e836 100644
--- a/src/bioimage_cpp/distance/__init__.py
+++ b/src/bioimage_cpp/distance/__init__.py
@@ -1,8 +1,13 @@
 """Distance transforms."""
 
-from ._distance import distance_transform, vector_difference_transform
+from ._distance import (
+    distance_transform,
+    non_maximum_distance_suppression,
+    vector_difference_transform,
+)
 
 __all__ = [
     "distance_transform",
+    "non_maximum_distance_suppression",
     "vector_difference_transform",
 ]
diff --git a/src/bioimage_cpp/distance/_distance.py b/src/bioimage_cpp/distance/_distance.py
index 1890162..9a42afe 100644
--- a/src/bioimage_cpp/distance/_distance.py
+++ b/src/bioimage_cpp/distance/_distance.py
@@ -215,3 +215,80 @@ def vector_difference_transform(
         return_vectors=True,
         number_of_threads=number_of_threads,
     )
+
+
+_NMS_DISPATCH = {
+    np.dtype(np.int64): _core._non_maximum_distance_suppression_int64,
+    np.dtype(np.uint64): _core._non_maximum_distance_suppression_uint64,
+    np.dtype(np.int32): _core._non_maximum_distance_suppression_int32,
+    np.dtype(np.uint32): _core._non_maximum_distance_suppression_uint32,
+}
+
+
+def non_maximum_distance_suppression(
+    distance_map: np.ndarray,
+    points: np.ndarray,
+    number_of_threads: int = 1,
+) -> np.ndarray:
+    """Filter candidate points by non-maximum suppression on a distance map.
+
+    For each input point ``p_i`` with distance value ``d_i =
+    distance_map[p_i]``, keep the point with the largest ``distance_map``
+    value among all points within Euclidean distance ``d_i`` of ``p_i``
+    (including ``p_i`` itself). The unique set of such "dominant" points is
+    returned, ordered by ascending input index. This mirrors
+    ``nifty.filters.nonMaximumDistanceSuppression``.
+
+    Parameters
+    ----------
+    distance_map
+        Float array of any ndim ``D``. Coerced to C-contiguous ``float32`` if
+        a different float dtype or layout is supplied.
+    points
+        Integer array of shape ``(N, D)``; each row is a coordinate into
+        ``distance_map`` in NumPy axis order. Supported dtypes:
+        ``int64``, ``uint64``, ``int32``, ``uint32``.
+    number_of_threads
+        Reserved for future parallelization; currently single-threaded.
+
+    Returns
+    -------
+    np.ndarray
+        Filtered subset of ``points`` with shape ``(K, D)`` and the same
+        dtype as ``points``. ``K <= N``.
+
+    Notes
+    -----
+    Uses an ``O(N^2)`` pairwise distance matrix internally; suitable for ``N``
+    up to roughly 30k points. For larger candidate sets, threshold the
+    distance map more aggressively before calling.
+    """
+    function = "non_maximum_distance_suppression"
+
+    distance_map = np.ascontiguousarray(distance_map, dtype=np.float32)
+    if distance_map.ndim < 1:
+        raise ValueError(
+            f"{function}: distance_map must have ndim >= 1, got ndim={distance_map.ndim}"
+        )
+
+    points = np.ascontiguousarray(points)
+    if points.ndim != 2:
+        raise ValueError(
+            f"{function}: points must have ndim == 2, got ndim={points.ndim}"
+        )
+    if points.shape[1] != distance_map.ndim:
+        raise ValueError(
+            f"{function}: points.shape[1] must equal distance_map.ndim "
+            f"({distance_map.ndim}), got points.shape[1]={points.shape[1]}"
+        )
+
+    dispatch = _NMS_DISPATCH.get(points.dtype)
+    if dispatch is None:
+        supported = ", ".join(str(dt) for dt in ("int64", "uint64", "int32", "uint32"))
+        raise TypeError(
+            f"{function}: points must have one of dtypes [{supported}], "
+            f"got dtype={points.dtype}"
+        )
+
+    n_threads = _normalize_threads(number_of_threads, function)
+    return dispatch(distance_map, points, n_threads)
diff --git a/tests/distance/test_non_maximum_distance_suppression.py b/tests/distance/test_non_maximum_distance_suppression.py
new file mode 100644
index 0000000..a152064
--- /dev/null
+++ b/tests/distance/test_non_maximum_distance_suppression.py
@@ -0,0 +1,144 @@
+"""Tests for bioimage_cpp.distance.non_maximum_distance_suppression."""
+
+import numpy as np
+import pytest
+
+import bioimage_cpp as bic
+
+nms = bic.distance.non_maximum_distance_suppression
+
+
+def test_empty_points_returns_empty():
+    dm = np.ones((10, 10), dtype=np.float32)
+    out = nms(dm, np.zeros((0, 2), dtype=np.int64))
+    assert out.shape == (0, 2)
+    assert out.dtype == np.int64
+
+
+def test_single_point_returns_itself():
+    dm = np.zeros((11, 11), dtype=np.float32)
+    dm[5, 5] = 3.0
+    out = nms(dm, np.array([[5, 5]], dtype=np.int64))
+    assert out.tolist() == [[5, 5]]
+
+
+def test_two_close_points_keeps_higher_value():
+    # Both points sit within each other's dynamic neighborhood; only the
+    # one with the larger distance value survives.
+    dm = np.zeros((11, 11), dtype=np.float32)
+    dm[5, 5] = 5.0
+    dm[5, 6] = 4.0
+    pts = np.array([[5, 5], [5, 6]], dtype=np.int64)
+    out = nms(dm, pts)
+    assert out.tolist() == [[5, 5]]
+
+
+def test_two_far_points_both_survive():
+    dm = np.zeros((20, 20), dtype=np.float32)
+    dm[2, 2] = 1.0
+    dm[15, 15] = 1.0
+    pts = np.array([[2, 2], [15, 15]], dtype=np.int64)
+    out = nms(dm, pts)
+    # Far apart relative to their radius of 1.0 -> both kept, original order.
+    assert out.tolist() == [[2, 2], [15, 15]]
+
+
+def test_zero_radius_point_keeps_itself():
+    # A point whose distance value is 0 has an empty neighborhood except for
+    # itself, so it is always retained.
+    dm = np.zeros((10, 10), dtype=np.float32)
+    dm[3, 3] = 4.0  # high-value neighbor nearby
+    pts = np.array([[3, 4], [3, 3]], dtype=np.int64)  # first has value 0
+    out = nms(dm, pts)
+    out_set = {tuple(row) for row in out.tolist()}
+    assert (3, 4) in out_set  # kept because its own radius is 0
+    assert (3, 3) in out_set  # the dominant peak
+
+
+@pytest.mark.parametrize("shape", [(20, 20), (10, 12, 14)])
+def test_subset_and_includes_global_max(shape):
+    scipy_ndi = pytest.importorskip("scipy.ndimage")
+    rng = np.random.default_rng(0)
+    mask = rng.random(shape) > 0.2
+    dm = scipy_ndi.distance_transform_edt(mask).astype(np.float32)
+    coords = np.argwhere(dm > 1.5).astype(np.int64)
+    if len(coords) == 0:
+        pytest.skip("no candidate points for this random mask")
+
+    out = nms(dm, coords)
+    assert out.ndim == 2
+    assert out.shape[1] == len(shape)
+    assert out.shape[0] <= coords.shape[0]
+
+    # Every output point must be one of the input points.
+    in_set = {tuple(row) for row in coords.tolist()}
+    for row in out.tolist():
+        assert tuple(row) in in_set
+
+    # The global maximum of the distance map is always its own best point.
+    gmax = np.unravel_index(int(np.argmax(dm)), dm.shape)
+    if dm[gmax] > 1.5:
+        assert list(gmax) in out.tolist()
+
+
+@pytest.mark.parametrize("dtype", [np.int64, np.uint64, np.int32, np.uint32])
+def test_dtype_dispatch_equivalent(dtype):
+    dm = np.zeros((12, 12), dtype=np.float32)
+    dm[3, 3] = 5.0
+    dm[3, 4] = 4.0
+    dm[9, 9] = 2.0
+    pts = np.array([[3, 3], [3, 4], [9, 9]], dtype=dtype)
+    out = nms(dm, pts)
+    assert out.dtype == np.dtype(dtype)
+    # (3,3) suppresses (3,4); (9,9) is far and survives.
+    assert out.tolist() == [[3, 3], [9, 9]]
+
+
+def test_distance_map_float64_is_coerced():
+    dm = np.zeros((11, 11), dtype=np.float64)
+    dm[5, 5] = 5.0
+    dm[5, 6] = 4.0
+    pts = np.array([[5, 5], [5, 6]], dtype=np.int64)
+    out = nms(dm, pts)
+    assert out.tolist() == [[5, 5]]
+
+
+def test_deterministic():
+    scipy_ndi = pytest.importorskip("scipy.ndimage")
+    rng = np.random.default_rng(7)
+    mask = rng.random((40, 40)) > 0.25
+    dm = scipy_ndi.distance_transform_edt(mask).astype(np.float32)
+    coords = np.argwhere(dm > 1.0).astype(np.int64)
+    a = nms(dm, coords)
+    b = nms(dm, coords)
+    assert np.array_equal(a, b)
+
+
+def test_invalid_points_ndim():
+    dm = np.ones((10, 10), dtype=np.float32)
+    with pytest.raises(ValueError):
+        nms(dm, np.array([1, 2, 3], dtype=np.int64))
+
+
+def test_invalid_points_axis_length():
+    dm = np.ones((10, 10), dtype=np.float32)
+    with pytest.raises(ValueError):
+        nms(dm, np.array([[1, 2, 3]], dtype=np.int64))
+
+
+def test_invalid_dtype():
+    dm = np.ones((10, 10), dtype=np.float32)
+    with pytest.raises(TypeError):
+        nms(dm, np.array([[1.0, 2.0]], dtype=np.float32))
+
+
+def test_out_of_bounds_coordinate_raises():
+    dm = np.ones((10, 10), dtype=np.float32)
+    with pytest.raises((ValueError, RuntimeError)):
+        nms(dm, np.array([[10, 0]], dtype=np.int64))
+
+
+def test_negative_coordinate_raises():
+    dm = np.ones((10, 10), dtype=np.float32)
+    with pytest.raises((ValueError, RuntimeError)):
+        nms(dm, np.array([[-1, 0]], dtype=np.int64))