perf(features): parallelize ISS keypoints with bounded CPU workers

CHANGELOG.md

@@ -21,6 +21,12 @@ removed no sooner than the next major (see `docs/API_STABILITY.md`).
 
 ### Changed
 
+- ISS keypoint detection shares one KD-tree across bounded scoped CPU workers
+  for covariance saliency and local-maximum selection under `feature-iss`.
+  Saliency completes before selection; small/one-worker clouds stay serial.
+  Equal maxima, selected point order, schema, attributes and metadata remain
+  unchanged. Parallel selection uses one temporary byte per point, without a
+  new dependency or public signature.
 - Point-cloud boundary detection shares one KD-tree across bounded scoped CPU
   workers under `feature-boundary`. Workers reuse angular scratch and write
   disjoint mask bands directly; small/one-worker clouds stay serial. Boundary

crates/spatialrust-features/src/iss.rs

@@ -10,7 +10,7 @@ use spatialrust_core::{
     HasPositions3, PointBuffer, PointBufferSet, PointCloud, SpatialError, SpatialResult,
 };
 use spatialrust_math::{symmetric_eigen3, Mat3};
-use spatialrust_search::{KdTree, RadiusSearchIndex};
+use spatialrust_search::{parallel_worker_count, KdTree, RadiusSearchIndex};
 
 /// Configuration for [`IssKeypointDetector`].
 #[derive(Clone, Copy, Debug, PartialEq)]
@@ -76,53 +76,87 @@ impl IssKeypointDetector {
     }
 
     /// Detects keypoints and returns their indices and a keypoint sub-cloud.
