Nms (#46)

* Add lifted nhood construction functionality

* Add NMS functionality

Author: constantinpape

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.

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()

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