From d4148f39ea0d4ce34ab913c9873af98a746718a6 Mon Sep 17 00:00:00 2001 From: rsasaki0109 Date: Tue, 13 Oct 2026 10:37:34 +0900 Subject: [PATCH] Parallelize exact KD-tree construction and resident GICP preparation --- CHANGELOG.md | 8 + README.md | 8 + .../src/gicp_reduction.rs | 29 ++- .../src/gicp_tree_build_tests.rs | 184 ++++++++++++++++ crates/spatialrust-search/src/kdtree.rs | 11 + crates/spatialrust-search/src/kdtree_build.rs | 146 +++++++++++++ .../src/kdtree_build_benchmark.rs | 43 ++++ .../src/kdtree_build_tests.rs | 127 +++++++++++ crates/spatialrust-search/src/lib.rs | 3 + crates/spatialrust/src/lib.rs | 2 +- .../tests/parallel_kdtree_build_smoke.rs | 17 ++ docs/ARCHITECTURE.md | 10 + docs/FEATURE_MATRIX.md | 10 + docs/ROADMAP.md | 1 + notes/2026-10-13_parallel_kdtree_build.md | 202 ++++++++++++++++++ 15 files changed, 797 insertions(+), 4 deletions(-) create mode 100644 crates/spatialrust-registration/src/gicp_tree_build_tests.rs create mode 100644 crates/spatialrust-search/src/kdtree_build.rs create mode 100644 crates/spatialrust-search/src/kdtree_build_benchmark.rs create mode 100644 crates/spatialrust-search/src/kdtree_build_tests.rs create mode 100644 crates/spatialrust/tests/parallel_kdtree_build_smoke.rs create mode 100644 notes/2026-10-13_parallel_kdtree_build.md diff --git a/CHANGELOG.md b/CHANGELOG.md index e94216ba..f92ece19 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,14 @@ removed no sooner than the next major (see `docs/API_STABILITY.md`). ### Changed +- Add explicit bounded parallel KD-tree construction under search `parallel` + and facade `search-parallel`. Disjoint preorder node/point ranges retain exact + canonical partition order and search ties, without unsafe, per-worker trees + or a node merge. Native large inputs use available-CPU-bounded workers; + small/one-worker/WASM builds remain serial. Resident GPU GICP opts into the + same preparation with unchanged equations/transfer receipts. Ordinary tree + constructors and CPU GICP retain serial placement. + - Add explicit resident grid GICP covariance preparation under `register-gicp-gpu`. Reuse original GPU positions and canonical radius-grid arithmetic, eliminating source KD-tree preparation, per-point neighborhood diff --git a/README.md b/README.md index bf840467..fd5c2fc5 100644 --- a/README.md +++ b/README.md @@ -690,6 +690,14 @@ the existing cell/device limits; choose the bounded radius method explicitly for sparse extents. Only 256-byte equations are downloaded per evaluation. See [resident grid evidence](notes/2026-10-13_gpu_gicp_resident_grid_covariances.md). +`KdTree::from_slices_parallel` / `from_point_cloud_parallel` under +`spatialrust-search/parallel` (facade `search-parallel`) build large trees with +bounded scoped CPU workers while retaining exact point order, node layout and +search ties. Below 32,768 points, one-worker hosts and WASM keep serial builds; +ordinary constructors remain serial. Resident GPU GICP uses this preparation +without changing device data, equations or transfer receipts. See +[parallel tree evidence](notes/2026-10-13_parallel_kdtree_build.md). + Point-to-plane ICP also exposes `PointToPlaneIcp::align_with_gpu_reduction(&runtime, &source, &target)` under `register-icp-point-to-plane-gpu`. Transforms, nearest matching and plane diff --git a/crates/spatialrust-registration/src/gicp_reduction.rs b/crates/spatialrust-registration/src/gicp_reduction.rs index 604e3e29..35b6c655 100644 --- a/crates/spatialrust-registration/src/gicp_reduction.rs +++ b/crates/spatialrust-registration/src/gicp_reduction.rs @@ -70,6 +70,8 @@ enum CovariancePreparation { ResidentGrid(u32), } +type TreeBuilder = fn(&[f32], &[f32], &[f32]) -> KdTree; + impl GicpRegistration { /// Retains source/target/covariances on GPU for transform, nearest matching, /// rotated covariance inverse, Hessian/gradient and Euclidean fitness passes. @@ -163,6 +165,23 @@ impl GicpRegistration { source: &PointCloud, target: &PointCloud, preparation: CovariancePreparation, + ) -> SpatialResult { + self.align_reduction_with_tree_builder( + runtime, + source, + target, + preparation, + KdTree::from_slices_parallel, + ) + } + + fn align_reduction_with_tree_builder( + &self, + runtime: &WgpuRuntime, + source: &PointCloud, + target: &PointCloud, + preparation: CovariancePreparation, + build_tree: TreeBuilder, ) -> SpatialResult { self.validate_inputs(source, target)?; let c = self.config; @@ -208,7 +227,7 @@ impl GicpRegistration { if sx.iter().chain(sy).chain(sz).chain(tx).chain(ty).chain(tz).any(|v| !v.is_finite()) { return Err(SpatialError::InvalidArgument("GPU GICP requires finite XYZ".to_owned())); } - let tree = KdTree::from_slices(tx, ty, tz); + let tree = build_tree(tx, ty, tz); let mut transfers = TransferStats::default(); let mut covariance_batches = 0; let input = if let CovariancePreparation::ResidentGrid(capacity) = preparation { @@ -222,7 +241,7 @@ impl GicpRegistration { epsilon, )? } else if !matches!(preparation, CovariancePreparation::Legacy) { - let source_tree = KdTree::from_slices(sx, sy, sz); + let source_tree = build_tree(sx, sy, sz); covariance_batches = 2; match preparation { CovariancePreparation::ResidentKnn => { @@ -268,7 +287,7 @@ impl GicpRegistration { let tree = if let Some(tree) = tree { tree } else { - built = KdTree::from_slices(x, y, z); + built = build_tree(x, y, z); &built }; Ok(covariances(x, y, z, tree, c.k_neighbors, c.epsilon) @@ -355,3 +374,7 @@ impl GicpRegistration { }) } } + +#[cfg(all(test, not(target_arch = "wasm32")))] +#[path = "gicp_tree_build_tests.rs"] +mod tree_build_tests; diff --git a/crates/spatialrust-registration/src/gicp_tree_build_tests.rs b/crates/spatialrust-registration/src/gicp_tree_build_tests.rs new file mode 100644 index 00000000..40cf82ef --- /dev/null +++ b/crates/spatialrust-registration/src/gicp_tree_build_tests.rs @@ -0,0 +1,184 @@ +use super::*; +use crate::{transform_point_cloud, GicpConfig}; +use spatialrust_core::PointCloudBuilder; +use spatialrust_math::{Quat, Vec3}; + +fn corner(n: usize, scale: f32) -> PointCloud { + let mut b = PointCloudBuilder::xyz(); + for i in 0..n { + let a = (i % 17) as f32 * 0.125 * scale; + let c = (i / 17) as f32 * 0.125 * scale; + for p in [ + [a, c, 0.0], + [a, 0.0, c + 0.0625 * scale], + [0.0, a + 0.0625 * scale, c + 0.0625 * scale], + ] { + b.push_point(p).unwrap(); + } + } + b.build().unwrap() +} + +#[test] +#[ignore = "requires a native wgpu adapter; run explicitly with --include-ignored"] +fn parallel_tree_preparation_preserves_complete_gicp_diagnostics_exactly() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let mut checked = 0; + for n in [17, 63, 64, 65, 257, 10923] { + for scale in [1e-4_f32, 1.0, 1e4] { + let target = corner(n, scale); + let source = transform_point_cloud( + &target, + Isometry3::new( + Quat::from_axis_angle(Vec3::new(0.0, 0.0, 1.0), 0.002), + Vec3::new(0.01 * scale, -0.005 * scale, 0.008 * scale), + ), + ) + .unwrap(); + for mode in 0..4 { + for stop in 0..3 { + let registration = GicpRegistration::new(GicpConfig { + max_iterations: 3, + max_correspondence_distance: 0.25 * scale, + k_neighbors: 20, + covariance_radius: if mode >= 2 { Some(0.25 * scale) } else { None }, + transformation_epsilon: if stop == 1 { scale.max(1.0) as f64 } else { 0.0 }, + fitness_epsilon: if stop == 2 { f64::from(scale).powi(2) } else { 0.0 }, + ..Default::default() + }); + let prep = match mode { + 0 => CovariancePreparation::Legacy, + 1 => CovariancePreparation::ResidentKnn, + 2 => CovariancePreparation::ResidentRadius(128), + _ => CovariancePreparation::ResidentGrid(128), + }; + let before = registration + .align_reduction_with_tree_builder( + &runtime, + &source, + &target, + prep, + KdTree::from_slices, + ) + .unwrap(); + let actual = registration + .align_reduction_impl(&runtime, &source, &target, prep) + .unwrap(); + assert_eq!(actual, before, "n={n} scale={scale} mode={mode} stop={stop}"); + checked += source.len(); + } + } + } + } + eprintln!( + "PARALLEL_GICP_TREE_BUILD exact_diagnostics_conditions=216 compared_points={checked}" + ); +} + +fn fixture(n: usize, profile: usize) -> PointCloud { + let mut b = PointCloudBuilder::xyz(); + for i in 0..n { + let j = if profile == 2 { i / 2 } else { i }; + b.push_point([(j % 32) as f32 * 0.5, (j / 32 % 32) as f32 * 0.5, (j / 1024) as f32 * 0.5]) + .unwrap(); + } + b.build().unwrap() +} + +#[test] +#[ignore = "completed GICP calls with serial/parallel tree preparation; run this exact test alone in release mode"] +fn completed_parallel_tree_grid_gicp_calls() { + use std::time::Instant; + let runtime = WgpuRuntime::new_headless().unwrap(); + let info = runtime.adapter_info(); + println!( + "PARALLEL_GICP_TREE_ENV backend={} device_type={} available_workers={}", + info.backend, + info.device_type, + std::thread::available_parallelism().map_or(1, |n| n.get()) + ); + for n in [4097, 32768, 65537, 262145] { + for profile in 0..3 { + let cloud = fixture(n, profile); + let target = fixture(n + if profile == 2 { 17 } else { 0 }, profile); + let source = transform_point_cloud( + &cloud, + Isometry3::new(Quat::::identity(), Vec3::new(0.03125, -0.015625, 0.0078125)), + ) + .unwrap(); + let registration = GicpRegistration::new(GicpConfig { + max_iterations: 3, + max_correspondence_distance: 0.25, + covariance_radius: Some([0.51, 0.8, 1.05][profile]), + transformation_epsilon: 0.0, + fitness_epsilon: 0.0, + ..Default::default() + }); + let capacity = [32, 64, 128][profile]; + let prep = CovariancePreparation::ResidentGrid(capacity); + let expected = registration + .align_reduction_with_tree_builder( + &runtime, + &source, + &target, + prep, + KdTree::from_slices, + ) + .unwrap(); + for _ in 0..3 { + assert_eq!( + registration + .align_reduction_with_tree_builder( + &runtime, + &source, + &target, + prep, + KdTree::from_slices + ) + .unwrap(), + expected + ); + assert_eq!( + registration + .align_with_gpu_resident_grid_covariances( + &runtime, &source, &target, capacity + ) + .unwrap(), + expected + ); + } + let (mut serial_ms, mut parallel_ms) = (vec![], vec![]); + for round in 0..20 { + for offset in 0..2 { + let mode = (round + offset) % 2; + let start = Instant::now(); + let result = if mode == 0 { + registration + .align_reduction_with_tree_builder( + &runtime, + &source, + &target, + prep, + KdTree::from_slices, + ) + .unwrap() + } else { + registration + .align_with_gpu_resident_grid_covariances( + &runtime, &source, &target, capacity, + ) + .unwrap() + }; + let ms = start.elapsed().as_secs_f64() * 1000.0; + if mode == 0 { + serial_ms.push(ms); + } else { + parallel_ms.push(ms); + } + assert_eq!(result, expected); + } + } + println!("PARALLEL_GICP_TREE_BENCH n={n} profile={profile} iterations={} h2d={} d2h={} fitness={} serial_ms={serial_ms:?} parallel_ms={parallel_ms:?}",expected.result.iterations,expected.transfers.host_to_device_bytes(),expected.transfers.device_to_host_bytes(),expected.result.fitness); + } + } +} diff --git a/crates/spatialrust-search/src/kdtree.rs b/crates/spatialrust-search/src/kdtree.rs index 1f4dfde2..e22f2e58 100644 --- a/crates/spatialrust-search/src/kdtree.rs +++ b/crates/spatialrust-search/src/kdtree.rs @@ -12,6 +12,13 @@ const SEARCH_STACK_SIZE: usize = 64; const AXIS_LEAF: u8 = u8::MAX; const INVALID_NODE: u32 = u32::MAX; +#[cfg(feature = "parallel")] +#[path = "kdtree_build.rs"] +mod parallel_build; + +#[cfg(feature = "parallel")] +pub use parallel_build::PARALLEL_KDTREE_BUILD_MIN_POINTS; + /// Cache-friendly KD-tree for 3D point clouds. #[derive(Clone, Debug)] pub struct KdTree { @@ -723,6 +730,10 @@ mod nearest_tests; #[path = "kdtree_nearest_benchmark.rs"] mod nearest_benchmark; +#[cfg(all(test, feature = "parallel", not(target_arch = "wasm32")))] +#[path = "kdtree_build_benchmark.rs"] +mod build_benchmark; + fn partition_by_axis( x: &[f32], y: &[f32], diff --git a/crates/spatialrust-search/src/kdtree_build.rs b/crates/spatialrust-search/src/kdtree_build.rs new file mode 100644 index 00000000..d78220f8 --- /dev/null +++ b/crates/spatialrust-search/src/kdtree_build.rs @@ -0,0 +1,146 @@ +//! Disjoint subtree construction with the unchanged partition/tie convention. +use super::*; + +/// Minimum points before the explicit KD-tree builder uses scoped CPU workers. +pub const PARALLEL_KDTREE_BUILD_MIN_POINTS: usize = 32_768; +#[cfg(not(target_arch = "wasm32"))] +const MIN_WORKER_POINTS: usize = PARALLEL_KDTREE_BUILD_MIN_POINTS / 2; + +impl KdTree { + /// Builds with bounded scoped CPU workers under the `parallel` feature. + /// Point permutation, split bits, preorder node layout and query ties match + /// `from_slices` exactly. Native workers per call are bounded by available + /// parallelism and one per 16,384 points; inputs below 32,768 points or one + /// useful worker keep the serial constructor. WASM also stays serial. + /// Axes are copied once; one exact-sized node array is split into disjoint + /// subtree ranges without locks, unsafe, per-worker trees or a merge copy. + /// This is a per-call bound, not a process-wide thread/stack limit. + /// Panics on unequal axis lengths, or counts exceeding u32 original indices. + #[must_use] + pub fn from_slices_parallel(x: &[f32], y: &[f32], z: &[f32]) -> Self { + if x.len() < PARALLEL_KDTREE_BUILD_MIN_POINTS { + return Self::from_slices(x, y, z); + } + #[cfg(not(target_arch = "wasm32"))] + { + let workers = std::thread::available_parallelism().map_or(1, |n| n.get()); + from_slices_with_workers(x, y, z, workers) + } + #[cfg(target_arch = "wasm32")] + { + assert!(x.len() <= u32::MAX as usize, "KD-tree original indices must fit u32"); + Self::from_slices(x, y, z) + } + } + + /// Builds an explicit parallel tree from a cloud, retaining position order. + pub fn from_point_cloud_parallel(cloud: &PointCloud) -> SpatialResult { + let (x, y, z) = cloud.positions3()?; + Ok(Self::from_slices_parallel(x, y, z)) + } +} + +#[cfg(not(target_arch = "wasm32"))] +fn worker_count(points: usize, available: usize) -> usize { + if points < PARALLEL_KDTREE_BUILD_MIN_POINTS { + 1 + } else { + available.max(1).min(points / MIN_WORKER_POINTS) + } +} + +// All levels above a power-of-two frontier split. At that frontier, bucket +// sizes are floor(N/frontier) or ceil(N/frontier). With the power-of-two leaf +// size, either every bucket splits once, or only N-LEAF_SIZE*frontier do. +// This gives exact preorder subtree sizes in O(1), including odd boundaries. +#[cfg(not(target_arch = "wasm32"))] +fn node_count(points: usize) -> usize { + if points == 0 { + return 0; + } + if points <= LEAF_SIZE { + return 1; + } + debug_assert!(LEAF_SIZE.is_power_of_two()); + let quotient = points / LEAF_SIZE; + let frontier = 1usize << (usize::BITS - 1 - quotient.leading_zeros()); + if points / frontier == LEAF_SIZE { + 2 * frontier - 1 + 2 * (points - LEAF_SIZE * frontier) + } else { + 4 * frontier - 1 + } +} + +#[cfg(not(target_arch = "wasm32"))] +fn from_slices_with_workers(x: &[f32], y: &[f32], z: &[f32], available: usize) -> KdTree { + assert_eq!(x.len(), y.len()); + assert_eq!(x.len(), z.len()); + assert!(x.len() <= u32::MAX as usize, "KD-tree original indices must fit u32"); + let workers = worker_count(x.len(), available); + if workers == 1 { + return KdTree::from_slices(x, y, z); + } + let mut points_order: Vec = (0..x.len() as u32).collect(); + let empty = + KdNode { split: 0.0, left: INVALID_NODE, right: INVALID_NODE, start: 0, end: 0, axis: 0 }; + let mut nodes = vec![empty; node_count(x.len())]; + build_into([x, y, z], &mut points_order, &mut nodes, 0, 0, workers); + KdTree { x: x.to_vec(), y: y.to_vec(), z: z.to_vec(), points_order, nodes, root: 0 } +} + +#[cfg(not(target_arch = "wasm32"))] +fn build_into( + xyz: [&[f32]; 3], + order: &mut [u32], + nodes: &mut [KdNode], + point_offset: usize, + node_offset: usize, + workers: usize, +) { + let count = order.len(); + let (node, children) = nodes.split_first_mut().expect("nonempty subtree node range"); + node.start = point_offset as u32; + node.end = (point_offset + count) as u32; + if count <= LEAF_SIZE { + node.axis = AXIS_LEAF; + return; + } + let [x, y, z] = xyz; + let (axis, minimum) = select_axis(x, y, z, order, 0, count); + let mid = count / 2; + select_nth_by_axis(x, y, z, order, 0, count, (axis, minimum), mid); + node.axis = axis; + node.split = coordinate(x, y, z, order[mid], axis); + let left_count = node_count(mid); + let left_offset = node_offset + 1; + let right_offset = left_offset + left_count; + node.left = left_offset as u32; + node.right = right_offset as u32; + let (left_order, right_order) = order.split_at_mut(mid); + let (left_nodes, right_nodes) = children.split_at_mut(left_count); + if workers == 1 || count < PARALLEL_KDTREE_BUILD_MIN_POINTS { + build_into(xyz, left_order, left_nodes, point_offset, left_offset, 1); + build_into(xyz, right_order, right_nodes, point_offset + mid, right_offset, 1); + } else { + let left_workers = workers / 2; + let right_workers = workers - left_workers; + std::thread::scope(|scope| { + let left = scope.spawn(move || { + build_into(xyz, left_order, left_nodes, point_offset, left_offset, left_workers) + }); + build_into( + xyz, + right_order, + right_nodes, + point_offset + mid, + right_offset, + right_workers, + ); + left.join().expect("bounded KD-tree build worker panicked"); + }); + } +} + +#[cfg(all(test, not(target_arch = "wasm32")))] +#[path = "kdtree_build_tests.rs"] +mod tests; diff --git a/crates/spatialrust-search/src/kdtree_build_benchmark.rs b/crates/spatialrust-search/src/kdtree_build_benchmark.rs new file mode 100644 index 00000000..29c4bab5 --- /dev/null +++ b/crates/spatialrust-search/src/kdtree_build_benchmark.rs @@ -0,0 +1,43 @@ +use super::{ + partition_tests::{fixture, snapshot}, + *, +}; +use std::time::Instant; + +#[test] +#[ignore = "completed KD-tree construction benchmark; run this exact test alone in release mode"] +fn completed_parallel_kdtree_builds() { + let available = std::thread::available_parallelism().map_or(1, |n| n.get()); + println!("PARALLEL_KDTREE_BUILD_ENV available_workers={available} threshold={PARALLEL_KDTREE_BUILD_MIN_POINTS}"); + for n in [4097, 32768, 65537, 262145] { + for kind in [0, 3, 8] { + let [x, y, z] = fixture(n, kind); + let expected = snapshot(&KdTree::from_slices(&x, &y, &z)); + for _ in 0..3 { + assert_eq!(snapshot(&KdTree::from_slices(&x, &y, &z)), expected); + assert_eq!(snapshot(&KdTree::from_slices_parallel(&x, &y, &z)), expected); + } + let (mut serial_ms, mut parallel_ms) = (vec![], vec![]); + for round in 0..20 { + for offset in 0..2 { + let mode = (round + offset) % 2; + let start = Instant::now(); + let tree = if mode == 0 { + KdTree::from_slices(&x, &y, &z) + } else { + KdTree::from_slices_parallel(&x, &y, &z) + }; + let ms = start.elapsed().as_secs_f64() * 1000.0; + if mode == 0 { + serial_ms.push(ms); + } else { + parallel_ms.push(ms); + } + assert_eq!(snapshot(&tree), expected); + std::hint::black_box(&tree); + } + } + println!("PARALLEL_KDTREE_BUILD_BENCH n={n} kind={kind} serial_ms={serial_ms:?} parallel_ms={parallel_ms:?}"); + } + } +} diff --git a/crates/spatialrust-search/src/kdtree_build_tests.rs b/crates/spatialrust-search/src/kdtree_build_tests.rs new file mode 100644 index 00000000..4163535c --- /dev/null +++ b/crates/spatialrust-search/src/kdtree_build_tests.rs @@ -0,0 +1,127 @@ +use super::*; +use crate::kdtree::{ + partition_tests::{fixture, neighbor_words, snapshot}, + previous, +}; + +#[test] +fn node_ranges_and_worker_policy_are_exact_and_bounded() { + let mut reference = vec![0usize; 131100]; + for n in 1..reference.len() { + reference[n] = if n <= LEAF_SIZE { 1 } else { 1 + reference[n / 2] + reference[n - n / 2] }; + assert_eq!(node_count(n), reference[n], "n={n}"); + } + for n in [0, 1, 16384, 32767, 32768, 32769, 65537, 131072, usize::MAX] { + for available in [0, 1, 2, 6, usize::MAX] { + let workers = worker_count(n, available); + assert!(workers >= 1 && workers <= available.max(1)); + if n < 32768 { + assert_eq!(workers, 1); + } else { + assert!(workers <= n / 16384); + } + } + } + // The original-index limit also keeps every computed node ID below u32. + assert!(node_count(u32::MAX as usize) < u32::MAX as usize); +} + +fn check_queries(actual: &KdTree, serial: &KdTree) -> usize { + let mut checked = 0; + for q in [[0.0; 3], [1.0; 3], [0.25, 0.5, 0.0], [f32::NAN, 0.0, 0.0], [f32::INFINITY, 0.0, 0.0]] + { + for k in [0, 1, 3, 17, 65] { + assert_eq!( + neighbor_words(actual.nearest_k(q[0], q[1], q[2], k)), + neighbor_words(serial.nearest_k(q[0], q[1], q[2], k)) + ); + checked += 1; + } + for r in [0.0, 0.1, 2.0, f32::INFINITY, f32::NAN] { + assert_eq!( + neighbor_words(actual.radius_search(q[0], q[1], q[2], r)), + neighbor_words(serial.radius_search(q[0], q[1], q[2], r)) + ); + checked += 1; + } + for bound in [-1.0, 0.0, 1.0, f32::INFINITY, f32::NAN] { + let bits = |v: Option| v.map(|n| (n.index, n.distance_squared.to_bits())); + assert_eq!( + bits(actual.nearest_one_within(q[0], q[1], q[2], bound)), + bits(serial.nearest_one_within(q[0], q[1], q[2], bound)) + ); + checked += 1; + } + } + checked +} + +#[test] +fn parallel_trees_match_frozen_layout_words_coordinates_and_queries() { + let mut layouts = 0; + let mut queries = 0; + for n in [0, 1, 15, 16, 17, 31, 32, 33, 257, 1025, 32767, 32768, 32769, 65537, 131073] { + for kind in 0..9 { + let [x, y, z] = fixture(n, kind); + let serial = KdTree::from_slices(&x, &y, &z); + let expected = snapshot(&serial); + if n <= 1025 { + assert_eq!( + expected, + previous::snapshot(&previous::KdTree::from_slices(&x, &y, &z)) + ); + } + for workers in [1, 2, 6] { + let actual = from_slices_with_workers(&x, &y, &z, workers); + assert_eq!(snapshot(&actual), expected, "n={n} kind={kind} workers={workers}"); + for (a, b) in + [actual.x.as_slice(), &actual.y, &actual.z].into_iter().zip([&x, &y, &z]) + { + assert_eq!( + a.iter().map(|v| v.to_bits()).collect::>(), + b.iter().map(|v| v.to_bits()).collect::>() + ); + } + queries += check_queries(&actual, &serial); + layouts += 1; + } + } + } + eprintln!("PARALLEL_KDTREE_BUILD exact_layouts={layouts} exact_query_results={queries}"); +} + +#[test] +fn public_cloud_builder_preserves_inputs_and_concurrent_calls() { + let [x, y, z] = fixture(65537, 8); + let before = [x.clone(), y.clone(), z.clone()]; + let expected = snapshot(&KdTree::from_slices(&x, &y, &z)); + std::thread::scope(|scope| { + let jobs: Vec<_> = (0..4) + .map(|_| { + scope.spawn(|| { + for _ in 0..2 { + assert_eq!(snapshot(&KdTree::from_slices_parallel(&x, &y, &z)), expected); + } + }) + }) + .collect(); + for job in jobs { + job.join().unwrap(); + } + }); + assert_eq!([x, y, z], before); + let mut b = spatialrust_core::PointCloudBuilder::xyz(); + for p in [[0.0; 3], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]] { + b.push_point(p).unwrap(); + } + let cloud = b.build().unwrap(); + assert_eq!( + snapshot(&KdTree::from_point_cloud_parallel(&cloud).unwrap()), + snapshot(&KdTree::from_point_cloud(&cloud).unwrap()) + ); + assert!(std::panic::catch_unwind(|| KdTree::from_slices_parallel(&[0.0], &[], &[0.0])).is_err()); + // Independent geometry expectation, including inclusive bound/self ties. + let t = KdTree::from_slices_parallel(&[0.0, 1.0, 1.0], &[0.0; 3], &[0.0; 3]); + assert_eq!(t.nearest_one_within(0.5, 0.0, 0.0, 0.25).unwrap().index, 0); + assert_eq!(t.radius_search(0.0, 0.0, 0.0, 1.0).len(), 3); +} diff --git a/crates/spatialrust-search/src/lib.rs b/crates/spatialrust-search/src/lib.rs index 06a401b2..b5f3c212 100644 --- a/crates/spatialrust-search/src/lib.rs +++ b/crates/spatialrust-search/src/lib.rs @@ -45,6 +45,9 @@ pub use parallel_chunked::{ #[cfg(feature = "search-kdtree")] pub use kdtree::{KdTree, KdTreeLayout, KdTreeNode}; +#[cfg(all(feature = "search-kdtree", feature = "parallel"))] +pub use kdtree::PARALLEL_KDTREE_BUILD_MIN_POINTS; + #[cfg(feature = "search-graph")] pub use graph::{knn_graph, radius_graph, NeighborGraph}; diff --git a/crates/spatialrust/src/lib.rs b/crates/spatialrust/src/lib.rs index df93d16b..561992ac 100644 --- a/crates/spatialrust/src/lib.rs +++ b/crates/spatialrust/src/lib.rs @@ -138,7 +138,7 @@ pub use spatialrust_search::{ pub use spatialrust_search::{ nearest_k_spatial_tensor_parallel, nearest_k_spatial_tensor_parallel_into, radius_search_spatial_tensor_parallel, radius_search_spatial_tensor_parallel_into, - PARALLEL_CHUNK_QUERY_MIN_POINTS, + PARALLEL_CHUNK_QUERY_MIN_POINTS, PARALLEL_KDTREE_BUILD_MIN_POINTS, }; #[cfg(feature = "search-graph")] diff --git a/crates/spatialrust/tests/parallel_kdtree_build_smoke.rs b/crates/spatialrust/tests/parallel_kdtree_build_smoke.rs new file mode 100644 index 00000000..2ed5e536 --- /dev/null +++ b/crates/spatialrust/tests/parallel_kdtree_build_smoke.rs @@ -0,0 +1,17 @@ +#![cfg(feature = "search-parallel")] +#[test] +fn facade_exposes_exact_parallel_tree_construction() { + use spatialrust::{KdTree, PointCloudBuilder, PARALLEL_KDTREE_BUILD_MIN_POINTS}; + let mut b = PointCloudBuilder::xyz(); + for i in 0..PARALLEL_KDTREE_BUILD_MIN_POINTS + 1 { + b.push_point([(i % 127) as f32, (i / 127) as f32, 0.0]).unwrap(); + } + let cloud = b.build().unwrap(); + let serial = KdTree::from_point_cloud(&cloud).unwrap(); + let parallel = KdTree::from_point_cloud_parallel(&cloud).unwrap(); + assert_eq!(serial.layout().point_order(), parallel.layout().point_order()); + assert_eq!( + serial.layout().nodes().collect::>(), + parallel.layout().nodes().collect::>() + ); +} diff --git a/docs/ARCHITECTURE.md b/docs/ARCHITECTURE.md index 4f04e1ce..b7ba4a49 100644 --- a/docs/ARCHITECTURE.md +++ b/docs/ARCHITECTURE.md @@ -93,6 +93,16 @@ bytes per grid offset entry. CPU grid construction/target tree/6x6 solve remain; dense-grid cell limits reject without automatically selecting another route. See [resident grid evidence](../notes/2026-10-13_gpu_gicp_resident_grid_covariances.md). +The optional CPU search `parallel` feature exposes explicit KD-tree construction +with bounded scoped subtree workers. One exact-sized preorder node array and +the point permutation are split into disjoint ranges, with no unsafe, locks, +per-worker trees or node merge. Canonical axis/partition code preserves all +layout/query bits. Ordinary constructors stay serial; small/one-worker/WASM +parallel calls use them directly. Resident GPU GICP opts into this host +preparation without altering device staging or equations. No production +dependency or core model change is introduced. See +[parallel tree evidence](../notes/2026-10-13_parallel_kdtree_build.md). + `PointToPlaneIcp::align_with_gpu_reduction` under the optional `register-icp-point-to-plane-gpu` feature retains original positions/tree and raw target normals. Five passes transform originals, perform the existing diff --git a/docs/FEATURE_MATRIX.md b/docs/FEATURE_MATRIX.md index 34982664..0845db2c 100644 --- a/docs/FEATURE_MATRIX.md +++ b/docs/FEATURE_MATRIX.md @@ -431,6 +431,16 @@ workspace because its build requires a Python toolchain. work. Accepted f32 bounds, the exact cell cap and positive-infinite-radius one-cell behavior remain compatible. See [grid validation regression receipt](../notes/2026-10-12_validated_grid_bounds.md). +- `spatialrust-search/parallel` adds explicit `KdTree::from_slices_parallel` + and `from_point_cloud_parallel`, also exposed under facade `search-parallel`. + Native subtrees use disjoint point/node ranges in one preallocated preorder + array, reusing canonical axis/partition code. Layout, split bits and query + ties remain exact. From 32,768 points, workers are bounded by available CPU + parallelism and one per 16,384 points; small/one-worker/WASM calls stay serial. + Ordinary constructors retain serial placement. Resident GICP GPU methods + opt into this host preparation; GPU math and receipts are unchanged. This + is a per-call worker bound, not a global thread/stack budget. See + [parallel tree measurements](../notes/2026-10-13_parallel_kdtree_build.md). - Shared CPU `build_grid` uses stable reverse scattering and one cell-offset table, eliminating the auxiliary cursor clone. Logical live heap payload drops by `4 * (cells + 1)` bytes; returned CSR order/sentinels and GPU uploads diff --git a/docs/ROADMAP.md b/docs/ROADMAP.md index 97eb3a12..a1e9ac8f 100644 --- a/docs/ROADMAP.md +++ b/docs/ROADMAP.md @@ -1079,6 +1079,7 @@ v1.3 admission gates. | Temporary storage | Complete | Private submitted scratch leases plus byte-bounded buffer/texture retention with explicit controls and stats | repeated-frame allocation/retained-byte measurements, CPU/GPU parity with stale scratch and live outputs, concurrent submissions, abandoned encoders, resolution/size churn and zero/shrinking budgets; `notes/2026-10-11_bounded_gpu_temporary_storage.md` | | CPU contrast enhancement | Complete | Bounded parallel CLAHE tile/row work and equalization lookup application with shared horizontal geometry | Exact 1/2/6-worker ROI, metadata, odd-grid, clipping, rounding and degenerate-image parity; checked grid allocation; before/after and one-worker Criterion comparison; `notes/2026-10-11_parallel_contrast_enhancement.md` | | Bounded CPU point queries | Complete | Group original radius/k-NN query chunks into balanced, available-CPU-bounded scoped workers with worker-sized dispatch metadata | Unsorted serial/append parity, exact callback boundaries, deterministic 1/2/6-worker and tiny-chunk limits, minimal feature/WASM checks and completed-work comparison; `notes/2026-10-11_bounded_parallel_point_queries.md` | +| Bounded CPU KD-tree construction | Complete | Explicit parallel constructor, exact preallocated preorder ranges and unchanged partitions; resident GPU GICP host preparation | Frozen/layout/query-bit checks with 1/2/6 workers, 216 exact GICP diagnostics conditions, feature/facade/WASM checks and completed tree/registration comparisons; `notes/2026-10-13_parallel_kdtree_build.md` | | GPU image histogram | Complete | Exact two-pass logical-channel histograms with bounded private partials, resident bin ownership, explicit native readback and optional shared timestamps | Packed/strided/ROI/channel and hot-bin CPU parity, stale scratch/live output/concurrency ownership, limit/overflow checks, GPU regressions and completed-work/transfer measurements; `notes/2026-10-11_resident_gpu_image_histograms.md` | | GPU histogram CDF | Complete | Work-efficient inclusive 256-bin integer scan into an independently owned resident result, with explicit readback and one-pass profiling | Exact u32::MAX counts, all logical channels/ROI and grayscale chains, source recycling/live output/concurrency ownership, cumulative receipts and completed-work placement comparisons; `notes/2026-10-11_resident_gpu_image_cdf.md` | | GPU histogram equalization | Complete | Resident grayscale histogram/CDF/lookup/application with CPU-compatible binary64 rounding and bounded temporary ownership | Constant-byte/metadata/ROI and RGB-gray chain parity, 2,309,120 independent LUT comparisons through maximum counters, pool reuse/concurrency, five-pass timing and completed-work placement measurements; `notes/2026-10-11_resident_gpu_histogram_equalization.md` | diff --git a/notes/2026-10-13_parallel_kdtree_build.md b/notes/2026-10-13_parallel_kdtree_build.md new file mode 100644 index 00000000..a4ebd0d0 --- /dev/null +++ b/notes/2026-10-13_parallel_kdtree_build.md @@ -0,0 +1,202 @@ +# Bounded parallel KD-tree construction — 2026-10-13 + +`KdTree::from_slices_parallel` / `from_point_cloud_parallel` under search +`parallel` (facade `search-parallel`) construct large CPU trees with bounded +scoped subtree workers. GPU GICP's resident reduction/covariance methods use +this host preparation; ordinary tree constructors and ordinary CPU GICP keep +serial placement. GPU algorithms, input staging, equations and receipts are +unchanged. This is CPU parallelism supporting GPU preparation, not device-side +KD-tree construction. + +## Structure, exactness and worker boundaries + +Canonical axis selection and minimum-pivot partition code are reused. Absolute +preorder node/point ranges are calculated before sibling work. One exact-sized +node array and the point permutation are split into disjoint mutable slices; +there are no per-worker trees, locks, unsafe or node-merging copies. Each +subtree scans/partitions its points in the original order. Local range indices +translate exactly to old global pivot indices, retaining all tie permutations. +Original coordinate arrays are copied once after construction. + +For the existing power-of-two leaf size 16, let F be the highest power of two +not exceeding floor(N/16). Levels above F fully split. If floor(N/F)=16, +node count is 2*F-1+2*(N-16*F); otherwise it is 4*F-1. Counts <=16 use one +node, empty inputs zero. This O(1) recurrence gives exact preorder sibling +ranges, including odd boundaries. The planning formula is checked against an +independent recursive-count table for 131,099 sizes and the u32 index boundary. + +From **32,768 points**, workers are bounded by available CPU parallelism and +floor(N/16,384). Root selection/partition remains serial, then child budgets +split between a scoped helper and the caller. Every spawning subtree has at +least 32,768 points, so children have at least 16,384. Concurrent helpers are +bounded by the supplied budget; completed calls join them. Below the threshold, +one useful worker or WASM keeps the original serial constructor. Tiny public +calls do not query CPU availability. Bounds apply per call, not process-wide; +multiple callers can still run concurrently. Unequal axes and >u32 original +indices panic, as documented for the explicit constructors. + +Resident GICP's internal tree-builder boundary permits a serial control for +verification/measurement. Production selects the explicit parallel constructor +for target and required source trees. The post-preparation pose/evaluation loop +is unchanged. Both preparation paths produce the same tree layout and upload +bytes; no new device copies or downloads are introduced. No production +thread-pool dependency, core model or GPU shader change is added. + +## Correctness evidence + +- **405** exact full layouts with forced 1/2/6 worker budgets span empty, + 1/15/16/17/31/32/33/257/1025, threshold/odd 32767/32768/32769, 65537 and + 131073 points. Nine morphologies include duplicates/outliers/planes/lines, + NaN payloads/infinities/signed zeros and distributed finite coordinates. + Point order, all split/node words and copied coordinate bits match serial. + Small cases also match the frozen pre-optimization serial builder. +- **30,375** ordered k-NN/radius/gated-nearest results match bit for bit, + including nonfinite queries, gates and ties. Independent midpoint-tie and + inclusive-radius expectations, immutable source arrays, malformed axes, + cloud construction and four concurrent public callers are checked. +- **216** completed native GICP conditions compare the previous serial tree + preparation and new public path: six sizes, three scales, four covariance + modes (legacy, resident k-NN, bounded radius, grid) and three stopping modes. + **1,230,012** compared source points retain exact complete diagnostics: + poses, fitness, histories, counts, stopping and transfer receipts. The + largest fixture crosses the parallel threshold. GPU arithmetic is unchanged. + +## Completed measurements + +Two isolated release surveys per family, three warmups and 20 alternating +paired samples per condition. Construction measures complete APIs, including +node/point allocation, planning, axis/partition work, coordinate copies and +thread creation/joins, excluding output drop and validation snapshots. Kinds +0/3/8 are coincident points, curved patches and distributed finite points. + +Registration measures complete resident-grid GICP calls using the same +implementation/configuration except serial versus parallel tree construction. +Runtime initialization/fixture creation is excluded; CPU tree/grid preparation, +uploads, GPU covariance/equation work, waits/readback, CPU 6x6 solves/composition +and final fitness are included. Sizes 4097/32768/65537/262145 use fixed 0.5 +spacing on a 32x32 grid; profiles 0/1 have radii 0.51/0.8 and capacities 32/64. +Profile 2 duplicates grid positions, adds 17 targets, radius 1.05 and capacity +128. Binary translation `(0.03125,-0.015625,0.0078125)`, distance gate 0.25, +epsilon 0.001, three updates and zero stopping thresholds are fixed. All warm +and timed outputs match their serial references exactly, including receipts. +No other local builds/tests/statistical work runs during surveys. + +Native CPU availability is 12 workers; policy budgets are 1/2/4/12 for the four +sizes. Registration uses Vulkan/DiscreteGpu without recording the adapter +model. These are consecutive local surveys, not cross-device/day evidence. + +### construction survey 1 + +| N | Kind/profile | Serial ms | Parallel ms | Serial/parallel [95%] | +| ---: | ---: | ---: | ---: | --- | +| 4,097 | 0 | 0.219 | 0.223 | 0.99 [0.98–1.00] | +| 4,097 | 3 | 0.521 | 0.527 | 1.00 [0.98–1.01] | +| 4,097 | 8 | 0.567 | 0.569 | 1.00 [0.99–1.01] | +| 32,768 | 0 | 2.241 | 1.565 | 1.43 [1.38–1.47] | +| 32,768 | 3 | 4.033 | 2.620 | 1.47 [1.41–1.54] | +| 32,768 | 8 | 7.091 | 4.350 | 1.60 [1.50–1.69] | +| 65,537 | 0 | 4.862 | 2.184 | 2.15 [2.07–2.22] | +| 65,537 | 3 | 8.692 | 3.459 | 2.55 [2.45–2.64] | +| 65,537 | 8 | 16.529 | 6.060 | 2.48 [2.28–2.68] | +| 262,145 | 0 | 27.686 | 10.520 | 2.72 [2.62–2.82] | +| 262,145 | 3 | 50.530 | 14.552 | 3.32 [3.02–3.64] | +| 262,145 | 8 | 91.752 | 24.697 | 3.37 [3.01–3.74] | + +### construction survey 2 + +| N | Kind/profile | Serial ms | Parallel ms | Serial/parallel [95%] | +| ---: | ---: | ---: | ---: | --- | +| 4,097 | 0 | 0.246 | 0.245 | 1.00 [0.99–1.01] | +| 4,097 | 3 | 0.513 | 0.499 | 1.01 [0.99–1.04] | +| 4,097 | 8 | 0.879 | 0.881 | 0.99 [0.97–1.01] | +| 32,768 | 0 | 2.697 | 1.891 | 1.38 [1.26–1.52] | +| 32,768 | 3 | 5.111 | 3.124 | 1.63 [1.53–1.73] | +| 32,768 | 8 | 8.616 | 5.160 | 1.63 [1.55–1.71] | +| 65,537 | 0 | 6.348 | 3.514 | 1.88 [1.80–1.96] | +| 65,537 | 3 | 9.575 | 4.169 | 2.35 [2.21–2.50] | +| 65,537 | 8 | 18.554 | 7.287 | 2.44 [2.30–2.60] | +| 262,145 | 0 | 26.271 | 9.427 | 2.75 [2.64–2.87] | +| 262,145 | 3 | 43.020 | 11.913 | 3.55 [3.36–3.76] | +| 262,145 | 8 | 80.755 | 21.493 | 3.75 [3.56–3.95] | + +### gicp survey 1 + +| N | Kind/profile | Serial ms | Parallel ms | Serial/parallel [95%] | +| ---: | ---: | ---: | ---: | --- | +| 4,097 | 0 | 6.424 | 6.494 | 1.00 [0.97–1.04] | +| 4,097 | 1 | 6.465 | 6.343 | 1.01 [0.97–1.05] | +| 4,097 | 2 | 7.156 | 7.096 | 1.04 [0.98–1.10] | +| 32,768 | 0 | 21.479 | 19.466 | 1.11 [1.04–1.19] | +| 32,768 | 1 | 21.806 | 19.504 | 1.10 [1.06–1.15] | +| 32,768 | 2 | 23.211 | 19.070 | 1.17 [1.10–1.26] | +| 65,537 | 0 | 33.916 | 28.439 | 1.21 [1.09–1.32] | +| 65,537 | 1 | 44.369 | 30.647 | 1.27 [1.18–1.37] | +| 65,537 | 2 | 41.717 | 33.129 | 1.25 [1.17–1.34] | +| 262,145 | 0 | 139.662 | 95.646 | 1.40 [1.28–1.55] | +| 262,145 | 1 | 117.238 | 84.482 | 1.36 [1.26–1.46] | +| 262,145 | 2 | 138.657 | 90.196 | 1.48 [1.39–1.56] | + +### gicp survey 2 + +| N | Kind/profile | Serial ms | Parallel ms | Serial/parallel [95%] | +| ---: | ---: | ---: | ---: | --- | +| 4,097 | 0 | 7.088 | 7.099 | 0.99 [0.95–1.03] | +| 4,097 | 1 | 6.475 | 6.384 | 1.01 [0.98–1.06] | +| 4,097 | 2 | 6.545 | 6.570 | 0.99 [0.95–1.04] | +| 32,768 | 0 | 16.184 | 13.855 | 1.18 [1.14–1.22] | +| 32,768 | 1 | 17.245 | 14.891 | 1.10 [1.00–1.19] | +| 32,768 | 2 | 19.592 | 17.434 | 1.11 [1.04–1.20] | +| 65,537 | 0 | 31.721 | 25.296 | 1.16 [1.07–1.25] | +| 65,537 | 1 | 29.991 | 26.272 | 1.17 [1.09–1.25] | +| 65,537 | 2 | 31.590 | 24.876 | 1.24 [1.13–1.37] | +| 262,145 | 0 | 110.475 | 75.764 | 1.39 [1.26–1.53] | +| 262,145 | 1 | 114.617 | 76.459 | 1.48 [1.37–1.60] | +| 262,145 | 2 | 121.355 | 81.975 | 1.45 [1.36–1.55] | + +Times are medians in milliseconds; ratios are geometric means of paired +serial/parallel latency. >1 favors parallel. Within-survey 95% intervals use +10,000 paired log-ratio bootstrap draws, seed +20261013+24*(family_is_gicp)+12*(survey-1)+condition_index, sorted endpoints +249/9749. All samples, including adverse results, are retained. Intervals do +not establish portable gains, cold GPU/runtime latency or a global crossover. + +At 65,537 points, complete construction improves **1.88–2.55x**, and complete +GPU GICP calls **1.16–1.27x**. At 262,145 points, construction improves +**2.72–3.75x**, complete GICP **1.36–1.48x**. The corresponding within-survey +intervals exceed one in both repeats. Threshold-size construction also improves; +registration improvements there are smaller, and one interval lies near one. +4,097-point calls use serial construction in both cases; their near-one ratios +are measurement variability and are not parallel speedup claims. + +All complete registration diagnostics remain exact, not merely poses within a +tolerance. H2D receipts retain the resident-grid formula +20*S+36*T+32*target_nodes+4*(source_cells+target_cells+2); D2H remains 1,024 +bytes for three updates. The benefit is CPU preparation time. It does not +remove CPU grid construction, target-tree preparation or the 6x6 solver. +No extra GPU phase is claimed. + +## Final verification and delivery + +- Search all-feature suite: **58 passed, 6 ignored**; registration all-feature + library: **35 passed, 3 ignored**. The explicit native 216-condition exact + GICP comparison and all four exact release benchmark runs pass separately. + New minimal facade smoke passes. +- All-feature/all-target search/registration Clippy, minimal search parallel + without KD-tree, CPU-only GICP Clippy, WASM search/optional GPU registration + compile and workspace formatting pass. WASM is compile checked, not browser + executed; its constructor remains serial. +- Serial construction, axis selection and partition source are unchanged + against parent 713ba20. CPU GICP and the common resident update loop remain + unchanged. No GPU shader, core model, dependency or unsafe changes are added. +- Raw timing/interval results are retained in + `C:\Users\rsasa\Workspace\SpatialRust\target\automation\parallel-kdtree-build-measurements.json`. + Complete construction logs: + `C:\Users\rsasa\Workspace\SpatialRust\target\automation\parallel-kdtree-build-construction-1.log` + and + `C:\Users\rsasa\Workspace\SpatialRust\target\automation\parallel-kdtree-build-construction-2.log`. + Complete GICP logs: + `C:\Users\rsasa\Workspace\SpatialRust\target\automation\parallel-kdtree-build-gicp-1.log` + and + `C:\Users\rsasa\Workspace\SpatialRust\target\automation\parallel-kdtree-build-gicp-2.log`. + Publication files/PR text omit the local adapter model; scans match zero. + Existing 124 named CI checks remain required before merge.