+    ///
+    /// Shares one KD-tree across bounded scoped CPU workers. Saliency completes
+    /// before local-maximum selection reads it. Small/one-worker calls run on
+    /// the caller; neighborhood arithmetic, ties and point order are preserved.
     pub fn detect(&self, input: &PointCloud) -> SpatialResult<IssKeypointResult> {
+        self.detect_with_workers(input, parallel_worker_count(input.len()))
+    }
+
+    fn detect_with_workers(
+        &self,
+        input: &PointCloud,
+        workers: usize,
+    ) -> SpatialResult<IssKeypointResult> {
         if self.config.salient_radius <= 0.0 || self.config.non_max_radius <= 0.0 {
             return Err(SpatialError::InvalidArgument("ISS radii must be positive".to_owned()));
         }
 
         let (x, y, z) = input.positions3()?;
         let len = input.len();
         let tree = KdTree::from_slices(x, y, z);
+        let workers = workers.max(1).min(len.max(1));
 
         // Saliency = smallest eigenvalue (lambda3) where the eigenvalue-ratio
         // tests pass; NaN marks a non-salient point.
         let mut saliency = vec![f32::NAN; len];
-        for i in 0..len {
-            let neighbors = tree.radius_search(x[i], y[i], z[i], self.config.salient_radius);
-            if neighbors.len() < self.config.min_neighbors {
-                continue;
-            }
-            let Some(eigenvalues) = neighborhood_eigenvalues(x, y, z, &neighbors) else {
-                continue;
-            };
-            // Ascending order: l3 <= l2 <= l1.
-            let (l3, l2, l1) = (eigenvalues[0], eigenvalues[1], eigenvalues[2]);
-            if l1 <= 0.0 || l2 <= 0.0 {
-                continue;
-            }
-            if (l2 / l1) < f64::from(self.config.gamma_21)
-                && (l3 / l2) < f64::from(self.config.gamma_32)
-            {
-                saliency[i] = l3 as f32;
+        fill_bands(&mut saliency, workers, |start, band| {
+            for (offset, value) in band.iter_mut().enumerate() {
+                let i = start + offset;
+                let neighbors = tree.radius_search(x[i], y[i], z[i], self.config.salient_radius);
+                if neighbors.len() < self.config.min_neighbors {
+                    continue;
+                }
+                let Some(eigenvalues) = neighborhood_eigenvalues(x, y, z, &neighbors) else {
+                    continue;
+                };
+                // Ascending order: l3 <= l2 <= l1.
+                let (l3, l2, l1) = (eigenvalues[0], eigenvalues[1], eigenvalues[2]);
+                if l1 <= 0.0 || l2 <= 0.0 {
+                    continue;
+                }
+                if (l2 / l1) < f64::from(self.config.gamma_21)
+                    && (l3 / l2) < f64::from(self.config.gamma_32)
+                {
+                    *value = l3 as f32;
+                }
             }
-        }
+        });
 
-        // Non-maximum suppression: keep points whose saliency is the strict
-        // maximum within `non_max_radius`.
-        let mut indices = Vec::new();
-        for i in 0..len {
+        // The completed saliency array is immutable throughout this phase.
+        // Equal maxima all survive, retaining the original >= tie behavior.
+        let is_keypoint = |i: usize| {
             let s = saliency[i];
             if s.is_nan() {
-                continue;
+                return false;
             }
             let neighbors = tree.radius_search(x[i], y[i], z[i], self.config.non_max_radius);
-            let is_local_max = neighbors.iter().all(|n| {
+            neighbors.iter().all(|n| {
                 let other = saliency[n.index];
                 n.index == i || other.is_nan() || s >= other
+            })
+        };
+        let mut indices = Vec::new();
+        if workers == 1 {
+            for i in 0..len {
+                if is_keypoint(i) {
+                    indices.push(i);
+                }
+            }
+        } else {
+            // One temporary byte per point avoids worker-local index vectors
+            // and a second result-sized index allocation during ordered joins.
+            let mut selected = vec![false; len];
+            fill_bands(&mut selected, workers, |start, band| {
+                for (offset, value) in band.iter_mut().enumerate() {
+                    *value = is_keypoint(start + offset);
+                }
             });
-            if is_local_max {
-                indices.push(i);
+            for (i, selected) in selected.into_iter().enumerate() {
+                if selected {
+                    indices.push(i);
+                }
             }
         }
 
@@ -131,6 +165,34 @@ impl IssKeypointDetector {
     }
 }
 
+/// Fills disjoint output bands, joining all workers before the next phase.
+fn fill_bands<T: Send, F: Fn(usize, &mut [T]) + Sync>(output: &mut [T], workers: usize, fill: F) {
+    if workers == 1 {
+        fill(0, output);
+        return;
+    }
+    let chunk_size = output.len().div_ceil(workers);
+    std::thread::scope(|scope| {
+        for (band, chunk) in output.chunks_mut(chunk_size).enumerate() {
+            let fill = &fill;
+            scope.spawn(move || fill(band * chunk_size, chunk));
+        }
+    });
+}
+
+#[cfg(all(test, not(target_arch = "wasm32")))]
+#[path = "../tests/fixtures/iss_serial.rs"]
+#[allow(dead_code)]
+mod previous;
+
+#[cfg(all(test, not(target_arch = "wasm32")))]
+#[path = "iss_parallel_tests.rs"]
+mod parallel_tests;
+
+#[cfg(all(test, not(target_arch = "wasm32")))]
+#[path = "iss_parallel_benchmark.rs"]
+mod parallel_benchmark;
+
 /// Eigenvalues (ascending) of the neighborhood covariance about its centroid.
 fn neighborhood_eigenvalues(
     x: &[f32],

crates/spatialrust-features/src/iss_parallel_benchmark.rs

@@ -0,0 +1,54 @@
+use super::{
+    parallel_tests::{assert_result, fixture, previous_result},
+    *,
+};
+use std::{hint::black_box, time::Instant};
+
+#[test]
+#[ignore = "completed ISS detection comparison; explicit native invocation only"]
+fn paired_previous_serial_parallel_iss() {
+    let available = std::thread::available_parallelism().map_or(1, |n| n.get());
+    eprintln!("available_parallelism={available}; 2 warmups + 20 rotating rounds; default shared staging policy");
+    eprintln!("Timers include validation, KD-tree build, saliency allocation/covariance, local-maximum selection, worker spawn/join, ordered indices, schema/attribute gathering and scratch/tree destruction; exclude fixture creation, result verification and returned result drop.");
+    eprintln!("points,kind,radius,workers,round,previous_ms,serial_ms,parallel_ms");
+    for points in [512usize, 16384, 32768, 65536] {
+        for kind in 0..3 {
+            let radius = if kind == 0 { 0.11 } else { 0.16 };
+            let config = IssKeypointConfig::with_radii(radius, radius * 0.75);
+            let cloud = fixture(points, kind);
+            let expected = previous_result(&cloud, config);
+            let detector = IssKeypointDetector::new(config);
+            let workers = parallel_worker_count(points);
+            let run = |mode| {
+                if mode == 0 {
+                    let previous = previous_result(&cloud, config);
+                    IssKeypointResult { indices: previous.indices, keypoints: previous.keypoints }
+                } else if mode == 1 {
+                    detector.detect_with_workers(&cloud, 1).unwrap()
+                } else {
+                    detector.detect(&cloud).unwrap()
+                }
+            };
+            for _ in 0..2 {
+                for mode in 0..3 {
+                    assert_result(&run(mode), &expected);
+                }
+            }
+            for round in 0..20 {
+                let mut ms = [0.0; 3];
+                for offset in 0..3 {
+                    let mode = (round + offset) % 3;
+                    let start = Instant::now();
+                    let actual = run(mode);
+                    ms[mode] = start.elapsed().as_secs_f64() * 1000.0;
+                    assert_result(&actual, &expected);
+                    black_box(&actual);
+                }
+                eprintln!(
+                    "{points},{kind},{radius},{workers},{round},{:.6},{:.6},{:.6}",
+                    ms[0], ms[1], ms[2]
+                );
+            }
+        }
+    }
+}