From f087794f358c2800a59a0c3576dd38ed6b598483 Mon Sep 17 00:00:00 2001 From: rsasaki0109 Date: Mon, 12 Oct 2026 23:33:43 +0900 Subject: [PATCH] feat(gpu): compute resident FPFH descriptors with explicit transfers --- .github/workflows/ci.yml | 4 + CHANGELOG.md | 7 + README.md | 12 +- crates/spatialrust-gpu/src/kernels/fpfh.rs | 433 ++++++++++++++++++ .../spatialrust-gpu/src/kernels/fpfh_tests.rs | 119 +++++ crates/spatialrust-gpu/src/kernels/mod.rs | 2 + crates/spatialrust-gpu/src/lib.rs | 3 +- crates/spatialrust-gpu/src/runtime.rs | 2 + crates/spatialrust-gpu/src/shaders/fpfh.wgsl | 108 +++++ crates/spatialrust-registration/Cargo.toml | 1 + crates/spatialrust-registration/src/fpfh.rs | 37 ++ crates/spatialrust-registration/src/lib.rs | 3 + .../tests/gpu_fpfh.rs | 216 +++++++++ .../tests/gpu_fpfh_benchmark.rs | 72 +++ crates/spatialrust/Cargo.toml | 1 + crates/spatialrust/tests/gpu_fpfh_smoke.rs | 36 ++ docs/ARCHITECTURE.md | 12 +- docs/FEATURE_MATRIX.md | 19 +- docs/ROADMAP.md | 5 +- notes/2026-10-12_gpu_fpfh_descriptors.md | 167 +++++++ 20 files changed, 1251 insertions(+), 8 deletions(-) create mode 100644 crates/spatialrust-gpu/src/kernels/fpfh.rs create mode 100644 crates/spatialrust-gpu/src/kernels/fpfh_tests.rs create mode 100644 crates/spatialrust-gpu/src/shaders/fpfh.wgsl create mode 100644 crates/spatialrust-registration/tests/gpu_fpfh.rs create mode 100644 crates/spatialrust-registration/tests/gpu_fpfh_benchmark.rs create mode 100644 crates/spatialrust/tests/gpu_fpfh_smoke.rs create mode 100644 notes/2026-10-12_gpu_fpfh_descriptors.md diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index f9084ad7..a808f6bf 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -182,6 +182,10 @@ jobs: args: -p spatialrust-registration --features register-fpfh - name: spatialrust-registration-fpfh-meta args: -p spatialrust --features register-fpfh + - name: spatialrust-registration-fpfh-gpu + args: -p spatialrust-registration --features register-fpfh-gpu + - name: spatialrust-registration-fpfh-gpu-meta + args: -p spatialrust --features register-fpfh-gpu - name: spatialrust-metrics args: -p spatialrust-metrics --features metrics-distance - name: spatialrust-metrics-meta diff --git a/CHANGELOG.md b/CHANGELOG.md index 05a2b930..c5ce8d39 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,13 @@ removed no sooner than the next major (see `docs/API_STABILITY.md`). ### Changed +- Add opt-in `register-fpfh-gpu` descriptor computation: reusable radius + neighborhoods feed GPU SPFH and distance-weighted FPFH passes, with explicit + sample/tree/query upload and descriptor/status readback receipts. Input and + neighborhood ownership are checked; incomplete/malformed rows invalidate + the entire device result and return errors at download. CPU FPFH/RANSAC + defaults remain unchanged; GPU angular/bin rounding is documented. + - Add explicitly uploaded reusable GPU query batches and resident KD-tree nearest-k/radius neighborhoods under `gpu-wgpu`. k=0..=64 retains CPU selection and stable distance order; radius rows retain traversal order and diff --git a/README.md b/README.md index e4314cac..ca4f7389 100644 --- a/README.md +++ b/README.md @@ -653,8 +653,16 @@ GPU k-NN/radius search can retain its results for later device passes with `GpuQueryBatch::upload` and `GpuKdTree::{nearest_k_resident,radius_neighbors_resident}` (`gpu-wgpu`). Radius capacity overflow is counted explicitly and complete readback rejects it. -This provides the neighborhood foundation for feature compute; FPFH/ISS/MLS -remain CPU algorithms. See [contracts and measurements](notes/2026-10-12_gpu_bounded_neighborhoods.md). +This provides the neighborhood foundation for feature compute. See +[contracts and measurements](notes/2026-10-12_gpu_bounded_neighborhoods.md). + +FPFH descriptors have an explicit GPU path under `register-fpfh-gpu`: +`fpfh_descriptors_gpu(&runtime, &cloud, radius, neighbor_capacity)` runs radius +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; ISS/MLS remain CPU. +See [FPFH measurements](notes/2026-10-12_gpu_fpfh_descriptors.md). ### Library diff --git a/crates/spatialrust-gpu/src/kernels/fpfh.rs b/crates/spatialrust-gpu/src/kernels/fpfh.rs new file mode 100644 index 00000000..1b18c970 --- /dev/null +++ b/crates/spatialrust-gpu/src/kernels/fpfh.rs @@ -0,0 +1,433 @@ +use std::sync::Arc; + +use bytemuck::{Pod, Zeroable}; +use spatialrust_core::{SpatialError, SpatialResult, TransferDirection, TransferStats}; +use spatialrust_search::KdTree; +use wgpu::util::DeviceExt; + +use crate::{scratch::ScratchBuffers, GpuKdTree, GpuNeighborhoods, GpuQueryBatch, WgpuRuntime}; + +/// Three eleven-bin angular histograms per FPFH descriptor. +pub const GPU_FPFH_DESCRIPTOR_LEN: usize = 33; + +#[repr(C)] +#[derive(Clone, Copy, Pod, Zeroable)] +struct Sample { + position: [f32; 4], + normal: [f32; 4], +} + +#[repr(C)] +#[derive(Clone, Copy, Pod, Zeroable)] +struct Params { + count: u32, + capacity: u32, + padding: [u32; 2], +} + +fn invalid() -> SpatialError { + SpatialError::InvalidArgument("GPU FPFH requires matching input/runtime ownership, finite XYZ/unit normals, a positive finite radius square and supported buffers/dispatch".to_owned()) +} + +fn plan(points: u64, limits: &wgpu::Limits) -> SpatialResult<(u32, u64)> { + let bytes = points.checked_mul(132).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 + || 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 < 5 + || limits.max_bindings_per_bind_group < 6 + || limits.max_bind_groups == 0 + || limits.max_uniform_buffers_per_shader_stage == 0 + || limits.max_uniform_buffer_binding_size < 16 + || bound < 32 + || bytes + .checked_add(16) + .map_or(true, |v| v > limits.max_buffer_size || v > isize::MAX as u64) + { + return Err(invalid()); + } + Ok((points as u32, bytes)) +} + +/// Explicitly uploaded FPFH positions/unit normals and their GPU KD-tree/query +/// batch. Original point order is retained; CPU tree construction is external. +/// All subsequent neighborhood and descriptor computation stays on GPU. +pub struct GpuFpfhInput { + samples: Option, + tree: GpuKdTree, + queries: GpuQueryBatch, + points: u32, + identity: Arc<()>, + input_identity: Arc<()>, + transfers: TransferStats, +} + +impl GpuFpfhInput { + /// Uploads a CPU KD-tree and matching original-order finite unit normals. + /// Payload is 64 bytes per point plus 32 per tree node (samples, original + /// queries and KD-tree points/nodes). No readback is performed. Coordinate + /// extents must have a finite f32 squared length, as in CPU FPFH. + pub fn upload( + runtime: &WgpuRuntime, + tree: &KdTree, + normals: &[[f32; 3]], + ) -> SpatialResult { + let layout = tree.layout(); + let (x, y, z) = layout.positions(); + let (points, _) = plan(x.len() as u64, &runtime.device().limits())?; + if normals.len() != x.len() { + return Err(SpatialError::BufferLengthMismatch { + expected: x.len(), + found: normals.len(), + }); + } + if x.iter().chain(y).chain(z).any(|v| !v.is_finite()) + || normals.iter().any(|n| { + let norm = n[0] * n[0] + n[1] * n[1] + n[2] * n[2]; + n.iter().any(|v| !v.is_finite()) || !norm.is_finite() || (norm - 1.0).abs() > 1e-3 + }) + { + return Err(invalid()); + } + if points > 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(), + ); + let samples = if points == 0 { + None + } else { + let samples: Vec<_> = positions + .iter() + .zip(normals) + .map(|(p, n)| Sample { + position: [p[0], p[1], p[2], 0.0], + normal: [n[0], n[1], n[2], 0.0], + }) + .collect(); + transfers.record(TransferDirection::HostToDevice, u64::from(points) * 32); + Some(runtime.device().create_buffer_init(&wgpu::util::BufferInitDescriptor { + label: Some("fpfh-samples"), + contents: bytemuck::cast_slice(&samples), + usage: wgpu::BufferUsages::STORAGE, + })) + }; + Ok(Self { + samples, + tree, + queries, + points, + identity: runtime.profiling_identity(), + input_identity: Arc::new(()), + transfers, + }) + } + + /// Original point/descriptor count. + #[must_use] + pub const fn point_count(&self) -> u32 { + self.points + } + + /// Payload bytes explicitly uploaded at construction. + #[must_use] + pub const fn upload_transfers(&self) -> TransferStats { + self.transfers + } + + /// Builds reusable GPU radius neighborhoods for this input, including self. + /// No readback or host wait occurs. All matches are counted; capacity + /// overflow is checked on GPU during descriptor computation and at download. + pub fn radius_neighborhoods_resident( + &self, + runtime: &WgpuRuntime, + radius: f32, + capacity: u32, + ) -> SpatialResult { + if !Arc::ptr_eq(&self.identity, &runtime.profiling_identity()) + || !radius.is_finite() + || radius <= 0.0 + || !(radius * radius).is_finite() + || radius * radius == 0.0 + { + return Err(invalid()); + } + let neighborhoods = + self.tree.radius_neighbors_resident(runtime, &self.queries, radius, capacity)?; + Ok(GpuFpfhNeighborhoods { + neighborhoods, + input_identity: Arc::clone(&self.input_identity), + radius, + }) + } + + /// Computes GPU SPFH and distance-weighted FPFH using resident neighborhoods + /// from this exact input. Both passes share one submission; results remain + /// on device. Overflow/malformed rows poison the entire output to zero and + /// set status; download rejects that output, and GPU consumers must check + /// status before using it. No automatic CPU fallback or payload copy occurs. + /// GPU angular/bin rounding may differ from CPU near boundaries. + pub fn descriptors_resident( + &self, + runtime: &WgpuRuntime, + neighbors: &GpuFpfhNeighborhoods, + ) -> SpatialResult { + if !Arc::ptr_eq(&self.identity, &runtime.profiling_identity()) + || !Arc::ptr_eq(&self.input_identity, &neighbors.input_identity) + { + return Err(invalid()); + } + let (_, bytes) = plan(u64::from(self.points), &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 output = GpuFpfhDescriptors { + buffer: create("fpfh-descriptors", bytes.max(4)), + status: create("fpfh-status", 16), + points: self.points, + bytes, + identity: Arc::clone(&self.identity), + }; + if self.points == 0 { + return Ok(output); + } + let cache = runtime.fpfh_pipelines.get_or_init(|| FpfhPipelines::new(runtime.device())); + let mut scratch = ScratchBuffers::new(runtime); + let spfh = scratch.acquire("fpfh-spfh", bytes); + let params = Params { + count: self.points, + capacity: neighbors.neighborhoods.capacity(), + padding: [0; 2], + }; + let uniform = runtime.device().create_buffer_init(&wgpu::util::BufferInitDescriptor { + label: Some("fpfh-params"), + contents: bytemuck::bytes_of(¶ms), + usage: wgpu::BufferUsages::UNIFORM, + }); + let buffers = [ + &uniform, + self.samples.as_ref().unwrap(), + neighbors.neighborhoods.buffer(), + &spfh, + &output.buffer, + &output.status, + ]; + let entries: Vec<_> = buffers + .iter() + .enumerate() + .map(|(i, b)| wgpu::BindGroupEntry { + binding: i as u32, + resource: b.as_entire_binding(), + }) + .collect(); + let group = runtime.device().create_bind_group(&wgpu::BindGroupDescriptor { + label: Some("fpfh-bindings"), + layout: &cache.layout, + entries: &entries, + }); + let mut encoder = runtime.device().create_command_encoder(&Default::default()); + for pipeline in [&cache.spfh, &cache.fpfh] { + let mut pass = encoder.begin_compute_pass(&wgpu::ComputePassDescriptor { + label: Some("fpfh-compute"), + timestamp_writes: None, + }); + pass.set_pipeline(pipeline); + pass.set_bind_group(0, &group, &[]); + pass.dispatch_workgroups(self.points.div_ceil(64), 1, 1); + } + scratch.submit_and_recycle(encoder); + Ok(output) + } +} + +/// Radius rows tied to one immutable FPFH input; usable by both device phases. +pub struct GpuFpfhNeighborhoods { + neighborhoods: GpuNeighborhoods, + input_identity: Arc<()>, + radius: f32, +} + +impl GpuFpfhNeighborhoods { + /// The underlying radius rows/counts. Treat storage as immutable. + #[must_use] + pub fn neighborhoods(&self) -> &GpuNeighborhoods { + &self.neighborhoods + } + /// Radius used for this neighborhood graph. + #[must_use] + pub const fn radius(&self) -> f32 { + self.radius + } +} + +/// Owned resident FPFH descriptors with GPU validity status. Nonempty outputs +/// have 33 f32 values per original point. Custom GPU consumers must reject +/// status.bad_rows>0 before accepting any descriptor; failed outputs are zero. +pub struct GpuFpfhDescriptors { + buffer: wgpu::Buffer, + status: wgpu::Buffer, + points: u32, + bytes: u64, + identity: Arc<()>, +} + +impl GpuFpfhDescriptors { + /// Flat original-order 33-f32 descriptor rows; empty results have 4 unused bytes. + #[must_use] + pub fn buffer(&self) -> &wgpu::Buffer { + &self.buffer + } + /// Four u32 status words: bad rows, last bad row plus one, maximum required + /// radius capacity, malformed rows. All zero means valid descriptors. + #[must_use] + pub fn status_buffer(&self) -> &wgpu::Buffer { + &self.status + } + /// Number of descriptor rows. + #[must_use] + pub const fn point_count(&self) -> u32 { + self.points + } + /// Explicit descriptor/status readback bytes (zero for empty input). + #[must_use] + pub const fn readback_bytes(&self) -> u64 { + if self.points == 0 { + 0 + } else { + self.bytes + 16 + } + } + /// Explicitly downloads descriptors and GPU status in one submission. Rejects + /// incomplete/malformed neighborhoods before decoding descriptors. Transfers + /// include 132 bytes per point plus 16 status bytes; no neighborhood/SPFH + /// readback occurs. This waits for completion and propagates GPU errors. + pub fn download( + &self, + runtime: &WgpuRuntime, + ) -> SpatialResult<(Vec<[f32; 33]>, TransferStats)> { + if !Arc::ptr_eq(&self.identity, &runtime.profiling_identity()) { + return Err(invalid()); + } + if self.points == 0 { + return Ok((Vec::new(), TransferStats::default())); + } + let staging = runtime.device().create_buffer(&wgpu::BufferDescriptor { + label: Some("fpfh-readback"), + size: self.bytes + 16, + 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, self.bytes); + encoder.copy_buffer_to_buffer(&self.status, 0, &staging, self.bytes, 16); + 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 status: &[u32] = bytemuck::cast_slice(&mapped[self.bytes as usize..]); + if status[0] > 0 { + let error=SpatialError::InvalidArgument(format!("GPU FPFH rejected {} incomplete/malformed neighborhood rows; last row {}, required capacity {}, malformed rows {}",status[0],status[1].saturating_sub(1),status[2],status[3])); + drop(mapped); + staging.unmap(); + return Err(error); + } + let descriptors = bytemuck::cast_slice::(&mapped[..self.bytes as usize]) + .chunks_exact(GPU_FPFH_DESCRIPTOR_LEN) + .map(|row| row.try_into().expect("FPFH rows contain exactly 33 values")) + .collect(); + drop(mapped); + staging.unmap(); + let mut transfers = TransferStats::default(); + transfers.record(TransferDirection::DeviceToHost, self.bytes + 16); + Ok((descriptors, transfers)) + } +} + +pub(crate) struct FpfhPipelines { + layout: wgpu::BindGroupLayout, + spfh: wgpu::ComputePipeline, + fpfh: wgpu::ComputePipeline, +} +impl FpfhPipelines { + fn new(device: &wgpu::Device) -> Self { + let entries: Vec<_> = (0..6) + .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("fpfh-layout"), + entries: &entries, + }); + let pipeline_layout = device.create_pipeline_layout(&wgpu::PipelineLayoutDescriptor { + label: Some("fpfh-layout"), + bind_group_layouts: &[&layout], + push_constant_ranges: &[], + }); + let shader = device.create_shader_module(wgpu::ShaderModuleDescriptor { + label: Some("fpfh-shader"), + source: wgpu::ShaderSource::Wgsl(include_str!("../shaders/fpfh.wgsl").into()), + }); + let create = |entry| { + device.create_compute_pipeline(&wgpu::ComputePipelineDescriptor { + label: Some("fpfh-pipeline"), + layout: Some(&pipeline_layout), + module: &shader, + entry_point: Some(entry), + compilation_options: Default::default(), + cache: None, + }) + }; + Self { spfh: create("build_spfh"), fpfh: create("build_fpfh"), layout } + } +} + +#[cfg(all(test, not(target_arch = "wasm32")))] +#[path = "fpfh_tests.rs"] +mod tests; diff --git a/crates/spatialrust-gpu/src/kernels/fpfh_tests.rs b/crates/spatialrust-gpu/src/kernels/fpfh_tests.rs new file mode 100644 index 00000000..f34f07cf --- /dev/null +++ b/crates/spatialrust-gpu/src/kernels/fpfh_tests.rs @@ -0,0 +1,119 @@ +use super::*; + +fn read_words(runtime: &WgpuRuntime, buffer: &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(buffer, 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 +} + +#[test] +fn fpfh_plan_rejects_buffer_dispatch_and_binding_boundaries() { + let limits = wgpu::Limits::default(); + assert!(plan(u64::MAX, &limits).is_err()); + assert!(plan(u64::from(u32::MAX) + 1, &limits).is_err()); + let max = u64::from(limits.max_storage_buffer_binding_size) / 132; + assert!(plan(max, &limits).is_ok()); + assert!(plan(max + 1, &limits).is_err()); + let mut restricted = limits.clone(); + restricted.max_compute_workgroups_per_dimension = 1; + assert!(plan(64, &restricted).is_ok()); + assert!(plan(65, &restricted).is_err()); + for field in 0..5 { + restricted = limits.clone(); + match field { + 0 => restricted.max_compute_workgroup_size_x = 63, + 1 => restricted.max_storage_buffers_per_shader_stage = 4, + 2 => restricted.max_bindings_per_bind_group = 5, + 3 => restricted.max_uniform_buffer_binding_size = 15, + _ => restricted.max_buffer_size = 31, + } + assert!(plan(1, &restricted).is_err()); + } +} + +#[test] +fn fpfh_upload_validates_normals_extent_and_empty_receipts() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let tree = KdTree::from_slices(&[0.0], &[0.0], &[0.0]); + for normal in [[0.0; 3], [0.0, 0.0, 2.0], [f32::NAN, 0.0, 1.0], [f32::INFINITY, 0.0, 1.0]] { + assert!(GpuFpfhInput::upload(&runtime, &tree, &[normal]).is_err()); + } + assert!(GpuFpfhInput::upload(&runtime, &tree, &[]).is_err()); + let huge = KdTree::from_slices(&[f32::MAX, -f32::MAX], &[0.0; 2], &[0.0; 2]); + assert!(GpuFpfhInput::upload(&runtime, &huge, &[[0.0, 0.0, 1.0]; 2]).is_err()); + let empty = GpuFpfhInput::upload(&runtime, &KdTree::from_slices(&[], &[], &[]), &[]).unwrap(); + assert_eq!(empty.upload_transfers(), TransferStats::default()); + let neighbors = empty.radius_neighborhoods_resident(&runtime, 1.0, 0).unwrap(); + let output = empty.descriptors_resident(&runtime, &neighbors).unwrap(); + assert_eq!(output.readback_bytes(), 0); + assert_eq!(output.download(&runtime).unwrap().1, TransferStats::default()); +} + +#[test] +fn overflow_and_malformed_device_neighbors_poison_all_descriptors() { + 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 tree = KdTree::from_slices(&x, &zero, &zero); + let input = GpuFpfhInput::upload(&runtime, &tree, &[[0.0, 0.0, 1.0]; 65]).unwrap(); + let short = input.radius_neighborhoods_resident(&runtime, 100.0, 3).unwrap(); + let bad = input.descriptors_resident(&runtime, &short).unwrap(); + assert_eq!(read_words(&runtime, bad.status_buffer(), 16), [65, 65, 65, 0]); + assert!(read_words(&runtime, bad.buffer(), 65 * 132).iter().all(|&v| v == 0)); + assert!(bad.download(&runtime).is_err()); + let full = input.radius_neighborhoods_resident(&runtime, 0.2, 3).unwrap(); + 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: full.neighborhoods().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())); + let poisoned = input.descriptors_resident(&runtime, &full).unwrap(); + assert_eq!(read_words(&runtime, poisoned.status_buffer(), 16), [1, 1, 2, 1]); + assert!(read_words(&runtime, poisoned.buffer(), 65 * 132).iter().all(|&v| v == 0)); + assert!(poisoned.download(&runtime).is_err()); + let fresh = input.radius_neighborhoods_resident(&runtime, 0.2, 3).unwrap(); + let valid = input.descriptors_resident(&runtime, &fresh).unwrap(); + assert_eq!(read_words(&runtime, valid.status_buffer(), 16), [0; 4]); + assert!(valid.download(&runtime).unwrap().0.iter().flatten().any(|&v| v > 0.0)); +} diff --git a/crates/spatialrust-gpu/src/kernels/mod.rs b/crates/spatialrust-gpu/src/kernels/mod.rs index 5d35415a..5e8c3d22 100644 --- a/crates/spatialrust-gpu/src/kernels/mod.rs +++ b/crates/spatialrust-gpu/src/kernels/mod.rs @@ -50,7 +50,9 @@ mod voxel_segments_gpu; mod voxel_gather; pub(crate) mod cluster_union; +pub(crate) mod fpfh; pub(crate) mod nearest; +pub use fpfh::{GpuFpfhDescriptors, GpuFpfhInput, GpuFpfhNeighborhoods, GPU_FPFH_DESCRIPTOR_LEN}; pub use nearest::{ GpuKdTree, GpuNeighborhoodKind, GpuNeighborhoodReadback, GpuNeighborhoods, GpuQueryBatch, MAX_GPU_K_NEIGHBORS, diff --git a/crates/spatialrust-gpu/src/lib.rs b/crates/spatialrust-gpu/src/lib.rs index 10ae6957..737ff219 100644 --- a/crates/spatialrust-gpu/src/lib.rs +++ b/crates/spatialrust-gpu/src/lib.rs @@ -52,7 +52,8 @@ pub use device::{GpuDevice, WgpuDevice}; #[cfg(feature = "gpu-wgpu")] pub use kernels::{ - GpuNeighborhoodKind, GpuNeighborhoodReadback, GpuNeighborhoods, GpuQueryBatch, + GpuFpfhDescriptors, GpuFpfhInput, GpuFpfhNeighborhoods, GpuNeighborhoodKind, + GpuNeighborhoodReadback, GpuNeighborhoods, GpuQueryBatch, GPU_FPFH_DESCRIPTOR_LEN, MAX_GPU_K_NEIGHBORS, }; diff --git a/crates/spatialrust-gpu/src/runtime.rs b/crates/spatialrust-gpu/src/runtime.rs index 4de0a6bc..317c1d51 100644 --- a/crates/spatialrust-gpu/src/runtime.rs +++ b/crates/spatialrust-gpu/src/runtime.rs @@ -25,6 +25,7 @@ pub struct WgpuRuntime { pub(crate) cluster_pipelines: OnceLock, pub(crate) nearest_pipeline: OnceLock, pub(crate) neighborhood_pipelines: OnceLock, + pub(crate) fpfh_pipelines: OnceLock, pub(crate) point_estimator_pipelines: OnceLock, max_gather_channels: u32, @@ -388,6 +389,7 @@ impl WgpuRuntime { cluster_pipelines: OnceLock::new(), nearest_pipeline: OnceLock::new(), neighborhood_pipelines: OnceLock::new(), + fpfh_pipelines: OnceLock::new(), point_estimator_pipelines: OnceLock::new(), max_gather_channels, upload_pool: GpuBufferPool::default(), diff --git a/crates/spatialrust-gpu/src/shaders/fpfh.wgsl b/crates/spatialrust-gpu/src/shaders/fpfh.wgsl new file mode 100644 index 00000000..8b425792 --- /dev/null +++ b/crates/spatialrust-gpu/src/shaders/fpfh.wgsl @@ -0,0 +1,108 @@ +struct Params { count: u32, capacity: u32, pad0: u32, pad1: u32 } +struct Sample { position: vec4, normal: vec4 } +struct Status { + bad_rows: atomic, last_bad_plus_one: atomic, + required_capacity: atomic, invalid_rows: atomic, +} +@group(0) @binding(0) var params: Params; +@group(0) @binding(1) var samples: array; +@group(0) @binding(2) var neighbors: array; +@group(0) @binding(3) var spfh: array; +@group(0) @binding(4) var descriptors: array; +@group(0) @binding(5) var status: Status; + +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 bin(value: f32, lo: f32, hi: f32) -> u32 { + return u32(clamp((value-lo)/(hi-lo),0.0,0.999999)*11.0); +} +fn angle(y: f32, x: f32) -> f32 { + // Explicit signed-zero handling avoids atan2(0,0) backend differences. + if y == 0.0 { + let sign = bitcast(y) & 0x80000000u; + if (bitcast(x) & 0x80000000u) != 0u { + return bitcast(bitcast(3.1415927410125732) | sign); + } + return bitcast(sign); + } + return atan2(y,x); +} +fn reject(i: u32, total: u32, invalid: bool) { + atomicAdd(&status.bad_rows,1u); + atomicMax(&status.last_bad_plus_one,i+1u); + atomicMax(&status.required_capacity,total); + if invalid { atomicAdd(&status.invalid_rows,1u); } +} + +@compute @workgroup_size(64) +fn build_spfh(@builtin(global_invocation_id) id: vec3) { + let i=id.x; + if i>=params.count { return; } + let base=i*(2u+2u*params.capacity); + let stored=neighbors[base]; let total=neighbors[base+1u]; + for (var b=0u;b<33u;b+=1u) { spfh[i*33u+b]=0.0; } + if stored>params.capacity || totalparams.count { + reject(i,total,true);return; + } + if stored!=total { reject(i,total,false);return; } + var histogram: array; + let p=samples[i].position.xyz; let np=samples[i].normal.xyz; + for (var n=0u;n=params.count { reject(i,total,true);return; } + if j==i { continue; } + let diff=samples[j].position.xyz-p; + let distance=sqrt(dot3(diff,diff)); + if distance<1e-9 { continue; } + let pq=diff/distance; + let raw_v=cross3(pq,np); + let v_length=sqrt(dot3(raw_v,raw_v)); + if v_length<1e-9 { continue; } + let v=raw_v/v_length; + let w=cross3(np,v); + let nq=samples[j].normal.xyz; + let alpha=dot3(v,nq);let phi=dot3(np,pq); + let theta=angle(dot3(w,nq),dot3(np,nq)); + histogram[bin(alpha,-1.0,1.0)]+=1.0; + histogram[11u+bin(phi,-1.0,1.0)]+=1.0; + histogram[22u+bin(theta,-3.1415927410125732,3.1415927410125732)]+=1.0; + } + for (var sub=0u;sub<3u;sub+=1u) { + var sum=0.0; + for (var b=0u;b<11u;b+=1u) { sum+=histogram[sub*11u+b]; } + for (var b=0u;b<11u;b+=1u) { + let value=histogram[sub*11u+b]; + if sum>0.0 { spfh[i*33u+sub*11u+b]=value/sum*100.0; } + } + } +} + +@compute @workgroup_size(64) +fn build_fpfh(@builtin(global_invocation_id) id: vec3) { + let i=id.x; + if i>=params.count { return; } + if atomicLoad(&status.bad_rows)>0u { + for (var b=0u;b<33u;b+=1u) { descriptors[i*33u+b]=0.0; } + return; + } + let base=i*(2u+2u*params.capacity); + let stored=neighbors[base]; + var acc: array; + var sum_weight=0.0; + for (var n=0u;n(neighbors[base+3u+2u*n])); + if distance<1e-9 { continue; } + let weight=1.0/distance; + sum_weight+=weight; + for (var b=0u;b<33u;b+=1u) { acc[b]+=weight*spfh[j*33u+b]; } + } + for (var b=0u;b<33u;b+=1u) { + var contribution=0.0; + if sum_weight>0.0 { contribution=acc[b]/sum_weight; } + descriptors[i*33u+b]=spfh[i*33u+b]+contribution; + } +} diff --git a/crates/spatialrust-registration/Cargo.toml b/crates/spatialrust-registration/Cargo.toml index d7ccb576..09fa4594 100644 --- a/crates/spatialrust-registration/Cargo.toml +++ b/crates/spatialrust-registration/Cargo.toml @@ -17,6 +17,7 @@ register-gicp = ["dep:spatialrust-search"] register-gicp-gpu = ["register-gicp", "dep:spatialrust-gpu"] register-ndt = [] register-fpfh = ["dep:spatialrust-search"] +register-fpfh-gpu = ["register-fpfh", "dep:spatialrust-gpu"] [dependencies] spatialrust-core = { workspace = true } diff --git a/crates/spatialrust-registration/src/fpfh.rs b/crates/spatialrust-registration/src/fpfh.rs index 4ec446d6..4677208a 100644 --- a/crates/spatialrust-registration/src/fpfh.rs +++ b/crates/spatialrust-registration/src/fpfh.rs @@ -63,6 +63,43 @@ fn fpfh_descriptors_with_workers( Ok(compute_fpfh_with_workers(&set, feature_radius, workers)) } +/// Explicit GPU FPFH result and payload transfers; ordinary FPFH stays CPU. +#[cfg(feature = "register-fpfh-gpu")] +#[derive(Clone, Debug, PartialEq)] +pub struct GpuFpfhResult { + /// Original-order GPU-computed 33-bin descriptor rows. + pub descriptors: Vec, + /// Target/query/normal uploads plus descriptor/status download bytes. + pub transfers: spatialrust_core::TransferStats, +} + +/// Explicitly uploads a cloud, runs GPU radius search, SPFH and weighted FPFH, +/// then downloads descriptors/status. No automatic CPU fallback occurs. +/// `neighbor_capacity` includes self and must fit ALL radius matches; overflow +/// returns an error rather than descriptors computed from truncated rows. +/// Input validation matches CPU FPFH; GPU angular/bin rounding may differ near +/// boundaries. Descriptor matching and RANSAC registration remain CPU. +#[cfg(feature = "register-fpfh-gpu")] +pub fn fpfh_descriptors_gpu( + runtime: &spatialrust_gpu::WgpuRuntime, + cloud: &PointCloud, + feature_radius: f32, + neighbor_capacity: u32, +) -> SpatialResult { + validate_distance(feature_radius, "feature_radius")?; + let set = PointSet::from_cloud(cloud)?; + let normals: Vec<_> = set.normals.iter().map(|n| [n.x, n.y, n.z]).collect(); + let input = spatialrust_gpu::GpuFpfhInput::upload(runtime, &set.tree, &normals)?; + let neighborhoods = + input.radius_neighborhoods_resident(runtime, feature_radius, neighbor_capacity)?; + let output = input.descriptors_resident(runtime, &neighborhoods)?; + let (descriptors, download) = output.download(runtime)?; + let mut transfers = input.upload_transfers(); + transfers + .record(spatialrust_core::TransferDirection::DeviceToHost, download.device_to_host_bytes()); + Ok(GpuFpfhResult { descriptors, transfers }) +} + /// Configuration for [`FpfhRansacRegistration`]. #[derive(Clone, Copy, Debug, PartialEq)] pub struct FpfhRansacConfig { diff --git a/crates/spatialrust-registration/src/lib.rs b/crates/spatialrust-registration/src/lib.rs index 4ce4e483..656a721b 100644 --- a/crates/spatialrust-registration/src/lib.rs +++ b/crates/spatialrust-registration/src/lib.rs @@ -50,3 +50,6 @@ pub use icp::GpuIcpDiagnostics; #[cfg(feature = "register-gicp-gpu")] pub use gicp::GpuGicpDiagnostics; + +#[cfg(feature = "register-fpfh-gpu")] +pub use fpfh::{fpfh_descriptors_gpu, GpuFpfhResult}; diff --git a/crates/spatialrust-registration/tests/gpu_fpfh.rs b/crates/spatialrust-registration/tests/gpu_fpfh.rs new file mode 100644 index 00000000..bc7d5156 --- /dev/null +++ b/crates/spatialrust-registration/tests/gpu_fpfh.rs @@ -0,0 +1,216 @@ +#![cfg(all(feature = "register-fpfh-gpu", not(target_arch = "wasm32")))] + +use spatialrust_core::{ + DType, FieldSemantic, PointBuffer, PointBufferSet, PointCloud, PointField, PointSchema, + SpatialMetadata, +}; +use spatialrust_gpu::{GpuFpfhInput, WgpuRuntime}; +use spatialrust_registration::{fpfh_descriptors, fpfh_descriptors_gpu}; + +pub fn fixture(n: usize, kind: usize) -> PointCloud { + let mut columns: [Vec; 6] = std::array::from_fn(|_| Vec::new()); + for i in 0..n { + let j = if kind == 2 { i / 2 } else { i }; + let x = (j % 32) as f32 * 0.125 - 2.0; + let y = (j / 32) as f32 * 0.125 - 1.0; + let z = if kind == 0 { 0.0 } else { 0.02 * (x * x + y * y) }; + let nx = if kind == 0 { 0.0 } else { -0.04 * x }; + let ny = if kind == 0 { 0.0 } else { -0.04 * y }; + let len = (nx * nx + ny * ny + 1.0).sqrt(); + let normal = if kind == 3 { + match j % 6 { + 0 => [0.0, 0.0, 1.0], + 1 => [0.0, 1.0, 0.0], + 2 => [1.0, 0.0, 0.0], + 3 => [-1.0, 0.0, 0.0], + 4 => [0.0, -1.0, 0.0], + _ => [0.0, 0.0, -1.0], + } + } else { + [nx / len, ny / len, 1.0 / len] + }; + for (axis, value) in columns.iter_mut().zip([x, y, z, normal[0], normal[1], normal[2]]) { + axis.push(value); + } + } + cloud_from_columns(columns) +} + +pub fn cloud_from_columns(columns: [Vec; 6]) -> PointCloud { + let mut schema = PointSchema::new(); + let mut buffers = PointBufferSet::new(); + for ((name, semantic), values) in [ + ("px", FieldSemantic::PositionX), + ("py", FieldSemantic::PositionY), + ("pz", FieldSemantic::PositionZ), + ("nx", FieldSemantic::NormalX), + ("ny", FieldSemantic::NormalY), + ("nz", FieldSemantic::NormalZ), + ] + .into_iter() + .zip(columns) + { + schema = schema.with_field(PointField::scalar(name, semantic, DType::F32)); + buffers.insert(name, PointBuffer::F32(values)); + } + PointCloud::try_from_parts(schema, buffers, SpatialMetadata::default()).unwrap() +} + +pub fn compare(actual: &[[f32; 33]], expected: &[[f32; 33]], tolerance: f32) -> f32 { + assert_eq!(actual.len(), expected.len()); + let mut max = 0.0_f32; + for (i, (a, b)) in actual.iter().zip(expected).enumerate() { + for (bin, (&a, &b)) in a.iter().zip(b).enumerate() { + assert!(a.is_finite()); + let error = (a - b).abs(); + max = max.max(error); + assert!(error <= tolerance, "row={i} bin={bin} gpu={a} cpu={b} error={error}"); + } + } + max +} + +#[test] +#[ignore = "requires a native wgpu adapter; run explicitly with --include-ignored"] +fn gpu_fpfh_matches_cpu_curved_planar_duplicate_and_partial_clouds() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let mut words = 0; + let mut max_error = 0.0_f32; + for n in [0, 1, 2, 17, 63, 64, 65, 257, 4097] { + for kind in 0..4 { + let cloud = fixture(n, kind); + for radius in [0.16, 0.35] { + let expected = fpfh_descriptors(&cloud, radius).unwrap(); + let actual = fpfh_descriptors_gpu(&runtime, &cloud, radius, 128).unwrap(); + max_error = max_error.max(compare(&actual.descriptors, &expected, 0.002)); + words += n * 33; + 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 * 64 + nodes * 32) as u64 } + ); + assert_eq!( + actual.transfers.device_to_host_bytes(), + if n == 0 { 0 } else { (n * 132 + 16) as u64 } + ); + } + } + } + eprintln!( + "GPU_FPFH verified {words} descriptor values; maximum absolute CPU difference {max_error}" + ); +} + +#[test] +#[ignore = "requires a native wgpu adapter; run explicitly with --include-ignored"] +fn gpu_fpfh_has_independent_planar_and_degenerate_results() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let pair = fixture(2, 0); + let output = fpfh_descriptors_gpu(&runtime, &pair, 0.16, 2).unwrap(); + for descriptor in output.descriptors { + for (bin, value) in descriptor.into_iter().enumerate() { + assert_eq!(value, if bin % 11 == 5 { 200.0 } else { 0.0 }); + } + } + for cloud in [ + fixture(1, 0), + cloud_from_columns([ + vec![0.0; 2], + vec![0.0; 2], + vec![0.0; 2], + vec![0.0; 2], + vec![0.0; 2], + vec![1.0; 2], + ]), + cloud_from_columns([ + vec![0.0; 2], + vec![0.0; 2], + vec![0.0, 0.125], + vec![0.0; 2], + vec![0.0; 2], + vec![1.0; 2], + ]), + ] { + assert!(fpfh_descriptors_gpu(&runtime, &cloud, 0.2, 2) + .unwrap() + .descriptors + .iter() + .flatten() + .all(|&v| v == 0.0)); + } +} + +#[test] +#[ignore = "requires a native wgpu adapter; run explicitly with --include-ignored"] +fn gpu_fpfh_keeps_cpu_validation_and_rejects_capacity_overflow() { + 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_eq!( + fpfh_descriptors(&cloud, radius).unwrap_err().to_string(), + fpfh_descriptors_gpu(&runtime, &cloud, radius, 128).unwrap_err().to_string() + ); + } + for (field, value) in [("px", f32::NAN), ("ny", f32::INFINITY), ("nz", 2.0)] { + let (schema, mut buffers, metadata) = cloud.clone().into_parts(); + let mut values = vec![0.0; 65]; + values[0] = value; + buffers.insert(field, PointBuffer::F32(values)); + let malformed = PointCloud::try_from_parts(schema, buffers, metadata).unwrap(); + assert_eq!( + fpfh_descriptors(&malformed, 0.2).unwrap_err().to_string(), + fpfh_descriptors_gpu(&runtime, &malformed, 0.2, 128).unwrap_err().to_string() + ); + } + for capacity in [0, 1, 3] { + let error = fpfh_descriptors_gpu(&runtime, &cloud, 0.35, capacity).unwrap_err().to_string(); + assert!(error.contains("incomplete/malformed"), "{error}"); + } + let full = fpfh_descriptors_gpu(&runtime, &cloud, 0.35, 128).unwrap(); + compare(&full.descriptors, &fpfh_descriptors(&cloud, 0.35).unwrap(), 0.002); +} + +#[test] +#[ignore = "requires a native wgpu adapter; run explicitly with --include-ignored"] +fn gpu_fpfh_reuses_resident_neighborhoods_with_owned_outputs_and_runtime_guards() { + use spatialrust_core::{HasNormals3, HasPositions3}; + let runtime = WgpuRuntime::new_headless().unwrap(); + let cloud = fixture(257, 1); + let (x, y, z) = cloud.positions3().unwrap(); + let (nx, ny, nz) = cloud.normals3().unwrap(); + let tree = spatialrust_search::KdTree::from_slices(x, y, z); + let normals: Vec<_> = (0..nx.len()).map(|i| [nx[i], ny[i], nz[i]]).collect(); + let input = GpuFpfhInput::upload(&runtime, &tree, &normals).unwrap(); + let neighbors = input.radius_neighborhoods_resident(&runtime, 0.35, 128).unwrap(); + let expected = fpfh_descriptors(&cloud, 0.35).unwrap(); + let live = input.descriptors_resident(&runtime, &neighbors).unwrap(); + std::thread::scope(|scope| { + for _ in 0..4 { + let (runtime, input, neighbors, expected) = (&runtime, &input, &neighbors, &expected); + scope.spawn(move || { + for _ in 0..2 { + compare( + &input + .descriptors_resident(runtime, neighbors) + .unwrap() + .download(runtime) + .unwrap() + .0, + expected, + 0.002, + ); + } + }); + } + }); + let other = WgpuRuntime::new_headless().unwrap(); + assert!(input.descriptors_resident(&other, &neighbors).is_err()); + assert!(live.download(&other).is_err()); + let second = GpuFpfhInput::upload(&runtime, &tree, &normals).unwrap(); + assert!(second.descriptors_resident(&runtime, &neighbors).is_err()); + drop(input); + drop(neighbors); + drop(tree); + compare(&live.download(&runtime).unwrap().0, &expected, 0.002); +} diff --git a/crates/spatialrust-registration/tests/gpu_fpfh_benchmark.rs b/crates/spatialrust-registration/tests/gpu_fpfh_benchmark.rs new file mode 100644 index 00000000..e0b72da3 --- /dev/null +++ b/crates/spatialrust-registration/tests/gpu_fpfh_benchmark.rs @@ -0,0 +1,72 @@ +#![cfg(all(feature = "register-fpfh-gpu", not(target_arch = "wasm32")))] + +#[path = "gpu_fpfh.rs"] +mod reference; + +use spatialrust_core::HasNormals3; +use spatialrust_gpu::{GpuFpfhInput, WgpuRuntime}; +use spatialrust_registration::{fpfh_descriptors, fpfh_descriptors_gpu}; +use std::time::Instant; + +#[test] +#[ignore = "completed FPFH benchmark; run this exact test alone in release mode"] +fn completed_fpfh_calls() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let info = runtime.adapter_info(); + println!("GPU_FPFH_ENV backend={} device_type={}", info.backend, info.device_type); + for n in [257, 4097, 16384, 65537] { + for kind in [0, 1, 3] { + let radius = 0.35; + let capacity = 64; + let cloud = reference::fixture(n, kind); + let expected = fpfh_descriptors(&cloud, radius).unwrap(); + let tree = spatialrust_search::KdTree::from_point_cloud(&cloud).unwrap(); + let (nx, ny, nz) = cloud.normals3().unwrap(); + let normals: Vec<_> = (0..n).map(|i| [nx[i], ny[i], nz[i]]).collect(); + let input = GpuFpfhInput::upload(&runtime, &tree, &normals).unwrap(); + let neighborhoods = + input.radius_neighborhoods_resident(&runtime, radius, capacity).unwrap(); + let cpu = || fpfh_descriptors(&cloud, radius).unwrap(); + let gpu = + || fpfh_descriptors_gpu(&runtime, &cloud, radius, capacity).unwrap().descriptors; + for _ in 0..3 { + reference::compare(&cpu(), &expected, 0.002); + reference::compare(&gpu(), &expected, 0.002); + reference::compare( + &input + .descriptors_resident(&runtime, &neighborhoods) + .unwrap() + .download(&runtime) + .unwrap() + .0, + &expected, + 0.002, + ); + } + let (mut cpu_ms, mut gpu_ms, mut resident_ms) = (Vec::new(), Vec::new(), Vec::new()); + 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.descriptors_resident(&runtime, &neighborhoods).unwrap(); + runtime.wait_idle(); + resident_ms.push(start.elapsed().as_secs_f64() * 1000.0); + reference::compare(&output.download(&runtime).unwrap().0, &expected, 0.002); + } else { + let rows = if stage == 0 { cpu() } else { gpu() }; + let elapsed = start.elapsed().as_secs_f64() * 1000.0; + if stage == 0 { + cpu_ms.push(elapsed) + } else { + gpu_ms.push(elapsed) + }; + reference::compare(&rows, &expected, 0.002); + } + } + } + let workers = spatialrust_search::parallel_worker_count(n); + println!("GPU_FPFH_BENCH n={n} kind={kind} workers={workers} cpu_ms={cpu_ms:?} gpu_ms={gpu_ms:?} resident_ms={resident_ms:?}"); + } + } +} diff --git a/crates/spatialrust/Cargo.toml b/crates/spatialrust/Cargo.toml index cc41d919..02503e83 100644 --- a/crates/spatialrust/Cargo.toml +++ b/crates/spatialrust/Cargo.toml @@ -106,6 +106,7 @@ register-fpfh = [ "spatialrust-registration/register-fpfh", "spatialrust-search/search-kdtree", ] +register-fpfh-gpu = ["register-fpfh", "spatialrust-registration/register-fpfh-gpu", "gpu-wgpu"] mvp = [ "io-pcd", "io-ply", diff --git a/crates/spatialrust/tests/gpu_fpfh_smoke.rs b/crates/spatialrust/tests/gpu_fpfh_smoke.rs new file mode 100644 index 00000000..b659668c --- /dev/null +++ b/crates/spatialrust/tests/gpu_fpfh_smoke.rs @@ -0,0 +1,36 @@ +#![cfg(all(feature = "register-fpfh-gpu", not(target_arch = "wasm32")))] +use spatialrust::core::{DType, FieldSemantic, PointCloudBuilder, PointField, PointSchema}; +use spatialrust::gpu::WgpuRuntime; +use spatialrust::registration::{fpfh_descriptors, fpfh_descriptors_gpu}; + +#[test] +#[ignore = "requires a native wgpu adapter; run explicitly with --include-ignored"] +fn facade_gpu_fpfh_preserves_planar_descriptors_and_records_transfers() { + let mut schema = PointSchema::new(); + for (name, semantic) in [ + ("x", FieldSemantic::PositionX), + ("y", FieldSemantic::PositionY), + ("z", FieldSemantic::PositionZ), + ("nx", FieldSemantic::NormalX), + ("ny", FieldSemantic::NormalY), + ("nz", FieldSemantic::NormalZ), + ] { + schema = schema.with_field(PointField::scalar(name, semantic, DType::F32)); + } + let mut builder = PointCloudBuilder::new(schema); + for i in 0..257 { + builder + .push_point([(i % 17) as f32 * 0.125, (i / 17) as f32 * 0.125, 0.0, 0.0, 0.0, 1.0]) + .unwrap(); + } + let cloud = builder.build().unwrap(); + let runtime = WgpuRuntime::new_headless().unwrap(); + let expected = fpfh_descriptors(&cloud, 0.35).unwrap(); + let actual = fpfh_descriptors_gpu(&runtime, &cloud, 0.35, 64).unwrap(); + assert!(actual.transfers.host_to_device_bytes() >= 257 * 64); + assert_eq!(actual.transfers.device_to_host_bytes(), 257 * 132 + 16); + for (a, b) in actual.descriptors.iter().flatten().zip(expected.iter().flatten()) { + assert!((a - b).abs() <= 0.002); + } + assert!(fpfh_descriptors_gpu(&runtime, &cloud, 0.35, 1).is_err()); +} diff --git a/docs/ARCHITECTURE.md b/docs/ARCHITECTURE.md index 6bc167b0..d41d827f 100644 --- a/docs/ARCHITECTURE.md +++ b/docs/ARCHITECTURE.md @@ -30,7 +30,17 @@ 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; FPFH/ISS/MLS device compute follows. +not change feature algorithm placement; ISS/MLS device compute follows. + +FPFH device computation now consumes this radius-row representation through +`GpuFpfhInput` and neighborhoods tied to that exact uploaded input. Explicit +input upload retains original positions/unit normals plus target/query buffers. +SPFH normalization and ordered distance-weighted aggregation run on GPU in two +passes sharing one submission. GPU status detects capacity overflow or malformed +rows before downstream feature use; invalid results are zeroed and explicit +download rejects them. No neighborhood/SPFH readback is needed. The optional +`register-fpfh-gpu` wrapper retains CPU validation/tree construction and exposes +complete transfer receipts; CPU descriptor matching and RANSAC are unchanged. | Area | Decision | | --- | --- | diff --git a/docs/FEATURE_MATRIX.md b/docs/FEATURE_MATRIX.md index 451d67b8..bfa18541 100644 --- a/docs/FEATURE_MATRIX.md +++ b/docs/FEATURE_MATRIX.md @@ -6,6 +6,20 @@ should normally depend on `spatialrust` and select only the profiles they need. ## User-facing profiles +`register-fpfh-gpu` exposes +`fpfh_descriptors_gpu(runtime, cloud, radius, neighbor_capacity)` alongside the +ordinary CPU descriptor API. CPU validation and KD-tree construction precede +explicit GPU input upload; radius search, SPFH and weighted FPFH run on device. +Capacity includes self and must fit every radius match; incomplete/malformed +rows set GPU status, zero all descriptors and reject at readback. No automatic +CPU fallback or placement change occurs. `GpuFpfhInput` under `gpu-wgpu` also +exposes input-bound reusable `GpuFpfhNeighborhoods` and resident descriptors +for later device consumers, which must check the GPU status buffer. The two +feature passes share a submission without neighborhood/SPFH readback. Output +transfer is 132 bytes per point plus 16 status bytes; GPU angular rounding can +change bins near boundaries. Descriptor matching/RANSAC remain CPU. See +[FPFH measurements](../notes/2026-10-12_gpu_fpfh_descriptors.md). + `gpu-wgpu` exposes `GpuQueryBatch::upload` and `GpuKdTree::{nearest_k_resident,radius_neighbors_resident}` for shared GPU neighborhoods. Target/query uploads are explicit and reusable; search submits @@ -16,8 +30,8 @@ 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. Feature algorithms remain separate CPU -APIs until their device compute is added. Native correctness and completed/ +and preserve query/target semantics. GPU FPFH consumes 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). @@ -76,6 +90,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 and wgpu correspondence search | Explicit GPU target upload and per-batch query/readback; CPU pose solver | | `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 | | `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 eb9f49ec..a85e9e3d 100644 --- a/docs/ROADMAP.md +++ b/docs/ROADMAP.md @@ -1128,6 +1128,7 @@ 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; feature consumers planned | 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/ISS/MLS compute | Planned | Device feature/smoothing phases on the shared neighborhood representation | CPU/math references, attribute/schema preservation, numeric tolerances where reductions change, transfer-inclusive benefit | +| Feature neighborhoods | Resident API complete; FPFH consumer 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/MLS compute | Planned | Device keypoint/smoothing phases on the shared neighborhood representation | CPU/math references, attribute/schema preservation, numeric tolerances where reductions change, transfer-inclusive benefit | | 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-12_gpu_fpfh_descriptors.md b/notes/2026-10-12_gpu_fpfh_descriptors.md new file mode 100644 index 00000000..0f266f45 --- /dev/null +++ b/notes/2026-10-12_gpu_fpfh_descriptors.md @@ -0,0 +1,167 @@ +# Resident GPU FPFH descriptors + +Date: 2026-10-12. Named slice: GPU radius neighborhoods, SPFH and weighted +FPFH descriptors, following the resident neighborhood API in PR #166. + +## Execution and transfer contract + +The optional `register-fpfh-gpu` feature exposes +`fpfh_descriptors_gpu(runtime, cloud, radius, neighbor_capacity)` in registration +and the facade. The caller chooses GPU execution explicitly. This wrapper uses +the existing CPU input validation, semantic position/normal access and KD-tree +construction, then uploads, computes and downloads. Ordinary FPFH, descriptor +matching, RANSAC/global registration, ISS and MLS retain their current CPU paths. + +`GpuFpfhInput::upload` uploads original-order XYZ/unit normals together with a +GPU KD-tree and reusable original-order query batch. Positions and normals must +be finite, normals must satisfy the existing unit-length tolerance, and global +coordinate extents must have a finite f32 squared length. Counts, host sizes, +device buffer sizes, dispatch limits and bindings are checked before allocation. +The descriptor shader requires five storage bindings, one uniform binding and +64 invocations per workgroup. No new dependency, core API or unsafe code is added. + +`radius_neighborhoods_resident` computes radius rows once without readback or a +completion wait. Its typed result is tied to that exact uploaded input; another +input's neighborhoods are rejected even on the same runtime. Radius capacity +includes self and must accommodate ALL matches. An incomplete row is an error, +not an approximate descriptor computed from a truncated prefix. Radius must be +positive with a finite, nonzero f32 square. + +`descriptors_resident` consumes those rows in two ordered compute passes in one +submission. Each original point builds three normalized eleven-bin SPFH +histograms. The next pass adds ordered inverse-distance-weighted neighboring +SPFH values. Self, coincident points and degenerate Darboux frames follow the +existing CPU exclusion rules. Signed-zero `atan2` cases are handled explicitly. +SPFH storage is private scratch, completely written before reuse; returned +descriptor/status buffers are owned independently and never pooled. Calls can +reuse uploaded inputs/neighborhoods, retain earlier outputs and run concurrently. + +GPU status contains four u32 words: bad row count, last bad row plus one, +maximum required radius capacity, and malformed row count. SPFH validates row +counts and target indices before accessing samples. If any row is incomplete or +malformed, the second pass zeros ALL descriptor rows. Custom GPU consumers must +inspect status before accepting results and treat exposed storage as immutable +application data. Explicit download rejects bad status before decoding rows. + +For N points and M CPU-tree nodes, input upload is `64*N + 32*M` H2D payload +bytes: 32*N samples, 16*N queries and 16*N tree positions plus 32*M nodes. +Resident neighborhood/descriptor execution transfers zero payload bytes. +Explicit download copies `132*N + 16` D2H bytes in one submission, including +33-f32 descriptors and status, and waits for completion. Empty input has zero +payload transfers. Uniform controls are excluded from payload receipts. There +is no neighborhood/SPFH readback or automatic CPU fallback. + +GPU angular arithmetic, `atan2` approximation and bin selection may differ from +CPU near histogram boundaries. This API does not promise bitwise parity or a +universal small numerical bound across inputs/backends. The observed error below +applies to the selected fixtures. WASM validation is compilation only; browser +execution of the synchronous download path was not tested. + +## Correctness and regression evidence + +- Three new GPU primitive tests cover limits, invalid inputs/empty receipts, + overflow and an externally poisoned target index. Exact status and complete + zeroed output are checked, followed by successful scratch reuse. +- Four explicitly run adapter integration tests compare 1,205,424 descriptor + values across 0/1/2/17/63/64/65/257/4097 points, two radii, planar/curved/ + coincident-pair fixtures and alternating signed axis normals. Absolute + tolerance is 0.002; maximum observed CPU difference is 0.000045776367. +- Independent planar-pair descriptors have exactly 200 in each central bin; + isolated/coincident/parallel-normal degeneracies are zero. Tests also verify + CPU validation errors, capacity errors, exact transfer receipts, concurrent + calls, retained outputs and runtime/input ownership guards. +- All-feature GPU regression tests: 199 passed, 23 ignored. Registration library + regression tests: 35 passed, one ignored. The facade GPU integration test + passed. After the WASM readback portability fix, the four adapter integration + tests and facade test passed again. +- All-feature GPU/registration Clippy with warnings denied, CPU-only registration + Clippy, minimal GPU compilation, WASM GPU-FPFH compilation and workspace + formatting passed. Source comparison confirms existing CPU FPFH math and + registration code are unchanged. CI adds isolated registration/facade feature + profiles; hosted CI compiles adapter tests rather than proving GPU execution. + +## Completed-call measurements + +Release builds on a native Vulkan discrete-device backend. Adapter model names +are deliberately omitted. Three warmups for each stage precede 20 paired rounds; +CPU/full-GPU/resident stage order rotates each round. Two surveys ran separately +without another local build, test or statistical analysis during measurement. +All samples are retained, including adverse and variable results. + +The CPU baseline is the current public `fpfh_descriptors`, including its existing +bounded CPU parallelism: four workers at 65,537 points, one at the smaller sizes. +Full GPU calls include CPU validation/lifting/tree construction, all input +uploads, GPU radius search, SPFH, weighted FPFH and descriptor/status download +with completion wait. Resident timings include descriptor allocation, both +descriptor passes, submission and `wait_idle`; they exclude CPU preparation, +input upload, neighborhood generation and final download. Correctness checks +and dropping returned descriptor values happen outside timers. Runtime/pipeline +creation is warmed outside timing. Resident columns are not complete-call +speedups, GPU timestamp measurements or global-registration timings. + +Fixtures use radius 0.35 and capacity 64: P is a planar grid with common +Z +normals; C is a curved grid with analytic normals; A uses the curved positions +with alternating signed axis normals. Each fixture uses 32 columns and spacing +0.125. Reported times are medians in milliseconds. Speedup is the geometric +mean of the 20 paired CPU/full-GPU ratios, with a percentile 95% bootstrap +interval from 10,000 resamples of paired log ratios (deterministic seeds +20261012 plus the zero-based record index). It is not a ratio of medians and +does not establish deployment performance or a crossover threshold. + +| Survey | Points | Fixture | CPU workers | CPU ms | Full GPU ms | Resident ms | Paired full-call speedup (95% interval) | +| --- | ---: | :---: | ---: | ---: | ---: | ---: | --- | +| 1 | 257 | P | 1 | 0.565 | 1.570 | 0.239 | 0.39x (0.35-0.42) | +| 1 | 257 | C | 1 | 0.633 | 1.265 | 0.256 | 0.50x (0.49-0.52) | +| 1 | 257 | A | 1 | 0.704 | 1.284 | 0.258 | 0.56x (0.53-0.58) | +| 1 | 4,097 | P | 1 | 13.434 | 2.906 | 0.430 | 4.74x (4.55-4.94) | +| 1 | 4,097 | C | 1 | 13.649 | 3.139 | 0.435 | 4.33x (4.08-4.60) | +| 1 | 4,097 | A | 1 | 13.712 | 3.057 | 0.423 | 4.81x (4.49-5.16) | +| 1 | 16,384 | P | 1 | 52.217 | 9.765 | 4.880 | 6.85x (5.85-8.07) | +| 1 | 16,384 | C | 1 | 47.275 | 25.295 | 13.840 | 1.88x (1.76-2.02) | +| 1 | 16,384 | A | 1 | 50.727 | 25.147 | 13.767 | 2.01x (1.89-2.13) | +| 1 | 65,537 | P | 4 | 87.486 | 29.268 | 4.425 | 2.43x (2.04-2.90) | +| 1 | 65,537 | C | 4 | 62.074 | 23.937 | 3.032 | 2.44x (2.11-2.76) | +| 1 | 65,537 | A | 4 | 62.031 | 45.835 | 22.891 | 1.73x (1.49-2.02) | +| 2 | 257 | P | 1 | 0.586 | 1.113 | 0.224 | 0.52x (0.51-0.53) | +| 2 | 257 | C | 1 | 0.636 | 1.266 | 0.275 | 0.50x (0.49-0.52) | +| 2 | 257 | A | 1 | 0.647 | 1.102 | 0.224 | 0.57x (0.56-0.59) | +| 2 | 4,097 | P | 1 | 11.668 | 2.575 | 0.380 | 4.70x (4.50-4.92) | +| 2 | 4,097 | C | 1 | 11.257 | 2.466 | 0.368 | 4.62x (4.47-4.78) | +| 2 | 4,097 | A | 1 | 11.844 | 2.601 | 0.373 | 4.69x (4.44-4.96) | +| 2 | 16,384 | P | 1 | 48.882 | 4.940 | 1.154 | 10.00x (9.59-10.43) | +| 2 | 16,384 | C | 1 | 50.195 | 24.807 | 13.671 | 2.56x (2.29-2.87) | +| 2 | 16,384 | A | 1 | 46.026 | 25.883 | 13.921 | 1.91x (1.79-2.04) | +| 2 | 65,537 | P | 4 | 88.117 | 23.325 | 4.438 | 3.31x (2.85-3.77) | +| 2 | 65,537 | C | 4 | 71.337 | 50.672 | 28.314 | 1.61x (1.39-1.91) | +| 2 | 65,537 | A | 4 | 64.949 | 21.837 | 3.024 | 2.86x (2.59-3.08) | + +Timings vary substantially between surveys and stages, and resident results +are not monotonic in point count. The data does not establish the cause of +that variation. Small fixtures can favor CPU; warmed resident results must not +be used to imply transfer-inclusive gains on arbitrary inputs. + +## Reproduction + +From `C:/Users/rsasa/Workspace/SpatialRust`, run adapter correctness explicitly: + +```powershell +cargo test -p spatialrust-registration --no-default-features --features register-fpfh-gpu --test gpu_fpfh -- --include-ignored --nocapture +cargo test -p spatialrust --no-default-features --features register-fpfh-gpu --test gpu_fpfh_smoke -- --include-ignored +``` + +Run this exact benchmark filter alone twice. The filter excludes the included +reference module's four adapter tests from timing runs: + +```powershell +cargo test -p spatialrust-registration --no-default-features --features register-fpfh-gpu --test gpu_fpfh_benchmark --release completed_fpfh_calls -- --include-ignored --nocapture +``` + +Local untracked receipts: + +- `C:/Users/rsasa/Workspace/SpatialRust/target/automation/gpu-fpfh-final-integration-tests.log` +- `C:/Users/rsasa/Workspace/SpatialRust/target/automation/gpu-fpfh-facade-tests.log` +- `C:/Users/rsasa/Workspace/SpatialRust/target/automation/gpu-fpfh-native-tests.log` +- `C:/Users/rsasa/Workspace/SpatialRust/target/automation/gpu-fpfh-registration-tests.log` +- `C:/Users/rsasa/Workspace/SpatialRust/target/automation/gpu-fpfh-benchmark-1.log` +- `C:/Users/rsasa/Workspace/SpatialRust/target/automation/gpu-fpfh-benchmark-2.log` +- `C:/Users/rsasa/Workspace/SpatialRust/target/automation/gpu-fpfh-measurements.json`