From 2ba7b6575f8234450dbbc54ed325e32158aa65e1 Mon Sep 17 00:00:00 2001 From: rsasaki0109 Date: Mon, 12 Oct 2026 11:36:15 +0900 Subject: [PATCH] perf(metrics): stage bounded parallel cloud distance queries --- CHANGELOG.md | 6 + crates/spatialrust-metrics/src/distance.rs | 91 +++++++- .../src/distance_parallel_benchmark.rs | 60 +++++ .../src/distance_parallel_tests.rs | 206 ++++++++++++++++++ .../tests/fixtures/distance_serial.rs | 81 +++++++ .../tests/distance_parallel_smoke.rs | 23 ++ docs/ARCHITECTURE.md | 13 ++ docs/FEATURE_MATRIX.md | 7 + docs/ROADMAP.md | 1 + notes/2026-10-12_parallel_cloud_distances.md | 194 +++++++++++++++++ 10 files changed, 677 insertions(+), 5 deletions(-) create mode 100644 crates/spatialrust-metrics/src/distance_parallel_benchmark.rs create mode 100644 crates/spatialrust-metrics/src/distance_parallel_tests.rs create mode 100644 crates/spatialrust-metrics/tests/fixtures/distance_serial.rs create mode 100644 crates/spatialrust/tests/distance_parallel_smoke.rs create mode 100644 notes/2026-10-12_parallel_cloud_distances.md diff --git a/CHANGELOG.md b/CHANGELOG.md index d82c0470..539239d5 100644 --- a/CHANGELOG.md +++ b/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. diff --git a/crates/spatialrust-metrics/src/distance.rs b/crates/spatialrust-metrics/src/distance.rs index e2c4cd31..0d162166 100644 --- a/crates/spatialrust-metrics/src/distance.rs +++ b/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 { + 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 { 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 (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]| { + // 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 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::*; diff --git a/crates/spatialrust-metrics/src/distance_parallel_benchmark.rs b/crates/spatialrust-metrics/src/distance_parallel_benchmark.rs new file mode 100644 index 00000000..432be974 --- /dev/null +++ b/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::>()); + 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] + ); + } + } + } +} diff --git a/crates/spatialrust-metrics/src/distance_parallel_tests.rs b/crates/spatialrust-metrics/src/distance_parallel_tests.rs new file mode 100644 index 00000000..43f1e376 --- /dev/null +++ b/crates/spatialrust-metrics/src/distance_parallel_tests.rs @@ -0,0 +1,206 @@ +use super::*; +use spatialrust_core::{ + DType, FieldSemantic, PointBuffer, PointBufferSet, PointField, PointSchema, SpatialMetadata, + StandardSchemas, Timestamp, +}; + +fn from_points(points: impl IntoIterator) -> PointCloud { + let mut columns: [Vec; 3] = std::array::from_fn(|_| Vec::new()); + for p in points { + for (column, value) in columns.iter_mut().zip(p) { + column.push(value); + } + } + let mut buffers = PointBufferSet::new(); + for (name, values) in ["x", "y", "z"].into_iter().zip(columns) { + buffers.insert(name, PointBuffer::F32(values)); + } + PointCloud::try_from_parts(StandardSchemas::point_xyz(), buffers, Default::default()).unwrap() +} + +fn fixture(len: usize, kind: usize, shifted: bool) -> PointCloud { + from_points((0..len).map(|i| { + let j = if kind == 2 { i / 2 } else { i }; + let x = (j % 128) as f32 * 0.05 - 3.2; + let y = (j / 128) as f32 * 0.05 - 1.6; + let z = if kind == 0 { + 0.0 + } else { + 0.015 * (x * x + y * y) + ((j * 37 % 101) as f32 / 100.0 - 0.5) * 0.08 + }; + if shifted { + [x + 0.017, y - 0.011, z + 0.023] + } else { + [x, y, z] + } + })) +} + +pub(super) fn fixture_pair(len: usize, kind: usize) -> (PointCloud, PointCloud) { + let target_len = if kind == 2 { len.div_ceil(2) } else { len }; + (fixture(len, kind, false), fixture(target_len, kind.min(1), true)) +} + +pub(super) fn words(value: CloudDistances) -> [u64; 6] { + [ + value.mean_sq_a_to_b, + value.mean_sq_b_to_a, + value.max_a_to_b, + value.max_b_to_a, + value.chamfer(), + value.hausdorff(), + ] + .map(f64::to_bits) +} + +pub(super) fn previous_words(a: &PointCloud, b: &PointCloud) -> [u64; 6] { + let p = previous::cloud_distances(a, b).unwrap(); + [p.mean_sq_a_to_b, p.mean_sq_b_to_a, p.max_a_to_b, p.max_b_to_a, p.chamfer(), p.hausdorff()] + .map(f64::to_bits) +} + +#[test] +fn frozen_metric_words_match_across_staging_windows_and_asymmetric_clouds() { + let mut comparisons = 0; + for len in [1usize, 17, 4095, 4096, 4103, 32767, 32768, 32771, 65537] { + for kind in 0..3 { + let (a, b) = fixture_pair(len, kind); + let expected = previous_words(&a, &b); + for workers in [1, 2, 6] { + assert_eq!( + words(cloud_distances_with_workers(&a, &b, workers, workers).unwrap()), + expected, + "len={len}, kind={kind}, workers={workers}" + ); + comparisons += 6; + } + assert_eq!(words(cloud_distances(&a, &b).unwrap()), expected); + comparisons += 6; + if len == 17 || len == 32771 { + assert_eq!(chamfer_distance(&a, &b).unwrap().to_bits(), expected[4]); + assert_eq!(hausdorff_distance(&a, &b).unwrap().to_bits(), expected[5]); + comparisons += 2; + } + } + } + eprintln!("Verified {comparisons} exact f64 metric words against frozen serial distances."); +} + +#[test] +fn independent_asymmetric_statistics_and_ties_are_preserved() { + let a = from_points([[0.0, 0.0, 0.0], [2.0, 0.0, 0.0]]); + let b = from_points([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [4.0, 0.0, 0.0]]); + let expected = [0.5, 5.0 / 3.0, 1.0, 2.0, 0.5 + 5.0 / 3.0, 2.0].map(f64::to_bits); + for (wa, wb) in [(0, 0), (1, 6), (2, 1), (6, 6), (usize::MAX, usize::MAX)] { + assert_eq!(words(cloud_distances_with_workers(&a, &b, wa, wb).unwrap()), expected); + } + let a = from_points([[-1.0, 0.0, 0.0], [1.0, 0.0, 0.0], [1.0, 0.0, 0.0]]); + let b = from_points([[0.0, 0.0, 0.0]]); + assert_eq!( + words(cloud_distances_with_workers(&a, &b, 2, 1).unwrap()), + [1.0, 1.0, 1.0, 1.0, 2.0, 1.0].map(f64::to_bits) + ); +} + +#[test] +fn ordered_accumulation_survives_multiple_windows_and_one_point_tail() { + let len = 65537; + let a = from_points( + (0..len).map(|i| [if i == 0 { 1e8 } else { (i - 1) as f32 / 65536.0 }, 0.0, 0.0]), + ); + let b = from_points([[0.0, 0.0, 0.0]]); + let first_sq = f64::from(1e8_f32 * 1e8_f32); + let ordered_mean = first_sq / len as f64; + // Every later squared distance is below 1 and rounds away at this + // magnitude. Distinct small coordinates avoid stressing the unrelated + // KD-tree builder's equal-coordinate partitioning in this math test. + let small_sum: f64 = (0..(len - 1)) + .map(|i| { + let x = i as f32 / 65536.0; + f64::from(x * x) + }) + .sum(); + let regrouped_mean = (first_sq + small_sum) / len as f64; + assert_ne!(ordered_mean.to_bits(), regrouped_mean.to_bits()); + let expected = previous_words(&a, &b); + assert_eq!(expected[0], ordered_mean.to_bits()); + for (wa, wb) in [(1, 1), (2, 1), (6, 2)] { + assert_eq!(words(cloud_distances_with_workers(&a, &b, wa, wb).unwrap()), expected); + } + assert_eq!(words(cloud_distances(&a, &b).unwrap()), expected); +} + +#[test] +fn nonfinite_distances_missing_neighbors_and_empty_errors_match_serial() { + let zero = from_points([[0.0, 0.0, 0.0]]); + for points in [ + vec![[f32::from_bits(0x7fc00037), 0.0, 0.0]; 33], + vec![[f32::INFINITY, 0.0, 0.0]; 33], + vec![[f32::MAX, f32::MAX, f32::MAX]; 33], + vec![[-0.0, 0.0, -0.0]; 33], + (0..33).map(|i| [if i % 2 == 0 { f32::NAN } else { i as f32 }, 0.0, 0.0]).collect(), + ] { + let cloud = from_points(points); + for (a, b) in [(&cloud, &zero), (&zero, &cloud)] { + let expected = previous_words(a, b); + for workers in [1, 2, 6] { + assert_eq!( + words(cloud_distances_with_workers(a, b, workers, workers).unwrap()), + expected + ); + } + } + } + let nan = from_points([[f32::NAN, 0.0, 0.0]]); + assert!(cloud_distances_with_workers(&nan, &zero, 2, 1).unwrap().mean_sq_a_to_b.is_nan()); + let no_target = KdTree::from_slices(&[], &[], &[]); + assert_eq!(directed(&[0.0, f32::NAN, 1.0], &[0.0; 3], &[0.0; 3], &no_target, 2), (0.0, 0.0)); + let empty = from_points([]); + for (a, b) in [(&empty, &zero), (&zero, &empty), (&empty, &empty)] { + let expected = previous::cloud_distances(a, b).unwrap_err().to_string(); + for (wa, wb) in [(0, 0), (1, 1), (2, 6), (usize::MAX, usize::MAX)] { + assert_eq!( + cloud_distances_with_workers(a, b, wa, wb).unwrap_err().to_string(), + expected + ); + } + assert_eq!(chamfer_distance(a, b).unwrap_err().to_string(), expected); + assert_eq!(hausdorff_distance(a, b).unwrap_err().to_string(), expected); + } +} + +#[test] +fn renamed_xyz_attributes_metadata_and_concurrent_calls_are_preserved() { + let (source, b) = fixture_pair(4103, 1); + let mut schema = PointSchema::new(); + let mut buffers = PointBufferSet::new(); + for (field, name) in source.schema().fields().iter().zip(["px", "py", "pz"]) { + schema = schema.with_field(PointField::scalar(name, field.semantic, field.dtype)); + buffers.insert(name, source.field(&field.name).unwrap().clone()); + } + let id = PointBuffer::U32((0..source.len()).map(|i| i as u32).collect()); + schema = schema.with_field(PointField::scalar("id", FieldSemantic::Unknown, DType::U32)); + buffers.insert("id", id.clone()); + let a = PointCloud::try_from_parts( + schema, + buffers, + SpatialMetadata::new("sensor", Timestamp::from_nanos(37)), + ) + .unwrap(); + let input_before = a.clone(); + let expected = previous_words(&a, &b); + std::thread::scope(|scope| { + let calls: Vec<_> = [(2, 6), (6, 2)] + .into_iter() + .map(|(wa, wb)| { + let (a, b) = (&a, &b); + scope.spawn(move || cloud_distances_with_workers(a, b, wa, wb).unwrap()) + }) + .collect(); + for call in calls { + assert_eq!(words(call.join().unwrap()), expected); + } + }); + assert_eq!(a, input_before); + assert_eq!(a.field("id").unwrap(), &id); +} diff --git a/crates/spatialrust-metrics/tests/fixtures/distance_serial.rs b/crates/spatialrust-metrics/tests/fixtures/distance_serial.rs new file mode 100644 index 00000000..c3efe4df --- /dev/null +++ b/crates/spatialrust-metrics/tests/fixtures/distance_serial.rs @@ -0,0 +1,81 @@ +//! Chamfer and Hausdorff distances between two point clouds. + +use spatialrust_core::{HasPositions3, PointCloud, SpatialError, SpatialResult}; +use spatialrust_search::{KdTree, NearestNeighborIndex}; + +/// A bundle of directed and symmetric cloud-to-cloud distance statistics. +/// +/// "Directed" means nearest-neighbor distances from one cloud to the other; the +/// metrics are not symmetric on their own, so both directions are reported. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct CloudDistances { + /// Mean squared NN distance from `a` to `b`. + pub mean_sq_a_to_b: f64, + /// Mean squared NN distance from `b` to `a`. + pub mean_sq_b_to_a: f64, + /// Largest NN distance from `a` to `b` (directed Hausdorff). + pub max_a_to_b: f64, + /// Largest NN distance from `b` to `a` (directed Hausdorff). + pub max_b_to_a: f64, +} + +impl CloudDistances { + /// Symmetric Chamfer distance: the sum of the two mean squared NN distances. + #[must_use] + pub fn chamfer(&self) -> f64 { + self.mean_sq_a_to_b + self.mean_sq_b_to_a + } + + /// Symmetric Hausdorff distance: the larger of the two directed maxima. + #[must_use] + pub fn hausdorff(&self) -> f64 { + self.max_a_to_b.max(self.max_b_to_a) + } +} + +/// Computes directed NN distance statistics between `a` and `b` in both +/// directions with a single KD-tree per side. +pub fn cloud_distances(a: &PointCloud, b: &PointCloud) -> SpatialResult { + if a.is_empty() || b.is_empty() { + return Err(SpatialError::InvalidArgument("both clouds must be non-empty".to_owned())); + } + let (ax, ay, az) = a.positions3()?; + let (bx, by, bz) = b.positions3()?; + + 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); + + 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) { + let mut sum_sq = 0.0_f64; + let mut max = 0.0_f64; + for i in 0..x.len() { + if let Some(neighbor) = target.nearest_one(x[i], y[i], z[i]) { + let d_sq = f64::from(neighbor.distance_squared); + sum_sq += d_sq; + let d = d_sq.sqrt(); + if d > max { + max = d; + } + } + } + (sum_sq / x.len() as f64, max) +} + +/// Symmetric Chamfer distance between two clouds (sum of mean squared NN +/// distances in both directions). Lower is better; zero for identical clouds. +pub fn chamfer_distance(a: &PointCloud, b: &PointCloud) -> SpatialResult { + Ok(cloud_distances(a, b)?.chamfer()) +} + +/// Symmetric Hausdorff distance between two clouds (the largest NN distance in +/// either direction). Captures the worst-case discrepancy. +pub fn hausdorff_distance(a: &PointCloud, b: &PointCloud) -> SpatialResult { + Ok(cloud_distances(a, b)?.hausdorff()) +} diff --git a/crates/spatialrust/tests/distance_parallel_smoke.rs b/crates/spatialrust/tests/distance_parallel_smoke.rs new file mode 100644 index 00000000..f4b810b9 --- /dev/null +++ b/crates/spatialrust/tests/distance_parallel_smoke.rs @@ -0,0 +1,23 @@ +#![cfg(feature = "metrics-distance")] + +use spatialrust::{chamfer_distance, cloud_distances, hausdorff_distance, PointCloudBuilder}; + +#[test] +fn facade_metrics_preserve_asymmetric_directed_statistics() { + let mut a = PointCloudBuilder::xyz(); + for x in [0.0, 2.0] { + a.push_point([x, 0.0, 0.0]).unwrap(); + } + let mut b = PointCloudBuilder::xyz(); + for x in [0.0, 1.0, 4.0] { + b.push_point([x, 0.0, 0.0]).unwrap(); + } + let (a, b) = (a.build().unwrap(), b.build().unwrap()); + let result = cloud_distances(&a, &b).unwrap(); + assert_eq!(result.mean_sq_a_to_b.to_bits(), 0.5_f64.to_bits()); + assert_eq!(result.mean_sq_b_to_a.to_bits(), (5.0_f64 / 3.0).to_bits()); + assert_eq!(result.max_a_to_b, 1.0); + assert_eq!(result.max_b_to_a, 2.0); + assert_eq!(chamfer_distance(&a, &b).unwrap(), 0.5 + 5.0 / 3.0); + assert_eq!(hausdorff_distance(&a, &b).unwrap(), 2.0); +} diff --git a/docs/ARCHITECTURE.md b/docs/ARCHITECTURE.md index 6e28d66b..e55de435 100644 --- a/docs/ARCHITECTURE.md +++ b/docs/ARCHITECTURE.md @@ -324,6 +324,19 @@ centroid/covariance arithmetic, eigenvalue ratios and the existing >= tie rule are preserved. Per-query neighbor vectors scale with radius support; worker caps remain per call rather than a scratch-byte or process-wide CPU budget. No floating-point reduction spans workers and no device transfer is introduced. +Cloud distances under `spatialrust-metrics/metrics-distance` build one KD-tree +per cloud and independently bound each direction's query workers with the same +point-staging policy. Parallel queries fill disjoint slices of a reusable window +of at most `workers * DEFAULT_SPATIAL_TENSOR_CHUNK_SIZE` optional f32 distances. +Each window joins before its values are converted to f64, added and square-rooted +in original point order; sums are never regrouped across workers or windows. +Missing neighbors remain distinct from returned NaNs. One-worker directions use +the unchanged serial loop and allocate no score window. The two directions run +sequentially, so their windows are not live together. Scratch scales with worker +budget rather than total query count; KD-tree storage remains input-sized, and +worker caps are per call rather than a process-wide CPU limit. Each active band +reuses one one-neighbor result vector with the same KD-tree traversal and fresh +per-query accumulator, avoiding a heap allocation for every parallel query. `spatialrust-image-io` depends on storage, never the reverse; standard codecs are additive, while TIFF and OpenEXR remain independently gated. Calibration datasets, robust solver controls, and residual receipts live in diff --git a/docs/FEATURE_MATRIX.md b/docs/FEATURE_MATRIX.md index 797bde79..6d63f1bc 100644 --- a/docs/FEATURE_MATRIX.md +++ b/docs/FEATURE_MATRIX.md @@ -218,6 +218,13 @@ workspace because its build requires a Python toolchain. collection. Equal maxima, selected point order, schema, all attribute types and metadata are preserved, without another dependency or device transfer. See [ISS detection measurements](../notes/2026-10-12_parallel_iss_keypoints.md). +- `metrics-distance` shares each target KD-tree across bounded CPU query workers + for Chamfer/Hausdorff statistics. A reusable optional-distance window holds at + most one 16,384-point chunk per worker; ordered f64 accumulation preserves + exact statistics without worker-local sums. Each band reuses its one-neighbor + query buffer. Small/one-worker directions keep + the original serial loop. Neither dependencies nor device transfers change. + See [cloud distance measurements](../notes/2026-10-12_parallel_cloud_distances.md). - `GpuStageProfiler` (under `gpu-wgpu`) records named pass timestamps through `run_gpu_vision_chain_profiled` (`gpu-image`) and `run_aoso_voxel_normal_frame_profiled` (`gpu-aoso-staging`). Create an optional diff --git a/docs/ROADMAP.md b/docs/ROADMAP.md index dd6d8844..e6cecae6 100644 --- a/docs/ROADMAP.md +++ b/docs/ROADMAP.md @@ -1097,6 +1097,7 @@ v1.3 admission gates. | Bounded CPU FPFH descriptors | Complete | Shared KD-tree, scoped SPFH/weighting phases and direct disjoint descriptor writes using existing point staging | 18,402,648 frozen exact descriptor words with 1/2/6 workers and public staging, independent planar/degenerate cases, validation/concurrency/seeded alignment, feature/facade/WASM checks and twelve completed-descriptor-API comparisons; `notes/2026-10-12_parallel_fpfh_descriptors.md` | | Bounded CPU boundary detection | Complete | Shared KD-tree, reusable worker angle scratch and direct disjoint mask writes using existing point staging | 1,081,992 frozen exact decisions with 1/2/6 workers and public staging, independent complete rim/degenerate/normal cases, all attribute types and metadata, validation/concurrency, feature/facade/WASM checks and twelve completed-mask-API comparisons; `notes/2026-10-12_parallel_boundary_detection.md` | | Bounded CPU ISS keypoints | Complete | Shared KD-tree, scoped covariance/selection phases and ordered selection with existing point staging | 1,081,992 frozen exact selection decisions and selected clouds with 1/2/6 workers and public staging, independent covariance/tied maxima/degenerate support, all attribute types and metadata, validation/concurrency, feature/facade/WASM checks and twelve completed-detection-API comparisons; `notes/2026-10-12_parallel_iss_keypoints.md` | +| Bounded CPU cloud distances | Complete | Shared target trees, reusable worker-bounded score windows and original-order f64 accumulation in both directions | 660 frozen exact metric words with 1/2/6 workers and public staging, independent asymmetric metrics/ties, multi-window order-sensitive sums, NaN/None/infinity and empty behavior, renamed schema/concurrency, feature/facade/WASM checks and twelve completed-distance-API comparisons; `notes/2026-10-12_parallel_cloud_distances.md` | | GPU spatial execution guards | Complete | Pre-allocation runtime/dimension/device/profiler checks for nearest resize and Sobel, including checked resize coordinate products | CPU/frozen full-RGBA parity, invalid-limit/capacity and ROI/concurrent/live-output/pool checks, image/point regressions and feature/WASM builds; `notes/2026-10-12_gpu_spatial_dispatch_validation.md` | | Bounded GPU gray dispatch | Complete | Balanced multi-row flat workgroups beyond the single-axis limit, independently lazy kernels and pre-allocation arithmetic/device/profiler validation | All 16,777,216 RGB values, 120 forced grid mappings, maximum-counter device probes, ownership/capacity/ROI/concurrency/reuse checks, native regressions and resident timing; `notes/2026-10-12_bounded_gpu_gray_dispatch.md` | | GPU AI execution guards | Complete | Pre-allocation image/index/buffer/dispatch and profiler validation while preserving the original CHW shader | 23,857,560 frozen output words across 192 cases, ROI/concurrency/live tensor reuse, negative-limit/capacity checks, 4096x4096 six-stage chain and feature/native/WASM validation; `notes/2026-10-12_validated_gpu_ai_packing.md` | diff --git a/notes/2026-10-12_parallel_cloud_distances.md b/notes/2026-10-12_parallel_cloud_distances.md new file mode 100644 index 00000000..39f358e7 --- /dev/null +++ b/notes/2026-10-12_parallel_cloud_distances.md @@ -0,0 +1,194 @@ +# Bounded CPU cloud distance metrics + +Date: 2026-10-12. This named slice follows +[PR #149](https://github.com/rsasaki0109/SpatialRust/pull/149), source +`e86cef4f3114dbd213f04e6243af4a898f931e79`. + +## Result and execution boundary + +Chamfer/Hausdorff statistics previously queried and accumulated every nearest +distance on the caller. Under `metrics-distance`, `cloud_distances` now shares +each target KD-tree across bounded scoped CPU query workers. It still builds +exactly one tree per side and computes a-to-b before b-to-a. Public Chamfer and +Hausdorff convenience functions use the same path; signatures and metric +definitions are unchanged. + +The implementation is +`C:/Users/rsasa/Workspace/SpatialRust/crates/spatialrust-metrics/src/distance.rs`. +Each direction independently reuses `spatialrust-search::parallel_worker_count`: +one worker below 4,096 points; otherwise at most available CPU parallelism and +`max(1, floor(query_points / 16,384))` useful chunks. With default chunking, +fewer than 32,768 query points retain serial execution. Worker injection is +private and used by tests. The one-worker directed loop matches the parent +verbatim, including its original query allocation behavior. + +Parallel workers fill disjoint slices of an optional f32 distance window. Its +capacity is at most `min(query_points, workers * 16,384)` entries; the product +uses saturating arithmetic before limiting it to the query length. Each window +joins completely, then its values are converted to f64 and accumulated in +original point order. The sum is not reset between windows and no worker-local +sums are combined. Division, square roots, directed maximum comparisons and +symmetric statistics preserve their original arithmetic. Missing neighbors stay +distinct from a returned NaN. Empty inputs retain the original error. + +Each active band reuses a one-neighbor vector through `KdTree::nearest_k_into` +with k=1. The existing `nearest_one` uses that same traversal via `nearest_k`; +the reusable API clears its result and initializes a fresh query accumulator. +This avoids a heap allocation per parallel query without changing backend +traversal, tie handling or distance math. KD-tree construction is unchanged. + +An optional distance occupies eight bytes on the measurement host. Parallel +score scratch is therefore at most 128 KiB per worker there, or 1.5 MiB at the +host's maximum twelve-worker budget, independent of total query count above +that window size. The two directions' windows are not live together. Each +worker also holds its one-neighbor vector and existing search state. The two +KD-trees remain input-sized. These per-call budgets do not establish a +process-wide CPU limit; concurrent callers have independent budgets. No +all-query neighborhood table, new dependency, global mutable state or device +transfer is introduced. + +## Correctness and feature validation + +The frozen parent module is +`C:/Users/rsasa/Workspace/SpatialRust/crates/spatialrust-metrics/tests/fixtures/distance_serial.rs`. +It matches source e86cef4 after line-ending/trailing-newline normalization, with +only the original test module omitted. It is compiled solely in native test +builds and retains the complete independent serial public metric APIs. + +Five tests in +`C:/Users/rsasa/Workspace/SpatialRust/crates/spatialrust-metrics/src/distance_parallel_tests.rs` +cover: + +- **660 exact f64 output words** across the four directed statistics, symmetric + Chamfer/Hausdorff methods and public wrappers. Inputs have + 1/17/4,095/4,096/4,103/32,767/32,768/32,771/65,537 points, planar/noisy curved/ + duplicate geometry and asymmetric target counts, with forced 1/2/6 workers + and public staging. This includes multi-window and partial-band boundaries. +- Independent asymmetric means/maxima for two points versus three, plus tied + nearest distances and duplicate points. Different worker counts per side, + zero requests and oversized requests preserve those exact expectations. +- A 65,537-query order-sensitive case with a first squared distance around + 1e16 and distinct later distances below one. Each later addition rounds away + in the original sum, while regrouping them changes the result. Forced two + workers cross two full windows and a one-point tail; public staging crosses + a full four-worker window and a one-point tail. Both retain frozen words. +- NaN payloads, infinities, squared-distance overflow and signed zero retain + frozen outputs in both directions. Returned NaNs remain distinct from + missing-neighbor rows on a private empty-tree query. Empty public inputs + retain frozen error strings through all three APIs. +- Renamed XYZ semantics, an ID attribute and metadata remain untouched during + concurrent calls with different per-side worker budgets. + +The ordered-sum fixture initially used tens of thousands of identical small +coordinates. Its pre-existing KD-tree equal-coordinate partitioning dominated +the test. Distinct small coordinates retain the same numerical stress and +window boundaries while keeping this correctness test focused on accumulation; +KD-tree partition optimization is outside this slice. + +Distance-feature tests pass: eight native tests and one ignored benchmark. One +facade smoke passes. All-target strict crate Clippy and facade smoke Clippy +pass. Disabling the sole metrics feature compiles without its optional search +dependency. Crate and facade WASM compilation passes; browser execution or +threading performance is not established. Formatting and diff checks pass. +No workflow or check names change; all 114 named CI checks gate merge. + +```text +cargo test -p spatialrust-metrics --no-default-features --features metrics-distance --lib -- --test-threads=1 --nocapture +cargo clippy -p spatialrust-metrics --no-default-features --features metrics-distance --all-targets -- -D warnings +cargo check -p spatialrust-metrics --no-default-features +cargo test -p spatialrust --no-default-features --features metrics-distance --test distance_parallel_smoke +cargo clippy -p spatialrust --no-default-features --features metrics-distance --test distance_parallel_smoke -- -D warnings +cargo check -p spatialrust-metrics --target wasm32-unknown-unknown --no-default-features --features metrics-distance +cargo check -p spatialrust --target wasm32-unknown-unknown --no-default-features --features metrics-distance +cargo fmt --all --check +git diff --check +``` + +## Completed distance-API measurements + +Release native Windows, Rust/Cargo 1.98.1. Twelve conditions combine +512/16,384/32,768/65,536 source points with planar, noisy curved and duplicate +curved geometry. The first two geometries have equal target counts; duplicate +sources use a unique noisy curved target with half as many points. The target +is translated by (0.017, -0.011, 0.023). This host reports twelve available CPU +workers. Default budgets are 1/1 below 32,768 points, 2/2 then 4/4 for equal +sides, and 2/1 then 4/2 for the large asymmetric conditions. + +The ignored diagnostic is +`C:/Users/rsasa/Workspace/SpatialRust/crates/spatialrust-metrics/src/distance_parallel_benchmark.rs`: + +```text +cargo test -p spatialrust-metrics --no-default-features --features metrics-distance --release --lib paired_previous_serial_parallel_distances -- --ignored --nocapture --test-threads=1 +``` + +Each condition runs two warmups and twenty rotating rounds of frozen serial +`cloud_distances`, new forced-one-worker execution and the current public API. +Each run verifies all six statistic words for all 792 completed warmup/measured +results. Timers include validation, both tree builds, both query directions, +window allocation/reuse, worker spawn/join, ordered f64 accumulation and +tree/scratch destruction. Fixture construction and exact-word verification are +outside the timer. Returned statistics have no heap allocation to drop. +No other tests or builds launched by this task overlap either run. Desktop +activity and other callers' CPU usage are not controlled. + +Table times are per-variant medians. Ratios are medians of twenty paired old/new +times, with 10,000 paired median bootstrap resamples (seed 20261012, sorted +indices 249 and 9749). Above one favors the new path. Intervals describe each +run's samples rather than portability or repeat-run guarantees. Cross-run raw +times are not a controlled comparison of the two staging implementations. + +### Initial staging with per-query result allocation + +The initial parallel implementation called `nearest_one` for every row. All +six large-condition paired medians exceeded one, but the 32,768 versus +16,384-point asymmetric case did not establish improvement. This survey is +retained in full. It motivated reusing the existing one-neighbor query buffer +inside each band; the final implementation was then retested and remeasured. + +| Points A | Geometry | Points B | Workers A/B | Old / forced one / public ms | Old / forced one (95% CI) | Old / public (95% CI) | +| ---: | --- | ---: | --- | ---: | ---: | ---: | +| 512 | plane | 512 | 1/1 | 0.247 / 0.246 / 0.246 | 1.001 (0.999–1.004) | 1.004 (0.992–1.007) | +| 512 | noisy curved | 512 | 1/1 | 0.305 / 0.301 / 0.299 | 1.033 (0.995–1.089) | 1.002 (0.986–1.058) | +| 512 | duplicate / unique curved | 256 | 1/1 | 0.255 / 0.294 / 0.257 | 0.983 (0.940–0.998) | 1.005 (0.994–1.016) | +| 16,384 | plane | 16,384 | 1/1 | 14.274 / 13.800 / 14.125 | 1.016 (0.991–1.038) | 1.008 (0.991–1.034) | +| 16,384 | noisy curved | 16,384 | 1/1 | 12.463 / 12.375 / 12.367 | 1.016 (0.992–1.024) | 1.008 (0.992–1.025) | +| 16,384 | duplicate / unique curved | 8,192 | 1/1 | 9.219 / 8.908 / 8.850 | 1.028 (1.006–1.043) | 1.021 (0.999–1.036) | +| 32,768 | plane | 32,768 | 2/2 | 25.596 / 25.075 / 20.888 | 1.009 (0.989–1.030) | 1.218 (1.129–1.368) | +| 32,768 | noisy curved | 32,768 | 2/2 | 25.434 / 25.479 / 23.395 | 1.006 (0.978–1.020) | 1.099 (1.065–1.205) | +| 32,768 | duplicate / unique curved | 16,384 | 2/1 | 18.400 / 18.457 / 18.205 | 0.996 (0.989–1.013) | 1.014 (0.990–1.095) | +| 65,536 | plane | 65,536 | 4/4 | 53.005 / 54.227 / 33.735 | 0.995 (0.974–1.002) | 1.567 (1.537–1.630) | +| 65,536 | noisy curved | 65,536 | 4/4 | 53.464 / 53.495 / 32.933 | 0.994 (0.962–1.036) | 1.641 (1.593–1.699) | +| 65,536 | duplicate / unique curved | 32,768 | 4/2 | 38.205 / 38.278 / 28.593 | 1.001 (0.989–1.018) | 1.373 (1.329–1.433) | + +### Final staging with per-band query-buffer reuse + +| Points A | Geometry | Points B | Workers A/B | Old / forced one / public ms | Old / forced one (95% CI) | Old / public (95% CI) | +| ---: | --- | ---: | --- | ---: | ---: | ---: | +| 512 | plane | 512 | 1/1 | 0.265 / 0.261 / 0.262 | 1.011 (1.004–1.105) | 1.011 (0.979–1.105) | +| 512 | noisy curved | 512 | 1/1 | 0.285 / 0.281 / 0.287 | 1.007 (0.928–1.221) | 1.008 (0.955–1.107) | +| 512 | duplicate / unique curved | 256 | 1/1 | 0.216 / 0.217 / 0.202 | 1.003 (0.938–1.082) | 1.000 (0.991–1.079) | +| 16,384 | plane | 16,384 | 1/1 | 12.376 / 12.382 / 12.528 | 0.996 (0.987–1.011) | 0.996 (0.979–1.004) | +| 16,384 | noisy curved | 16,384 | 1/1 | 12.335 / 12.457 / 12.376 | 0.995 (0.980–1.009) | 0.994 (0.987–1.003) | +| 16,384 | duplicate / unique curved | 8,192 | 1/1 | 9.020 / 8.901 / 9.009 | 1.006 (0.997–1.020) | 0.995 (0.976–1.011) | +| 32,768 | plane | 32,768 | 2/2 | 25.387 / 25.299 / 18.322 | 1.006 (0.996–1.021) | 1.390 (1.302–1.607) | +| 32,768 | noisy curved | 32,768 | 2/2 | 25.420 / 25.158 / 19.657 | 1.009 (0.988–1.032) | 1.305 (1.256–1.386) | +| 32,768 | duplicate / unique curved | 16,384 | 2/1 | 18.472 / 18.253 / 16.306 | 1.013 (1.003–1.023) | 1.138 (1.118–1.163) | +| 65,536 | plane | 65,536 | 4/4 | 55.003 / 57.278 / 31.419 | 1.004 (0.984–1.032) | 1.804 (1.763–1.978) | +| 65,536 | noisy curved | 65,536 | 4/4 | 53.961 / 53.467 / 29.345 | 1.006 (0.995–1.038) | 1.849 (1.817–1.897) | +| 65,536 | duplicate / unique curved | 32,768 | 4/2 | 37.726 / 37.602 / 24.518 | 1.001 (0.994–1.020) | 1.548 (1.523–1.683) | + +In the final run all six parallel public conditions favor staging, with paired +median speedups 1.138–1.849 and interval lower bounds above one. All six final +one-worker public intervals include one. The initial 512/256-point forced-one +control's adverse ratio 0.983 (0.940–0.998) is retained above. No small-input +neutrality, fixed speedup or controlled cross-run improvement is claimed. +Serial KD-tree construction and ordered accumulation still limit scaling; the +existing conservative worker policy is retained. The measurements cover the +complete host-returning metric API, not GPU execution or a registration pipeline. + +Raw samples, summaries, analysis and the initial implementation snapshot remain +private ignored files under +`C:/Users/rsasa/Workspace/SpatialRust/target/automation/parallel-distances-*`. +Published measurements describe execution conditions without a local adapter +model name. GPU metrics, parallel tree construction and parallel reassociation +of distance sums are outside this slice.