perf(features): parallelize boundary detection with bounded CPU workers

CHANGELOG.md

@@ -21,6 +21,11 @@ removed no sooner than the next major (see `docs/API_STABILITY.md`).
 
 ### Changed
 
+- 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
+  decisions, selected point order, schema, attributes and metadata are preserved
+  without a new dependency or public signature.
 - GPU CHW image packing validates output indices, buffer/binding/dispatch limits
   and profiler ownership/capacity before work allocation or pipeline setup.
   The original f32 normalization shader, CHW order, one-pass receipts and

crates/spatialrust-features/src/boundary.rs

@@ -11,7 +11,7 @@ use spatialrust_core::{
     SpatialResult,
 };
 use spatialrust_math::Vec3;
-use spatialrust_search::{KdTree, RadiusSearchIndex};
+use spatialrust_search::{parallel_worker_count, KdTree, RadiusSearchIndex};
 
 /// Configuration for [`BoundaryDetector`].
 #[derive(Clone, Copy, Debug, PartialEq)]
@@ -68,7 +68,19 @@ impl BoundaryDetector {
     }
 
     /// Computes the per-point boundary mask (`true` = boundary).
+    ///
+    /// Shares one KD-tree across bounded scoped CPU workers using the existing
+    /// point-staging policy. Small clouds and one-worker calls run on the caller.
+    /// Point order, neighborhood order and per-point arithmetic are preserved.
     pub fn boundary_mask(&self, input: &PointCloud) -> SpatialResult<Vec<bool>> {
+        self.boundary_mask_with_workers(input, parallel_worker_count(input.len()))
+    }
+
+    fn boundary_mask_with_workers(
+        &self,
+        input: &PointCloud,
+        workers: usize,
+    ) -> SpatialResult<Vec<bool>> {
         if self.config.search_radius <= 0.0 || self.config.search_radius.is_nan() {
             return Err(SpatialError::InvalidArgument("search_radius must be positive".to_owned()));
         }
@@ -82,29 +94,46 @@ impl BoundaryDetector {
         let tree = KdTree::from_slices(x, y, z);
 
         let mut mask = vec![false; len];
-        let mut angles: Vec<f32> = Vec::new();
-        for i in 0..len {
-            let p = Vec3::new(x[i], y[i], z[i]);
-            let normal = Vec3::new(nx[i], ny[i], nz[i]);
-            let (u, v) = tangent_basis(normal);
-
-            let neighbors = tree.radius_search(p.x, p.y, p.z, self.config.search_radius);
-            angles.clear();
-            for neighbor in &neighbors {
-                let j = neighbor.index;
-                if j == i {
-                    continue;
+        let fill = |start: usize, band: &mut [bool]| {
+            // Reuse angular scratch within each worker, without retaining an
+            // all-point neighbor table or allocating a worker result vector.
+            let mut angles: Vec<f32> = Vec::new();
+            for (offset, is_boundary) in band.iter_mut().enumerate() {
+                let i = start + offset;
+                let p = Vec3::new(x[i], y[i], z[i]);
+                let normal = Vec3::new(nx[i], ny[i], nz[i]);
+                let (u, v) = tangent_basis(normal);
+
+                let neighbors = tree.radius_search(p.x, p.y, p.z, self.config.search_radius);
+                angles.clear();
+                for neighbor in &neighbors {
+                    let j = neighbor.index;
+                    if j == i {
+                        continue;
+                    }
+                    let d = Vec3::new(x[j] - p.x, y[j] - p.y, z[j] - p.z);
+                    let (du, dv) = (d.dot(u), d.dot(v));
+                    if du * du + dv * dv > 1e-12 {
+                        angles.push(dv.atan2(du));
+                    }
                 }
-                let d = Vec3::new(x[j] - p.x, y[j] - p.y, z[j] - p.z);
-                let (du, dv) = (d.dot(u), d.dot(v));
-                if du * du + dv * dv > 1e-12 {
-                    angles.push(dv.atan2(du));
+                if angles.len() < self.config.min_neighbors {
+                    continue;
                 }
+                *is_boundary = max_angle_gap(&mut angles) > self.config.angle_threshold;
             }
-            if angles.len() < self.config.min_neighbors {
-                continue;
-            }
-            mask[i] = max_angle_gap(&mut angles) > self.config.angle_threshold;
+        };
+        let workers = workers.max(1).min(len);
+        if workers == 1 {
+            fill(0, &mut mask);
+        } else {
+            let chunk_size = len.div_ceil(workers);
+            std::thread::scope(|scope| {
+                for (band, chunk) in mask.chunks_mut(chunk_size).enumerate() {
+                    let fill = &fill;
+                    scope.spawn(move || fill(band * chunk_size, chunk));
+                }
+            });
         }
         Ok(mask)
     }
@@ -119,6 +148,19 @@ impl BoundaryDetector {
     }
 }
 
+#[cfg(all(test, not(target_arch = "wasm32")))]
+#[path = "../tests/fixtures/boundary_serial.rs"]
+#[allow(dead_code)]
+mod previous;
+
+#[cfg(all(test, not(target_arch = "wasm32")))]
+#[path = "boundary_parallel_tests.rs"]
+mod parallel_tests;
+
+#[cfg(all(test, not(target_arch = "wasm32")))]
+#[path = "boundary_parallel_benchmark.rs"]
+mod parallel_benchmark;
+
 /// Builds an orthonormal basis `(u, v)` spanning the plane perpendicular to `n`.
 fn tangent_basis(n: Vec3<f32>) -> (Vec3<f32>, Vec3<f32>) {
     let helper = if n.x.abs() < 0.9 { Vec3::new(1.0, 0.0, 0.0) } else { Vec3::new(0.0, 1.0, 0.0) };

crates/spatialrust-features/src/boundary_parallel_benchmark.rs

@@ -0,0 +1,52 @@
+use super::{
+    parallel_tests::{fixture, previous_mask},
+    *,
+};
+use std::{hint::black_box, time::Instant};
+
+#[test]
+#[ignore = "completed boundary API comparison; explicit native invocation only"]
+fn paired_previous_serial_parallel_boundary() {
+    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, mask allocation, angle-gap evaluation, worker spawn/join and scratch/tree destruction; exclude fixture creation, mask verification and returned mask 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 config = BoundaryConfig {
+                search_radius: if kind == 0 { 0.11 } else { 0.16 },
+                min_neighbors: 3,
+                ..Default::default()
+            };
+            let cloud = fixture(points, kind);
+            let expected = previous_mask(&cloud, config);
+            let detector = BoundaryDetector::new(config);
+            let workers = parallel_worker_count(points);
+            let run = |mode| match mode {
+                0 => previous_mask(&cloud, config),
+                1 => detector.boundary_mask_with_workers(&cloud, 1).unwrap(),
+                _ => detector.boundary_mask(&cloud).unwrap(),
+            };
+            for _ in 0..2 {
+                for mode in 0..3 {
+                    assert_eq!(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_eq!(actual, expected);
+                    black_box(&actual);
+                }
+                eprintln!(
+                    "{points},{kind},{},{workers},{round},{:.6},{:.6},{:.6}",
+                    config.search_radius, ms[0], ms[1], ms[2]
+                );
+            }
+        }
+    }
+}