diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 73d86219..04c2ca95 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -122,6 +122,10 @@ jobs: args: -p spatialrust-filtering --features filter-mls - name: spatialrust-filter-mls-meta args: -p spatialrust --features filter-mls + - name: spatialrust-filter-mls-gpu + args: -p spatialrust-filtering --features filter-mls-gpu + - name: spatialrust-filter-mls-gpu-meta + args: -p spatialrust --features filter-mls-gpu - name: spatialrust-filter-fps args: -p spatialrust-filtering --features filter-fps - name: spatialrust-filter-fps-meta diff --git a/CHANGELOG.md b/CHANGELOG.md index b20de50d..518cf41d 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 opt-in `filter-mls-gpu` smoothing for constant, linear and quadratic + surfaces. Reusable input-bound radius rows feed GPU weighted frames and + polynomial fits in three passes; explicit upload/download receipts and + per-point unchanged/plane/polynomial outcomes accompany the result. CPU tree + construction and schema/attribute reconstruction remain. Normalized f32 + Cholesky uses a different rank test from CPU f64; nearly singular fits can + produce large coordinate differences, documented with adverse fixtures. + Incomplete/malformed rows or numerical errors invalidate all device output. + - Add opt-in `feature-iss-gpu` keypoint detection: shared uploaded positions feed GPU radius graphs, compensated covariance, eigenvalue saliency and local maxima. Explicit payload receipts cover upload and flag/status download; diff --git a/README.md b/README.md index d64a36b0..a416559a 100644 --- a/README.md +++ b/README.md @@ -661,7 +661,7 @@ FPFH descriptors have an explicit GPU path under `register-fpfh-gpu`: search, SPFH and weighted FPFH on device and reports all payload transfers. Resident APIs reuse uploaded input and neighborhoods. Capacity overflow returns an error; GPU rounding can change angular bins near boundaries. CPU matching/ -RANSAC and ordinary descriptor defaults are retained; MLS remains CPU. +RANSAC and ordinary descriptor defaults are retained. See [FPFH measurements](notes/2026-10-12_gpu_fpfh_descriptors.md). ISS keypoints have an explicit `feature-iss-gpu` path: @@ -672,6 +672,16 @@ remain host stages. Capacity includes self; incomplete graphs reject. GPU f32 arithmetic can change selection near thresholds and NMS ties compared with CPU f64. See [ISS validation and measurements](notes/2026-10-13_gpu_iss_keypoints.md). +MLS smoothing has an explicit `filter-mls-gpu` path: +`smoother.filter_with_gpu(&runtime, &cloud, neighbor_capacity)` runs radius search, +weighted local frames and order-0/1/2 surface fits on GPU. The result preserves +schema, metadata and attributes and reports transfers, per-point modes and +unchanged/plane/polynomial counts. Resident APIs reuse uploaded input and radius +rows. CPU tree construction and output reconstruction remain host stages. +GPU f32 rank tests can choose a different fallback from CPU f64, with large +differences on nearly singular fits. See +[MLS precision and measurements](notes/2026-10-13_gpu_mls_smoothing.md). + ### Library Load or save by file extension: diff --git a/crates/spatialrust-filtering/Cargo.toml b/crates/spatialrust-filtering/Cargo.toml index 60209aa0..ca311671 100644 --- a/crates/spatialrust-filtering/Cargo.toml +++ b/crates/spatialrust-filtering/Cargo.toml @@ -16,6 +16,7 @@ filter-outlier = ["dep:spatialrust-search"] filter-crop = [] filter-fps = [] filter-mls = ["dep:spatialrust-search"] +filter-mls-gpu = ["filter-mls", "dep:spatialrust-gpu", "spatialrust-gpu/gpu-wgpu"] [dependencies] spatialrust-core = { workspace = true } diff --git a/crates/spatialrust-filtering/src/lib.rs b/crates/spatialrust-filtering/src/lib.rs index b068a676..1604d0c7 100644 --- a/crates/spatialrust-filtering/src/lib.rs +++ b/crates/spatialrust-filtering/src/lib.rs @@ -31,6 +31,9 @@ pub use fps::{FarthestPointSampling, FarthestPointSamplingConfig}; #[cfg(feature = "filter-mls")] pub use mls::{MlsConfig, MlsSmoothing}; +#[cfg(feature = "filter-mls-gpu")] +pub use mls::{GpuMlsFilterResult, GpuMlsPositionResult}; + #[cfg(feature = "filter-outlier")] pub use outlier::{ RadiusOutlierConfig, RadiusOutlierRemoval, StatisticalOutlierConfig, StatisticalOutlierRemoval, diff --git a/crates/spatialrust-filtering/src/mls.rs b/crates/spatialrust-filtering/src/mls.rs index ad8f9bf9..4d76c0c4 100644 --- a/crates/spatialrust-filtering/src/mls.rs +++ b/crates/spatialrust-filtering/src/mls.rs @@ -40,6 +40,34 @@ impl MlsConfig { } } +/// Explicit GPU MLS positions, device outcomes and payload receipt. +#[cfg(feature = "filter-mls-gpu")] +#[derive(Clone, Debug, PartialEq)] +pub struct GpuMlsPositionResult { + /// Original-order explicitly downloaded GPU coordinates. + pub positions: Vec>, + /// Per-point device outcomes aligned with positions. + pub modes: Vec, + /// Device fitting/fallback counts. + pub diagnostics: spatialrust_gpu::GpuMlsDiagnostics, + /// Tree/query upload and padded position/status download payload. + pub transfers: spatialrust_core::TransferStats, +} + +/// Explicit GPU MLS filtered cloud, outcome counts and payload receipt. +#[cfg(feature = "filter-mls-gpu")] +#[derive(Clone, Debug, PartialEq)] +pub struct GpuMlsFilterResult { + /// Smoothed cloud preserving schema, metadata and non-position fields. + pub cloud: PointCloud, + /// Per-point device outcomes aligned with cloud rows. + pub modes: Vec, + /// Device fitting/fallback counts. + pub diagnostics: spatialrust_gpu::GpuMlsDiagnostics, + /// Input upload and position/status download payload. + pub transfers: spatialrust_core::TransferStats, +} + /// Moving Least Squares smoothing filter. #[derive(Clone, Copy, Debug, PartialEq)] pub struct MlsSmoothing { @@ -69,6 +97,82 @@ impl MlsSmoothing { self.smoothed_positions_with_workers(input, parallel_worker_count(input.len())) } + /// Explicitly uploads input, computes GPU radius neighborhoods, weighted + /// frames and order-0/1/2 polynomial fits, and downloads positions/status. + /// CPU tree construction remains; no automatic CPU fallback. Capacity includes + /// self and must fit all radius matches. GPU configuration requires finite + /// XYZ/extents and a positive finite nonzero f32 radius square, stricter than + /// CPU MLS. GPU normalized f32 math/rank tests can change coordinates and + /// fit/fallback decisions compared with CPU f64. Numerical underflow errors + /// propagate; insufficient support and singular fits use device fallbacks. + #[cfg(feature = "filter-mls-gpu")] + pub fn smoothed_positions_with_gpu( + &self, + runtime: &spatialrust_gpu::WgpuRuntime, + input: &PointCloud, + neighbor_capacity: u32, + ) -> SpatialResult { + let radius = self.config.search_radius; + if !radius.is_finite() + || radius <= 0.0 + || !(radius * radius).is_finite() + || radius * radius == 0.0 + || self.config.polynomial_order > 2 + { + return Err(SpatialError::InvalidArgument( + "GPU MLS requires a positive finite nonzero radius square and order 0..=2" + .to_owned(), + )); + } + let (x, y, z) = input.positions3()?; + if x.iter().chain(y).chain(z).any(|v| !v.is_finite()) { + return Err(SpatialError::InvalidArgument("GPU MLS requires finite XYZ".to_owned())); + } + let tree = KdTree::from_slices(x, y, z); + let device_input = spatialrust_gpu::GpuMlsInput::upload(runtime, &tree)?; + let rows = device_input.neighborhoods_resident(runtime, radius, neighbor_capacity)?; + let output = device_input.smooth_resident( + runtime, + &rows, + self.config.polynomial_order, + u32::try_from(self.config.min_neighbors).unwrap_or(u32::MAX), + )?; + let readback = output.download(runtime)?; + let mut transfers = device_input.upload_transfers(); + transfers.record( + spatialrust_core::TransferDirection::DeviceToHost, + readback.transfers.device_to_host_bytes(), + ); + let positions = + readback.positions.into_iter().map(|p| Vec3::new(p[0], p[1], p[2])).collect(); + Ok(GpuMlsPositionResult { + positions, + modes: readback.modes, + diagnostics: readback.diagnostics, + transfers, + }) + } + + /// Explicit GPU MLS filtering with CPU schema/attribute reconstruction. The + /// ordinary filter trait retains CPU execution. See the position API for + /// capacity, precision and numeric-error contracts; all payload is reported. + #[cfg(feature = "filter-mls-gpu")] + pub fn filter_with_gpu( + &self, + runtime: &spatialrust_gpu::WgpuRuntime, + input: &PointCloud, + neighbor_capacity: u32, + ) -> SpatialResult { + let result = self.smoothed_positions_with_gpu(runtime, input, neighbor_capacity)?; + let cloud = build_output(input, &result.positions)?; + Ok(GpuMlsFilterResult { + cloud, + modes: result.modes, + diagnostics: result.diagnostics, + transfers: result.transfers, + }) + } + fn smoothed_positions_with_workers( &self, input: &PointCloud, diff --git a/crates/spatialrust-filtering/tests/gpu_mls.rs b/crates/spatialrust-filtering/tests/gpu_mls.rs new file mode 100644 index 00000000..5f1319be --- /dev/null +++ b/crates/spatialrust-filtering/tests/gpu_mls.rs @@ -0,0 +1,344 @@ +#![cfg(all(feature = "filter-mls-gpu", not(target_arch = "wasm32")))] + +use spatialrust_core::{ + DType, FieldSemantic, HasPositions3, PointBuffer, PointBufferSet, PointCloud, PointField, + PointSchema, SpatialMetadata, StandardSchemas, Timestamp, +}; +use spatialrust_filtering::{GpuMlsPositionResult, MlsConfig, MlsSmoothing}; +use spatialrust_gpu::{GpuMlsInput, GpuMlsMode, WgpuRuntime}; +use spatialrust_math::Vec3; + +pub fn fixture(n: usize, kind: usize) -> PointCloud { + fixture_with_width(n, kind, if kind == 4 { 128 } else { 32 }) +} +pub fn fixture_with_width(n: usize, kind: usize, width: usize) -> PointCloud { + let mut columns: [Vec; 3] = std::array::from_fn(|_| Vec::with_capacity(n)); + for i in 0..n { + let j = if kind == 2 { i / 2 } else { i }; + let x = (j % width) as f32 * 0.05 - 0.8; + let y = (j / width) as f32 * 0.05 - 0.4; + let noise = ((j * 37 % 101) as f32 / 100.0 - 0.5) * 0.02; + let z = match kind { + 0 => noise, + 1 | 2 | 4 => 0.015 * (x * x + y * y) + noise, + _ => 0.0, + }; + for (column, value) in columns.iter_mut().zip([x, y, z]) { + column.push(value); + } + } + from_columns(columns) +} +pub fn from_columns(columns: [Vec; 3]) -> PointCloud { + let mut buffers = PointBufferSet::new(); + for (name, values) in ["x", "y", "z"].into_iter().zip(columns) { + buffers.insert(name, PointBuffer::F32(values)); + } + PointCloud::try_from_parts(StandardSchemas::point_xyz(), buffers, Default::default()).unwrap() +} +pub fn compare(actual: &[Vec3], expected: &[Vec3], radius: f32) -> f32 { + assert_eq!(actual.len(), expected.len()); + let mut maximum = 0.0_f32; + for (i, (a, b)) in actual.iter().zip(expected).enumerate() { + for (axis, (a, b)) in [a.x, a.y, a.z].into_iter().zip([b.x, b.y, b.z]).enumerate() { + assert!(a.is_finite()); + let error = (a - b).abs(); + maximum = maximum.max(error); + let ulp = f32::from_bits(b.abs().to_bits() + 1) - b.abs(); + let bound = 5e-5 + radius * 1e-3 + 8.0 * ulp; + assert!( + error <= bound, + "point={i} axis={axis} gpu={a} cpu={b} error={error} bound={bound}" + ); + } + } + maximum +} + +// Validate the coordinates of each reported outcome against independent CPU +// math for that outcome; retain error against the requested CPU order separately. +pub fn audit( + cloud: &PointCloud, + config: MlsConfig, + actual: &GpuMlsPositionResult, + requested: &[Vec3], +) -> (f32, f32) { + let plane = MlsSmoothing::new(MlsConfig { polynomial_order: 0, ..config }) + .smoothed_positions(cloud) + .unwrap(); + let (x, y, z) = cloud.positions3().unwrap(); + let mut reference = Vec::with_capacity(cloud.len()); + let mut requested_error = 0.0_f32; + for (i, mode) in actual.modes.iter().enumerate() { + reference.push(match mode { + GpuMlsMode::Unchanged => Vec3::new(x[i], y[i], z[i]), + GpuMlsMode::PlaneFallback => plane[i], + GpuMlsMode::Polynomial => requested[i], + }); + if *mode == GpuMlsMode::Unchanged { + let p = actual.positions[i]; + assert_eq!( + [p.x.to_bits(), p.y.to_bits(), p.z.to_bits()], + [x[i].to_bits(), y[i].to_bits(), z[i].to_bits()] + ); + } + let a = actual.positions[i]; + let b = requested[i]; + requested_error = + requested_error.max((a.x - b.x).abs().max((a.y - b.y).abs()).max((a.z - b.z).abs())); + } + (compare(&actual.positions, &reference, config.search_radius), requested_error) +} + +#[test] +#[ignore = "requires a native wgpu adapter; run explicitly with --include-ignored"] +fn gpu_mls_compares_all_orders_coordinates_modes_and_receipts_with_cpu() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let mut coordinates = 0; + let mut maximum = 0.0_f32; + let mut requested_maximum = 0.0_f32; + for n in [0, 1, 2, 17, 63, 64, 65, 257, 4097] { + for kind in 0..5 { + let cloud = fixture(n, kind); + for order in 0..=2 { + let config = + MlsConfig { search_radius: 0.16, polynomial_order: order, min_neighbors: 6 }; + let smoother = MlsSmoothing::new(config); + let expected = smoother.smoothed_positions(&cloud).unwrap(); + let actual = smoother.smoothed_positions_with_gpu(&runtime, &cloud, 128).unwrap(); + let (difference, requested_error) = audit(&cloud, config, &actual, &expected); + requested_maximum = requested_maximum.max(requested_error); + if n == 4097 && matches!(kind, 0..=2) && order == 2 { + assert!( + actual.diagnostics.polynomial as usize >= n * 95 / 100, + "well-supported surfaces must use quadratic fits: {:?}", + actual.diagnostics + ); + } + maximum = maximum.max(difference); + coordinates += n * 3; + assert_eq!(actual.modes.len(), n); + for (mode, count) in [ + (GpuMlsMode::Unchanged, actual.diagnostics.unchanged), + (GpuMlsMode::PlaneFallback, actual.diagnostics.plane_fallback), + (GpuMlsMode::Polynomial, actual.diagnostics.polynomial), + ] { + assert_eq!(actual.modes.iter().filter(|&&m| m == mode).count(), count as usize); + } + let tree = spatialrust_search::KdTree::from_point_cloud(&cloud).unwrap(); + let nodes = tree.layout().nodes().len(); + assert_eq!( + actual.transfers.host_to_device_bytes(), + if n == 0 { 0 } else { (n * 32 + nodes * 32) as u64 } + ); + assert_eq!( + actual.transfers.device_to_host_bytes(), + if n == 0 { 0 } else { (n * 16 + 32) as u64 } + ); + eprintln!("GPU_MLS_ACCURACY n={n} kind={kind} order={order} max_mode_reference_error={difference} max_requested_cpu_error={requested_error} unchanged={} plane={} polynomial={}",actual.diagnostics.unchanged,actual.diagnostics.plane_fallback,actual.diagnostics.polynomial); + } + } + } + eprintln!("GPU_MLS compared {coordinates} coordinates across all orders; maximum mode-reference difference {maximum}; maximum requested-order CPU difference {requested_maximum}"); +} + +#[test] +#[ignore = "requires a native wgpu adapter; run explicitly with --include-ignored"] +fn gpu_mls_preserves_no_support_and_validates_capacity_configuration_and_extent() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let cloud = fixture(65, 1); + for radius in [0.0, -1.0, f32::NAN, f32::INFINITY, f32::MAX, f32::MIN_POSITIVE] { + assert!(MlsSmoothing::new(MlsConfig::with_radius(radius)) + .smoothed_positions_with_gpu(&runtime, &cloud, 128) + .is_err()); + } + assert!(MlsSmoothing::new(MlsConfig { polynomial_order: 3, ..Default::default() }) + .smoothed_positions_with_gpu(&runtime, &cloud, 128) + .is_err()); + let smoother = MlsSmoothing::new(MlsConfig::with_radius(0.16)); + for capacity in [0, 1, 3] { + assert!(smoother.smoothed_positions_with_gpu(&runtime, &cloud, capacity).is_err()); + } + for columns in [ + [vec![f32::NAN], vec![0.0], vec![0.0]], + [vec![f32::MAX, -f32::MAX], vec![0.0; 2], vec![0.0; 2]], + ] { + assert!(smoother + .smoothed_positions_with_gpu(&runtime, &from_columns(columns), 128) + .is_err()); + } + let input = from_columns([vec![-0.0, 1.0, 2.0], vec![0.0; 3], vec![-0.0; 3]]); + for order in 0..=2 { + let result = MlsSmoothing::new(MlsConfig { + search_radius: 3.0, + polynomial_order: order, + min_neighbors: usize::MAX, + }) + .smoothed_positions_with_gpu(&runtime, &input, 3) + .unwrap(); + assert_eq!(result.diagnostics.unchanged, 3); + let (x, y, z) = input.positions3().unwrap(); + for (i, p) in result.positions.iter().enumerate() { + assert_eq!( + [p.x.to_bits(), p.y.to_bits(), p.z.to_bits()], + [x[i].to_bits(), y[i].to_bits(), z[i].to_bits()] + ); + } + } + let translated = from_columns([ + (0..25).map(|i| 1e6 + (i % 5) as f32 * 0.25).collect(), + (0..25).map(|i| -1e6 + (i / 5) as f32 * 0.25).collect(), + vec![7.0; 25], + ]); + let s = MlsSmoothing::new(MlsConfig::with_radius(1.0)); + compare( + &s.smoothed_positions_with_gpu(&runtime, &translated, 25).unwrap().positions, + &s.smoothed_positions(&translated).unwrap(), + 1.0, + ); + let tiny = from_columns([vec![0.0, 1e-22, 2e-22], vec![0.0; 3], vec![0.0; 3]]); + assert!(MlsSmoothing::new(MlsConfig { + search_radius: 1.0, + min_neighbors: 0, + ..Default::default() + }) + .smoothed_positions_with_gpu(&runtime, &tiny, 3) + .is_err()); +} + +#[test] +#[ignore = "requires a native wgpu adapter; run explicitly with --include-ignored"] +fn gpu_mls_reduces_plane_noise_and_preserves_all_attribute_bits() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let mut points = Vec::new(); + let mut seed = 12345_u64; + for x in 0..25 { + for y in 0..25 { + seed = seed.wrapping_mul(6_364_136_223_846_793_005).wrapping_add(1); + let noise = (((seed >> 32) as u32 as f32 / u32::MAX as f32) * 2.0 - 1.0) * 0.02; + points.push([x as f32 * 0.05, y as f32 * 0.05, noise]); + } + } + let mut schema = PointSchema::new(); + let mut buffers = PointBufferSet::new(); + for (axis, (name, semantic)) in [ + ("px", FieldSemantic::PositionX), + ("py", FieldSemantic::PositionY), + ("pz", FieldSemantic::PositionZ), + ] + .into_iter() + .enumerate() + { + schema = schema.with_field(PointField::scalar(name, semantic, DType::F32)); + buffers.insert(name, PointBuffer::F32(points.iter().map(|p| p[axis]).collect())); + } + for (name, dtype, buffer) in [ + ( + "float", + DType::F32, + PointBuffer::F32( + (0..625) + .map(|i| f32::from_bits(if i % 2 == 0 { 0x7fc00037 } else { 0x80000000 })) + .collect(), + ), + ), + ( + "double", + DType::F64, + PointBuffer::F64( + (0..625) + .map(|i| { + f64::from_bits(if i % 2 == 0 { + 0x7ff8000000000037 + } else { + 0x8000000000000000 + }) + }) + .collect(), + ), + ), + ("byte", DType::U8, PointBuffer::U8((0..625).map(|i| i as u8).collect())), + ("short", DType::U16, PointBuffer::U16((0..625).map(|i| i as u16).collect())), + ("index", DType::U32, PointBuffer::U32((0..625).map(|i| i as u32).collect())), + ("signed", DType::I32, PointBuffer::I32((0..625).map(|i| -i).collect())), + ] { + schema = schema.with_field(PointField::scalar(name, FieldSemantic::Unknown, dtype)); + buffers.insert(name, buffer); + } + let cloud = PointCloud::try_from_parts( + schema, + buffers, + SpatialMetadata::new("map", Timestamp::from_nanos(37)), + ) + .unwrap(); + let smoother = + MlsSmoothing::new(MlsConfig { search_radius: 0.2, polynomial_order: 1, min_neighbors: 6 }); + let output = smoother.filter_with_gpu(&runtime, &cloud, 128).unwrap(); + assert_eq!(output.cloud.schema(), cloud.schema()); + assert_eq!(output.cloud.metadata(), cloud.metadata()); + let (_, _, z) = output.cloud.positions3().unwrap(); + let interior: Vec<_> = + (0..625).filter(|&i| (4..21).contains(&(i / 25)) && (4..21).contains(&(i % 25))).collect(); + let before = (interior.iter().map(|&i| points[i][2].powi(2)).sum::() + / interior.len() as f32) + .sqrt(); + let after = + (interior.iter().map(|&i| z[i].powi(2)).sum::() / interior.len() as f32).sqrt(); + assert!(after < before * 0.6, "noise RMS before={before} after={after}"); + for field in cloud.schema().fields().iter().filter(|f| f.semantic == FieldSemantic::Unknown) { + match (cloud.field(&field.name).unwrap(), output.cloud.field(&field.name).unwrap()) { + (PointBuffer::F32(a), PointBuffer::F32(b)) => assert_eq!( + a.iter().map(|v| v.to_bits()).collect::>(), + b.iter().map(|v| v.to_bits()).collect::>() + ), + (PointBuffer::F64(a), PointBuffer::F64(b)) => assert_eq!( + a.iter().map(|v| v.to_bits()).collect::>(), + b.iter().map(|v| v.to_bits()).collect::>() + ), + (a, b) => assert_eq!(a, b), + } + } + eprintln!("GPU_MLS plane RMS before={before} after={after}"); +} + +#[test] +#[ignore = "requires a native wgpu adapter; run explicitly with --include-ignored"] +fn gpu_mls_retains_outputs_and_reuses_scratch_with_concurrent_runtime_guards() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let other = WgpuRuntime::new_headless().unwrap(); + let cloud = fixture(257, 1); + let tree = spatialrust_search::KdTree::from_point_cloud(&cloud).unwrap(); + let input = GpuMlsInput::upload(&runtime, &tree).unwrap(); + let different = GpuMlsInput::upload(&runtime, &tree).unwrap(); + let rows = input.neighborhoods_resident(&runtime, 0.16, 128).unwrap(); + assert!(input.neighborhoods_resident(&other, 0.16, 128).is_err()); + assert!(input.smooth_resident(&other, &rows, 2, 6).is_err()); + assert!(different.smooth_resident(&runtime, &rows, 2, 6).is_err()); + let retained = input.smooth_resident(&runtime, &rows, 2, 6).unwrap(); + let expected = retained.download(&runtime).unwrap(); + std::thread::scope(|scope| { + let calls: Vec<_> = (0..4) + .map(|_| { + scope.spawn(|| { + for _ in 0..2 { + let output = input + .smooth_resident(&runtime, &rows, 2, 6) + .unwrap() + .download(&runtime) + .unwrap(); + assert_eq!(output, expected); + } + }) + }) + .collect(); + for call in calls { + call.join().unwrap(); + } + }); + drop(rows); + drop(input); + drop(tree); + assert_eq!(retained.download(&runtime).unwrap(), expected); + assert!(retained.download(&other).is_err()); + assert!(runtime.scratch_buffer_pool_stats().cached_resources >= 2); +} diff --git a/crates/spatialrust-filtering/tests/gpu_mls_benchmark.rs b/crates/spatialrust-filtering/tests/gpu_mls_benchmark.rs new file mode 100644 index 00000000..ae997e5c --- /dev/null +++ b/crates/spatialrust-filtering/tests/gpu_mls_benchmark.rs @@ -0,0 +1,89 @@ +#![cfg(all(feature = "filter-mls-gpu", not(target_arch = "wasm32")))] + +#[path = "gpu_mls.rs"] +mod reference; + +use spatialrust_core::HasPositions3; +use spatialrust_filtering::{GpuMlsPositionResult, MlsConfig, MlsSmoothing, PointCloudFilter}; +use spatialrust_gpu::{GpuMlsInput, WgpuRuntime}; +use spatialrust_math::Vec3; +use std::time::Instant; + +#[test] +#[ignore = "completed MLS benchmark; run this exact test alone in release mode"] +fn completed_mls_calls() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let info = runtime.adapter_info(); + println!("GPU_MLS_ENV backend={} device_type={}", info.backend, info.device_type); + for n in [257, 4097, 16384, 65537] { + for kind in 0..3 { + let order = if kind == 0 { 1 } else { 2 }; + let capacity = if kind == 2 { 96 } else { 64 }; + let config = + MlsConfig { search_radius: 0.16, polynomial_order: order, min_neighbors: 6 }; + let cloud = reference::fixture_with_width(n, kind, 128); + let smoother = MlsSmoothing::new(config); + let cpu_expected = smoother.filter(&cloud).unwrap(); + let cpu_positions = smoother.smoothed_positions(&cloud).unwrap(); + let gpu_expected = smoother.filter_with_gpu(&runtime, &cloud, capacity).unwrap(); + let (x, y, z) = gpu_expected.cloud.positions3().unwrap(); + let positions: Vec<_> = + x.iter().zip(y).zip(z).map(|((&x, &y), &z)| Vec3::new(x, y, z)).collect(); + let audited = GpuMlsPositionResult { + positions: positions.clone(), + modes: gpu_expected.modes.clone(), + diagnostics: gpu_expected.diagnostics, + transfers: gpu_expected.transfers, + }; + let (mode_reference_error, requested_cpu_error) = + reference::audit(&cloud, config, &audited, &cpu_positions); + let tree = spatialrust_search::KdTree::from_point_cloud(&cloud).unwrap(); + let input = GpuMlsInput::upload(&runtime, &tree).unwrap(); + let rows = input.neighborhoods_resident(&runtime, 0.16, capacity).unwrap(); + let resident_expected = input + .smooth_resident(&runtime, &rows, order, 6) + .unwrap() + .download(&runtime) + .unwrap(); + for _ in 0..3 { + assert_eq!(smoother.filter(&cloud).unwrap(), cpu_expected); + assert_eq!( + smoother.filter_with_gpu(&runtime, &cloud, capacity).unwrap(), + gpu_expected + ); + assert_eq!( + input + .smooth_resident(&runtime, &rows, order, 6) + .unwrap() + .download(&runtime) + .unwrap(), + resident_expected + ); + } + let (mut cpu_ms, mut gpu_ms, mut resident_ms) = (vec![], vec![], vec![]); + for round in 0..20 { + for offset in 0..3 { + let stage = (round + offset) % 3; + let start = Instant::now(); + if stage == 2 { + let output = input.smooth_resident(&runtime, &rows, order, 6).unwrap(); + runtime.wait_idle(); + resident_ms.push(start.elapsed().as_secs_f64() * 1000.0); + assert_eq!(output.download(&runtime).unwrap(), resident_expected); + } else if stage == 0 { + let output = smoother.filter(&cloud).unwrap(); + cpu_ms.push(start.elapsed().as_secs_f64() * 1000.0); + assert_eq!(output, cpu_expected); + } else { + let output = smoother.filter_with_gpu(&runtime, &cloud, capacity).unwrap(); + gpu_ms.push(start.elapsed().as_secs_f64() * 1000.0); + assert_eq!(output, gpu_expected); + } + } + } + let workers = spatialrust_search::parallel_worker_count(n); + let d = gpu_expected.diagnostics; + println!("GPU_MLS_BENCH n={n} kind={kind} order={order} workers={workers} capacity={capacity} unchanged={} plane={} polynomial={} max_mode_reference_error={mode_reference_error} max_requested_cpu_error={requested_cpu_error} cpu_ms={cpu_ms:?} gpu_ms={gpu_ms:?} resident_ms={resident_ms:?}",d.unchanged,d.plane_fallback,d.polynomial); + } + } +} diff --git a/crates/spatialrust-gpu/src/kernels/mls.rs b/crates/spatialrust-gpu/src/kernels/mls.rs new file mode 100644 index 00000000..486d7c56 --- /dev/null +++ b/crates/spatialrust-gpu/src/kernels/mls.rs @@ -0,0 +1,473 @@ +use crate::{scratch::ScratchBuffers, GpuKdTree, GpuNeighborhoods, GpuQueryBatch, WgpuRuntime}; +use bytemuck::{Pod, Zeroable}; +use spatialrust_core::{SpatialError, SpatialResult, TransferDirection, TransferStats}; +use spatialrust_search::KdTree; +use std::sync::Arc; +use wgpu::util::DeviceExt; + +#[repr(C)] +#[derive(Clone, Copy, Pod, Zeroable)] +struct Params { + count: u32, + capacity: u32, + min_neighbors: u32, + padding: u32, + radius_squared: f32, + weight_scale_squared: f32, + unused: [f32; 2], +} + +fn invalid() -> SpatialError { + SpatialError::InvalidArgument("GPU MLS requires matching input/runtime ownership, finite XYZ/extents, positive finite nonzero radius squares, order 0..=2 and supported buffers/dispatch".to_owned()) +} + +fn valid_radius(radius: f32) -> bool { + radius.is_finite() && radius > 0.0 && (radius * radius).is_finite() && radius * radius > 0.0 +} + +fn plan(points: u64, limits: &wgpu::Limits) -> SpatialResult { + let bytes = points.checked_mul(192).ok_or_else(invalid)?; + let bound = limits.max_buffer_size.min(u64::from(limits.max_storage_buffer_binding_size)); + if points > u64::from(u32::MAX) + || bytes > bound + || bytes > isize::MAX as u64 + || bound < 192 + || limits.max_compute_workgroup_size_x < 64 + || limits.max_compute_workgroup_size_y == 0 + || limits.max_compute_workgroup_size_z == 0 + || limits.max_compute_invocations_per_workgroup < 64 + || points.div_ceil(64) > u64::from(limits.max_compute_workgroups_per_dimension) + || limits.max_storage_buffers_per_shader_stage < 6 + || limits.max_bindings_per_bind_group < 7 + || limits.max_bind_groups == 0 + || limits.max_uniform_buffers_per_shader_stage == 0 + || limits.max_uniform_buffer_binding_size < 32 + || points + .checked_mul(16) + .and_then(|n| n.checked_add(32)) + .map_or(true, |n| n > limits.max_buffer_size || n > isize::MAX as u64) + { + return Err(invalid()); + } + Ok(points as u32) +} + +/// Explicitly uploaded MLS positions/tree/queries, reusable on one runtime. +/// CPU tree construction is external; subsequent search and surface fitting run on GPU. +pub struct GpuMlsInput { + tree: GpuKdTree, + queries: GpuQueryBatch, + count: u32, + identity: Arc<()>, + input_identity: Arc<()>, + transfers: TransferStats, +} + +impl GpuMlsInput { + /// Uploads a finite CPU tree and original-order queries. The query buffer is + /// also the covariance position buffer: 32 bytes per point plus 32 per tree + /// node, with no duplicate sample upload. Finite f32 extent squares required. + pub fn upload(runtime: &WgpuRuntime, tree: &KdTree) -> SpatialResult { + let layout = tree.layout(); + let (x, y, z) = layout.positions(); + let count = plan(x.len() as u64, &runtime.device().limits())?; + if x.iter().chain(y).chain(z).any(|v| !v.is_finite()) { + return Err(invalid()); + } + if count > 0 { + let extent: f64 = [x, y, z] + .iter() + .map(|axis| { + let lo = axis.iter().copied().fold(f32::INFINITY, f32::min); + let hi = axis.iter().copied().fold(f32::NEG_INFINITY, f32::max); + (f64::from(hi) - f64::from(lo)).powi(2) + }) + .sum(); + if !(extent as f32).is_finite() { + return Err(invalid()); + } + } + let positions: Vec<_> = (0..x.len()).map(|i| [x[i], y[i], z[i]]).collect(); + let tree = GpuKdTree::upload(runtime, tree)?; + let queries = GpuQueryBatch::upload(runtime, &positions)?; + let mut transfers = tree.upload_transfers(); + transfers.record( + TransferDirection::HostToDevice, + queries.upload_transfers().host_to_device_bytes(), + ); + Ok(Self { + tree, + queries, + count, + identity: runtime.profiling_identity(), + input_identity: Arc::new(()), + transfers, + }) + } + + /// Original point count. + #[must_use] + pub const fn point_count(&self) -> u32 { + self.count + } + + /// Payload explicitly uploaded at construction (zero for empty input). + #[must_use] + pub const fn upload_transfers(&self) -> TransferStats { + self.transfers + } + + /// Builds reusable radius rows tied to this input, including self. Every + /// match is counted; fitting rejects overflow. No readback or completion wait. + pub fn neighborhoods_resident( + &self, + runtime: &WgpuRuntime, + radius: f32, + capacity: u32, + ) -> SpatialResult { + if !Arc::ptr_eq(&self.identity, &runtime.profiling_identity()) || !valid_radius(radius) { + return Err(invalid()); + } + let rows = self.tree.radius_neighbors_resident(runtime, &self.queries, radius, capacity)?; + Ok(GpuMlsNeighborhoods { rows, radius, input_identity: Arc::clone(&self.input_identity) }) + } + + /// Fits Gaussian-weighted local frames and constant/linear/quadratic surfaces + /// on GPU, excluding the query itself after the minimum-support test. Sparse + /// points stay unchanged; insufficient polynomial support/ill-conditioned + /// solves use the plane on device. The normalized f32 solve uses a relative + /// Cholesky rank threshold of 1e-6; CPU f64 uses Gaussian elimination with + /// absolute pivot threshold 1e-12. + /// Coordinates and fit/fallback decisions can differ near rank/eigenvalue + /// degeneracy. Unsupported extent underflow is a numerical error. Three passes + /// share one submission, without payload copies/host wait. Invalid rows or + /// numerics zero ALL output; GPU consumers must check status. + pub fn smooth_resident( + &self, + runtime: &WgpuRuntime, + neighbors: &GpuMlsNeighborhoods, + polynomial_order: u8, + min_neighbors: u32, + ) -> SpatialResult { + if !Arc::ptr_eq(&self.identity, &runtime.profiling_identity()) + || !Arc::ptr_eq(&self.input_identity, &neighbors.input_identity) + || polynomial_order > 2 + { + return Err(invalid()); + } + plan(u64::from(self.count), &runtime.device().limits())?; + let create = |label, size| { + runtime.device().create_buffer(&wgpu::BufferDescriptor { + label: Some(label), + size, + usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_SRC, + mapped_at_creation: false, + }) + }; + let result = GpuMlsPositions { + buffer: create("mls-positions", (u64::from(self.count) * 16).max(16)), + status: create("mls-status", 32), + count: self.count, + identity: Arc::clone(&self.identity), + }; + if self.count == 0 { + return Ok(result); + } + let cache = runtime.mls_pipelines.get_or_init(|| MlsPipelines::new(runtime.device())); + let params = Params { + count: self.count, + capacity: neighbors.rows.capacity(), + min_neighbors, + padding: 0, + radius_squared: neighbors.radius * neighbors.radius, + weight_scale_squared: (neighbors.radius * 0.5).powi(2).max(1e-12), + unused: [0.0; 2], + }; + let uniform = runtime.device().create_buffer_init(&wgpu::util::BufferInitDescriptor { + label: Some("mls-params"), + contents: bytemuck::bytes_of(¶ms), + usage: wgpu::BufferUsages::UNIFORM, + }); + let mut scratch = ScratchBuffers::new(runtime); + let frames = scratch.acquire("mls-frames", u64::from(self.count) * 64); + let systems = scratch.acquire("mls-systems", u64::from(self.count) * 192); + let buffers = [ + &uniform, + self.queries.buffer().unwrap(), + neighbors.rows.buffer(), + &frames, + &systems, + &result.buffer, + &result.status, + ]; + let entries: Vec<_> = buffers + .iter() + .enumerate() + .map(|(binding, buffer)| wgpu::BindGroupEntry { + binding: binding as u32, + resource: buffer.as_entire_binding(), + }) + .collect(); + let group = runtime.device().create_bind_group(&wgpu::BindGroupDescriptor { + label: Some("mls-bindings"), + layout: &cache.layout, + entries: &entries, + }); + let mut encoder = runtime.device().create_command_encoder(&Default::default()); + for pipeline in [&cache.frames, &cache.fits[polynomial_order as usize], &cache.finalize] { + let mut pass = encoder.begin_compute_pass(&wgpu::ComputePassDescriptor { + label: Some("mls-compute"), + timestamp_writes: None, + }); + pass.set_pipeline(pipeline); + pass.set_bind_group(0, &group, &[]); + pass.dispatch_workgroups(self.count.div_ceil(64), 1, 1); + } + scratch.submit_and_recycle(encoder); + Ok(result) + } +} + +/// Radius rows tied to one immutable GPU MLS input. Capacity includes self. +pub struct GpuMlsNeighborhoods { + rows: GpuNeighborhoods, + radius: f32, + input_identity: Arc<()>, +} +impl GpuMlsNeighborhoods { + /// Underlying rows/counts for device consumers. Treat storage as immutable. + #[must_use] + pub fn rows(&self) -> &GpuNeighborhoods { + &self.rows + } +} + +/// Device counts of MLS outcomes; valid only when the first five status words are zero. +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)] +pub struct GpuMlsDiagnostics { + /// Points retained without fitting (sparse/zero support or zero extent). + pub unchanged: u32, + /// Points using the plane after insufficient polynomial support or rank failure. + pub plane_fallback: u32, + /// Points using the requested constant/linear/quadratic polynomial fit. + pub polynomial: u32, +} + +/// Per-point outcome of device MLS processing. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum GpuMlsMode { + /// Original coordinates retained without fitting. + Unchanged, + /// Plane projection after insufficient support or failed rank test. + PlaneFallback, + /// Requested constant/linear/quadratic fit evaluated on device. + Polynomial, +} + +/// Explicit position/status download with outcome counts and payload accounting. +#[derive(Clone, Debug, PartialEq)] +pub struct GpuMlsReadback { + /// Original-order smoothed XYZ. + pub positions: Vec<[f32; 3]>, + /// Per-point outcomes aligned with positions; no additional device transfer. + pub modes: Vec, + /// Device fit/fallback counts. + pub diagnostics: GpuMlsDiagnostics, + /// Payload copied during this explicit readback. + pub transfers: TransferStats, +} + +/// Owned resident MLS positions/modes and global status. No buffers are pooled; +/// outputs survive input drops and other calls. GPU consumers must check status. +pub struct GpuMlsPositions { + buffer: wgpu::Buffer, + status: wgpu::Buffer, + count: u32, + identity: Arc<()>, +} +impl GpuMlsPositions { + /// Original-order vec4: XYZ and mode (0 unchanged, 1 plane, 2 polynomial). + /// Invalid results are zeroed; empty input has 16 unused bytes. + #[must_use] + pub fn buffer(&self) -> &wgpu::Buffer { + &self.buffer + } + /// Eight u32 words: bad rows, required capacity, malformed rows, last bad row + /// plus one, numerical errors, unchanged/plane/polynomial counts. First five + /// words must ALL be zero; remaining counts are diagnostics for valid output. + #[must_use] + pub fn status_buffer(&self) -> &wgpu::Buffer { + &self.status + } + /// Original point/output count. + #[must_use] + pub const fn point_count(&self) -> u32 { + self.count + } + /// Position/mode and status readback bytes including padding; empty zero. + #[must_use] + pub const fn readback_bytes(&self) -> u64 { + if self.count == 0 { + 0 + } else { + self.count as u64 * 16 + 32 + } + } + /// Explicitly downloads positions/modes and status in one submission, waits + /// for completion and rejects invalid status/nonfinite rows. Payload is + /// 16 bytes per point plus 32 status bytes, empty input zero. No graph readback. + pub fn download(&self, runtime: &WgpuRuntime) -> SpatialResult { + if !Arc::ptr_eq(&self.identity, &runtime.profiling_identity()) { + return Err(invalid()); + } + if self.count == 0 { + return Ok(GpuMlsReadback { + positions: vec![], + modes: vec![], + diagnostics: Default::default(), + transfers: Default::default(), + }); + } + let bytes = u64::from(self.count) * 16; + let staging = runtime.device().create_buffer(&wgpu::BufferDescriptor { + label: Some("mls-readback"), + size: bytes + 32, + usage: wgpu::BufferUsages::COPY_DST | wgpu::BufferUsages::MAP_READ, + mapped_at_creation: false, + }); + let mut encoder = runtime.device().create_command_encoder(&Default::default()); + encoder.copy_buffer_to_buffer(&self.buffer, 0, &staging, 0, bytes); + encoder.copy_buffer_to_buffer(&self.status, 0, &staging, bytes, 32); + runtime.queue().submit(Some(encoder.finish())); + let slice = staging.slice(..); + let (sender, receiver) = std::sync::mpsc::channel(); + slice.map_async(wgpu::MapMode::Read, move |r| { + let _ = sender.send(r); + }); + runtime.device().poll(wgpu::Maintain::Wait); + receiver + .recv() + .map_err(|e| SpatialError::InvalidArgument(e.to_string()))? + .map_err(|e| SpatialError::InvalidArgument(e.to_string()))?; + let mapped = slice.get_mapped_range(); + let words: &[u32] = bytemuck::cast_slice(&mapped); + let status = &words[bytes as usize / 4..]; + let decoded = (|| { + if status[..5].iter().any(|&v| v != 0) + || u64::from(status[5]) + u64::from(status[6]) + u64::from(status[7]) + != u64::from(self.count) + { + return Err(SpatialError::InvalidArgument(format!("GPU MLS rejected incomplete/malformed neighborhoods or numerical output; status={status:?}"))); + } + let mut positions = Vec::with_capacity(self.count as usize); + let mut modes = Vec::with_capacity(self.count as usize); + let mut counts = [0u32; 3]; + for row in words[..bytes as usize / 4].chunks_exact(4) { + let values: [f32; 4] = std::array::from_fn(|i| f32::from_bits(row[i])); + if values.iter().any(|v| !v.is_finite()) || ![0.0, 1.0, 2.0].contains(&values[3]) { + return Err(invalid()); + } + counts[values[3] as usize] += 1; + modes.push(match values[3] as u32 { + 0 => GpuMlsMode::Unchanged, + 1 => GpuMlsMode::PlaneFallback, + _ => GpuMlsMode::Polynomial, + }); + positions.push([values[0], values[1], values[2]]); + } + if counts != status[5..] { + return Err(invalid()); + } + Ok(( + positions, + modes, + GpuMlsDiagnostics { + unchanged: status[5], + plane_fallback: status[6], + polynomial: status[7], + }, + )) + })(); + drop(mapped); + staging.unmap(); + let (positions, modes, diagnostics) = decoded?; + let mut transfers = TransferStats::default(); + transfers.record(TransferDirection::DeviceToHost, bytes + 32); + Ok(GpuMlsReadback { positions, modes, diagnostics, transfers }) + } +} + +pub(crate) struct MlsPipelines { + layout: wgpu::BindGroupLayout, + frames: wgpu::ComputePipeline, + fits: [wgpu::ComputePipeline; 3], + finalize: wgpu::ComputePipeline, +} +impl MlsPipelines { + fn new(device: &wgpu::Device) -> Self { + let entries: Vec<_> = (0..7) + .map(|binding| wgpu::BindGroupLayoutEntry { + binding, + visibility: wgpu::ShaderStages::COMPUTE, + ty: wgpu::BindingType::Buffer { + ty: if binding == 0 { + wgpu::BufferBindingType::Uniform + } else { + wgpu::BufferBindingType::Storage { read_only: binding < 3 } + }, + has_dynamic_offset: false, + min_binding_size: None, + }, + count: None, + }) + .collect(); + let layout = device.create_bind_group_layout(&wgpu::BindGroupLayoutDescriptor { + label: Some("mls-layout"), + entries: &entries, + }); + let pipeline_layout = device.create_pipeline_layout(&wgpu::PipelineLayoutDescriptor { + label: Some("mls-layout"), + bind_group_layouts: &[&layout], + push_constant_ranges: &[], + }); + let common = include_str!("../shaders/mls_common.wgsl"); + let frame_shader = device.create_shader_module(wgpu::ShaderModuleDescriptor { + label: Some("mls-frame-shader"), + source: wgpu::ShaderSource::Wgsl( + format!("{common}{}", include_str!("../shaders/mls_frames.wgsl")).into(), + ), + }); + let create = |shader: &wgpu::ShaderModule, entry| { + device.create_compute_pipeline(&wgpu::ComputePipelineDescriptor { + label: Some("mls-pipeline"), + layout: Some(&pipeline_layout), + module: shader, + entry_point: Some(entry), + compilation_options: Default::default(), + cache: None, + }) + }; + let frames = create(&frame_shader, "prepare_frames"); + let fits = std::array::from_fn(|order| { + let terms = [1, 3, 6][order]; + let template = if order == 0 { + include_str!("../shaders/mls_constant.wgsl").to_owned() + } else { + include_str!("../shaders/mls_fit.wgsl").replace( + "const FIT_TERMS: u32 = 6u;", + &format!("const FIT_TERMS: u32 = {terms}u;"), + ) + }; + let shader = device.create_shader_module(wgpu::ShaderModuleDescriptor { + label: Some("mls-fit-shader"), + source: wgpu::ShaderSource::Wgsl(format!("{common}{template}").into()), + }); + create(&shader, "fit_surface") + }); + Self { frames, fits, finalize: create(&frame_shader, "finalize"), layout } + } +} + +#[cfg(all(test, not(target_arch = "wasm32")))] +#[path = "mls_tests.rs"] +mod tests; diff --git a/crates/spatialrust-gpu/src/kernels/mls_tests.rs b/crates/spatialrust-gpu/src/kernels/mls_tests.rs new file mode 100644 index 00000000..2d571c14 --- /dev/null +++ b/crates/spatialrust-gpu/src/kernels/mls_tests.rs @@ -0,0 +1,167 @@ +use super::*; + +fn read_words(runtime: &WgpuRuntime, source: &wgpu::Buffer, bytes: u64) -> Vec { + let staging = runtime.device().create_buffer(&wgpu::BufferDescriptor { + label: None, + size: bytes, + usage: wgpu::BufferUsages::MAP_READ | wgpu::BufferUsages::COPY_DST, + mapped_at_creation: false, + }); + let mut encoder = runtime.device().create_command_encoder(&Default::default()); + encoder.copy_buffer_to_buffer(source, 0, &staging, 0, bytes); + runtime.queue().submit(Some(encoder.finish())); + let slice = staging.slice(..); + let (sender, receiver) = std::sync::mpsc::channel(); + slice.map_async(wgpu::MapMode::Read, move |r| { + let _ = sender.send(r); + }); + runtime.device().poll(wgpu::Maintain::Wait); + receiver.recv().unwrap().unwrap(); + let mapped = slice.get_mapped_range(); + let words = bytemuck::cast_slice(&mapped).to_vec(); + drop(mapped); + staging.unmap(); + words +} + +fn poison(runtime: &WgpuRuntime, buffer: &wgpu::Buffer) { + let shader=runtime.device().create_shader_module(wgpu::ShaderModuleDescriptor { + label:None,source:wgpu::ShaderSource::Wgsl("@group(0) @binding(0) var words:array; @compute @workgroup_size(1) fn poison(){words[2]=0xffffffffu;}".into()), + }); + let pipeline = runtime.device().create_compute_pipeline(&wgpu::ComputePipelineDescriptor { + label: None, + layout: None, + module: &shader, + entry_point: Some("poison"), + compilation_options: Default::default(), + cache: None, + }); + 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 { + 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())); +} + +#[test] +fn mls_plan_checks_counts_storage_bindings_and_dispatch() { + let limits = wgpu::Limits::default(); + assert!(plan(u64::MAX, &limits).is_err()); + assert!(plan(u64::from(u32::MAX) + 1, &limits).is_err()); + let maximum = (u64::from(limits.max_storage_buffer_binding_size) / 192) + .min(u64::from(limits.max_compute_workgroups_per_dimension) * 64); + assert!(plan(maximum, &limits).is_ok()); + assert!(plan(maximum + 1, &limits).is_err()); + for case in 0..6 { + let mut restricted = limits.clone(); + match case { + 0 => restricted.max_compute_workgroup_size_x = 63, + 1 => restricted.max_storage_buffers_per_shader_stage = 5, + 2 => restricted.max_bindings_per_bind_group = 6, + 3 => restricted.max_uniform_buffer_binding_size = 31, + 4 => restricted.max_buffer_size = 191, + _ => restricted.max_compute_workgroups_per_dimension = 0, + } + assert!(plan(1, &restricted).is_err()); + } +} + +#[test] +fn mls_upload_parameters_and_empty_receipts_are_checked() { + let runtime = WgpuRuntime::new_headless().unwrap(); + for x in [vec![f32::NAN], vec![f32::MAX, -f32::MAX]] { + let zero = vec![0.0; x.len()]; + assert!(GpuMlsInput::upload(&runtime, &KdTree::from_slices(&x, &zero, &zero)).is_err()); + } + let input = GpuMlsInput::upload(&runtime, &KdTree::from_slices(&[], &[], &[])).unwrap(); + assert_eq!(input.upload_transfers(), TransferStats::default()); + for radius in [0.0, -1.0, f32::NAN, f32::INFINITY, f32::MAX, f32::MIN_POSITIVE] { + assert!(input.neighborhoods_resident(&runtime, radius, 0).is_err()); + } + let rows = input.neighborhoods_resident(&runtime, 1.0, 0).unwrap(); + assert!(input.smooth_resident(&runtime, &rows, 3, 0).is_err()); + for order in 0..=2 { + let output = input.smooth_resident(&runtime, &rows, order, 0).unwrap(); + assert_eq!(output.readback_bytes(), 0); + let data = output.download(&runtime).unwrap(); + assert!(data.positions.is_empty()); + assert_eq!(data.diagnostics, GpuMlsDiagnostics::default()); + assert_eq!(data.transfers, TransferStats::default()); + } +} + +#[test] +fn overflow_poison_and_extent_underflow_zero_every_output() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let x: Vec<_> = (0..65).map(|i| i as f32 * 0.125).collect(); + let zero = vec![0.0; 65]; + let input = GpuMlsInput::upload(&runtime, &KdTree::from_slices(&x, &zero, &zero)).unwrap(); + let rows = input.neighborhoods_resident(&runtime, 100.0, 3).unwrap(); + let bad = input.smooth_resident(&runtime, &rows, 2, 0).unwrap(); + assert_eq!(read_words(&runtime, bad.status_buffer(), 32), [65, 65, 0, 65, 0, 0, 0, 0]); + assert!(read_words(&runtime, bad.buffer(), 65 * 16).iter().all(|&v| v == 0)); + assert!(bad.download(&runtime).is_err()); + let rows = input.neighborhoods_resident(&runtime, 0.2, 3).unwrap(); + poison(&runtime, rows.rows().buffer()); + let bad = input.smooth_resident(&runtime, &rows, 2, 6).unwrap(); + assert_eq!(&read_words(&runtime, bad.status_buffer(), 32)[..5], [1, 2, 1, 1, 0]); + assert!(read_words(&runtime, bad.buffer(), 65 * 16).iter().all(|&v| v == 0)); + assert!(bad.download(&runtime).is_err()); + let tiny = GpuMlsInput::upload( + &runtime, + &KdTree::from_slices(&[0.0, 1e-22, 2e-22], &[0.0; 3], &[0.0; 3]), + ) + .unwrap(); + let rows = tiny.neighborhoods_resident(&runtime, 1.0, 3).unwrap(); + let bad = tiny.smooth_resident(&runtime, &rows, 1, 0).unwrap(); + assert_eq!(read_words(&runtime, bad.status_buffer(), 32)[4], 3); + assert!(read_words(&runtime, bad.buffer(), 3 * 16).iter().all(|&v| v == 0)); + assert!(bad.download(&runtime).is_err()); + let rows = input.neighborhoods_resident(&runtime, 0.2, 3).unwrap(); + let good = input.smooth_resident(&runtime, &rows, 2, 6).unwrap().download(&runtime).unwrap(); + assert_eq!(good.diagnostics.unchanged, 65); + assert_eq!(good.transfers.device_to_host_bytes(), 65 * 16 + 32); +} + +#[test] +fn analytic_quadratic_fit_removes_the_query_spike_and_preserves_orders() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let (mut x, mut y, mut z) = (vec![], vec![], vec![]); + for a in -2..=2 { + for b in -2..=2 { + let px = a as f32 * 0.05; + let py = b as f32 * 0.05; + x.push(px); + y.push(py); + z.push(if a == 0 && b == 0 { 0.02 } else { 0.5 * (px * px + py * py) }); + } + } + let input = GpuMlsInput::upload(&runtime, &KdTree::from_slices(&x, &y, &z)).unwrap(); + let rows = input.neighborhoods_resident(&runtime, 0.3, 25).unwrap(); + let linear = input.smooth_resident(&runtime, &rows, 1, 6).unwrap().download(&runtime).unwrap(); + let quad = input.smooth_resident(&runtime, &rows, 2, 6).unwrap().download(&runtime).unwrap(); + assert!(linear.positions[12][2] > 0.002); + assert!(quad.positions[12][2].abs() < 1e-6, "quadratic center {:?}", quad.positions[12]); + assert!(quad.diagnostics.polynomial > 0); + assert_eq!( + quad.diagnostics.unchanged + quad.diagnostics.plane_fallback + quad.diagnostics.polynomial, + 25 + ); + let sparse = + input.smooth_resident(&runtime, &rows, 2, u32::MAX).unwrap().download(&runtime).unwrap(); + assert_eq!(sparse.diagnostics.unchanged, 25); + for (i, p) in sparse.positions.iter().enumerate() { + assert_eq!(p.map(f32::to_bits), [x[i].to_bits(), y[i].to_bits(), z[i].to_bits()]); + } +} diff --git a/crates/spatialrust-gpu/src/kernels/mod.rs b/crates/spatialrust-gpu/src/kernels/mod.rs index 5903f53a..00992bca 100644 --- a/crates/spatialrust-gpu/src/kernels/mod.rs +++ b/crates/spatialrust-gpu/src/kernels/mod.rs @@ -52,7 +52,12 @@ mod voxel_gather; pub(crate) mod cluster_union; pub(crate) mod fpfh; pub(crate) mod iss; +pub(crate) mod mls; pub use iss::{GpuIssInput, GpuIssKeypoints, GpuIssNeighborhoods}; +pub use mls::{ + GpuMlsDiagnostics, GpuMlsInput, GpuMlsMode, GpuMlsNeighborhoods, GpuMlsPositions, + GpuMlsReadback, +}; pub(crate) mod nearest; pub use fpfh::{GpuFpfhDescriptors, GpuFpfhInput, GpuFpfhNeighborhoods, GPU_FPFH_DESCRIPTOR_LEN}; pub use nearest::{ diff --git a/crates/spatialrust-gpu/src/lib.rs b/crates/spatialrust-gpu/src/lib.rs index 5346610f..315b2501 100644 --- a/crates/spatialrust-gpu/src/lib.rs +++ b/crates/spatialrust-gpu/src/lib.rs @@ -53,8 +53,9 @@ pub use device::{GpuDevice, WgpuDevice}; #[cfg(feature = "gpu-wgpu")] pub use kernels::{ GpuFpfhDescriptors, GpuFpfhInput, GpuFpfhNeighborhoods, GpuIssInput, GpuIssKeypoints, - GpuIssNeighborhoods, GpuNeighborhoodKind, GpuNeighborhoodReadback, GpuNeighborhoods, - GpuQueryBatch, GPU_FPFH_DESCRIPTOR_LEN, MAX_GPU_K_NEIGHBORS, + GpuIssNeighborhoods, GpuMlsDiagnostics, GpuMlsInput, GpuMlsMode, GpuMlsNeighborhoods, + GpuMlsPositions, GpuMlsReadback, GpuNeighborhoodKind, GpuNeighborhoodReadback, + GpuNeighborhoods, GpuQueryBatch, GPU_FPFH_DESCRIPTOR_LEN, MAX_GPU_K_NEIGHBORS, }; #[cfg(feature = "gpu-wgpu")] diff --git a/crates/spatialrust-gpu/src/runtime.rs b/crates/spatialrust-gpu/src/runtime.rs index 14d96186..887f636e 100644 --- a/crates/spatialrust-gpu/src/runtime.rs +++ b/crates/spatialrust-gpu/src/runtime.rs @@ -27,6 +27,7 @@ pub struct WgpuRuntime { pub(crate) neighborhood_pipelines: OnceLock, pub(crate) fpfh_pipelines: OnceLock, pub(crate) iss_pipelines: OnceLock, + pub(crate) mls_pipelines: OnceLock, pub(crate) point_estimator_pipelines: OnceLock, max_gather_channels: u32, @@ -392,6 +393,7 @@ impl WgpuRuntime { neighborhood_pipelines: OnceLock::new(), fpfh_pipelines: OnceLock::new(), iss_pipelines: OnceLock::new(), + mls_pipelines: OnceLock::new(), point_estimator_pipelines: OnceLock::new(), max_gather_channels, upload_pool: GpuBufferPool::default(), diff --git a/crates/spatialrust-gpu/src/shaders/mls_common.wgsl b/crates/spatialrust-gpu/src/shaders/mls_common.wgsl new file mode 100644 index 00000000..a79605e2 --- /dev/null +++ b/crates/spatialrust-gpu/src/shaders/mls_common.wgsl @@ -0,0 +1,45 @@ +struct Params { + count: u32, capacity: u32, min_neighbors: u32, padding: u32, + radius_squared: f32, weight_scale_squared: f32, pad0: f32, pad1: f32, +} +struct Status { + bad_rows: atomic, required_capacity: atomic, malformed: atomic, + last_bad_plus_one: atomic, numerical: atomic, + unchanged: atomic, plane: atomic, polynomial: atomic, +} +@group(0) @binding(0) var params: Params; +@group(0) @binding(1) var positions: array>; +@group(0) @binding(2) var neighbors: array; +struct Frame { centroid: vec4, normal: vec4, u: vec4, v: vec4 } +@group(0) @binding(3) var frames: array; +@group(0) @binding(4) var systems: array; +@group(0) @binding(5) var output: array>; +@group(0) @binding(6) var status: Status; + +fn add(sum: vec2, value: f32) -> vec2 { + let s=sum.x+value;let v=s-sum.x; + let e=(sum.x-(s-v))+(value-v)+sum.y;let h=s+e; + return vec2(h,e-(h-s)); +} +fn dot3(a: vec3, b: vec3) -> f32 {return a.x*b.x+a.y*b.y+a.z*b.z;} +fn cross3(a: vec3, b: vec3) -> vec3 { + return vec3(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x); +} +fn unit(v: vec3) -> vec3 {return v/sqrt(dot3(v,v));} +fn finite(value: f32) -> bool {return abs(value)<=3.4028234663852886e38;} +fn reject(i: u32, total: u32, malformed: bool) { + atomicAdd(&status.bad_rows,1u);atomicMax(&status.required_capacity,total); + atomicMax(&status.last_bad_plus_one,i+1u); + if malformed {atomicAdd(&status.malformed,1u);} +} +fn numerical(i: u32) {atomicAdd(&status.numerical,1u);atomicMax(&status.last_bad_plus_one,i+1u);} +fn weight(base: u32, n: u32) -> f32 { + return exp(-bitcast(neighbors[base+3u+2u*n])/params.weight_scale_squared); +} +fn write_result(i: u32, value: vec3, mode: u32) { + if !finite(value.x) || !finite(value.y) || !finite(value.z) {numerical(i);return;} + output[i]=vec4(value,f32(mode)); + if mode==0u {atomicAdd(&status.unchanged,1u);} + else if mode==1u {atomicAdd(&status.plane,1u);} + else {atomicAdd(&status.polynomial,1u);} +} diff --git a/crates/spatialrust-gpu/src/shaders/mls_constant.wgsl b/crates/spatialrust-gpu/src/shaders/mls_constant.wgsl new file mode 100644 index 00000000..2ad80159 --- /dev/null +++ b/crates/spatialrust-gpu/src/shaders/mls_constant.wgsl @@ -0,0 +1,20 @@ +fn fit_mls_point(i: u32) { + if i>=params.count {return;} + if atomicLoad(&status.bad_rows)>0u || atomicLoad(&status.numerical)>0u {return;} + let scale=frames[i].centroid.w;if scale==0.0 {return;} + let centroid=frames[i].centroid.xyz;let normal=frames[i].normal.xyz;let u=frames[i].u.xyz;let v=frames[i].v.xyz; + let p=positions[i].xyz;let base=i*(2u+2u*params.capacity);let stored=neighbors[base]; + var support=0u; + for (var n=0u;n) {fit_mls_point(id.x);} diff --git a/crates/spatialrust-gpu/src/shaders/mls_fit.wgsl b/crates/spatialrust-gpu/src/shaders/mls_fit.wgsl new file mode 100644 index 00000000..9271f1d4 --- /dev/null +++ b/crates/spatialrust-gpu/src/shaders/mls_fit.wgsl @@ -0,0 +1,64 @@ +const FIT_TERMS: u32 = 6u; +fn basis_term(term: u32, s: f32, t: f32) -> f32 { + switch term {case 0u: {return 1.0;} case 1u: {return s;} case 2u: {return t;} + case 3u: {return s*s;} case 4u: {return s*t;} default: {return t*t;}} +} +fn fit_mls_point(i: u32) { + if i>=params.count {return;} + if atomicLoad(&status.bad_rows)>0u || atomicLoad(&status.numerical)>0u {return;} + let scale=frames[i].centroid.w;if scale==0.0 {return;} + let centroid=frames[i].centroid.xyz;let normal=frames[i].normal.xyz;let u=frames[i].u.xyz;let v=frames[i].v.xyz; + let p=positions[i].xyz;let base=i*(2u+2u*params.capacity);let stored=neighbors[base]; + var support=0u; + for (var n=0u;n) {fit_mls_point(id.x);} diff --git a/crates/spatialrust-gpu/src/shaders/mls_frames.wgsl b/crates/spatialrust-gpu/src/shaders/mls_frames.wgsl new file mode 100644 index 00000000..bfe8f366 --- /dev/null +++ b/crates/spatialrust-gpu/src/shaders/mls_frames.wgsl @@ -0,0 +1,79 @@ +fn smallest_normal(matrix: array,3>, epsilon: f32) -> vec3 { + var a=matrix; + var v=array,3>(vec3(1.0,0.0,0.0),vec3(0.0,1.0,0.0),vec3(0.0,0.0,1.0)); + for (var iteration=0u;iteration<32u;iteration+=1u) { + var p=0u;var q=1u;var off=abs(a[0][1]); + if abs(a[0][2])>off {p=0u;q=2u;off=abs(a[0][2]);} + if abs(a[1][2])>off {p=1u;q=2u;off=abs(a[1][2]);} + if off==0.0 || off(delta)&0x80000000u)!=0u); + let t=sign_delta*h/(abs(d)+sqrt(d*d+h*h));let c=1.0/sqrt(1.0+t*t);let s=t*c; + a[p][p]=c*c*app-2.0*s*c*apq+s*s*aqq; + a[q][q]=s*s*app+2.0*s*c*apq+c*c*aqq;a[p][q]=0.0;a[q][p]=0.0; + for (var r=0u;r<3u;r+=1u) { + if r!=p && r!=q { + let arp=a[r][p];let arq=a[r][q];a[r][p]=c*arp-s*arq;a[p][r]=a[r][p]; + a[r][q]=s*arp+c*arq;a[q][r]=a[r][q]; + } + let vrp=v[r][p];let vrq=v[r][q];v[r][p]=c*vrp-s*vrq;v[r][q]=s*vrp+c*vrq; + } + } + var index=0u; + if a[1][1](v[0][index],v[1][index],v[2][index])); +} + +fn prepare_frame(i: u32) { + if i>=params.count {return;} + output[i]=vec4(0.0); + frames[i]=Frame(vec4(0.0),vec4(0.0),vec4(0.0),vec4(0.0)); + let base=i*(2u+2u*params.capacity);let stored=neighbors[base];let total=neighbors[base+1u]; + if stored>params.capacity || totalparams.count {reject(i,total,true);return;} + if stored!=total {reject(i,total,false);return;} + var support=0u;var scale=0.0;let p=positions[i].xyz; + for (var n=0u;n(neighbors[base+3u+2u*n]); + if j>=params.count || !(distance>=0.0 && distance<=params.radius_squared) {reject(i,total,true);return;} + if j!=i { + support+=1u;let d=abs(positions[j].xyz-p);scale=max(scale,max(d.x,max(d.y,d.z))); + } + } + if stored(0.0);var mean: array,3>; + for (var n=0u;n; + for (var r=0u;r<3u;r+=1u) {offset[r]=(mean[r].x+mean[r].y)/sum_weight;} + let centroid=p+offset*scale; + var covariance: array,6>; + for (var n=0u;n(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 b=0u;b<6u;b+=1u) {covariance[b]=add(covariance[b],w*products[b]);} + } + var cov: array; + for (var b=0u;b<6u;b+=1u) {cov[b]=(covariance[b].x+covariance[b].y)/sum_weight;} + let normal=smallest_normal(array,3>(vec3(cov[0],cov[3],cov[4]),vec3(cov[3],cov[1],cov[5]),vec3(cov[4],cov[5],cov[2])),(1e-12/scale_squared)/sum_weight); + var helper=vec3(1.0,0.0,0.0); + if abs(normal.x)>=0.9 {helper=vec3(0.0,1.0,0.0);} + let u=unit(cross3(normal,helper));let v=cross3(normal,u); + frames[i]=Frame(vec4(centroid,scale),vec4(normal,0.0),vec4(u,0.0),vec4(v,0.0)); +} +@compute @workgroup_size(64) +fn prepare_frames(@builtin(global_invocation_id) id: vec3) {prepare_frame(id.x);} +@compute @workgroup_size(64) +fn finalize(@builtin(global_invocation_id) id: vec3) { + let i=id.x;if i>=params.count {return;} + if atomicLoad(&status.bad_rows)>0u || atomicLoad(&status.numerical)>0u {output[i]=vec4(0.0);} +} diff --git a/crates/spatialrust/Cargo.toml b/crates/spatialrust/Cargo.toml index c0252a7e..ab4193f2 100644 --- a/crates/spatialrust/Cargo.toml +++ b/crates/spatialrust/Cargo.toml @@ -44,6 +44,7 @@ filter-mls = [ "spatialrust-filtering/filter-mls", "spatialrust-search/search-kdtree", ] +filter-mls-gpu = ["filter-mls", "spatialrust-filtering/filter-mls-gpu", "gpu-wgpu"] feature-normal = ["spatialrust-features/feature-normal", "spatialrust-search/search-kdtree"] feature-iss = ["spatialrust-features/feature-iss", "spatialrust-search/search-kdtree"] feature-iss-gpu = ["feature-iss", "spatialrust-features/feature-iss-gpu", "gpu-wgpu"] diff --git a/crates/spatialrust/src/lib.rs b/crates/spatialrust/src/lib.rs index 01576562..6826e717 100644 --- a/crates/spatialrust/src/lib.rs +++ b/crates/spatialrust/src/lib.rs @@ -167,6 +167,9 @@ pub use spatialrust_filtering::{FarthestPointSampling, FarthestPointSamplingConf #[cfg(feature = "filter-mls")] pub use spatialrust_filtering::{MlsConfig, MlsSmoothing}; +#[cfg(feature = "filter-mls-gpu")] +pub use spatialrust_filtering::{GpuMlsFilterResult, GpuMlsPositionResult}; + #[cfg(feature = "feature-normal")] pub use spatialrust_features::{ orient_normal_towards_viewpoint, FeatureEstimator, KdTreeNeighborhood, NeighborhoodProvider, diff --git a/crates/spatialrust/tests/gpu_mls_smoke.rs b/crates/spatialrust/tests/gpu_mls_smoke.rs new file mode 100644 index 00000000..4849bcc8 --- /dev/null +++ b/crates/spatialrust/tests/gpu_mls_smoke.rs @@ -0,0 +1,30 @@ +#![cfg(all(feature = "filter-mls-gpu", not(target_arch = "wasm32")))] + +use spatialrust::{HasPositions3, MlsConfig, MlsSmoothing, PointCloudBuilder}; +use spatialrust_gpu::WgpuRuntime; + +#[test] +#[ignore = "requires a native wgpu adapter; run explicitly with --include-ignored"] +fn facade_gpu_mls_exposes_quadratic_surface_and_complete_receipts() { + let mut builder = PointCloudBuilder::xyz(); + for x in -2..=2 { + for y in -2..=2 { + let x = x as f32 * 0.05; + let y = y as f32 * 0.05; + let z = if x == 0.0 && y == 0.0 { 0.02 } else { 0.5 * (x * x + y * y) }; + builder.push_point([x, y, z]).unwrap(); + } + } + let cloud = builder.build().unwrap(); + let runtime = WgpuRuntime::new_headless().unwrap(); + let smoother = MlsSmoothing::new(MlsConfig::with_radius(0.3)); + let output = smoother.filter_with_gpu(&runtime, &cloud, 25).unwrap(); + let (_, _, z) = output.cloud.positions3().unwrap(); + assert!(z[12].abs() < 1e-6); + assert_eq!(output.cloud.schema(), cloud.schema()); + assert_eq!(output.cloud.metadata(), cloud.metadata()); + assert_eq!(output.modes.len(), 25); + assert!(output.transfers.host_to_device_bytes() >= 25 * 32); + assert_eq!(output.transfers.device_to_host_bytes(), 25 * 16 + 32); + assert!(smoother.filter_with_gpu(&runtime, &cloud, 1).is_err()); +} diff --git a/docs/ARCHITECTURE.md b/docs/ARCHITECTURE.md index 7e0ea80a..50cdb204 100644 --- a/docs/ARCHITECTURE.md +++ b/docs/ARCHITECTURE.md @@ -30,7 +30,7 @@ Nearest-k rows retain CPU stable ordering; radius rows retain traversal order and expose overflow. Complete CPU readback rejects incomplete neighborhoods; GPU consumers must inspect counts before using rows. Target/query payloads remain resident between passes. This adds no core/search dependency and does -not change feature algorithm placement; MLS device compute follows. +not change ordinary feature algorithm placement; explicit device consumers follow. FPFH device computation now consumes this radius-row representation through `GpuFpfhInput` and neighborhoods tied to that exact uploaded input. Explicit @@ -52,6 +52,23 @@ retaining CPU tree construction, ordered flag compaction and attribute gathering GPU f32 selection can differ near thresholds/ties from CPU f64; ordinary CPU behavior remains unchanged. Eigen diagnostics require separate explicit readback. +MLS uses the same input-bound radius representation through `GpuMlsInput`. +The uploaded query positions are shared with frame computation. Three device +passes in one submission validate complete neighborhoods, compute compensated +weighted frames/Jacobi eigenvectors, fit order-0/1/2 surfaces and finalize global +validity. Private leased scratch stores 64-byte frames and 192-byte fit systems +per point; owned original-order position/mode buffers survive input drops. +`filter-mls-gpu` keeps CPU tree construction and schema/attribute reconstruction +at explicit upload/download boundaries. Minimum support includes self, then +fitting excludes only the query index. Sparse points retain their coordinates; +insufficient polynomial support or rank failure uses a device plane fallback. +Normalized f32 Cholesky uses relative rank threshold 1e-6; CPU f64 Gaussian +elimination uses absolute 1e-12. Nearly singular fits can differ substantially; +per-point modes and counts expose outcomes, without a CPU equivalence promise. +Malformed/incomplete rows or evaluated local-extent underflow invalidate and +zero all output; consumers must inspect the first five status words. See +[MLS evidence](../notes/2026-10-13_gpu_mls_smoothing.md). + | Area | Decision | | --- | --- | | Core model | `SpatialTensor` + `SpatialIndex` + `SpatialAlgorithm` + `SpatialRuntime` | diff --git a/docs/FEATURE_MATRIX.md b/docs/FEATURE_MATRIX.md index 4278303c..ec46f653 100644 --- a/docs/FEATURE_MATRIX.md +++ b/docs/FEATURE_MATRIX.md @@ -6,6 +6,21 @@ should normally depend on `spatialrust` and select only the profiles they need. ## User-facing profiles +`filter-mls-gpu` exposes `MlsSmoothing::{smoothed_positions_with_gpu,filter_with_gpu}` +with a supplied runtime and self-inclusive neighborhood capacity. Radius search, +weighted local frames and constant/linear/quadratic fits run on GPU; CPU tree +construction and schema/attribute reconstruction remain. `GpuMlsInput` under +`gpu-wgpu` retains input-bound neighborhoods and owned position/mode results for +device consumers. H2D payload is 32 bytes per point plus 32 per tree node; +D2H is 16 bytes per point plus 32 status bytes, both zero when empty. No graph +readback or automatic CPU fallback occurs. Finite XYZ/extents and positive finite +nonzero radius squares are required; evaluated nonzero local extents with +subnormal/zero squares reject. Invalid rows/numerics zero all positions/modes. +Downloaded outcomes distinguish unchanged points, plane fallback and requested +polynomial fits. GPU f32 relative Cholesky rank tests differ from CPU f64 +absolute Gaussian pivots; nearly singular fits can have large coordinate +differences. See [MLS precision and measurements](../notes/2026-10-13_gpu_mls_smoothing.md). + `feature-iss-gpu` exposes `IssKeypointDetector::detect_with_gpu` with a supplied runtime and salient/non-max capacities. Both graphs, compensated local covariance, eigenvalue saliency and local maxima run on GPU; CPU tree @@ -44,7 +59,7 @@ CPU traversal order, counts all matches, and supports capacity=0 count-only search. `download` accounts for padded readback bytes; complete row access rejects overflow while `truncated_rows` explicitly returns a traversal prefix. Self/duplicate points are included. GPU feature consumers must check counts -and preserve query/target semantics. GPU FPFH/ISS consume radius rows as described +and preserve query/target semantics. GPU FPFH/ISS/MLS consume radius rows as described above; other feature algorithms retain their separate CPU APIs. Native correctness and completed/ resident measurements are in [the neighborhood note](../notes/2026-10-12_gpu_bounded_neighborhoods.md). @@ -106,6 +121,7 @@ before allocation; no dependency, public signature or image readback is added. | `register-gicp-gpu` | GICP, GPU radius covariance and explicit GPU correspondences | Caller-owned runtime and full payload receipt on explicit alignment; CPU 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 | | `pipeline-mvp-gpu` | `mvp` GPU stages plus `gpu-wgpu` | GPU-enabled MVP pipeline | | `gpu-aoso-staging` | AoSoA core packing plus GPU-resident frame APIs | Chained GPU execution with explicit readback | | `serde` | serde derives for core/math metadata and schemas | Configuration and metadata serialization | diff --git a/docs/ROADMAP.md b/docs/ROADMAP.md index 2726e2af..41a72c31 100644 --- a/docs/ROADMAP.md +++ b/docs/ROADMAP.md @@ -1128,8 +1128,8 @@ Dependencies and explicit transfers come before broader algorithm coverage. | --- | --- | --- | --- | | ICP correspondences | Complete | Reusable runtime-owned GPU KD-tree and explicit batched nearest uploads/readback; opt-in ICP update/fitness integration with CPU solver | CPU/brute nearest ties and gates, ownership/reuse/capacity, pose/trace/receipt checks and completed query/registration comparisons; `notes/2026-10-12_gpu_icp_correspondences.md` | | GICP correspondences | Complete | Reuse GPU target for unbounded matching followed by host gating; explicit radius covariance runtime and full receipt with CPU Hessian/pose solver | Pose/stopping/gating/receipt/unequal-cloud checks, error propagation, completed comparisons for CPU/radius-GPU covariance; `notes/2026-10-12_gpu_gicp_correspondences.md` | -| Feature neighborhoods | Resident API complete; FPFH/ISS consumers added | Reused target/query buffers, GPU k-NN and bounded radius rows with exact counts/overflow, explicit readback and raw storage for device consumers | CPU/brute membership/order, degenerate/tied/partial/ownership/capacity cases, device consumer test and completed/resident measurements; `notes/2026-10-12_gpu_bounded_neighborhoods.md` | +| Feature neighborhoods | Resident API complete; FPFH/ISS/MLS consumers added | Reused target/query buffers, GPU k-NN and bounded radius rows with exact counts/overflow, explicit readback and raw storage for device consumers | CPU/brute membership/order, degenerate/tied/partial/ownership/capacity cases, device consumer test and completed/resident measurements; `notes/2026-10-12_gpu_bounded_neighborhoods.md` | | FPFH descriptors | Complete | Input-bound reusable radius rows, GPU SPFH and weighted FPFH, explicit upload/readback with overflow status; CPU matching/RANSAC retained | CPU/analytic descriptor comparisons, input/ownership/overflow/error guards, native/facade/WASM builds and completed comparisons against current parallel CPU; `notes/2026-10-12_gpu_fpfh_descriptors.md` | | ISS keypoints | Complete | Shared original queries/positions, input-bound radius graphs, GPU covariance/eigenvalue saliency and local maxima; explicit flag/status download, CPU compaction/attributes | Analytic eigen/tie checks, f64 math/error-band and independent NMS references, schema/attribute bits, ownership/overflow/underflow guards and completed parallel-CPU comparisons; `notes/2026-10-13_gpu_iss_keypoints.md` | -| MLS compute | Planned | Device smoothing phases on the shared neighborhood representation | CPU/math references, attribute/schema preservation, numeric tolerances where reductions change, transfer-inclusive benefit | +| MLS compute | Complete | Input-bound reusable radius rows, GPU weighted frames and order-0/1/2 fits, explicit position/mode/status download; CPU tree and reconstruction | Analytic quadratic/noise checks, per-mode CPU references with raw requested-order errors retained, schema/attribute bits, ownership/overflow/underflow guards and completed parallel-CPU comparisons; `notes/2026-10-13_gpu_mls_smoothing.md` | | Registration updates | Planned | Explicit correspondence reduction/pose-update boundaries after GPU search is established | Pose/convergence/fitness references, numerical conditioning, source/target residency and end-to-end registration timing | diff --git a/notes/2026-10-13_gpu_mls_smoothing.md b/notes/2026-10-13_gpu_mls_smoothing.md new file mode 100644 index 00000000..84f16de4 --- /dev/null +++ b/notes/2026-10-13_gpu_mls_smoothing.md @@ -0,0 +1,221 @@ +# Resident GPU MLS smoothing + +Date: 2026-10-13. Named slice: constant, linear and quadratic Moving Least +Squares surface fits on GPU, following the resident radius-neighborhood API. + +## Execution and transfer contract + +`filter-mls-gpu` exposes +`MlsSmoothing::smoothed_positions_with_gpu(runtime, cloud, neighbor_capacity)` +and `filter_with_gpu` in filtering and the facade. Ordinary smoothing retains +its existing bounded CPU path. The explicit wrapper builds one CPU KD-tree, +uploads the tree and original queries, runs radius search and MLS on GPU, then +downloads positions/modes/status. Schema, metadata and all non-position fields +are reconstructed on CPU using the existing output builder. There is no +intermediate neighborhood/frame/system readback or automatic CPU fallback. + +`GpuMlsInput::upload` under `gpu-wgpu` owns reusable input on one runtime. Query +positions also serve as the frame samples, without another sample upload. +For N original points and M tree nodes, H2D payload is `32*N + 32*M` bytes: +16*N tree positions, 16*N queries and 32*M nodes. CPU tree construction is +external to this resident API. `neighborhoods_resident` submits radius search +without readback or host completion wait. Rows are tied to that exact input; +passing a different input's neighborhoods or a different runtime rejects. + +Capacity and minimum support include self. Every radius match must fit every +row, even rows below minimum support; truncated prefixes reject. After that +support test, fitting excludes only the query index, retaining coincident +points at other indices. Distances use an inclusive radius gate and Gaussian +weights `exp(-distance_squared / max((radius/2)^2, 1e-12))`. + +`smooth_resident` submits three ordered device passes in one submission: + +1. Validate counts, indices and distance gates; prepare weighted centroid, + covariance, smallest-eigenvector normal and orthonormal local frame. +2. Fit and evaluate the requested constant/linear/quadratic height field. +3. Finalize validity, zeroing ALL output if any row or numerical error occurred. + +Frames use query-relative coordinates and actual local extent for scaling, +compensated f32 weight/centroid/covariance sums, and max-pivot Jacobi with +32 iterations and the physical covariance stopping rule of 1e-12. Centroids +are reconstructed as absolute f32 XYZ before covariance/fitting. The constant +fit is a scalar weighted mean height. Linear and quadratic fits use normalized +bases `[1,s,t]` and `[1,s,t,s*s,s*t,t*t]`, direct f32 Gram/RHS accumulation and +in-place Cholesky. A relative Schur-complement threshold of +`1e-6 * largest_absolute_Gram_entry` selects the plane on rank failure. + +Sparse/zero-weight/zero-extent points retain exact original coordinates. +Insufficient polynomial support or failed rank tests produce plane projections +on device. No point is fitted by a CPU fallback. Input XYZ and global extent +squares must be finite; radius must be positive with a finite nonzero f32 +square. An evaluated nonzero local extent with a subnormal/zero f32 square +sets global numerical-error status. These checks are stricter than ordinary +CPU MLS. Huge usize support requirements map to u32::MAX. + +Per-runtime lazy pipelines share six storage bindings and one 32-byte uniform, +with 64 invocations per workgroup. Private scratch leases contain 64*N frame +bytes and 192*N system bytes; both are allocated for every nonempty order. +Consumed scratch rows are initialized before use. Leases remain exclusive +through submission, then return to the existing byte-bounded pool. Result +buffers are owned separately and never pooled. All sizes, u32 indexing, host +bounds, shader bindings and dispatch limits are checked before allocation. + +`GpuMlsPositions` owns original-order vec4 rows `(XYZ,mode)`, where mode +0 means unchanged, 1 plane fallback and 2 the requested polynomial. Empty +physical position storage has 16 unused bytes. Eight u32 status words contain +bad rows, required capacity, malformed rows, last bad row plus one, numerical +errors, unchanged count, plane count and polynomial count. The FIRST FIVE +words must all be zero. The last three are outcome diagnostics, not errors; +their sum must equal N. All XYZ/mode words are zeroed on global errors. Custom +GPU consumers must inspect status and treat exposed storage as immutable. + +`download` explicitly copies position/mode/status in one staging submission, +waits, validates status before decoding and verifies finite XYZ, valid modes +and matching counts. D2H payload is `16*N + 32` bytes. Empty H2D/D2H transfers +are both zero. Uniform controls and private scratch are excluded from payload +receipts. Native Vulkan execution is verified; WASM evidence is compilation +only, without a browser synchronous-download claim. No core dependency or +unsafe code is added. + +## Precision and correctness evidence + +GPU MLS is opt-in and DOES NOT promise CPU coordinate or fallback parity. +CPU MLS uses unnormalized f64 Gram/RHS sums and Gaussian elimination with an +absolute pivot threshold of 1e-12; the GPU route uses normalized f32 Cholesky +with the relative threshold above. Eigenvectors, near-rank fits and fallback +decisions can differ. Per-point modes/counts expose the actual GPU outcome. +Use the existing CPU API when its exact decisions are required. + +Four explicitly run adapter integration tests cover 0/1/2/17/63/64/65/257/4097 +points, five noisy/curved/coincident/planar/skinny fixture kinds and orders +0/1/2: 135 conditions and 205,470 coordinates. Each coordinate is checked +against independent existing CPU math for its REPORTED outcome: exact original +XYZ for unchanged, CPU order-0 weighted plane for plane fallback, and the CPU +requested order for polynomial fits. The selected-fixture coordinate bound is +`5e-5 + search_radius*1e-3 + 8*ulp(reference_coordinate)`; maximum observed error +against that mode reference is 0.0001719594. This is fixture evidence, not a +general numerical guarantee or a bound on requested-order CPU differences. + +The raw difference from the requested CPU order is also retained. Its maximum +is **0.97287786** on the 65-point curved/noisy fixture with quadratic order: +GPU chose 52 plane fallbacks and 13 polynomial fits. This adverse, nearly +singular case is not equivalent to CPU quadratic smoothing; the per-mode +bound does not erase that difference. For the 4,097-point well-supported +noisy, curved and coincident quadratic fixtures, all 4,097 points fit +polynomials and the raw CPU difference is at most 0.0000004768. + +An independent symmetric paraboloid with a 0.02 center spike distinguishes +linear and quadratic execution: linear center height exceeds 0.002 while +quadratic center height is below 1e-6. A deterministic noisy plane's height +RMS decreases from 0.01141399 to 0.0021848176 (about 81%). These tests establish +actual polynomial fitting and smoothing beyond the mode-conditioned comparisons. + +Other checks cover custom position names, exact schema/metadata preservation, +all scalar attribute types including NaN payloads/signed zero, exact unchanged +XYZ bits, full transfer receipts, invalid config/coordinates/extents, +self-inclusive overflow, poisoned indices, numerical underflow, translated +planes, empty inputs, and fresh output after invalid work. Four concurrent +callers reuse scratch while keeping earlier output alive; results remain +usable after input/tree/neighborhood drops and reject mismatched ownership. + +Local validation passed: GPU all-features tests 207 passed / 23 ignored, +filtering all-features library tests 40 passed / 1 ignored, four explicitly run +adapter integration tests, and one explicitly run facade GPU test. Clippy +passed for GPU/filtering all features and the CPU-only MLS configuration. +GPU no-default-features, filtering MLS-GPU wasm32-unknown-unknown and workspace +Rustfmt checks passed. Removing the added GPU result types/methods leaves CPU +MLS source exactly equal to baseline commit +`053e84fac38831cf33f43a0847e1a34d61b6cbd9`. + +## Completed-call measurements + +Release builds on a native Vulkan discrete-device backend. Three warmups per +stage precede 20 paired rounds with cyclic CPU/full-GPU/resident stage order. +Two surveys run separately with no other local build, test or statistical +analysis during timing. All samples are retained, including slow observations. +Runtime/pipeline creation is warmed outside timers. + +CPU timings use current public `MlsSmoothing::filter`, including its bounded +worker policy: four workers at 65,537 points and one at smaller selected sizes. +Full GPU timings use `filter_with_gpu`, including validation, CPU tree build, +tree/query upload, radius search, frames/fits/finalization, position/mode/status +download, completion wait and schema/attribute reconstruction. Result +verification/drop is outside timers; internal preparation/staging destruction +is included in the completed call. + +Resident timings include owned result allocation, private scratch acquisition, +all three MLS passes, submission and `wait_idle`. They exclude CPU preparation, +upload, neighborhood generation, download and cloud reconstruction. These are +prepared fitting-phase measurements, not GPU timestamps or full-call speedups. + +P is a noisy planar 128-column grid with spacing 0.05 and order 1; C adds +curved height and deterministic noise with order 2; D duplicates consecutive C +points and uses order 2. Search radius is 0.16 throughout, minimum support 6, +capacity 64 for P/C and 96 for D. Before timing, each GPU output is checked +against its mode reference and its raw requested-order CPU difference is +recorded. Repeated CPU and GPU calls match their respective references exactly. +Ratios on rows with different CPU/GPU fitting decisions do not imply +accuracy-equivalent output. + +Times below are medians in milliseconds. Ratios are geometric means of 20 +paired CPU/full-GPU timings. Percentile 95% bootstrap intervals use 10,000 +resamples of paired log ratios, seeds 20261013 plus zero-based record index, +and sorted entries 249/9749. These selected fixtures do not establish a +universal deployment crossover; condition geometry and scheduling matter. + +| Survey | Points | Fixture/order | Workers | Unchanged/plane/poly | Mode error | Requested CPU error | CPU ms | Full GPU ms | Resident ms | Paired full-call ratio (95% interval) | +| --- | ---: | :---: | ---: | ---: | ---: | ---: | ---: | ---: | ---: | --- | +| 1 | 257 | P/1 | 1 | 0/0/257 | 2.98e-08 | 2.98e-08 | 1.757 | 34.259 | 6.090 | 0.09x (0.06-0.14) | +| 1 | 257 | C/2 | 1 | 0/138/119 | 0.000172 | 0.00676 | 3.034 | 5.726 | 2.249 | 0.44x (0.35-0.55) | +| 1 | 257 | D/2 | 1 | 0/35/222 | 4.42e-07 | 4.42e-07 | 2.356 | 4.909 | 1.986 | 0.52x (0.46-0.60) | +| 1 | 4,097 | P/1 | 1 | 0/0/4097 | 1.19e-07 | 1.19e-07 | 77.003 | 15.063 | 4.116 | 5.31x (4.84-5.77) | +| 1 | 4,097 | C/2 | 1 | 0/0/4097 | 4.77e-07 | 4.77e-07 | 94.565 | 25.771 | 14.504 | 3.84x (3.51-4.18) | +| 1 | 4,097 | D/2 | 1 | 0/0/4097 | 2.38e-07 | 2.38e-07 | 80.714 | 15.875 | 11.662 | 5.13x (4.81-5.49) | +| 1 | 16,384 | P/1 | 1 | 0/0/16384 | 5.96e-08 | 5.96e-08 | 304.223 | 24.767 | 9.304 | 9.43x (7.68-11.22) | +| 1 | 16,384 | C/2 | 1 | 0/0/16384 | 4.77e-07 | 4.77e-07 | 279.677 | 33.858 | 25.175 | 7.78x (6.94-8.68) | +| 1 | 16,384 | D/2 | 1 | 0/0/16384 | 4.77e-07 | 4.77e-07 | 605.405 | 68.258 | 52.525 | 8.57x (7.91-9.20) | +| 1 | 65,537 | P/1 | 4 | 0/0/65537 | 4.77e-07 | 4.77e-07 | 209.168 | 58.353 | 38.299 | 3.45x (3.13-3.71) | +| 1 | 65,537 | C/2 | 4 | 0/1/65536 | 1.91e-06 | 1.91e-06 | 210.302 | 111.005 | 93.274 | 1.92x (1.83-2.02) | +| 1 | 65,537 | D/2 | 4 | 0/0/65537 | 9.54e-07 | 9.54e-07 | 371.191 | 226.562 | 204.646 | 1.61x (1.48-1.76) | +| 2 | 257 | P/1 | 1 | 0/0/257 | 2.98e-08 | 2.98e-08 | 0.756 | 1.435 | 0.259 | 0.56x (0.53-0.60) | +| 2 | 257 | C/2 | 1 | 0/138/119 | 0.000172 | 0.00676 | 0.907 | 1.231 | 0.288 | 0.72x (0.70-0.75) | +| 2 | 257 | D/2 | 1 | 0/35/222 | 4.42e-07 | 4.42e-07 | 1.066 | 1.626 | 0.352 | 0.67x (0.65-0.68) | +| 2 | 4,097 | P/1 | 1 | 0/0/4097 | 1.19e-07 | 1.19e-07 | 27.877 | 3.067 | 0.720 | 9.08x (8.48-9.81) | +| 2 | 4,097 | C/2 | 1 | 0/0/4097 | 4.77e-07 | 4.77e-07 | 38.621 | 3.780 | 1.559 | 10.05x (8.90-11.36) | +| 2 | 4,097 | D/2 | 1 | 0/0/4097 | 2.38e-07 | 2.38e-07 | 50.834 | 11.088 | 7.999 | 4.58x (4.21-5.01) | +| 2 | 16,384 | P/1 | 1 | 0/0/16384 | 5.96e-08 | 5.96e-08 | 301.913 | 44.647 | 10.117 | 7.79x (6.26-9.68) | +| 2 | 16,384 | C/2 | 1 | 0/0/16384 | 4.77e-07 | 4.77e-07 | 387.395 | 42.101 | 26.952 | 8.45x (7.24-9.38) | +| 2 | 16,384 | D/2 | 1 | 0/0/16384 | 4.77e-07 | 4.77e-07 | 294.013 | 56.855 | 50.388 | 5.01x (4.74-5.29) | +| 2 | 65,537 | P/1 | 4 | 0/0/65537 | 4.77e-07 | 4.77e-07 | 203.909 | 59.630 | 38.248 | 3.43x (3.26-3.59) | +| 2 | 65,537 | C/2 | 4 | 0/1/65536 | 1.91e-06 | 1.91e-06 | 224.250 | 112.656 | 93.240 | 2.00x (1.92-2.09) | +| 2 | 65,537 | D/2 | 4 | 0/0/65537 | 9.54e-07 | 9.54e-07 | 429.181 | 226.458 | 204.603 | 1.94x (1.73-2.21) | + +All 257-point full GPU calls are slower in both surveys (paired ratios +0.09-0.72x). At 4,097 and 16,384 points, selected ratios are 3.84-10.05x +and 5.01-9.43x respectively. Against four CPU workers at 65,537 points, +ratios are 1.61-3.45x, with raw requested-order coordinate differences no +larger than 0.00000190735 and at most one GPU plane fallback. These measured +gains support an explicit GPU option on these fixtures, without changing +placement defaults. Substantial run-to-run timing variation is retained in +the table, especially for small calls and the first survey; no samples or +conditions are excluded. The 257-point C/2 row retains its 0.006760955 raw +CPU difference and 138 plane fallbacks, independently of timing. + +Reproduction from `C:/Users/rsasa/Workspace/SpatialRust`: + +```powershell +cargo test -p spatialrust-gpu --all-features +cargo test -p spatialrust-filtering --all-features --test gpu_mls -- --include-ignored --nocapture +cargo test -p spatialrust --no-default-features --features filter-mls-gpu --test gpu_mls_smoke -- --include-ignored +cargo test -p spatialrust-filtering --all-features --release --test gpu_mls_benchmark completed_mls_calls -- --include-ignored --nocapture --exact --test-threads=1 +``` + +The benchmark's exact filter is required: its imported reference module also +contains adapter integration tests. Benchmark runs are isolated from local +builds, tests and statistics, and intentionally ignored in normal CI. Explicit +adapter tests fail if no native adapter exists instead of silently skipping. + +Local validation logs, raw timing surveys and measurement JSON are under +`C:/Users/rsasa/Workspace/SpatialRust/target/automation/gpu-mls-*.log` and +`C:/Users/rsasa/Workspace/SpatialRust/target/automation/gpu-mls-measurements.json`. +These ignored local artifacts do not name the local adapter.