From 6a48d8bad83351d0fb4ac3995223ee65232024cd Mon Sep 17 00:00:00 2001 From: rsasaki0109 Date: Tue, 13 Oct 2026 09:20:53 +0900 Subject: [PATCH 1/2] Add GPU resident radius GICP covariance preparation --- CHANGELOG.md | 9 + README.md | 9 + .../src/kernels/gicp_covariance_tests.rs | 162 ++++++----- .../src/kernels/gicp_covariances.rs | 124 +++++++-- .../kernels/gicp_radius_covariance_tests.rs | 258 ++++++++++++++++++ .../src/kernels/gicp_reduction.rs | 3 +- .../src/shaders/gicp_covariances.wgsl | 16 +- .../src/gicp_reduction.rs | 77 +++++- .../gpu_gicp_resident_radius_covariances.rs | 205 ++++++++++++++ ...p_resident_radius_covariances_benchmark.rs | 171 ++++++++++++ ..._gicp_resident_radius_covariances_smoke.rs | 25 ++ docs/ARCHITECTURE.md | 13 + docs/FEATURE_MATRIX.md | 16 +- docs/ROADMAP.md | 3 +- ...13_gpu_gicp_resident_radius_covariances.md | 246 +++++++++++++++++ 15 files changed, 1238 insertions(+), 99 deletions(-) create mode 100644 crates/spatialrust-gpu/src/kernels/gicp_radius_covariance_tests.rs create mode 100644 crates/spatialrust-registration/tests/gpu_gicp_resident_radius_covariances.rs create mode 100644 crates/spatialrust-registration/tests/gpu_gicp_resident_radius_covariances_benchmark.rs create mode 100644 crates/spatialrust/tests/gpu_gicp_resident_radius_covariances_smoke.rs create mode 100644 notes/2026-10-13_gpu_gicp_resident_radius_covariances.md diff --git a/CHANGELOG.md b/CHANGELOG.md index d9c680d..fa07e78 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,15 @@ removed no sooner than the next major (see `docs/API_STABILITY.md`). ### Changed +- Add explicit resident radius covariance preparation for GICP under + `register-gicp-gpu`. Inclusive bounded GPU neighborhoods, centered/scaled + compensated central sums and plane eigenvectors feed resident equations + without covariance readback/re-upload. Caller capacity includes self and + duplicates; overflow rejects all output without truncation. Sparse rows + preserve epsilon*I. Radius must be positive finite with normal f32 square, + epsilon=1e-4..=1. CPU trees/6x6 solves remain; legacy routes retain behavior. + Device f32 accumulation and Jacobi iteration can differ from grid results. + - Add explicit resident k-NN covariance preparation for GICP under `register-gicp-gpu`. GPU bounded neighborhoods, centered/scaled sample covariance and plane eigenvectors feed the existing resident equation loop diff --git a/README.md b/README.md index 7422489..fcb7bd8 100644 --- a/README.md +++ b/README.md @@ -672,6 +672,15 @@ readback or re-upload occurs. This path requires covariance_radius=None, k=3..=64 and epsilon=1e-4..=1; CPU builds trees and solves 6x6 equations. See [resident covariance measurements](notes/2026-10-13_gpu_gicp_resident_covariances.md). +`GicpRegistration::align_with_gpu_resident_radius_covariances(&runtime, &source, +&target, neighbor_capacity)` prepares radius neighborhoods and covariances +without CPU readback/re-upload. Set `covariance_radius=Some(radius)` and choose +a positive capacity including self/duplicates; overflow rejects the result. +Radius must have a finite normal f32 square, epsilon must be in 1e-4..=1. +Sparse neighborhoods retain the grid route's epsilon*I fallback. CPU tree +construction and 6x6 solves remain; only 256-byte equations are downloaded. +See [resident radius evidence](notes/2026-10-13_gpu_gicp_resident_radius_covariances.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-gpu/src/kernels/gicp_covariance_tests.rs b/crates/spatialrust-gpu/src/kernels/gicp_covariance_tests.rs index 675819c..18f1f81 100644 --- a/crates/spatialrust-gpu/src/kernels/gicp_covariance_tests.rs +++ b/crates/spatialrust-gpu/src/kernels/gicp_covariance_tests.rs @@ -89,7 +89,7 @@ fn matrix(c: [f32; 6]) -> Mat3 { #[test] fn resident_covariance_plans_bound_counts_arity_epsilon_and_dispatch() { - assert_eq!(std::mem::size_of::(), 16); + assert_eq!(std::mem::size_of::(), 32); let limits = wgpu::Limits::default(); assert_eq!( covariance_plan(65537, 65554, 20, 0.001, &limits).unwrap(), @@ -296,76 +296,106 @@ fn malformed_resident_neighbor_rows_invalidate_both_covariance_sets() { let source = [[0.0; 3], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]]; let st = tree(&source); let tt = tree(&source[..2]); - for second in [false, true] { - for case in 0..3 { - let positions = GpuIcpInput::upload(&runtime, &source, &tt).unwrap(); - let source_tree = GpuKdTree::upload(&runtime, &st).unwrap(); - let sn = source_tree.nearest_k_resident(&runtime, &positions.source, 3).unwrap(); - let tn = positions.tree.nearest_k_resident(&runtime, &positions.target, 3).unwrap(); - let statement = match case { - 0 => "data[0]=0u;", - 1 => "data[2]=0xffffffffu;", - _ => "data[3]=0x7f800000u;", - }; - let shader=runtime.device().create_shader_module(wgpu::ShaderModuleDescriptor{label:None,source:wgpu::ShaderSource::Wgsl(format!("@group(0) @binding(0) var data:array; @compute @workgroup_size(1) fn main() {{{statement}}}").into())}); - let pipeline = - runtime.device().create_compute_pipeline(&wgpu::ComputePipelineDescriptor { + for radius_mode in [0, 1] { + for second in [false, true] { + for case in 0..if radius_mode == 0 { 3 } else { 6 } { + let positions = GpuIcpInput::upload(&runtime, &source, &tt).unwrap(); + let source_tree = GpuKdTree::upload(&runtime, &st).unwrap(); + let sn = if radius_mode == 0 { + source_tree.nearest_k_resident(&runtime, &positions.source, 3).unwrap() + } else { + source_tree + .radius_neighbors_resident(&runtime, &positions.source, 2.0, 3) + .unwrap() + }; + let tn = if radius_mode == 0 { + positions.tree.nearest_k_resident(&runtime, &positions.target, 3).unwrap() + } else { + positions + .tree + .radius_neighbors_resident(&runtime, &positions.target, 2.0, 3) + .unwrap() + }; + let statement = match case { + 0 => "data[0]=0u;", + 1 => "data[2]=0xffffffffu;", + 2 => "data[3]=0x7f800000u;", + 3 => "data[3]=0x40a00000u;", + 4 => "data[1]=0u;", + _ => "data[1]=3u;data[0]=3u;data[2]=0xffffffffu;", + }; + let shader=runtime.device().create_shader_module(wgpu::ShaderModuleDescriptor{label:None,source:wgpu::ShaderSource::Wgsl(format!("@group(0) @binding(0) var data:array; @compute @workgroup_size(1) fn main() {{{statement}}}").into())}); + let pipeline = + runtime.device().create_compute_pipeline(&wgpu::ComputePipelineDescriptor { + label: None, + layout: None, + module: &shader, + entry_point: Some("main"), + compilation_options: Default::default(), + cache: None, + }); + let buffer = if second { tn.buffer() } else { sn.buffer() }; + let group = runtime.device().create_bind_group(&wgpu::BindGroupDescriptor { label: None, - layout: None, - module: &shader, - entry_point: Some("main"), - compilation_options: Default::default(), - cache: None, + layout: &pipeline.get_bind_group_layout(0), + entries: &[wgpu::BindGroupEntry { + binding: 0, + resource: buffer.as_entire_binding(), + }], }); - let buffer = if second { tn.buffer() } else { sn.buffer() }; - let group = runtime.device().create_bind_group(&wgpu::BindGroupDescriptor { - label: None, - layout: &pipeline.get_bind_group_layout(0), - entries: &[wgpu::BindGroupEntry { - binding: 0, - resource: buffer.as_entire_binding(), - }], - }); - let mut encoder = runtime.device().create_command_encoder(&Default::default()); - { - let mut pass = encoder.begin_compute_pass(&wgpu::ComputePassDescriptor { + let mut encoder = runtime.device().create_command_encoder(&Default::default()); + { + let mut pass = encoder.begin_compute_pass(&wgpu::ComputePassDescriptor { + label: None, + timestamp_writes: None, + }); + pass.set_pipeline(&pipeline); + pass.set_bind_group(0, &group, &[]); + pass.dispatch_workgroups(1, 1, 1); + } + runtime.queue().submit(Some(encoder.finish())); + let covariances = runtime.device().create_buffer(&wgpu::BufferDescriptor { label: None, - timestamp_writes: None, + size: 6 * 32, + usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_SRC, + mapped_at_creation: false, }); - pass.set_pipeline(&pipeline); - pass.set_bind_group(0, &group, &[]); - pass.dispatch_workgroups(1, 1, 1); + encode_covariances( + &runtime, + &positions, + [&sn, &tn], + &covariances, + Params { + sources: 3, + targets: 2, + neighbors: 3, + epsilon: 0.001, + radius_squared: 4.0, + radius_mode, + padding: [0; 2], + }, + ); + let cov = words(&runtime, &covariances); + assert_eq!(cov[0], 1); + assert!(cov[1..].iter().all(|&v| v == 0)); + let transfers = positions.upload_transfers(); + let input = GpuGicpInput { + positions, + covariances, + weight_scale: 6000.0, + covariance_offset: 1, + transfers, + }; + let output = + input.equations_resident(&runtime, Isometry3::identity(), 2.0).unwrap(); + let w = words(&runtime, output.buffer()); + assert_eq!(w[0], 16); + assert!(w[1..].iter().all(|&v| v == 0)); + assert!(output.download(&runtime).is_err()); } - runtime.queue().submit(Some(encoder.finish())); - let covariances = runtime.device().create_buffer(&wgpu::BufferDescriptor { - label: None, - size: 6 * 32, - usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_SRC, - mapped_at_creation: false, - }); - encode_covariances( - &runtime, - &positions, - [&sn, &tn], - &covariances, - Params { sources: 3, targets: 2, neighbors: 3, epsilon: 0.001 }, - ); - let cov = words(&runtime, &covariances); - assert_eq!(cov[0], 1); - assert!(cov[1..].iter().all(|&v| v == 0)); - let transfers = positions.upload_transfers(); - let input = GpuGicpInput { - positions, - covariances, - weight_scale: 6000.0, - covariance_offset: 1, - transfers, - }; - let output = input.equations_resident(&runtime, Isometry3::identity(), 2.0).unwrap(); - let w = words(&runtime, output.buffer()); - assert_eq!(w[0], 16); - assert!(w[1..].iter().all(|&v| v == 0)); - assert!(output.download(&runtime).is_err()); } } } + +#[path = "gicp_radius_covariance_tests.rs"] +mod radius; diff --git a/crates/spatialrust-gpu/src/kernels/gicp_covariances.rs b/crates/spatialrust-gpu/src/kernels/gicp_covariances.rs index f84f843..4ac2664 100644 --- a/crates/spatialrust-gpu/src/kernels/gicp_covariances.rs +++ b/crates/spatialrust-gpu/src/kernels/gicp_covariances.rs @@ -12,6 +12,15 @@ struct Params { targets: u32, neighbors: u32, epsilon: f32, + radius_squared: f32, + radius_mode: u32, + padding: [u32; 2], +} + +#[derive(Clone, Copy)] +enum Preparation { + Knn(u32), + Radius { radius: f32, capacity: u32 }, } fn covariance_plan( @@ -21,17 +30,26 @@ fn covariance_plan( epsilon: f32, limits: &wgpu::Limits, ) -> SpatialResult { - if source == 0 - || target == 0 - || !(3..=MAX_GPU_K_NEIGHBORS).contains(&k) - || !(1e-4..=1.0).contains(&epsilon) - { + if !(3..=MAX_GPU_K_NEIGHBORS).contains(&k) { + return Err(invalid()); + } + bounded_covariance_plan(source, target, k, epsilon, limits) +} + +fn bounded_covariance_plan( + source: u64, + target: u64, + capacity: u32, + epsilon: f32, + limits: &wgpu::Limits, +) -> SpatialResult { + if source == 0 || target == 0 || capacity == 0 || !(1e-4..=1.0).contains(&epsilon) { return Err(invalid()); } let total = source.checked_add(target).and_then(|v| v.checked_add(1)).ok_or_else(invalid)?; plan(source, total, limits)?; let bound = limits.max_buffer_size.min(u64::from(limits.max_storage_buffer_binding_size)); - let stride = 8 + 8 * u64::from(k); + let stride = 8 + 8 * u64::from(capacity); for count in [source, target] { if count > u64::from(u32::MAX) || count.div_ceil(64) > u64::from(limits.max_compute_workgroups_per_dimension) @@ -50,6 +68,10 @@ fn covariance_plan( Ok(total * 32) } +fn valid_radius(radius: f32) -> bool { + radius.is_finite() && radius > 0.0 && (radius * radius).is_normal() +} + impl GpuGicpInput { /// Explicitly uploads finite source/target trees and original-order positions, /// then computes both k-NN plane-regularized covariances entirely on GPU. @@ -70,16 +92,58 @@ impl GpuGicpInput { target: &KdTree, k: u32, epsilon: f32, + ) -> SpatialResult { + Self::upload_prepared(runtime, source, target, Preparation::Knn(k), epsilon) + } + + /// Uploads original positions/trees and prepares complete inclusive radius + /// neighborhoods and plane covariances on GPU, without payload readback. + /// `capacity` is a positive per-point bound, including self and duplicates. + /// Overflow invalidates both covariance sets and every equation download; + /// rows are never silently truncated. Radius must be finite/positive with + /// a normal f32 square; epsilon must be in 1e-4..=1. + /// Fewer than three neighbors use epsilon*I, as in the existing grid route. + /// Other rows use centered/scaled compensated sums and 32 max-pivot Jacobi + /// steps with physical 1e-20 stopping on the unnormalized central sum. + /// f32 arithmetic, traversal order and repeated eigenvalues can differ + /// from the legacy grid's raw means and 16 three-pivot sweeps. + /// Transfers, ownership and CPU solve boundaries match the k-NN method. + pub fn upload_with_gpu_radius_covariances( + runtime: &WgpuRuntime, + source: &KdTree, + target: &KdTree, + radius: f32, + capacity: u32, + epsilon: f32, + ) -> SpatialResult { + if !valid_radius(radius) { + return Err(invalid()); + } + Self::upload_prepared( + runtime, + source, + target, + Preparation::Radius { radius, capacity }, + epsilon, + ) + } + + fn upload_prepared( + runtime: &WgpuRuntime, + source: &KdTree, + target: &KdTree, + preparation: Preparation, + epsilon: f32, ) -> SpatialResult { let (x, y, z) = source.layout().positions(); let targets = target.layout().positions().0.len(); - let bytes = covariance_plan( - x.len() as u64, - targets as u64, - k, - epsilon, - &runtime.device().limits(), - )?; + let (capacity, radius) = match preparation { + Preparation::Knn(k) => (k, None), + Preparation::Radius { radius, capacity } => (capacity, Some(radius)), + }; + let planner = if radius.is_some() { bounded_covariance_plan } else { covariance_plan }; + let bytes = + planner(x.len() as u64, targets as u64, capacity, epsilon, &runtime.device().limits())?; let bound = runtime .device() .limits() @@ -93,16 +157,42 @@ impl GpuGicpInput { let points: Vec<_> = (0..x.len()).map(|i| [x[i], y[i], z[i]]).collect(); let positions = GpuIcpInput::upload(runtime, &points, target)?; let source_tree = GpuKdTree::upload(runtime, source)?; - let source_neighbors = source_tree.nearest_k_resident(runtime, &positions.source, k)?; - let target_neighbors = positions.tree.nearest_k_resident(runtime, &positions.target, k)?; + let (source_neighbors, target_neighbors) = if let Some(radius) = radius { + ( + source_tree.radius_neighbors_resident( + runtime, + &positions.source, + radius, + capacity, + )?, + positions.tree.radius_neighbors_resident( + runtime, + &positions.target, + radius, + capacity, + )?, + ) + } else { + ( + source_tree.nearest_k_resident(runtime, &positions.source, capacity)?, + positions.tree.nearest_k_resident(runtime, &positions.target, capacity)?, + ) + }; let covariances = runtime.device().create_buffer(&wgpu::BufferDescriptor { label: Some("gicp-resident-covariances"), size: bytes, usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_SRC, mapped_at_creation: false, }); - let params = - Params { sources: x.len() as u32, targets: targets as u32, neighbors: k, epsilon }; + let params = Params { + sources: x.len() as u32, + targets: targets as u32, + neighbors: capacity, + epsilon, + radius_squared: radius.map_or(0.0, |r| r * r), + radius_mode: u32::from(radius.is_some()), + padding: [0; 2], + }; encode_covariances( runtime, &positions, diff --git a/crates/spatialrust-gpu/src/kernels/gicp_radius_covariance_tests.rs b/crates/spatialrust-gpu/src/kernels/gicp_radius_covariance_tests.rs new file mode 100644 index 0000000..f95e738 --- /dev/null +++ b/crates/spatialrust-gpu/src/kernels/gicp_radius_covariance_tests.rs @@ -0,0 +1,258 @@ +use super::*; + +fn radius_reference(points: &[[f32; 3]], i: usize, radius: f32) -> ([f64; 6], Mat3, usize) { + // Independent brute-force inclusive membership, including self/duplicates. + let neighbors: Vec<_> = points + .iter() + .filter(|p| { + let d: [f32; 3] = std::array::from_fn(|a| p[a] - points[i][a]); + d.iter().map(|v| v * v).sum::() <= radius * radius + }) + .collect(); + let zero = Mat3::from_rows([0.0; 3], [0.0; 3], [0.0; 3]); + if neighbors.len() < 3 { + return ( + [f64::from(0.001_f32), f64::from(0.001_f32), f64::from(0.001_f32), 0.0, 0.0, 0.0], + zero, + neighbors.len(), + ); + } + let anchor = points[i].map(f64::from); + let scale = neighbors + .iter() + .flat_map(|p| (0..3).map(|a| (f64::from(p[a]) - anchor[a]).abs())) + .fold(0.0_f64, f64::max); + if scale == 0.0 { + return ([0.001, 1.0, 1.0, 0.0, 0.0, 0.0], zero, neighbors.len()); + } + let mean: [f64; 3] = std::array::from_fn(|a| { + neighbors.iter().map(|p| (f64::from(p[a]) - anchor[a]) / scale).sum::() + / neighbors.len() as f64 + }); + let mut c = [[0.0; 3]; 3]; + for p in &neighbors { + let d: [f64; 3] = std::array::from_fn(|a| (f64::from(p[a]) - anchor[a]) / scale - mean[a]); + for a in 0..3 { + for b in 0..3 { + c[a][b] += d[a] * d[b]; + } + } + } + let raw = Mat3::from_rows(c[0], c[1], c[2]); + let e = symmetric_eigen3(raw); + let n: [f64; 3] = std::array::from_fn(|a| e.eigenvectors.m[a][0]); + let epsilon = f64::from(0.001_f32); + ( + [ + 1.0 + (epsilon - 1.0) * n[0] * n[0], + 1.0 + (epsilon - 1.0) * n[1] * n[1], + 1.0 + (epsilon - 1.0) * n[2] * n[2], + (epsilon - 1.0) * n[0] * n[1], + (epsilon - 1.0) * n[0] * n[2], + (epsilon - 1.0) * n[1] * n[2], + ], + raw, + neighbors.len(), + ) +} + +#[test] +fn resident_radius_covariances_match_brute_force_planes_and_spectral_invariants() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let mut checked = 0; + let mut maximum = 0.0_f64; + for n in [1, 2, 3, 63, 64, 65, 257, 1025] { + for scale in [1e-4_f32, 1.0, 1e4] { + let points: Vec<_> = (0..n) + .map(|i| { + let a = hash(i as u32 + 1); + let b = hash(i as u32 + 0x10203); + [(a + 4.0) * scale, (b - 3.0) * scale, (0.3 * a - 0.2 * b + 2.0) * scale] + }) + .collect(); + let tree = tree(&points); + for r in [0.01, 0.8, 20.0] { + let radius = r * scale; + let input = GpuGicpInput::upload_with_gpu_radius_covariances( + &runtime, &tree, &tree, radius, n as u32, 0.001, + ) + .unwrap(); + let cov = covariances(&runtime, &input); + assert_eq!(input.upload_transfers().device_to_host_bytes(), 0); + assert_eq!( + input.upload_transfers().host_to_device_bytes(), + (64 * (n + tree.layout().nodes().len())) as u64 + ); + for i in 0..n { + assert_eq!(cov[i], cov[n + i]); + let (expected, raw, count) = radius_reference(&points, i, radius); + let e = symmetric_eigen3(matrix(cov[i])); + let spectrum = symmetric_eigen3(raw).eigenvalues; + let weights = if count < 3 { [0.001; 3] } else { [0.001, 1.0, 1.0] }; + for (a, b) in e.eigenvalues.iter().zip(weights) { + assert!((a - b).abs() < 2e-6); + } + if count >= 3 { + let normal: [f64; 3] = std::array::from_fn(|a| e.eigenvectors.m[a][0]); + let energy = (0..3) + .flat_map(|a| (0..3).map(move |b| normal[a] * raw.m[a][b] * normal[b])) + .sum::(); + assert!( + energy <= spectrum[0] + 1e-10 + 3e-5 * spectrum[2].abs(), + "n={n} r={r} scale={scale} i={i} energy={energy}" + ); + } + if count < 3 + || spectrum[1] - spectrum[0] > (spectrum[2].abs() * 1e-4).max(1e-11) + { + for a in 0..6 { + let error = (f64::from(cov[i][a]) - expected[a]).abs(); + maximum = maximum.max(error); + assert!(error < 3e-4, "n={n} r={r} scale={scale} i={i} error={error}"); + } + } + checked += 1; + } + } + } + } + eprintln!("GPU_RESIDENT_RADIUS_COVARIANCE checked={checked}; maximum centered f64 component error={maximum}"); +} + +#[test] +fn radius_boundaries_sparse_fallback_duplicates_overflow_and_numeric_errors() { + let runtime = WgpuRuntime::new_headless().unwrap(); + // Exactly-on-radius samples must participate; the next representable value does not. + let points = [ + [0.0; 3], + [0.5, 0.0, 0.0], + [0.0, 0.5, 0.0], + [f32::from_bits(0.5_f32.to_bits() + 1), 0.0, 0.0], + ]; + let t = tree(&points); + let input = + GpuGicpInput::upload_with_gpu_radius_covariances(&runtime, &t, &t, 0.5, 4, 0.001).unwrap(); + assert_eq!(covariances(&runtime, &input)[0], [1.0, 1.0, 0.001, 0.0, 0.0, 0.0]); + for epsilon in [1e-4, 0.001, 1.0] { + let t = tree(&[[0.0; 3], [0.25, 0.0, 0.0]]); + let input = + GpuGicpInput::upload_with_gpu_radius_covariances(&runtime, &t, &t, 0.5, 2, epsilon) + .unwrap(); + assert!(covariances(&runtime, &input) + .iter() + .all(|c| *c == [epsilon, epsilon, epsilon, 0.0, 0.0, 0.0])); + } + let duplicate = tree(&[[0.0; 3]; 65]); + let singleton = tree(&[[0.0; 3]]); + let good = GpuGicpInput::upload_with_gpu_radius_covariances( + &runtime, &duplicate, &singleton, 0.5, 65, 0.001, + ) + .unwrap(); + let c = covariances(&runtime, &good); + assert_eq!(c[0], [0.001, 1.0, 1.0, 0.0, 0.0, 0.0]); + assert_eq!(c[65], [0.001, 0.001, 0.001, 0.0, 0.0, 0.0]); + for second in [false, true] { + let (s, t) = if second { (&singleton, &duplicate) } else { (&duplicate, &singleton) }; + let bad = GpuGicpInput::upload_with_gpu_radius_covariances(&runtime, s, t, 0.5, 64, 0.001) + .unwrap(); + let w = words(&runtime, &bad.covariances); + assert_eq!(w[0], 8); + assert!(w[1..].iter().all(|v| *v == 0)); + let eq = bad.equations_resident(&runtime, Isometry3::identity(), 0.25).unwrap(); + let w = words(&runtime, eq.buffer()); + assert_eq!(w[0], 16); + assert!(w[1..].iter().all(|v| *v == 0)); + assert!(eq.download(&runtime).is_err()); + } + let tiny = tree(&[[0.0; 3], [1e-22, 0.0, 0.0], [0.0, 1e-22, 0.0]]); + let bad = GpuGicpInput::upload_with_gpu_radius_covariances(&runtime, &tiny, &t, 0.5, 4, 0.001) + .unwrap(); + let w = words(&runtime, &bad.covariances); + assert_eq!(w[0], 2); + assert!(w[1..].iter().all(|v| *v == 0)); + for r in [0.0, -1.0, f32::NAN, f32::INFINITY, 1e20, 1e-22] { + assert!(GpuGicpInput::upload_with_gpu_radius_covariances(&runtime, &t, &t, r, 4, 0.001) + .is_err()); + } + for capacity in [0, u32::MAX] { + assert!(GpuGicpInput::upload_with_gpu_radius_covariances( + &runtime, &t, &t, 0.5, capacity, 0.001 + ) + .is_err()); + } + assert!(GpuGicpInput::upload_with_gpu_radius_covariances( + &runtime, + &tree(&[]), + &t, + 0.5, + 4, + 0.001 + ) + .is_err()); + assert!(GpuGicpInput::upload_with_gpu_radius_covariances( + &runtime, + &tree(&[[f32::NAN, 0.0, 0.0]]), + &t, + 0.5, + 4, + 0.001 + ) + .is_err()); + assert_eq!(bounded_covariance_plan(1, 2, 128, 0.001, &wgpu::Limits::default()).unwrap(), 128); + assert!(bounded_covariance_plan(65537, 65537, 1024, 0.001, &wgpu::Limits::default()).is_err()); +} + +#[test] +fn radius_covariances_match_unique_legacy_grid_planes_and_feed_owned_equations() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let other = WgpuRuntime::new_headless().unwrap(); + let points: Vec<_> = (0..257) + .map(|i| { + let a = hash(i + 1); + let b = hash(i + 0x10203); + [a, b, 0.3 * a - 0.2 * b + 2.0] + }) + .collect(); + let t = tree(&points); + let input = + GpuGicpInput::upload_with_gpu_radius_covariances(&runtime, &t, &t, 20.0, 257, 0.001) + .unwrap(); + let c = covariances(&runtime, &input); + let (x, y, z) = t.layout().positions(); + let old = + crate::estimate_plane_covariances_grid_gpu_with_receipt(&runtime, x, y, z, 20.0, 0.001) + .unwrap() + .0; + for (a, b) in c[..257].iter().zip(&old) { + for axis in 0..6 { + assert!((a[axis] - b[axis]).abs() < 3e-5); + } + } + let uploaded = GpuGicpInput::upload(&runtime, &points, &c[..257], &t, &c[257..]).unwrap(); + let pose = Isometry3::new(Quat::::identity(), Vec3::new(0.01, -0.005, 0.002)); + let result = input.equations_resident(&runtime, pose, 0.25).unwrap(); + let actual = result.download(&runtime).unwrap(); + let expected = + uploaded.equations_resident(&runtime, pose, 0.25).unwrap().download(&runtime).unwrap(); + assert_eq!(actual.squared_error_sum, expected.squared_error_sum); + for a in 0..6 { + assert!( + (actual.gradient[a] * actual.weight_scale + - expected.gradient[a] * expected.weight_scale) + .abs() + < 2e-4 * (1.0 + (expected.gradient[a] * expected.weight_scale).abs()) + ); + for b in 0..6 { + assert!( + (actual.hessian[a][b] * actual.weight_scale + - expected.hessian[a][b] * expected.weight_scale) + .abs() + < 2e-4 * (1.0 + (expected.hessian[a][b] * expected.weight_scale).abs()) + ); + } + } + assert!(input.equations_resident(&other, pose, 0.25).is_err()); + drop(input); + assert_eq!(result.download(&runtime).unwrap(), actual); + assert!(result.download(&other).is_err()); +} diff --git a/crates/spatialrust-gpu/src/kernels/gicp_reduction.rs b/crates/spatialrust-gpu/src/kernels/gicp_reduction.rs index 5e3d6cb..1293b7f 100644 --- a/crates/spatialrust-gpu/src/kernels/gicp_reduction.rs +++ b/crates/spatialrust-gpu/src/kernels/gicp_reduction.rs @@ -365,7 +365,8 @@ impl GpuGicpNormalEquations { /// are Hessian upper triangle by row, 21..26 gradient and 27 Euclidean error /// normalized by distance squared. Status bits 1/2/4/8 indicate transform, /// match, nonfinite arithmetic and unsupported covariance conditioning; - /// bit 16 indicates invalid resident covariance preparation. + /// bit 16 indicates invalid resident covariance preparation, including + /// radius neighborhood capacity overflow. No truncated rows are consumed. /// Any nonzero status invalidates ALL payload, which is zeroed on device. /// GPU consumers must check status and preserve immutable application data. #[must_use] diff --git a/crates/spatialrust-gpu/src/shaders/gicp_covariances.wgsl b/crates/spatialrust-gpu/src/shaders/gicp_covariances.wgsl index d0717f9..7f7055e 100644 --- a/crates/spatialrust-gpu/src/shaders/gicp_covariances.wgsl +++ b/crates/spatialrust-gpu/src/shaders/gicp_covariances.wgsl @@ -1,4 +1,4 @@ -struct Params { sources:u32, targets:u32, neighbors:u32, epsilon:f32 } +struct Params { sources:u32, targets:u32, neighbors:u32, epsilon:f32, radius_squared:f32, radius_mode:u32, padding:vec2 } struct Covariance { diagonal:vec4, off:vec4 } struct Covariances { status:atomic, reserved:array, rows:array } @group(0) @binding(0) var params:Params; @@ -61,8 +61,13 @@ fn build_covariances(@builtin(global_invocation_id) id:vec3) { result.rows[index]=Covariance(vec4(0.0),vec4(0.0)); let second=index>=params.sources;let i=select(index,index-params.sources,second);let size=select(params.sources,params.targets,second); let base=i*(2u+2u*params.neighbors);let stored=word(base,second);let total=word(base+1u,second); - if stored!=min(size,params.neighbors)||total!=stored {atomicOr(&result.status,1u);return;} - for(var j=0u;j(word(base+3u+2u*j,second));if neighbor>=size||!(distance>=0.0)||!finite(distance) {atomicOr(&result.status,1u);return;}} + if params.radius_mode!=0u { + if total>params.neighbors {atomicOr(&result.status,8u);return;} + if total==0u||total>size||stored!=total {atomicOr(&result.status,1u);return;} + } else if stored!=min(size,params.neighbors)||total!=stored {atomicOr(&result.status,1u);return;} + for(var j=0u;j(word(base+3u+2u*j,second));if neighbor>=size||!(distance>=0.0)||!finite(distance)||(params.radius_mode!=0u&&distance>params.radius_squared) {atomicOr(&result.status,1u);return;}} + // Preserve the radius grid's sparse isotropic fallback, including self. + if params.radius_mode!=0u&&stored<3u {result.rows[index]=Covariance(vec4(vec3(params.epsilon),0.0),vec4(0.0));return;} // CPU's single-neighbor fallback uses identity covariance, regularized on X. if stored<2u {emit(index,vec3(1.0,0.0,0.0));return;} let origin=point(i,second);var scale=0.0; @@ -80,9 +85,10 @@ fn build_covariances(@builtin(global_invocation_id) id:vec3) { let products=array(d.x*d.x,d.y*d.y,d.z*d.z,d.x*d.y,d.x*d.z,d.y*d.z); for(var a=0u;a<6u;a+=1u) {sums[a]=add(sums[a],products[a]);} } - var cov:array;for(var a=0u;a<6u;a+=1u) {cov[a]=(sums[a].x+sums[a].y)/f32(stored-1u);if !finite(cov[a]) {atomicOr(&result.status,2u);return;}} + let divisor=select(f32(stored-1u),1.0,params.radius_mode!=0u); + var cov:array;for(var a=0u;a<6u;a+=1u) {cov[a]=(sums[a].x+sums[a].y)/divisor;if !finite(cov[a]) {atomicOr(&result.status,2u);return;}} let matrix=array,3>(vec3(cov[0],cov[3],cov[4]),vec3(cov[3],cov[1],cov[5]),vec3(cov[4],cov[5],cov[2])); - let n=normal(matrix,1e-12/scale_squared); + let n=normal(matrix,select(1e-12,1e-20,params.radius_mode!=0u)/scale_squared); if !finite3(n) {atomicOr(&result.status,2u);return;}emit(index,n); } diff --git a/crates/spatialrust-registration/src/gicp_reduction.rs b/crates/spatialrust-registration/src/gicp_reduction.rs index 3aae383..78001ea 100644 --- a/crates/spatialrust-registration/src/gicp_reduction.rs +++ b/crates/spatialrust-registration/src/gicp_reduction.rs @@ -62,6 +62,13 @@ fn fitness(equations: &GpuGicpEquations) -> f64 { } } +#[derive(Clone, Copy)] +enum CovariancePreparation { + Legacy, + ResidentKnn, + ResidentRadius(u32), +} + impl GicpRegistration { /// Retains source/target/covariances on GPU for transform, nearest matching, /// rotated covariance inverse, Hessian/gradient and Euclidean fitness passes. @@ -78,7 +85,7 @@ impl GicpRegistration { source: &PointCloud, target: &PointCloud, ) -> SpatialResult { - self.align_reduction_impl(runtime, source, target, false) + self.align_reduction_impl(runtime, source, target, CovariancePreparation::Legacy) } /// Computes k-NN neighborhoods and plane covariances on GPU, retaining the @@ -96,7 +103,33 @@ impl GicpRegistration { source: &PointCloud, target: &PointCloud, ) -> SpatialResult { - self.align_reduction_impl(runtime, source, target, true) + self.align_reduction_impl(runtime, source, target, CovariancePreparation::ResidentKnn) + } + + /// Computes complete inclusive radius neighborhoods and covariances on GPU + /// and retains them for the resident equation/update loop. Requires configured + /// covariance_radius=Some(finite positive radius) with a normal f32 square, + /// epsilon (as f32) in 1e-4..=1 and positive `neighbor_capacity` per point. + /// Capacity includes self/duplicates; overflow fails without truncation or + /// fallback. Fewer than three neighbors use epsilon*I as in the grid route. + /// Centered/scaled compensated f32 covariance and max-pivot Jacobi can differ + /// from grid accumulation/sweeps, especially at repeated eigenvalues/gates. + /// CPU builds both trees, solves 6x6 equations and composes/stops the pose. + /// Upload is 32*(S+T+source_nodes+target_nodes) bytes; downloads are only + /// 256 bytes per evaluation. No neighborhood/covariance payload crosses CPU. + pub fn align_with_gpu_resident_radius_covariances( + &self, + runtime: &WgpuRuntime, + source: &PointCloud, + target: &PointCloud, + neighbor_capacity: u32, + ) -> SpatialResult { + self.align_reduction_impl( + runtime, + source, + target, + CovariancePreparation::ResidentRadius(neighbor_capacity), + ) } fn align_reduction_impl( @@ -104,19 +137,28 @@ impl GicpRegistration { runtime: &WgpuRuntime, source: &PointCloud, target: &PointCloud, - resident_knn: bool, + preparation: CovariancePreparation, ) -> SpatialResult { self.validate_inputs(source, target)?; let c = self.config; let distance_squared = c.max_correspondence_distance * c.max_correspondence_distance; let epsilon = c.epsilon as f32; - if resident_knn + if matches!(preparation, CovariancePreparation::ResidentKnn) && (c.covariance_radius.is_some() || c.k_neighbors > spatialrust_gpu::MAX_GPU_K_NEIGHBORS as usize || !(1e-4..=1.0).contains(&epsilon)) { return Err(SpatialError::InvalidArgument("resident GPU GICP covariances require covariance_radius=None, k_neighbors=3..=64 and epsilon (as f32) in 1e-4..=1".to_owned())); } + if let CovariancePreparation::ResidentRadius(capacity) = preparation { + if capacity == 0 + || !(1e-4..=1.0).contains(&epsilon) + || c.covariance_radius + .map_or(true, |r| !r.is_finite() || r <= 0.0 || !(r * r).is_normal()) + { + return Err(SpatialError::InvalidArgument("resident GPU GICP radius covariances require a positive capacity including self, covariance_radius=Some(finite positive radius) with a normal f32 square and epsilon (as f32) in 1e-4..=1".to_owned())); + } + } if c.max_iterations == 0 || c.min_correspondences < 6 || !c.epsilon.is_finite() @@ -142,11 +184,32 @@ impl GicpRegistration { let tree = KdTree::from_slices(tx, ty, tz); let mut transfers = TransferStats::default(); let mut covariance_batches = 0; - let input = if resident_knn { - let k = u32::try_from(c.k_neighbors).map_err(|_| invalid())?; + let input = if !matches!(preparation, CovariancePreparation::Legacy) { let source_tree = KdTree::from_slices(sx, sy, sz); covariance_batches = 2; - GpuGicpInput::upload_with_gpu_covariances(runtime, &source_tree, &tree, k, epsilon)? + match preparation { + CovariancePreparation::ResidentKnn => { + let k = u32::try_from(c.k_neighbors).map_err(|_| invalid())?; + GpuGicpInput::upload_with_gpu_covariances( + runtime, + &source_tree, + &tree, + k, + epsilon, + )? + } + CovariancePreparation::ResidentRadius(capacity) => { + GpuGicpInput::upload_with_gpu_radius_covariances( + runtime, + &source_tree, + &tree, + c.covariance_radius.unwrap(), + capacity, + epsilon, + )? + } + CovariancePreparation::Legacy => unreachable!(), + } } else { let mut prepare = |x: &[f32], y: &[f32], diff --git a/crates/spatialrust-registration/tests/gpu_gicp_resident_radius_covariances.rs b/crates/spatialrust-registration/tests/gpu_gicp_resident_radius_covariances.rs new file mode 100644 index 0000000..6ba5a51 --- /dev/null +++ b/crates/spatialrust-registration/tests/gpu_gicp_resident_radius_covariances.rs @@ -0,0 +1,205 @@ +#![cfg(all(feature = "register-gicp-gpu", not(target_arch = "wasm32")))] +#[path = "gpu_gicp_reduction.rs"] +mod reference; +use spatialrust_core::{HasPositions3, PointCloudBuilder}; +use spatialrust_gpu::WgpuRuntime; +use spatialrust_math::{Isometry3, Quat, TransformPoint, Vec3}; +use spatialrust_registration::{transform_point_cloud, GicpConfig, GicpRegistration}; +use spatialrust_search::{KdTree, NearestNeighborIndex}; + +#[test] +#[ignore = "requires a native wgpu adapter; run explicitly with --include-ignored"] +fn resident_radius_gicp_matches_cpu_poses_stopping_fitness_and_transfer_receipts() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let mut maximum = 0.0_f32; + let mut compared = 0; + for n in [17, 63, 64, 65, 257, 1366] { + for scale in [1e-4, 1.0, 1e4] { + for radius in [0.075, 0.25, 0.5] { + for stop in 0..3 { + let target = reference::corner(n, scale); + let source = reference::shifted(&target, scale); + let saved_source = source.clone(); + let saved_target = target.clone(); + let mut config = + reference::registration(scale, Some(radius * scale), stop).config(); + config.k_neighbors = 6; + let registration = GicpRegistration::new(config); + let cpu = registration.align_with_diagnostics(&source, &target).unwrap(); + let previous = + registration.align_with_gpu_reduction(&runtime, &source, &target).unwrap(); + let actual = registration + .align_with_gpu_resident_radius_covariances(&runtime, &source, &target, 128) + .unwrap_or_else(|e| { + panic!("n={n} scale={scale} radius={radius} stop={stop}: {e}") + }); + maximum = maximum.max( + reference::compare_pose( + &source, + cpu.transform, + actual.result.transform, + scale, + ) / scale, + ); + reference::compare_pose( + &source, + previous.result.transform, + actual.result.transform, + scale, + ); + compared += source.len(); + assert_eq!(actual.result.iterations, cpu.iterations); + assert_eq!(actual.result.converged, cpu.converged); + assert_eq!(actual.covariance_batches, 2); + assert_eq!(actual.evaluation_batches, actual.result.iterations + 1); + assert_eq!(actual.history.len(), actual.result.iterations); + let (sx, sy, sz) = source.positions3().unwrap(); + let (tx, ty, tz) = target.positions3().unwrap(); + let st = KdTree::from_slices(sx, sy, sz); + let tt = KdTree::from_slices(tx, ty, tz); + assert_eq!( + actual.transfers.host_to_device_bytes(), + (32 * (source.len() + + target.len() + + st.layout().nodes().len() + + tt.layout().nodes().len())) as u64 + ); + assert_eq!( + actual.transfers.device_to_host_bytes(), + actual.evaluation_batches as u64 * 256 + ); + let mut error = 0.0; + for i in 0..source.len() { + let p = + actual.result.transform.transform_point(Vec3::new(sx[i], sy[i], sz[i])); + let j = tt.nearest_one(p.x, p.y, p.z).unwrap().index; + error += f64::from(p.x - tx[j]).powi(2) + + f64::from(p.y - ty[j]).powi(2) + + f64::from(p.z - tz[j]).powi(2); + } + error /= source.len() as f64; + assert!( + (actual.result.fitness - error).abs() + <= 1e-9 * f64::from(scale).powi(2) + 1e-4 * error + ); + assert!(actual.history.iter().all(|row| row.correspondences >= 6 + && row.evaluated_correspondences >= 6 + && row.fitness.is_finite() + && row.update_magnitude.is_finite())); + assert_eq!(source, saved_source); + assert_eq!(target, saved_target); + } + } + } + } + eprintln!("GPU_RESIDENT_RADIUS_GICP compared {compared} points in 162 conditions; maximum scale-normalized CPU pose-coordinate error {maximum}"); +} + +#[test] +#[ignore = "requires a native wgpu adapter; run explicitly with --include-ignored"] +fn resident_radius_gicp_rejects_unsupported_settings_and_preserves_concurrent_calls() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let target = reference::corner(85, 1.0); + let (x, y, z) = target.positions3().unwrap(); + let mut b = PointCloudBuilder::xyz(); + for i in (0..target.len()).step_by(2) { + b.push_point([x[i] + 0.015625, y[i], z[i] - 0.0078125]).unwrap(); + } + for i in 0..7 { + b.push_point([100.0 + i as f32, 100.0, 100.0]).unwrap(); + } + let source = b.build().unwrap(); + let registration = reference::registration(1.0, Some(0.25), 0); + let actual = registration + .align_with_gpu_resident_radius_covariances(&runtime, &source, &target, 128) + .unwrap(); + let cpu = registration.align_with_diagnostics(&source, &target).unwrap(); + reference::compare_pose(&source, cpu.transform, actual.result.transform, 1.0); + assert!(actual.history.iter().all(|row| row.correspondences as usize == source.len() - 7)); + std::thread::scope(|scope| { + let callers: Vec<_> = (0..4) + .map(|_| { + scope.spawn(|| { + assert_eq!( + registration + .align_with_gpu_resident_radius_covariances( + &runtime, &source, &target, 128 + ) + .unwrap(), + actual + ); + }) + }) + .collect(); + for c in callers { + c.join().unwrap(); + } + }); + for config in [ + GicpConfig { k_neighbors: 2, ..Default::default() }, + GicpConfig { epsilon: 1e-5, ..Default::default() }, + GicpConfig { epsilon: 1.1, ..Default::default() }, + GicpConfig { epsilon: f64::NAN, ..Default::default() }, + GicpConfig { max_iterations: 0, ..Default::default() }, + GicpConfig { min_correspondences: 5, ..Default::default() }, + GicpConfig { transformation_epsilon: -1.0, ..Default::default() }, + ] { + let config = GicpConfig { covariance_radius: Some(0.25), ..config }; + assert!(GicpRegistration::new(config) + .align_with_gpu_resident_radius_covariances(&runtime, &source, &target, 128) + .is_err()); + } + assert!(GicpRegistration::new(GicpConfig::default()) + .align_with_gpu_resident_radius_covariances(&runtime, &source, &target, 128) + .is_err()); + for radius in [0.0, -1.0, f32::NAN, f32::INFINITY, 1e20, 1e-22] { + assert!(GicpRegistration::new(GicpConfig { + covariance_radius: Some(radius), + ..Default::default() + }) + .align_with_gpu_resident_radius_covariances(&runtime, &source, &target, 128) + .is_err()); + } + for capacity in [0, 1, u32::MAX] { + assert!(registration + .align_with_gpu_resident_radius_covariances(&runtime, &source, &target, capacity) + .is_err()); + } + let far = transform_point_cloud( + &target, + Isometry3::new(Quat::::identity(), Vec3::new(100.0, 0.0, 0.0)), + ) + .unwrap(); + assert!(registration + .align_with_gpu_resident_radius_covariances(&runtime, &far, &target, 128) + .unwrap_err() + .to_string() + .contains("only 0 correspondences")); + let mut b = PointCloudBuilder::xyz(); + b.push_point([f32::INFINITY, 0.0, 0.0]).unwrap(); + let bad = b.build().unwrap(); + assert!(registration + .align_with_gpu_resident_radius_covariances(&runtime, &bad, &target, 128) + .is_err()); +} + +#[test] +#[ignore = "requires a native wgpu adapter; run explicitly with --include-ignored"] +fn resident_radius_gicp_supports_sparse_extent_beyond_dense_grid_capacity() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let cloud = reference::corner(1366, 1.0); + let registration = + GicpRegistration::new(GicpConfig { covariance_radius: Some(0.01), ..Default::default() }); + assert!(registration + .align_with_gpu_reduction(&runtime, &cloud, &cloud) + .unwrap_err() + .to_string() + .contains("cells")); + let actual = registration + .align_with_gpu_resident_radius_covariances(&runtime, &cloud, &cloud, 2) + .unwrap(); + assert!(actual.result.converged); + assert_eq!(actual.result.fitness, 0.0); + assert_eq!(actual.evaluation_batches, 2); + assert_eq!(actual.transfers.device_to_host_bytes(), 512); +} diff --git a/crates/spatialrust-registration/tests/gpu_gicp_resident_radius_covariances_benchmark.rs b/crates/spatialrust-registration/tests/gpu_gicp_resident_radius_covariances_benchmark.rs new file mode 100644 index 0000000..4529ad8 --- /dev/null +++ b/crates/spatialrust-registration/tests/gpu_gicp_resident_radius_covariances_benchmark.rs @@ -0,0 +1,171 @@ +#![cfg(all(feature = "register-gicp-gpu", not(target_arch = "wasm32")))] +#[path = "gpu_gicp_reduction.rs"] +mod reference; +use spatialrust_core::{HasPositions3, PointCloud, PointCloudBuilder}; +use spatialrust_gpu::{GpuGicpInput, WgpuRuntime}; +use spatialrust_math::{Isometry3, Quat, Vec3}; +use spatialrust_registration::{transform_point_cloud, GicpConfig, GicpRegistration}; +use std::time::Instant; + +fn fixture(n: usize, profile: usize) -> PointCloud { + let mut builder = PointCloudBuilder::xyz(); + for i in 0..n { + let j = if profile == 2 { i / 2 } else { i }; + builder + .push_point([ + (j % 32) as f32 * 0.5, + (j / 32 % 32) as f32 * 0.5, + (j / 1024) as f32 * 0.5, + ]) + .unwrap(); + } + builder.build().unwrap() +} + +#[test] +#[ignore = "completed resident radius covariance GICP benchmark; run this exact test alone in release mode"] +fn completed_resident_radius_covariance_gicp_calls() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let info = runtime.adapter_info(); + println!( + "GPU_GICP_RESIDENT_RADIUS_COVARIANCES_ENV backend={} device_type={}", + info.backend, info.device_type + ); + let limits = runtime.device().limits(); + println!("GPU_GICP_RESIDENT_RADIUS_COVARIANCES_LIMITS storage_bytes={} workgroups_per_dimension={} storage_bindings={}",limits.max_storage_buffer_binding_size,limits.max_compute_workgroups_per_dimension,limits.max_storage_buffers_per_shader_stage); + // `host` retains CPU matching/equations/solve; its radius covariances use + // the legacy shared GPU runtime. It is not a pure CPU radius benchmark. + for n in [257, 4097, 16384, 65537] { + 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, + k_neighbors: 20, + transformation_epsilon: 0.0, + fitness_epsilon: 0.0, + covariance_radius: Some([0.51, 0.8, 1.05][profile]), + ..Default::default() + }); + let host = registration.align_with_diagnostics(&source, &target).unwrap(); + let legacy = registration.align_with_gpu_reduction(&runtime, &source, &target).unwrap(); + let gpu = registration + .align_with_gpu_resident_radius_covariances( + &runtime, + &source, + &target, + [32, 64, 128][profile], + ) + .unwrap(); + assert_eq!(gpu.result.iterations, 3); + assert_eq!(gpu.evaluation_batches, 4); + assert_eq!(gpu.covariance_batches, 2); + let maximum = + reference::compare_pose(&source, host.transform, gpu.result.transform, 1.0); + reference::compare_pose(&source, legacy.result.transform, gpu.result.transform, 1.0); + let (sx, sy, sz) = source.positions3().unwrap(); + let (tx, ty, tz) = target.positions3().unwrap(); + let tree = spatialrust_search::KdTree::from_slices(tx, ty, tz); + let source_tree = spatialrust_search::KdTree::from_slices(sx, sy, sz); + let input = GpuGicpInput::upload_with_gpu_radius_covariances( + &runtime, + &source_tree, + &tree, + [0.51, 0.8, 1.05][profile], + [32, 64, 128][profile], + 0.001, + ) + .unwrap(); + assert_eq!( + gpu.transfers.host_to_device_bytes(), + (32 * (source.len() + + target.len() + + source_tree.layout().nodes().len() + + tree.layout().nodes().len())) as u64 + ); + assert_eq!(gpu.transfers.device_to_host_bytes(), 1024); + let expected = input + .equations_resident(&runtime, gpu.result.transform, 0.25) + .unwrap() + .download(&runtime) + .unwrap(); + assert_eq!(expected.correspondences as usize, source.len()); + for _ in 0..3 { + assert_eq!(registration.align_with_diagnostics(&source, &target).unwrap(), host); + assert_eq!( + registration.align_with_gpu_reduction(&runtime, &source, &target).unwrap(), + legacy + ); + assert_eq!( + registration + .align_with_gpu_resident_radius_covariances( + &runtime, + &source, + &target, + [32, 64, 128][profile] + ) + .unwrap(), + gpu + ); + assert_eq!( + input + .equations_resident(&runtime, gpu.result.transform, 0.25) + .unwrap() + .download(&runtime) + .unwrap(), + expected + ); + } + let (mut host_ms, mut legacy_ms, mut gpu_ms, mut resident_ms) = + (vec![], vec![], vec![], vec![]); + for round in 0..20 { + for offset in 0..4 { + let stage = (round + offset) % 4; + let start = Instant::now(); + match stage { + 0 => { + let output = + registration.align_with_diagnostics(&source, &target).unwrap(); + host_ms.push(start.elapsed().as_secs_f64() * 1000.0); + assert_eq!(output, host); + } + 1 => { + let output = registration + .align_with_gpu_reduction(&runtime, &source, &target) + .unwrap(); + legacy_ms.push(start.elapsed().as_secs_f64() * 1000.0); + assert_eq!(output, legacy); + } + 2 => { + let output = registration + .align_with_gpu_resident_radius_covariances( + &runtime, + &source, + &target, + [32, 64, 128][profile], + ) + .unwrap(); + gpu_ms.push(start.elapsed().as_secs_f64() * 1000.0); + assert_eq!(output, gpu); + } + _ => { + let output = input + .equations_resident(&runtime, gpu.result.transform, 0.25) + .unwrap(); + runtime.wait_idle(); + resident_ms.push(start.elapsed().as_secs_f64() * 1000.0); + assert_eq!(output.download(&runtime).unwrap(), expected); + } + } + } + } + println!("GPU_GICP_RESIDENT_RADIUS_COVARIANCES_BENCH n={n} profile={profile} iterations=3 max_pose_error={maximum} radius={} capacity={} host_fitness={} gpu_fitness={} matches={} legacy_h2d={} legacy_d2h={} gpu_h2d={} gpu_d2h={} host_ms={host_ms:?} legacy_ms={legacy_ms:?} gpu_ms={gpu_ms:?} resident_ms={resident_ms:?}",[0.51,0.8,1.05][profile],[32,64,128][profile],host.fitness,gpu.result.fitness,expected.correspondences,legacy.transfers.host_to_device_bytes(),legacy.transfers.device_to_host_bytes(),gpu.transfers.host_to_device_bytes(),gpu.transfers.device_to_host_bytes()); + } + } +} diff --git a/crates/spatialrust/tests/gpu_gicp_resident_radius_covariances_smoke.rs b/crates/spatialrust/tests/gpu_gicp_resident_radius_covariances_smoke.rs new file mode 100644 index 0000000..c77fec4 --- /dev/null +++ b/crates/spatialrust/tests/gpu_gicp_resident_radius_covariances_smoke.rs @@ -0,0 +1,25 @@ +#![cfg(all(feature = "register-gicp-gpu", not(target_arch = "wasm32")))] +#[test] +#[ignore = "requires a native wgpu adapter; run explicitly with --include-ignored"] +fn facade_exposes_resident_gicp_radius_covariance_preparation() { + use spatialrust::{ + gpu::WgpuRuntime, GicpConfig, GicpRegistration, GpuGicpReductionDiagnostics, + PointCloudBuilder, + }; + let mut b = PointCloudBuilder::xyz(); + for i in 0..65 { + b.push_point([(i % 5) as f32 * 0.25, (i / 5 % 5) as f32 * 0.25, (i / 25) as f32 * 0.25]) + .unwrap(); + } + let cloud = b.build().unwrap(); + let runtime = WgpuRuntime::new_headless().unwrap(); + let actual: GpuGicpReductionDiagnostics = + GicpRegistration::new(GicpConfig { covariance_radius: Some(0.5), ..Default::default() }) + .align_with_gpu_resident_radius_covariances(&runtime, &cloud, &cloud, 65) + .unwrap(); + assert!(actual.result.converged); + assert_eq!(actual.result.fitness, 0.0); + assert_eq!(actual.evaluation_batches, 2); + assert_eq!(actual.covariance_batches, 2); + assert_eq!(actual.transfers.device_to_host_bytes(), 512); +} diff --git a/docs/ARCHITECTURE.md b/docs/ARCHITECTURE.md index 82b6ea4..d582a3c 100644 --- a/docs/ARCHITECTURE.md +++ b/docs/ARCHITECTURE.md @@ -67,6 +67,19 @@ epsilon=1e-4..=1 path requires covariance_radius=None; f32 covariance and repeated eigenvalues can change poses/stopping relative to old raw products. See [resident covariance evidence](../notes/2026-10-13_gpu_gicp_resident_covariances.md). +`align_with_gpu_resident_radius_covariances` and +`GpuGicpInput::upload_with_gpu_radius_covariances` extend this preparation +boundary to complete inclusive radius rows with an explicit positive capacity. +Self/duplicates count toward capacity. Overflow zeros BOTH covariance sets +and propagates equation status 16; no truncated neighborhood is consumed. +Sparse rows (fewer than three samples) retain epsilon*I. Other rows use +centered/scaled compensated central sums and 32 max-pivot Jacobi iterations +with physical 1e-20 stopping, so grid accumulation/sweeps and degenerate +eigenvectors may differ. Finite positive radius with normal f32 square and +epsilon=1e-4..=1 are required. Upload/download ownership and receipts match +resident k-NN. CPU trees/6x6 updates remain; legacy preparation is unchanged. +See [resident radius evidence](../notes/2026-10-13_gpu_gicp_resident_radius_covariances.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 15df831..b0ee988 100644 --- a/docs/FEATURE_MATRIX.md +++ b/docs/FEATURE_MATRIX.md @@ -132,9 +132,21 @@ covariance_radius=None, plus existing reduction limits/valid configuration. Malformed neighborhoods or unsupported arithmetic zero all covariance rows; equation status bit 16 rejects their use without fallback. f32 eigenvectors, small/repeated spectra and CPU raw-product rounding can change results. -Radius covariance remains the existing explicit readback/re-upload route. +The existing reduction method retains its radius readback/re-upload route. See [resident covariance evidence](../notes/2026-10-13_gpu_gicp_resident_covariances.md). +`GicpRegistration::align_with_gpu_resident_radius_covariances` and +`GpuGicpInput::upload_with_gpu_radius_covariances` add complete inclusive radius +search, centered/scaled covariance and eigenvectors retained for registration. +Positive explicit capacity includes self/duplicates; overflow rejects all +equations without truncation/fallback. Requires configured positive finite +radius with normal f32 square and epsilon=1e-4..=1. Fewer than three samples +use epsilon*I; other rows use compensated central sums and 32 max-pivot Jacobi +steps at physical 1e-20 stopping. Arithmetic/traversal/repeated eigenvalues can +differ from grid raw means/sweeps. Upload is 32*(S+T+source_nodes+target_nodes), +download is only 256 bytes per evaluation; CPU trees and 6x6 solve remain. +See [resident radius evidence](../notes/2026-10-13_gpu_gicp_resident_radius_covariances.md). + `register-icp-point-to-plane-gpu` adds `PointToPlaneIcp::align_with_gpu_reduction` and diagnostic iteration history. `GpuPointToPlaneInput` under `gpu-wgpu` explicitly retains originals/tree and @@ -179,7 +191,7 @@ before allocation; no dependency, public signature or image readback is added. | `gpu-wgpu` | wgpu runtime and GPU kernels | Explicit GPU operation APIs | | `register-icp-gpu` | ICP GPU correspondence search and resident transforms/moment reductions | Explicit correspondence or untrimmed resident reduction API; CPU pose solver | | `register-icp-point-to-plane-gpu` | Resident transforms, nearest matching, plane Hessian/gradient/fitness reductions | Raw target normal weights; explicit receipts and CPU 6x6 pose solver | -| `register-gicp-gpu` | GICP correspondence/reductions, resident k-NN or explicit radius covariance | Caller-owned runtime and complete receipts; CPU 6x6 pose solver | +| `register-gicp-gpu` | GICP correspondence/reductions, resident k-NN/radius or legacy covariance preparation | Caller-owned runtime and complete receipts; CPU 6x6 pose solver | | `register-fpfh-gpu` | GPU radius search, SPFH and weighted FPFH | Explicit descriptor API and payload receipts; CPU matching/RANSAC | | `feature-iss-gpu` | GPU radius graphs, covariance/eigenvalue saliency and local maxima | Explicit detection/receipts; CPU tree, compaction and attribute gathering | | `filter-mls-gpu` | GPU radius search, weighted frames and order-0/1/2 surface fits | Explicit smoothing/receipts and per-point outcomes; CPU tree and attribute reconstruction | diff --git a/docs/ROADMAP.md b/docs/ROADMAP.md index 1e93071..41afae3 100644 --- a/docs/ROADMAP.md +++ b/docs/ROADMAP.md @@ -1136,4 +1136,5 @@ Dependencies and explicit transfers come before broader algorithm coverage. | GICP resident normal equations | Complete | Reusable positions/covariances, GPU inverse weighting and Hessian/gradient/Euclidean reductions; fixed 256-byte equations with CPU normalized solve/world conversion | Independent f64 inverse/Jacobian/full equations, 108 pose/scale/covariance/stopping cases, strict numeric/ownership/receipt guards and completed CPU/previous-GPU comparisons; `notes/2026-10-13_gpu_gicp_reduction.md` | | Point-to-plane resident normal equations | Complete | Original positions/tree/raw normals, GPU transforms/matching/plane Hessian/gradient/fitness; 256-byte equations with CPU normalized solve/world conversion | Independent raw-normal f64 equations, 108 pose/scale/weight/stopping cases, gating/ownership/schema/receipt guards and completed CPU comparisons; `notes/2026-10-13_gpu_point_to_plane_reduction.md` | | GICP resident k-NN covariances | Complete | Reused query/sample positions and target tree, bounded GPU k-NN/sample covariance/plane eigenvectors with status carried into equations; no covariance readback/re-upload | Independent f64/spectral/Rayleigh checks, 162 pose/scale/k/stopping cases, malformed/status/ownership/receipt guards and completed CPU/previous-GPU comparisons; `notes/2026-10-13_gpu_gicp_resident_covariances.md` | -| Remaining registration updates | Planned | Trimming, resident radius covariance preparation and small pose-update boundaries | Pose/convergence/fitness references, numerical conditioning, source/target residency and end-to-end registration timing | +| GICP resident radius covariances | Complete | Complete inclusive GPU radius neighborhoods with explicit capacity, compensated central sums/eigenvectors retained for equations; sparse epsilon*I and overflow rejection | Independent brute-force/f64/spectral references, 162 pose/scale/radius/stopping cases, sparse extent beyond grid capacity, receipts/ownership/status checks and completed API comparisons; `notes/2026-10-13_gpu_gicp_resident_radius_covariances.md` | +| Remaining registration updates | Planned | Trimming and small pose-update boundaries | Pose/convergence/fitness references, numerical conditioning, source/target residency and end-to-end registration timing | diff --git a/notes/2026-10-13_gpu_gicp_resident_radius_covariances.md b/notes/2026-10-13_gpu_gicp_resident_radius_covariances.md new file mode 100644 index 0000000..64e9e9f --- /dev/null +++ b/notes/2026-10-13_gpu_gicp_resident_radius_covariances.md @@ -0,0 +1,246 @@ +# Resident radius GICP covariance preparation — 2026-10-13 + +`GicpRegistration::align_with_gpu_resident_radius_covariances(runtime, source, +target, neighbor_capacity)` under `register-gicp-gpu` removes the explicit +radius covariance readback and re-upload from the previous GPU reduction path. +Both inclusive radius searches, covariance sums and plane eigenvectors remain +on the caller's GPU runtime and feed the existing resident equation loop. +CPU builds both KD-trees, solves normalized 6x6 equations, converts the world +step and composes/stops poses. Ordinary GICP and previous GPU methods retain +their covariance placement and fallback/error contracts. + +## Explicit preparation and transfer boundary + +`GpuGicpInput::upload_with_gpu_radius_covariances(runtime, source_tree, +target_tree, radius, capacity, epsilon)` reuses original source/target query +buffers as covariance samples and shares the target tree with registration. +The temporarily uploaded source tree serves preparation. Two existing radius +passes submit owned rows; a third submission builds/validates covariances then +zeros BOTH sets on any failure. No host wait or neighborhood/covariance +readback/re-upload is introduced. Submitted resources remain alive until GPU +completion; owned equations survive input drop and scratch reuse. + +Capacity is a positive per-point bound, including self and duplicate indices. +Radius membership uses the existing inclusive f32 squared-distance gate and +near-first tree traversal. Rows must have stored=total<=capacity and total<=N; +zero support, invalid indices, nonfinite/negative distances or distances beyond +the radius are rejected. Overflow never becomes a truncated covariance. +Capacity need not be <=64 (unlike k-NN); checked device storage/word/dispatch +limits constrain it. Neighborhood storage per dataset is N*(8+8*capacity) +bytes. This is a fixed-memory explicit contract, with no automatic host-driven +retry or placement choice. + +Covariance storage is 32*(S+T+1) bytes: a zero-reserved 32-byte status prefix +then source/target `[xx,yy,zz,xy,xz,yz,0,0]` rows in original order. Preparation +status 1/2/4/8 means malformed graph, unsupported arithmetic, invalid SPD/weight +bound or radius overflow. Every failure zeros all covariance rows and carries +equation status 16, zeroing all 256-byte equation payload before download rejects +it. Radius overflow details do not require a separate production readback. + +Upload is **32*(S+T+source_nodes+target_nodes)** payload bytes, including source +tree preparation. Only **256*(I+1)** bytes are downloaded for I updates; +`covariance_batches=2`. Device writes/control uniforms are excluded consistently +with existing receipts. The previous reduction method uploads +48*S+64*T+32*target_nodes plus each radius grid's 16*N+4*(cells+1), and downloads +32*(S+T) covariance bytes plus 256*(I+1) equation bytes. The new method removes +that covariance download and upload but adds source-tree staging; receipts +account for the complete tradeoff. + +## Math and limits + +With fewer than three neighbors, the device emits **epsilon*I**, preserving the +existing radius-grid sparse fallback. Coincident support of at least three +selects normal X, rebuilding diagonal `(epsilon,1,1)`. Other support uses the +query origin, actual maximum absolute local extent, compensated two-f32 means +and second-pass central products. Radius covariance uses unnormalized central +sums with physical absolute **1e-20** Jacobi stopping divided by extent squared. +The shared solver tracks eigenvectors for at most 32 maximum-pivot steps and +uses scaled hypot to avoid tangent overflow. The smallest normalized normal +rebuilds `I+(epsilon-1)*n*n^T`; packed Cholesky/inverse trace checks protect the +existing common weight normalization 6/epsilon. Sparse epsilon*I has inverse +trace 3/epsilon, inside that bound. + +Both trees must be nonempty/finite, capacity positive, epsilon as f32 in +1e-4..=1, radius finite/positive with normal finite f32 square. Nonzero local +extent squared must also be normal finite f32. Combined dispatches, node bytes, +neighborhood strides and u32 offsets are checked before allocation. Existing +GICP input/configuration/pose/conditioning restrictions still apply. + +The existing grid uses raw absolute f32 means, plain sums and 16 sweeps of +three rotations. This new explicit precision path uses centered/scaled sums, +a different traversal order and maximum-pivot iteration. It does not promise +bitwise grid equivalence. Eigenvalues with small/repeated gaps have no stable +unique normal; gates, poses, fitness and stopping can differ. The k-NN path +retains count-1 sample covariance and physical 1e-12 stopping, and its previous +f64 component error is unchanged. + +## Correctness evidence + +- Brute-force inclusive memberships, centered f64 covariance, SPD eigenweights + and normal Rayleigh energy cover 1/2/3/63/64/65/257/1025 points, scales + 1e-4/1/1e4 and radii 0.01/0.8/20 times scale: **13,320** source rows and + identical target rows. Unique-spectrum/sparse component error is at most + **1.5026412008e-5**; repeated spectra are checked without forcing an arbitrary + normal. Native legacy grid planes with unique normals match within 3e-5. +- Exact-radius/next-representable membership, sparse epsilon boundaries, + 65 duplicate samples at capacity 64 versus 65, source/target overflow and + tiny local extents exercise validity propagation and fully zero output. + Malformed source/target counts/indices/distances, including out-of-radius + finite distances, invalidate both covariance sets and all equations. +- Radius-derived covariances feed the existing uploaded-covariance equations + with matching Euclidean error and physically rescaled Hessian/gradient + tolerances. Runtime mismatches reject; outputs survive input drop. +- **162** registration conditions cover six cloud sizes, three scales, three + radii and three stopping configurations. Host-update/previous GPU pose, stopping, + fitness, immutable inputs and exact receipts agree within declared tolerance; + maximum scale-normalized host-update pose-coordinate error is **9.536743e-7** across + **148,392** compared points. Unequal clouds/outliers, four concurrent complete + alignments and invalid capacities/configuration are also checked. The host + update reference uses legacy shared-runtime GPU radius covariances; this + comparison does not imply pure CPU radius preparation. +- A 4,098-point sparse cloud with radius 0.01 exceeds the existing dense grid + cell cap but completes the resident-radius route at capacity two with zero + fitness and fixed-size readback. The new sparse representation avoids a + dense grid allocation. Minimal facade exposes the method. + +## Completed-call measurements + +The exact native release test is +`completed_resident_radius_covariance_gicp_calls` in +`C:\Users\rsasa\Workspace\SpatialRust\crates\spatialrust-registration\tests\gpu_gicp_resident_radius_covariances_benchmark.rs`. +Runtime construction is outside timing. Each condition has three warmups and +20 samples per stage, cyclically rotating host-update, previous GPU reduction, +new resident-radius alignment and prepared resident equations. Complete APIs +include validation, CPU trees/grid preparation, uploads, covariance computation, +all waits/readback, CPU solves/pose composition and final fitness. Prepared +equations include submission/completion for one evaluation on pre-uploaded, +precomputed input, excluding preparation/upload, readback and CPU solve. That +stage is recorded separately and is not an alignment latency claim. + +Four sizes (257/4097/16384/65537) use fixed 0.5 spacing, 32x32 grid coordinates. +Profiles 0/1 use unique points, radii 0.51/0.8 and capacities 32/64; profile 2 +duplicates each grid position, radius 1.05, capacity 128 and 17 extra targets. +Every source starts at binary translation `(0.03125,-0.015625,0.0078125)`. +Three updates, distance gate 0.25, epsilon 0.001 and zero stopping thresholds +are fixed. Every warm/timed output is checked against its warmed reference; +final correspondences must retain all source points. + +The **host-update** baseline uses ordinary CPU matching, Hessian/gradient, +fitness and pose solve, with legacy shared-runtime GPU radius covariances. +The **previous GPU** baseline uses explicit grid GPU covariance readback and +re-upload on the caller runtime. Neither baseline is labeled as pure CPU radius +preparation. Surveys run alone without simultaneous local builds, tests or +statistical processing. Backend, device type and relevant limits are recorded +without the local adapter model. + +### Survey 1 + +| N | Radius / capacity | Host updates ms | Previous GPU ms | Resident radius ms | Prepared equations ms | Host/new ratio [95%] | Previous/new ratio [95%] | +| ---: | --- | ---: | ---: | ---: | ---: | --- | --- | +| 257 | 0.51 / 32 | 3.195 | 9.315 | 6.597 | 0.478 | 0.49 [0.48–0.51] | 1.43 [1.39–1.47] | +| 257 | 0.8 / 64 | 3.265 | 9.489 | 6.669 | 0.692 | 0.50 [0.48–0.52] | 1.41 [1.34–1.47] | +| 257 | 1.05 / 128 | 4.146 | 12.016 | 8.560 | 0.616 | 0.49 [0.45–0.53] | 1.36 [1.25–1.45] | +| 4,097 | 0.51 / 32 | 11.646 | 11.338 | 8.067 | 0.819 | 1.54 [1.45–1.65] | 1.44 [1.37–1.51] | +| 4,097 | 0.8 / 64 | 13.040 | 11.757 | 9.262 | 0.835 | 1.39 [1.33–1.47] | 1.25 [1.21–1.28] | +| 4,097 | 1.05 / 128 | 12.774 | 12.282 | 10.604 | 0.777 | 1.18 [1.12–1.24] | 1.14 [1.11–1.17] | +| 16,384 | 0.51 / 32 | 64.579 | 21.636 | 18.509 | 1.729 | 3.48 [3.29–3.68] | 1.14 [1.09–1.19] | +| 16,384 | 0.8 / 64 | 56.786 | 18.855 | 18.636 | 1.792 | 3.02 [2.85–3.18] | 1.00 [0.96–1.04] | +| 16,384 | 1.05 / 128 | 41.507 | 18.955 | 21.746 | 1.730 | 1.96 [1.83–2.11] | 0.87 [0.84–0.90] | +| 65,537 | 0.51 / 32 | 211.292 | 49.043 | 51.044 | 3.529 | 4.12 [3.94–4.31] | 0.97 [0.93–1.01] | +| 65,537 | 0.8 / 64 | 214.940 | 50.190 | 58.759 | 3.645 | 3.70 [3.52–3.88] | 0.88 [0.84–0.93] | +| 65,537 | 1.05 / 128 | 200.584 | 52.286 | 75.540 | 3.552 | 2.57 [2.43–2.73] | 0.68 [0.65–0.72] | + +### Survey 2 + +| N | Radius / capacity | Host updates ms | Previous GPU ms | Resident radius ms | Prepared equations ms | Host/new ratio [95%] | Previous/new ratio [95%] | +| ---: | --- | ---: | ---: | ---: | ---: | --- | --- | +| 257 | 0.51 / 32 | 2.805 | 7.482 | 5.113 | 0.246 | 0.55 [0.53–0.56] | 1.45 [1.42–1.48] | +| 257 | 0.8 / 64 | 2.816 | 7.415 | 5.126 | 0.246 | 0.54 [0.52–0.55] | 1.42 [1.37–1.46] | +| 257 | 1.05 / 128 | 2.819 | 7.521 | 5.213 | 0.248 | 0.54 [0.52–0.55] | 1.45 [1.40–1.49] | +| 4,097 | 0.51 / 32 | 10.401 | 8.828 | 6.452 | 0.425 | 1.66 [1.58–1.74] | 1.39 [1.33–1.45] | +| 4,097 | 0.8 / 64 | 12.142 | 9.814 | 8.040 | 0.475 | 1.50 [1.41–1.59] | 1.24 [1.19–1.29] | +| 4,097 | 1.05 / 128 | 11.971 | 12.157 | 13.633 | 0.889 | 0.95 [0.88–1.02] | 0.90 [0.85–0.94] | +| 16,384 | 0.51 / 32 | 69.468 | 30.023 | 32.468 | 3.079 | 2.10 [1.98–2.21] | 0.91 [0.88–0.94] | +| 16,384 | 0.8 / 64 | 57.557 | 21.499 | 26.536 | 2.042 | 2.01 [1.88–2.16] | 0.77 [0.73–0.81] | +| 16,384 | 1.05 / 128 | 46.626 | 17.234 | 20.600 | 1.194 | 2.12 [1.91–2.33] | 0.80 [0.73–0.85] | +| 65,537 | 0.51 / 32 | 221.118 | 67.137 | 73.491 | 7.073 | 3.22 [2.94–3.58] | 0.93 [0.89–0.96] | +| 65,537 | 0.8 / 64 | 216.502 | 49.141 | 56.165 | 2.227 | 3.37 [2.90–3.89] | 0.83 [0.78–0.89] | +| 65,537 | 1.05 / 128 | 227.242 | 49.142 | 66.858 | 1.744 | 3.00 [2.62–3.37] | 0.68 [0.59–0.75] | + +Times are medians in milliseconds. Ratios are geometric means of paired +baseline/new latency; >1 is faster for the new path. The 95% intervals are +10,000 paired log-ratio bootstrap samples, seeded by +20261013+24*(survey-1)+2*condition_index+baseline_index (host=0, previous=1), +using sorted endpoints 249/9749. All raw samples, including adverse results, +are preserved. Intervals describe each survey, not backend/device portability, +cold-runtime latency or a universal crossover. Large within/between-survey +variance also appears in the separate prepared-equation stage. + +For 257 points, the new route is 1.36–1.45x faster than previous GPU reduction +but remains slower than host updates (host/new 0.49–0.55). At 4,097 points, +smaller radii improve on previous GPU by 1.24–1.44x; duplicated/larger support +changes from 1.14x in survey 1 to 0.90x in survey 2. At 16,384 points, gains do +not repeat across surveys and duplicated support is consistently slower. +At 65,537 points, host updates/new ratios are 2.57–4.12x, but previous GPU/new +ratios are only 0.68–0.97: **this route does not establish a large-cloud latency +improvement over the existing GPU implementation**. Source-tree staging and +fixed-capacity neighborhood allocation/computation are additional costs; +these completed-call measurements do not isolate their individual shares. +This is an opt-in residency/sparse-grid-capacity improvement, with measured +latency tradeoffs, and does not change default placement. + +All warmed and timed outputs matched their warmed references. Every final +prepared evaluation retained N correspondences. Maximum measured host/new +pose-coordinate error is **3.7252903e-9**; host Euclidean fitness is zero and +new fitness is at most **5.4620291e-19**. Relative error around zero is not a +useful precision measure. Both native surveys passed (37.43/39.98 seconds, +excluding compilation), using Vulkan/DiscreteGpu and reported storage limit +2,147,483,647 bytes, maximum dispatch dimension 65,535; required storage +bindings are eight. No local adapter model is recorded. + +| N | Radius | Previous H2D bytes | New H2D bytes | Previous D2H bytes | New D2H bytes | +| ---: | ---: | ---: | ---: | ---: | ---: | +| 257 | 0.51 | 40,056 | 18,560 | 17,472 | 1,024 | +| 257 | 0.8 | 39,032 | 18,560 | 17,472 | 1,024 | +| 257 | 1.05 | 40,632 | 20,064 | 18,016 | 1,024 | +| 4,097 | 0.51 | 637,144 | 295,040 | 263,232 | 1,024 | +| 4,097 | 0.8 | 615,992 | 295,040 | 263,232 | 1,024 | +| 4,097 | 1.05 | 610,640 | 296,672 | 263,776 | 1,024 | +| 16,384 | 0.51 | 2,540,128 | 1,179,584 | 1,049,600 | 1,024 | +| 16,384 | 0.8 | 2,456,808 | 1,179,584 | 1,049,600 | 1,024 | +| 16,384 | 1.05 | 2,434,456 | 1,181,216 | 1,050,144 | 1,024 | +| 65,537 | 0.51 | 10,183,856 | 4,718,720 | 4,195,392 | 1,024 | +| 65,537 | 0.8 | 9,830,712 | 4,718,720 | 4,195,392 | 1,024 | +| 65,537 | 1.05 | 9,730,760 | 4,720,352 | 4,195,936 | 1,024 | + +For 65,537 source points the previous radius-covariance D2H is +4,195,392–4,195,936 bytes; the new completed three-update call downloads +**1,024 bytes** (99.9756% less). H2D is 4,718,720–4,720,352 versus +9,730,760–10,183,856 bytes (51.49–53.67% less), including the new source tree. +The resident radius/neighborhood buffers still consume device memory; fewer +transfers do not imply faster complete calls or a smaller total memory footprint. + +## Final verification and delivery + +- GPU all-feature library: **200 passed, 20 ignored**; integration targets: + **27 passed, 3 ignored**, 227 passing total. Radius covariance tests and + expanded source/target malformed-row tests pass, as do previous k-NN and + independent f64 GICP equations. +- Registration all-feature library: **35 passed, 1 ignored**. Three new native + radius integration tests pass, plus previous k-NN's 162 conditions and both + previous reduction tests' 108 conditions. New minimal facade smoke passes. +- All-feature/all-target GPU/registration Clippy, CPU-only GICP Clippy, WASM + optional registration, minimal facade, benchmark-target Clippy and final + workspace formatting pass. WASM is compile checked, not browser executed. +- No core/math/search code, dependencies, unsafe or automatic placement rules + were added. CPU GICP and the common post-preparation pose/fitness/update loop + are unchanged against parent 5df2a62. +- Local measurement JSON: + `C:\Users\rsasa\Workspace\SpatialRust\target\automation\gpu-gicp-resident-radius-covariances-measurements.json`. + Raw surveys: + `C:\Users\rsasa\Workspace\SpatialRust\target\automation\gpu-gicp-resident-radius-covariances-bench-1.log` + and + `C:\Users\rsasa\Workspace\SpatialRust\target\automation\gpu-gicp-resident-radius-covariances-bench-2.log`. + Publication files/logs omit the local adapter model. Existing 124 named CI + checks remain required before merge, together with the pinned parent slice. + From 5688fc6ae76ada972d312c6498214ba7b3aa240f Mon Sep 17 00:00:00 2001 From: rsasaki0109 Date: Tue, 13 Oct 2026 09:21:09 +0900 Subject: [PATCH 2/2] Trim trailing measurement-note whitespace --- notes/2026-10-13_gpu_gicp_resident_radius_covariances.md | 1 - 1 file changed, 1 deletion(-) diff --git a/notes/2026-10-13_gpu_gicp_resident_radius_covariances.md b/notes/2026-10-13_gpu_gicp_resident_radius_covariances.md index 64e9e9f..dd4e00b 100644 --- a/notes/2026-10-13_gpu_gicp_resident_radius_covariances.md +++ b/notes/2026-10-13_gpu_gicp_resident_radius_covariances.md @@ -243,4 +243,3 @@ transfers do not imply faster complete calls or a smaller total memory footprint `C:\Users\rsasa\Workspace\SpatialRust\target\automation\gpu-gicp-resident-radius-covariances-bench-2.log`. Publication files/logs omit the local adapter model. Existing 124 named CI checks remain required before merge, together with the pinned parent slice. -