perf(metrics): bound parallel Chamfer and Hausdorff queries

CHANGELOG.md

@@ -21,6 +21,12 @@ removed no sooner than the next major (see `docs/API_STABILITY.md`).
 
 ### Changed
 
+- Chamfer/Hausdorff cloud metrics share each target KD-tree across bounded CPU
+  query workers under `metrics-distance`. A reusable score window holds at
+  most one default point chunk per worker, with one reused query buffer per
+  active band; f64 sums and maxima retain original
+  query order and exact output words. Small/one-worker sides keep the original
+  serial loop, with no new dependency, public signature or device transfer.
 - 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.

crates/spatialrust-metrics/src/distance.rs

@@ -1,7 +1,9 @@
 //! Chamfer and Hausdorff distances between two point clouds.
 
-use spatialrust_core::{HasPositions3, PointCloud, SpatialError, SpatialResult};
-use spatialrust_search::{KdTree, NearestNeighborIndex};
+use spatialrust_core::{
+    HasPositions3, PointCloud, SpatialError, SpatialResult, DEFAULT_SPATIAL_TENSOR_CHUNK_SIZE,
+};
+use spatialrust_search::{parallel_worker_count, KdTree, NearestNeighborIndex};
 
 /// A bundle of directed and symmetric cloud-to-cloud distance statistics.
 ///
@@ -35,7 +37,25 @@ impl CloudDistances {
 
 /// Computes directed NN distance statistics between `a` and `b` in both
 /// directions with a single KD-tree per side.
+///
+/// Large query sides share their target tree across bounded scoped CPU workers.
+/// A reusable window holds at most one default point chunk per worker. Distance
+/// accumulation retains input order; small/one-worker sides remain serial.
 pub fn cloud_distances(a: &PointCloud, b: &PointCloud) -> SpatialResult<CloudDistances> {
+    cloud_distances_with_workers(
+        a,
+        b,
+        parallel_worker_count(a.len()),
+        parallel_worker_count(b.len()),
+    )
+}
+
+fn cloud_distances_with_workers(
+    a: &PointCloud,
+    b: &PointCloud,
+    workers_a: usize,
+    workers_b: usize,
+) -> SpatialResult<CloudDistances> {
     if a.is_empty() || b.is_empty() {
         return Err(SpatialError::InvalidArgument("both clouds must be non-empty".to_owned()));
     }
@@ -45,14 +65,62 @@ pub fn cloud_distances(a: &PointCloud, b: &PointCloud) -> SpatialResult<CloudDis
     let tree_b = KdTree::from_slices(bx, by, bz);
     let tree_a = KdTree::from_slices(ax, ay, az);
 
-    let (mean_sq_a_to_b, max_a_to_b) = directed(ax, ay, az, &tree_b);
-    let (mean_sq_b_to_a, max_b_to_a) = directed(bx, by, bz, &tree_a);
+    let (mean_sq_a_to_b, max_a_to_b) = directed(ax, ay, az, &tree_b, workers_a);
+    let (mean_sq_b_to_a, max_b_to_a) = directed(bx, by, bz, &tree_a, workers_b);
 
     Ok(CloudDistances { mean_sq_a_to_b, mean_sq_b_to_a, max_a_to_b, max_b_to_a })
 }
 
 /// Mean squared and max NN distance from each query point to `target`.
-fn directed(x: &[f32], y: &[f32], z: &[f32], target: &KdTree) -> (f64, f64) {
+fn directed(x: &[f32], y: &[f32], z: &[f32], target: &KdTree, workers: usize) -> (f64, f64) {
+    let workers = workers.max(1).min(x.len().max(1));
+    if workers == 1 {
+        return directed_serial(x, y, z, target);
+    }
+    // Preserve Some(NaN) versus None, including the original sum/max behavior.
+    // Reuse a worker-bounded window instead of staging every query's result.
+    let window_size = x.len().min(workers.saturating_mul(DEFAULT_SPATIAL_TENSOR_CHUNK_SIZE));
+    let mut nearest = vec![None; window_size];
+    let mut sum_sq = 0.0_f64;
+    let mut max = 0.0_f64;
+    for start in (0..x.len()).step_by(window_size) {
+        let count = window_size.min(x.len() - start);
+        let active_workers = workers.min(count);
+        let fill = |offset: usize, band: &mut [Option<f32>]| {
+            // nearest_one uses this same one-neighbor traversal. Reuse its
+            // output buffer; nearest_k_into resets the query accumulator.
+            let mut neighbor = Vec::new();
+            for (i, distance) in band.iter_mut().enumerate() {
+                let i = start + offset + i;
+                target.nearest_k_into(x[i], y[i], z[i], 1, &mut neighbor);
+                *distance = neighbor.first().map(|n| n.distance_squared);
+            }
+        };
+        if active_workers == 1 {
+            fill(0, &mut nearest[..count]);
+        } else {
+            let chunk_size = count.div_ceil(active_workers);
+            std::thread::scope(|scope| {
+                for (band, chunk) in nearest[..count].chunks_mut(chunk_size).enumerate() {
+                    let fill = &fill;
+                    scope.spawn(move || fill(band * chunk_size, chunk));
+                }
+            });
+        }
+        // Never combine worker-local sums: f64 addition is order-sensitive.
+        for distance in nearest[..count].iter().flatten() {
+            let d_sq = f64::from(*distance);
+            sum_sq += d_sq;
+            let d = d_sq.sqrt();
+            if d > max {
+                max = d;
+            }
+        }
+    }
+    (sum_sq / x.len() as f64, max)
+}
+
+fn directed_serial(x: &[f32], y: &[f32], z: &[f32], target: &KdTree) -> (f64, f64) {
     let mut sum_sq = 0.0_f64;
     let mut max = 0.0_f64;
     for i in 0..x.len() {
@@ -80,6 +148,19 @@ pub fn hausdorff_distance(a: &PointCloud, b: &PointCloud) -> SpatialResult<f64>
     Ok(cloud_distances(a, b)?.hausdorff())
 }
 
+#[cfg(all(test, not(target_arch = "wasm32")))]
+#[path = "../tests/fixtures/distance_serial.rs"]
+#[allow(dead_code)]
+mod previous;
+
+#[cfg(all(test, not(target_arch = "wasm32")))]
+#[path = "distance_parallel_tests.rs"]
+mod parallel_tests;
+
+#[cfg(all(test, not(target_arch = "wasm32")))]
+#[path = "distance_parallel_benchmark.rs"]
+mod parallel_benchmark;
+
 #[cfg(test)]
 mod tests {
     use super::*;

crates/spatialrust-metrics/src/distance_parallel_benchmark.rs

@@ -0,0 +1,60 @@
+use super::{
+    parallel_tests::{fixture_pair, previous_words, words},
+    *,
+};
+use std::{hint::black_box, time::Instant};
+
+#[test]
+#[ignore = "completed cloud-distance comparison; explicit native invocation only"]
+fn paired_previous_serial_parallel_distances() {
+    let available = std::thread::available_parallelism().map_or(1, |n| n.get());
+    eprintln!("available_parallelism={available}; optional_distance_bytes={}; 2 warmups + 20 rotating rounds; default shared staging policy", std::mem::size_of::<Option<f32>>());
+    eprintln!("Timers include validation, both KD-tree builds, both directions, score-window allocation/reuse, worker spawn/join, ordered f64 accumulation and tree/scratch destruction; exclude fixtures and exact-word verification.");
+    eprintln!("points,kind,target_points,workers,round,previous_ms,serial_ms,parallel_ms");
+    for points in [512usize, 16384, 32768, 65536] {
+        for kind in 0..3 {
+            let (a, b) = fixture_pair(points, kind);
+            let expected = previous_words(&a, &b);
+            let wa = parallel_worker_count(a.len());
+            let wb = parallel_worker_count(b.len());
+            let run = |mode| {
+                if mode == 0 {
+                    let p = previous::cloud_distances(&a, &b).unwrap();
+                    CloudDistances {
+                        mean_sq_a_to_b: p.mean_sq_a_to_b,
+                        mean_sq_b_to_a: p.mean_sq_b_to_a,
+                        max_a_to_b: p.max_a_to_b,
+                        max_b_to_a: p.max_b_to_a,
+                    }
+                } else if mode == 1 {
+                    cloud_distances_with_workers(&a, &b, 1, 1).unwrap()
+                } else {
+                    cloud_distances(&a, &b).unwrap()
+                }
+            };
+            for _ in 0..2 {
+                for mode in 0..3 {
+                    assert_eq!(words(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!(words(actual), expected);
+                    black_box(actual);
+                }
+                eprintln!(
+                    "{points},{kind},{},{wa}/{wb},{round},{:.6},{:.6},{:.6}",
+                    b.len(),
+                    ms[0],
+                    ms[1],
+                    ms[2]
+                );
+            }
+        }
+    }
+}